16 #define MAX_FILES 20 // num max of file in the input file list
22 #define TREE_NAME "waveburst"
24 #define MIN_SNR_NET 5 // SNRnet inf
25 #define MAX_SNR_NET 100 // SNRnet sup
26 #define NSNR 20 // number of snr bins
28 #define MIN_NSNR 20 // minimum number of entries in the SNR bin
29 #define MAX_MEDIAN50 50 // max median50 to be included in the fit
30 #define MAX_MEDIAN90 100 // max median90 to be included in the fit
32 #define LINX // set Logx to false
46 if((ofName!=
"")&&(!ofName.Contains(
".png"))) {
47 cout <<
"Error Output File : " << ilist_file.Data() << endl;
48 cout <<
"Must have .png extension" << endl;
51 ofName.ReplaceAll(
".png",
".gif");
54 TString ptitle=
"Sky Localization : median error area vs SNR";
56 if((pType!=
"MEDIAN50")&&(pType!=
"MEDIAN90")&&(pType!=
"")) {
59 if(pType==
"") pType=
"MEDIAN50";
60 if(pType==
"MEDIAN50") ptitle=ptitle+
"50% error region";
61 if(pType==
"MEDIAN90") ptitle=ptitle+
"90% error region";
72 in.open(ilist_file.Data(),
ios::in);
73 if (!in.good()) {cout <<
"Error Opening File : " << ilist_file.Data() << endl;
exit(1);}
80 if (!in.good())
break;
81 if(str[0] !=
'#') size++;
83 cout <<
"size " << size << endl;
84 in.clear(ios::goodbit);
85 in.seekg(0, ios::beg);
86 if (size==0) {cout <<
"Error : File " << ilist_file.Data() <<
" is empty" << endl;
exit(1);}
95 in >> sfile >> sname >> stype >> scolors >> sltyle;
97 if(sfile[0]==
'#')
continue;
101 colors[NFILE]=scolors;
102 lstyle[NFILE]=sltyle;
103 cout << NFILE+1 <<
" " << ifile[NFILE] <<
" " << name[NFILE] <<
" "
104 << type[NFILE] <<
" " << colors[NFILE] <<
" " << lstyle[NFILE] << endl;
111 f1=
new TF1(
"medianfit",
medianfit,0.5,4,3);
112 f1->SetParNames(
"A",
"B",
"C");
123 for(
int n=0;
n<NFILE;
n++) {
124 cout <<
"Process File : " << ifile[
n].Data() << endl;
125 nevents[
n]=
GetMedian(ifile[
n], type[n], pType, snr[n], median[n], entries[n],
129 if(median[n][
i]>max_median) max_median=median[
n][
i];
136 gStyle->SetFrameBorderMode(0);
139 gStyle->SetTitleFont(72);
140 gStyle->SetMarkerColor(50);
141 gStyle->SetLineColor(kWhite);
142 gStyle->SetTitleW(0.98);
143 gStyle->SetTitleH(0.05);
144 gStyle->SetTitleY(0.98);
145 gStyle->SetFillColor(kWhite);
146 gStyle->SetLineColor(kWhite);
147 gStyle->SetTitleFont(12,
"D");
149 TCanvas *
canvas =
new TCanvas(
"median",
"median", 300,40, 600, 600);
151 canvas->ToggleEventStatus();
153 canvas->SetLogx(
false);
155 canvas->SetLogx(
true);
160 canvas->SetFillColor(kWhite);
169 for(
int n=0;
n<NFILE;
n++) {
173 gr[
n] =
new TGraphErrors(snr[
n].
size(),snr[
n].data,median[
n].data,esnr.
data,emedian.
data);
175 for(
int n=0;
n<NFILE;
n++) {
176 gr[
n]->SetLineWidth(0);
177 gr[
n]->SetLineColor(colors[
n]);
178 gr[
n]->SetLineStyle(1);
179 gr[
n]->SetMarkerColor(colors[n]);
180 gr[
n]->SetMarkerStyle(20);
181 gr[
n]->SetMarkerSize(1);
184 f1->SetLineColor(colors[n]);
185 f1->SetLineStyle(lstyle[n]);
188 gr[
n]->Fit(
"medianfit");
189 chi2[
n]=f1->GetChisquare();
190 A[
n]=f1->GetParameter(0);
191 B[
n]=f1->GetParameter(1);
192 C[
n]=f1->GetParameter(2);
196 TMultiGraph*
mg =
new TMultiGraph();
197 TString gTitle =
"Sky Localization "+pType+
" vs SNR";
198 if(pTitle!=
"") gTitle = gTitle+
" : "+pTitle;
199 mg->SetTitle(gTitle.Data());
200 for(
int n=0;
n<NFILE;
n++) mg->Add(gr[
n]);
206 mg->GetHistogram()->GetXaxis()->SetTitle(
"Injected SNRnet");
207 if(pType==
"MEDIAN50") mg->GetHistogram()->GetYaxis()->SetTitle(
"Median50 error angle (degrees)");
208 else mg->GetHistogram()->GetYaxis()->SetTitle(
"Median90 error angle (degrees)");
210 mg->GetHistogram()->GetXaxis()->SetLabelSize(0.04);
211 mg->GetHistogram()->GetYaxis()->SetLabelSize(0.04);
212 mg->GetHistogram()->GetXaxis()->SetTitleSize(0.04);
213 mg->GetHistogram()->GetYaxis()->SetTitleSize(0.04);
214 mg->GetHistogram()->GetXaxis()->SetTitleOffset(1.2);
215 mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.2);
216 mg->GetHistogram()->GetXaxis()->CenterTitle(
true);
217 mg->GetHistogram()->GetYaxis()->CenterTitle(
true);
218 mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
219 mg->GetHistogram()->GetXaxis()->SetTitleFont(42);
220 mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
221 mg->GetHistogram()->GetYaxis()->SetTitleFont(42);
227 double hleg = 0.8-NFILE*0.05;
228 leg =
new TLegend(0.5793173,hleg,0.9915385,0.8721805,NULL,
"brNDC");
230 leg->SetBorderSize(1);
231 leg->SetTextAlign(22);
232 leg->SetTextFont(12);
233 leg->SetLineColor(1);
234 leg->SetLineStyle(1);
235 leg->SetLineWidth(1);
236 leg->SetFillColor(0);
237 leg->SetFillStyle(1001);
238 leg->SetTextSize(0.05);
239 leg->SetLineColor(kBlack);
240 leg->SetFillColor(kWhite);
242 for(
int n=0;n<NFILE;n++) {
243 double median_snr15 = A[
n]+pow(
fabs(B[n])/15.,1)+pow(C[n]/15.,2);
246 sprintf(legLabel,
"SNR(15) = %2.1f",median_snr15);
247 leg->AddEntry(gr[n],legLabel,
"lp");
255 canvas->Print(gfileName);
257 pfileName.ReplaceAll(
".gif",
".png");
259 sprintf(cmd,
"convert %s %s",gfileName.Data(),pfileName.Data());
262 sprintf(cmd,
"rm %s",gfileName.Data());
272 double y = par[0]+pow(
fabs(par[1])/x[0],1)+pow(par[2]/x[0],2);
282 TFile *file0 = TFile::Open(ifName);
283 cout << ifName.Data() << endl;
284 if (file0==NULL) {cout <<
"Error Opening file" << endl;
exit(1);}
286 if (itree==NULL) {cout <<
"Error Opening tree" << endl;
exit(1);}
289 TList*
list = itree->GetUserInfo();
290 int nIFO=list->GetSize();
291 if (nIFO==0) {cout <<
"Error : no ifo present in the tree" << endl;
exit(1);}
292 for (
int n=0;
n<list->GetSize();
n++) {
302 sprintf(cut,
"abs(time[0]-time[%d])<%f && netcc[%d]>%f && rho[%d]>%f",
303 nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
304 if(T_vED>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && neted[0]/ecor<%f",tmp,T_vED);}
305 if(T_pen>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && penalty>%f",tmp,T_pen);}
306 if(type>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && type[1]==%d",tmp,type);}
307 if(T_ifar>0) {
strcpy(tmp,cut);
sprintf(cut,
"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}
321 double supsnr=fsnr+infsnr;
323 double supsnr=fsnr*infsnr;
328 sprintf(SNRnet,
"sqrt(%s)",tmp);
333 sprintf(icut,
"%s && %s>%f && %s<=%f",cut,SNRnet,infsnr,SNRnet,supsnr);
336 itree->Draw(
"erA[0]",icut,
"goff");
338 int size = itree->GetSelectedRows();
340 double* erA = itree->GetV1();
344 TMath::Sort(size,erA,index,
false);
346 entries[nsnr] =
size;
347 median[nsnr] = (pType==
"MEDIAN50") ? erA[index[
int(size*0.5)]] : erA[index[
int(size*0.9)]];
348 snr[nsnr] = (supsnr+infsnr)/2.;
349 if(median[nsnr]<max_median) {
350 cout <<
"Selected entries : " << entries[nsnr] <<
"\t SNRnet : "
351 << snr[nsnr] <<
"\t median : " << median[nsnr] << endl;
detectorParams getDetectorParams()
virtual size_t size() const
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<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t 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
double median(std::vector< double > &vec)
Complex Exp(double phase)
double fabs(const Complex &x)
strcpy(RunLabel, RUN_LABEL)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
virtual void resize(unsigned int)