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