31 #include "TMultiGraph.h" 41 #include "TRotation.h" 42 #include "Math/Vector3Dfwd.h" 43 #include "Math/Rotation3D.h" 50 #define COORDINATES "Geographic" 61 #define ANGULAR_DISTANCE 90 65 #define DISPLAY_WORLD_MAP 66 #define WORLD_MAP_DIR "$CWB_GWAT/data/" 78 if (!gROOT->GetClass(
"Polar3DVector")) gSystem->Load(
"libMathCore");
87 int nIFO = ifo.size();
91 if(ifo[
n]!=
"") pD[
n] =
new detector((
char*)ifo[
n].Data());
94 for(
int n=0;
n<
nIFO;
n++) cout <<
n <<
" " << pD[
n]->Name << endl;
112 h2->GetXaxis()->SetTitleSize(0.05);
113 h2->GetXaxis()->SetLabelSize(0.05);
114 h2->GetYaxis()->SetTitleSize(0.05);
115 h2->GetYaxis()->SetLabelSize(0.05);
116 h2->GetYaxis()->SetLabelFont(42);
117 h2->GetYaxis()->SetLabelFont(42);
118 h2->GetXaxis()->SetTitleFont(42);
119 h2->GetYaxis()->SetTitleFont(42);
120 h2->GetZaxis()->SetRangeUser(0,1.0);
133 sprintf(ofileName,
"%s/antpat_alignment.png",odir.Data());
134 cout <<
"Write : " << ofileName << endl;
135 gSM->
Print(ofileName);
143 sprintf(ofileName,
"%s/antpat_sensitivity.png",odir.Data());
144 cout <<
"Write : " << ofileName << endl;
145 gSM->
Print(ofileName);
151 cout <<
"order : " << gSM->
getOrder() << endl;
152 cout <<
"size : " << gSM->
size() << endl;
153 int apsize = gSM->
size();
154 double* ap =
new double[apsize];
155 for(
int i=0;i<apsize;i++) ap[i] = gSM->
get(i);
157 Int_t *iap =
new Int_t[apsize];
158 TMath::Sort(apsize,ap,iap,
true);
165 XYZVector
D1(pD[0]->Rv[0],pD[0]->Rv[1],pD[0]->Rv[2]);
166 XYZVector D2(pD[1]->Rv[0],pD[1]->Rv[1],pD[1]->Rv[2]);
169 XYZVector D12 = D1-D2;
170 TVector3 vD12(D12.X(),D12.Y(),D12.Z());
172 double th12 = r2d*vD12.Theta();
173 double ph12 = r2d*vD12.Phi();
174 cout <<
"coordinates D12 " << ph12 <<
" " << th12 << endl;
175 Polar3DVector v12(1, d2r*th12, d2r*ph12);
180 char wave_file_name[1024];
181 sprintf(wave_file_name,
"merge/wave_%s.%s.root",data_label.Data(),merge_label.Data());
183 TFile*
ifile = TFile::Open(wave_file_name);
184 if(ifile==NULL) {cout<<
"Error opening file : "<<wave_file_name<<endl;
exit(1);}
185 TTree*
itree = (TTree *) gROOT->FindObject(
"waveburst");
186 if(itree==NULL) {cout<<
"Error opening tree : "<<
"waveburst"<<endl;
exit(1);}
190 sprintf(cut,
"abs(time[0]-time[%d])<%f && netcc[%d]>%f && rho[%d]>%f",
191 nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
192 if(T_vED>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && neted[0]/ecor<%f",tmp,T_vED);}
193 if(T_pen>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && penalty>%f",tmp,T_pen);}
194 if(T_ifar>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}
196 itree->Draw(
"theta[0]:phi[0]:theta[1]:phi[1]",cut,
"goff");
197 int isize=itree->GetSelectedRows();
198 cout <<
"isize : " << isize << endl;
200 double* th0 = itree->GetV1();
201 double* ph0 = itree->GetV2();
202 double* th1 = itree->GetV3();
203 double* ph1 = itree->GetV4();
205 #ifdef DRAW_ANGULAR_DISTANCE 206 TH1F* h1 =
new TH1F(
"hist",
"hist",100,0,180);
215 double markerSize = isize>10000 ? 0.3 : 0.4;
217 for (
int i=0;i<
isize;i++) {
220 Polar3DVector v0(1, d2r*th0[i], d2r*ph0[i]);
221 Polar3DVector v1(1, d2r*th1[i], d2r*ph1[i]);
223 double dot = v0.Dot(v1);
224 double dOmega = r2d*TMath::ACos(dot);
227 double dot12 = v12.Dot(v1);
228 double dOmega12 = r2d*TMath::ACos(dot12);
233 #ifdef DRAW_ANGULAR_DISTANCE 237 if(binj) {ph=ph1[
i]; th=th1[
i];}
238 else {ph=ph0[
i]; th=th0[
i];}
241 gSM->
DrawMarker(ph, th, 20, markerSize, kBlack);
246 gSM->
DrawMarker(phi,theta, 20, markerSize, kBlack);
249 gEVT->
set(index,gEVT->
get(index)+1);
260 gNET->
DrawCircles(phi,theta,(Color_t)kWhite,1,1,
true);
265 #ifdef DRAW_ANGULAR_DISTANCE 266 gStyle->SetLineColor(kBlack);
273 sprintf(ofname,
"%s/inj_detected_antpat.png",odir.Data());
275 sprintf(ofname,
"%s/rec_detected_antpat.png",odir.Data());
277 cout <<
"Write : " << ofname << endl;
282 TCanvas*
canvas =
new TCanvas;
286 double*
x =
new double[apsize];
287 double*
y =
new double[apsize];
288 for(
int i=0;i<apsize;i++) {
290 x[
i]=pow(ap[iap[i]],4);
291 y[
i]=gEVT->
get(iap[i]);
295 for(
int i=1;i<apsize;i++) {
299 for(
int i=0;i<apsize;i++) {
303 TGraph*
gr =
new TGraph(apsize,x,y);
305 gr->SetLineColor(kRed);
306 gr->SetMarkerColor(kRed);
308 TH1F*
hist = gr->GetHistogram();
309 hist->SetStats(kFALSE);
312 hist->SetTitle(title);
314 hist->GetXaxis()->SetTitle(
"(F+^{2}+Fx^{2})^{2} Probability");
315 hist->GetYaxis()->SetTitle(
"Frequentist Probability");
316 hist->GetXaxis()->SetRangeUser(0,1);
317 hist->GetYaxis()->SetRangeUser(0,1);
detector * getifo(size_t n)
param: detector index
size_t add(detector *)
param: detector structure return number of detectors in the network
void DrawSitesArms(double mlength=600000., Color_t lcolor=kBlack, Size_t lwidth=1.0, Style_t lstyle=1)
void DrawAntennaPattern(int polarization=-1, int dpaletteId=0, bool btitle=true, int order=6)
void DrawCircles(double phi, double theta, double gps, Color_t lcolor=kBlack, Width_t lwidth=1, Style_t lstyle=1, bool labels=false)
void set(size_t i, double a)
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
void DrawSkyDistributionPRC(TString data_label, TString odir, TString merge_label, vector< TString > ifo, detectorParams *dP, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar, bool binj=true, TString polarization="TENSOR", bool save_antpat=false)
void DrawMarker(double phi, double theta, int marker, Size_t msize=1, Color_t tcolor=1)
void setPolarization(POLARIZATION polarization=TENSOR)
size_t getSkyIndex(double th, double ph)
param: theta param: phi
void DrawSitesShortLabel(Color_t tcolor=kBlack, Size_t tsize=0.052, Font_t tfont=32)
void DrawSites(Color_t mcolor=kBlack, Size_t msize=2.0, Style_t mstyle=20)
void SetWorldMapPath(TString worldMapPath)
void SetWorldMap(bool drawWorldMap=true)
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
strcpy(RunLabel, RUN_LABEL)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double get(size_t i)
param: sky index
void Print(TString pname)
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
double SpeedOfLightInVacuo()