Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DrawWRC.C
Go to the documentation of this file.
1 //
2 // Draw the Normalized Residual Errors & Reaconstructed SNRnet vs Injected SNRnet
3 // Note : this macro is used to generate the PRC report (cwb_report merge_label prc)
4 // Author : Gabriele Vedovato
5 
6 namespace DrawWRC_namespace { // used to avoid conflict with global variables
7 
8 #include <vector>
9 
10 void Plot(TString gtype, TString ptitle, TString xtitle, TString ytitle, int nIFO, float T_win,
11  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar);
12 
13 // --------------------------------------------------------
14 // Global variables
15 // --------------------------------------------------------
16 TCanvas* gCANVAS = NULL; // canvas object
17 TTree* gTRWAVE = NULL; // wave tree
18 
19 
20 void
21 DrawWRC(TString gtype, TString data_label, TString odir, TString merge_label, int nIFO, float T_win,
22  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {
23 
24  if(gtype!="nre" && gtype!="snr" && gtype!="ff" && gtype!="of" && gtype!="ofnre" && gtype!="offf" &&
25  gtype!="mch1" && gtype!="mch2" && gtype!="mch3" && gtype!="mch4" &&
26  gtype!="mch5" && gtype!="mch6" && gtype!="mch7" && gtype!="mch8" &&
27  gtype!="mch9" && gtype!="mch10" && gtype!="mch11" && gtype!="mch12") {
28  cout << "DrawWRC : Error - wrong input gtype (snr/nre/ff/of/ofnre/offf/mch1:12)" << endl;
29  gSystem->Exit(1);
30  }
31 
32  // ------------------------------------------------
33  // create canvas
34  // ------------------------------------------------
35  gCANVAS = new TCanvas("fom", "PRC", 300,40, 600, 600);
36  gCANVAS->Range(-19.4801,-9.25,-17.4775,83.25);
37  if(gtype=="ofnre" || gtype=="offf") gCANVAS->SetRightMargin(0.13);
38  gCANVAS->SetBorderSize(2);
39  gCANVAS->SetFrameFillColor(0);
40  gCANVAS->SetGridx();
41  gCANVAS->SetGridy();
42 
43  // ------------------------------------------------
44  // open wave/mdc trees
45  // ------------------------------------------------
46  // open wave file
47  char sim_file_name[1024];
48  sprintf(sim_file_name,"merge/wave_%s.%s.root",data_label.Data(),merge_label.Data());
49  TFile* fwave = TFile::Open(sim_file_name);
50  gTRWAVE = (TTree*) gROOT->FindObject("waveburst");
51 
52  TString ptitle;
53  TString xtitle;
54  TString ytitle;
55 
56  gCANVAS->SetLogx(true);
57  if(gtype=="snr") gCANVAS->SetLogy(true);
58  if(gtype=="snr") ptitle="Reconstructed SNRnet vs Injected SNRnet";
59  if(gtype=="nre") ptitle="Normalized Residual Energy (NRE)";
60  if(gtype=="ff") ptitle="Fitting Factor (FF)";
61  if(gtype=="of") ptitle="Overlap Factor (OF)";
62  if(gtype=="ofnre") {
63  ptitle="Overlap Factor vs Normalized Residual Energy (OFNRE)";
64  gCANVAS->SetLogx(false);
65  }
66  if(gtype=="offf") {
67  ptitle="Overlap Factor vs Fitting Factor";
68  gCANVAS->SetLogx(false);
69  }
70  if(gtype=="mch1") {
71  ptitle="Rec Chirp Mass vs Inj Chirp Mass (CM)";
72  gCANVAS->SetLogx(false);
73  }
74  if(gtype=="mch2") {
75  ptitle="Rec-Inj Chirp Mass vs Injected SNRnet";
76  gCANVAS->SetLogx(false);
77  }
78  if(gtype=="mch3") {
79  ptitle="Chirp Ellipticity vs Inj Chirp Mass";
80  gCANVAS->SetLogx(false);
81  }
82  if(gtype=="mch4") {
83  ptitle="Chirp Ellipticity vs Injected SNRnet";
84  gCANVAS->SetLogx(false);
85  }
86  if(gtype=="mch5") {
87  ptitle="Chirp Pixel Fraction vs Inj Chirp Mass";
88  gCANVAS->SetLogx(false);
89  }
90  if(gtype=="mch6") {
91  ptitle="Chirp Pixel Fraction vs Injected SNRnet";
92  gCANVAS->SetLogx(false);
93  }
94  if(gtype=="mch7") {
95  ptitle="Chirp Energy Fraction vs Inj Chirp Mass";
96  gCANVAS->SetLogx(false);
97  }
98  if(gtype=="mch8") {
99  ptitle="Chirp Energy Fraction vs Injected SNRnet";
100  gCANVAS->SetLogx(false);
101  }
102  if(gtype=="mch9") {
103  ptitle="Chirp Mass Error vs Inj Chirp Mass";
104  gCANVAS->SetLogx(false);
105  }
106  if(gtype=="mch10") {
107  ptitle="Chirp Mass Error vs Injected SNRnet";
108  gCANVAS->SetLogx(false);
109  }
110  if(gtype=="mch11") {
111  ptitle="PE Rec-Inj Chirp Mass vs Injected SNRnet";
112  gCANVAS->SetLogx(false);
113  }
114  if(gtype=="mch12") {
115  ptitle="PE Chirp Mass Error vs Injected SNRnet";
116  gCANVAS->SetLogx(false);
117  }
118  gStyle->SetOptStat(0);
119  //if(nIFO==2) xtitle = "sqrt(hrss[0]^2+hrss[1]^2)";
120  //if(nIFO==3) xtitle = "sqrt(hrss[0]^2+hrss[1]^2+hrss[2]^2)";
121  xtitle = "Injected SNRnet";
122  if(gtype=="snr") ytitle = "Reconstructed SNRnet";
123  if(gtype=="nre") ytitle = "NRE";
124  if(gtype=="ff") ytitle = "FF";
125  if(gtype=="of") ytitle = "OF";
126  if(gtype=="ofnre") {
127  xtitle = "1-NRE";
128  ytitle = "OF";
129  }
130  if(gtype=="offf") {
131  xtitle = "FF";
132  ytitle = "OF";
133  }
134  if(gtype=="mch1") {
135  xtitle = "Injected Chirp Mass";
136  ytitle = "Reconstructed Chirp Mass";
137  }
138  if(gtype=="mch2") ytitle = "Rec-Inj Chirp Mass";
139  if(gtype=="mch3") {
140  xtitle = "Injected Chirp Mass";
141  ytitle = "Chirp Ellipticity";
142  }
143  if(gtype=="mch4") ytitle = "Chirp Ellipticity";
144  if(gtype=="mch5") {
145  xtitle = "Injected Chirp Mass";
146  ytitle = "Chirp Pixel Fraction";
147  }
148  if(gtype=="mch6") ytitle = "Chirp Pixel Fraction";
149  if(gtype=="mch7") {
150  xtitle = "Injected Chirp Mass";
151  ytitle = "Chirp Energy Fraction";
152  }
153  if(gtype=="mch8") ytitle = "Chirp Energy Fraction";
154  if(gtype=="mch9") {
155  xtitle = "Injected Chirp Mass";
156  ytitle = "Chirp Mass Error";
157  }
158  if(gtype=="mch10") ytitle = "Chirp Mass Error";
159  if(gtype=="mch11") ytitle = "PE Rec-Inj Chirp Mass";
160  if(gtype=="mch12") ytitle = "PE Chirp Mass Error";
161  Plot(gtype, ptitle, xtitle, ytitle, nIFO, T_win, pp_inetcc, T_cor, pp_irho, T_cut, T_vED, T_pen, T_ifar);
162 
163  char ofname[1024];
164  if(gtype=="snr") sprintf(ofname,"%s/osnr_vs_isnr.gif",odir.Data());
165  if(gtype=="nre") sprintf(ofname,"%s/nre_vs_isnr.gif",odir.Data());
166  if(gtype=="ff") sprintf(ofname,"%s/ff_vs_isnr.gif",odir.Data());
167  if(gtype=="of") sprintf(ofname,"%s/of_vs_isnr.gif",odir.Data());
168  if(gtype=="ofnre") sprintf(ofname,"%s/of_vs_nre.gif",odir.Data());
169  if(gtype=="offf") sprintf(ofname,"%s/of_vs_ff.gif",odir.Data());
170  if(gtype=="mch1") sprintf(ofname,"%s/omch_vs_imch.gif",odir.Data());
171  if(gtype=="mch2") sprintf(ofname,"%s/dmch_vs_isnr.gif",odir.Data());
172  if(gtype=="mch3") sprintf(ofname,"%s/elch_vs_imch.gif",odir.Data());
173  if(gtype=="mch4") sprintf(ofname,"%s/elch_vs_isnr.gif",odir.Data());
174  if(gtype=="mch5") sprintf(ofname,"%s/pfch_vs_imch.gif",odir.Data());
175  if(gtype=="mch6") sprintf(ofname,"%s/pfch_vs_isnr.gif",odir.Data());
176  if(gtype=="mch7") sprintf(ofname,"%s/efch_vs_imch.gif",odir.Data());
177  if(gtype=="mch8") sprintf(ofname,"%s/efch_vs_isnr.gif",odir.Data());
178  if(gtype=="mch9") sprintf(ofname,"%s/emch_vs_imch.gif",odir.Data());
179  if(gtype=="mch10") sprintf(ofname,"%s/emch_vs_isnr.gif",odir.Data());
180  if(gtype=="mch11") sprintf(ofname,"%s/opemch_vs_isnr.gif",odir.Data());
181  if(gtype=="mch12") sprintf(ofname,"%s/epemch_vs_isnr.gif",odir.Data());
182  gCANVAS->Print(ofname);
183  char cmd[1024];
184  TString pfname(ofname);
185  pfname.ReplaceAll(".gif",".png");
186  sprintf(cmd,"convert %s %s",ofname,pfname.Data());
187  cout << cmd << endl;
188  gSystem->Exec(cmd);
189  sprintf(cmd,"rm %s",ofname);
190  cout << cmd << endl;
191  gSystem->Exec(cmd);
192  //exit(0);
193 
194 }
195 
196 void Plot(TString gtype, TString ptitle, TString xtitle, TString ytitle, int nIFO, float T_win,
197  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {
198 
199  if(gtype=="mch11" || gtype=="mch12") { // check if pe_mch leaf exist
200  TBranch* branch;
201  bool pe_mch_exists=false;
202  TIter next(gTRWAVE->GetListOfBranches());
203  while ((branch=(TBranch*)next())) {
204  if(TString("pe_mch").CompareTo(branch->GetName())==0) pe_mch_exists=true;
205  }
206  next.Reset();
207  if(!pe_mch_exists) return;
208  }
209 
210  char sel[256]="";
211  char tmp[256]="";
212  char cut[256]="";
213  char title[256]="";
214  char num[256]="";
215  char den[256]="";
216  char ffnum[256]="";
217  char iffden[256]="";
218  char offden[256]="";
219 
220  // DETECTED
221  sprintf(cut,"");
222  // fitting factor
223  strcpy(tmp,"ioSNR[0]");
224  for(int i=1;i<nIFO;i++) {sprintf(ffnum,"%s+ioSNR[%d]",tmp,i);strcpy(tmp,ffnum);}
225  strcpy(tmp,"iSNR[0]");
226  for(int i=1;i<nIFO;i++) {sprintf(iffden,"%s+iSNR[%d]",tmp,i);strcpy(tmp,iffden);}
227  strcpy(tmp,"oSNR[0]");
228  for(int i=1;i<nIFO;i++) {sprintf(offden,"%s+oSNR[%d]",tmp,i);strcpy(tmp,offden);}
229  // residual energy
230  strcpy(tmp,"(iSNR[0]+oSNR[0]-2*ioSNR[0])");
231  for(int i=1;i<nIFO;i++) {sprintf(num,"%s+(iSNR[%d]+oSNR[%d]-2*ioSNR[%d])",tmp,i,i,i);strcpy(tmp,num);}
232  strcpy(tmp,"iSNR[0]");
233  for(int i=1;i<nIFO;i++) {sprintf(den,"%s+iSNR[%d]",tmp,i);strcpy(tmp,den);}
234  if(gtype=="snr") sprintf(sel,"sqrt(likelihood):sqrt(%s)>>htemp",den);
235  if(gtype=="nre") sprintf(sel,"(%s)/(%s):sqrt(%s)>>htemp",num,den,den);
236  if(gtype=="ff") sprintf(sel,"(%s)/sqrt((%s)*(%s)):sqrt(%s)>>htemp",ffnum,iffden,iffden,den);
237  if(gtype=="of") sprintf(sel,"(%s)/sqrt((%s)*(%s)):sqrt(%s)>>htemp",ffnum,iffden,offden,den);
238  if(gtype=="ofnre") sprintf(sel,"(%s)/sqrt((%s)*(%s)):1-(%s)/(%s)>>htemp",ffnum,iffden,offden,num,den);
239  if(gtype=="offf") sprintf(sel,"(%s)/sqrt((%s)*(%s)):(%s)/sqrt((%s)*(%s))>>htemp",ffnum,iffden,offden,ffnum,iffden,iffden);
240  if(gtype=="mch1") sprintf(sel,"chirp[1]:chirp[0]>>htemp");
241  if(gtype=="mch2") sprintf(sel,"chirp[1]-chirp[0]:sqrt(%s)>>htemp",den);
242  if(gtype=="mch3") sprintf(sel,"chirp[3]:chirp[0]>>htemp");
243  if(gtype=="mch4") sprintf(sel,"chirp[3]:sqrt(%s)>>htemp",den);
244  if(gtype=="mch5") sprintf(sel,"chirp[4]:chirp[0]>>htemp");
245  if(gtype=="mch6") sprintf(sel,"chirp[4]:sqrt(%s)>>htemp",den);
246  if(gtype=="mch7") sprintf(sel,"chirp[5]:chirp[0]>>htemp");
247  if(gtype=="mch8") sprintf(sel,"chirp[5]:sqrt(%s)>>htemp",den);
248  if(gtype=="mch9") sprintf(sel,"chirp[2]:chirp[0]>>htemp");
249  if(gtype=="mch10") sprintf(sel,"chirp[2]:sqrt(%s)>>htemp",den);
250  if(gtype=="mch11") sprintf(sel,"pe_mch[0]-chirp[0]:sqrt(%s)>>htemp",den);
251  if(gtype=="mch12") sprintf(sel,"pe_mch[1]:sqrt(%s)>>htemp",den);
252  sprintf(cut,"abs(time[0]-time[%d])<%g && netcc[%d]>%g && rho[%d]>%g",
253  nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
254  if(T_vED>0) {strcpy(tmp,cut);sprintf(cut,"%s && neted[0]/ecor<%f",tmp,T_vED);}
255  if(T_pen>0) {strcpy(tmp,cut);sprintf(cut,"%s && penalty>%f",tmp,T_pen);}
256  if(T_ifar>0) {strcpy(tmp,cut);sprintf(cut,"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}
257  gTRWAVE->SetMarkerColor(kRed);
258  gTRWAVE->SetMarkerStyle(20);
259  gTRWAVE->SetMarkerSize(0.5);
260  gTRWAVE->SetLineColor(kRed);
261  gTRWAVE->SetFillColor(kRed);
262  if(gtype=="ofnre" || gtype=="offf") gTRWAVE->Draw(sel,cut,"colz");
263  else gTRWAVE->Draw(sel,cut);
264  int nwave = gTRWAVE->GetSelectedRows();
265  cout << "nwave : " << nwave << endl;
266  sprintf(title,"%s",cut);
267 
268  if(gtype=="nre") {
269  double* inre = gTRWAVE->GetV1();
270  for(int i=0;i<nwave;i++) if(inre[i]>1) inre[i]=1;
271  for(int i=0;i<nwave;i++) if(inre[i]<0.001) inre[i]=0.001;
272  }
273 
274  TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
275  sprintf(title,"%s : rec : %d",ptitle.Data(),nwave);
276  htemp->SetTitle(title);
277  htemp->GetXaxis()->SetTitle(xtitle);
278  htemp->GetYaxis()->SetTitle(ytitle);
279  htemp->GetXaxis()->SetMoreLogLabels();
280  if(gtype=="snr") htemp->GetYaxis()->SetMoreLogLabels();
281  if(gtype=="nre") htemp->GetYaxis()->SetRangeUser(0.0,1.0);
282  if(gtype=="ff") htemp->GetYaxis()->SetRangeUser(0.0,1.0);
283  if(gtype=="of") htemp->GetYaxis()->SetRangeUser(0.0,1.0);
284  if(gtype=="ofnre") {htemp->GetXaxis()->SetRangeUser(-4.0,1.0);htemp->GetYaxis()->SetRangeUser(-1.0,1.0);}
285  if(gtype=="offf") {htemp->GetXaxis()->SetRangeUser(-2.0,2.0);htemp->GetYaxis()->SetRangeUser(-1.0,1.0);}
286  htemp->GetXaxis()->SetTitleOffset(1.35);
287  htemp->GetYaxis()->SetTitleOffset(1.50);
288  htemp->GetXaxis()->CenterTitle(true);
289  htemp->GetYaxis()->CenterTitle(true);
290  htemp->GetXaxis()->SetLabelFont(42);
291  htemp->GetXaxis()->SetTitleFont(42);
292  htemp->GetYaxis()->SetLabelFont(42);
293  htemp->GetYaxis()->SetTitleFont(42);
294 
295  // draw reference line
296  if(gtype=="snr" || gtype=="mch1") {
297  double lmin=htemp->GetYaxis()->GetXmin();
298  double lmax=htemp->GetYaxis()->GetXmax();
299  if(htemp->GetXaxis()->GetXmin()>lmin) lmin=htemp->GetXaxis()->GetXmin();
300  if(htemp->GetXaxis()->GetXmax()<lmax) lmax=htemp->GetXaxis()->GetXmax();
301  TLine *line = new TLine(lmin,lmin,lmax,lmax);
302  line->SetLineColor(kBlue);
303  line->Draw();
304  } else {
305  }
306 
307  return;
308 }
309 
310 } // end namespace
311 
312 void
313 DrawWRC(TString gtype, TString data_label, TString odir, TString merge_label, int nIFO, float T_win,
314  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {
315 
316  DrawWRC_namespace::DrawWRC(gtype, data_label, odir, merge_label, nIFO, T_win,
317  pp_inetcc, T_cor, pp_irho, T_cut, T_vED, T_pen, T_ifar);
318  return;
319 }
cut
char cut[512]
Definition: check_slag_bkg.C:62
T_ifar
double T_ifar
Definition: cwb_pparameters.C:139
xtitle
char xtitle[1024]
Definition: cwb_epparameters.C:182
T_pen
double T_pen
Definition: cwb_pparameters.C:146
cmd
char cmd[1024]
Definition: cwb_condor_cleanup.C:20
nwave
int nwave
Definition: cbc_plots.C:820
TString
TString("c")
Definition: cwb_report_skymap.C:111
T_cor
double T_cor
Definition: cwb_pparameters.C:140
odir
char odir[1024]
Definition: cwb_mkhtml_header.C:148
pp_inetcc
i pp_inetcc
Definition: cwb_epparameters.C:84
i
i drho i
Definition: cwb_epparameters.C:85
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
data_label
char data_label[512]
Definition: test_config1.C:160
DrawWRC_namespace::gTRWAVE
TTree * gTRWAVE
Definition: DrawWRC.C:17
DrawWRC
void DrawWRC(TString gtype, TString data_label, TString odir, TString merge_label, int nIFO, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar)
Definition: DrawWRC.C:313
tmp
double * tmp
Definition: testWDM_5.C:31
next
TIter next(twave->GetListOfBranches())
merge_label
TString merge_label
Definition: cwb_mkhtml_header.C:49
DrawWRC_namespace::Plot
void Plot(TString gtype, TString ptitle, TString xtitle, TString ytitle, int nIFO, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar)
Definition: DrawWRC.C:196
num
int num
Definition: check_slag_bkg.C:71
sim_file_name
char sim_file_name[1024]
Definition: cwb_pparameters.C:166
DrawWRC_namespace::gCANVAS
TCanvas * gCANVAS
Definition: DrawWRC.C:16
title
char title[256]
Definition: SSeriesExample.C:1
line
TLine * line
Definition: compare_bkg.C:482
T_win
double T_win
Definition: cwb_pparameters.C:148
sel
TString sel("slag[1]:rho[1]")
branch
TBranch * branch
Definition: cwb_fix_live_slag_missed.C:42
strcpy
strcpy(RunLabel, RUN_LABEL)
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
T_cut
double T_cut
Definition: cwb_pparameters.C:141
DrawWRC_namespace
Definition: DrawWRC.C:6
T_vED
double T_vED
Definition: cwb_pparameters.C:145
DrawWRC_namespace::DrawWRC
void DrawWRC(TString gtype, TString data_label, TString odir, TString merge_label, int nIFO, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar)
Definition: DrawWRC.C:21
pp_irho
i drho pp_irho
Definition: cwb_epparameters.C:85