Logo Coherent WaveBurst  
Reference Guide
Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cbc_plots.C
Go to the documentation of this file.
1 
2 #define RHO_MIN 6.0
3 #define RHO_BIN 0.1
4 #define RHO_NBINS 1000
5 #define CONTOURS 20
6 
7 
8 #define LIV_FILE_NAME "live.txt"
9 
10 
11 
12 {
13 
14 
15  //###############################################################################################################################
16 // Palette
17 //###############################################################################################################################
18 
19 
20  #define NRGBs 6
21  #define NCont 99
22  gStyle->SetNumberContours(NCont);
23  double stops[NRGBs] = { 0.10, 0.25, 0.45, 0.60, 0.75, 1.00 };
24  double red[NRGBs] = { 0.00, 0.00, 0.00, 0.97, 0.97, 0.10 };
25  double green[NRGBs] = { 0.97, 0.30, 0.40, 0.97, 0.00, 0.00 };
26  double blue[NRGBs] = { 0.97, 0.97, 0.00, 0.00, 0.00, 0.00 };
27  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
28 
29  // -----------------------------------------------------------
30  // CBC
31  // -----------------------------------------------------------
32 
33  int NBINS_mass = (int)((MAX_MASS-MIN_MASS)/MASS_BIN);
34 
35 
36 #ifdef MAX_MASS2
37  float max_mass2 = MAX_MASS2;
38 #else
39  float max_mass2 = MAX_MASS;
40 #endif
41 
42 #ifdef MIN_MASS1
43  float min_mass1 = MIN_MASS1;
44 #else
45  float min_mass1 = MIN_MASS;
46 #endif
47 
48 #ifdef MAX_MASS1
49  float max_mass1 = MAX_MASS1;
50 #else
51  float max_mass1 = MAX_MASS;
52 #endif
53 
54 #ifdef MIN_MASS2
55  float min_mass2 = MIN_MASS2;
56 #else
57  float min_mass2 = MIN_MASS;
58 #endif
59 
60 
61 
62  int NBINS_mass1 = (int)((max_mass1-min_mass1)/MASS_BIN);
63  int NBINS_mass2 = (int)((max_mass2-min_mass2)/MASS_BIN);
64 
65 
66 
67  cout << "cbc_plots.C starts..."<<endl;
68 // cout << "Mass1 : ["<<MIN_MASS<<","<<MAX_MASS<<"] with "<<NBINS_mass<<" bins"<<endl;
69  cout << "Mass1 : ["<<min_mass1<<","<<max_mass1<<"] with "<<NBINS_mass1<<" bins"<<endl;
70  cout << "Mass2 : ["<<min_mass2<<","<<max_mass2<<"] with "<<NBINS_mass2<<" bins"<<endl;
71 
73  CWB::CBCTool cbcTool;
74 
75  TB.checkFile(gSystem->Getenv("CWB_ROOTLOGON_FILE"));
76  TB.checkFile(gSystem->Getenv("CWB_PARAMETERS_FILE"));
77  TB.checkFile(gSystem->Getenv("CWB_UPARAMETERS_FILE"));
78  TB.checkFile(gSystem->Getenv("CWB_PPARAMETERS_FILE"));
79  TB.checkFile(gSystem->Getenv("CWB_UPPARAMETERS_FILE"));
80  TB.checkFile(gSystem->Getenv("CWB_EPPARAMETERS_FILE"));
81 
82  TB.mkDir(netdir,true);
83 
84  gStyle->SetTitleFillColor(kWhite);
85  //FgStyle->SetLineColor(kWhite);
86  gStyle->SetNumberContours(256);
87  gStyle->SetCanvasColor(kWhite);
88  gStyle->SetStatBorderSize(1);
89  gStyle->SetOptStat(kFALSE);
90 
91  // remove the red box around canvas
92  gStyle->SetFrameBorderMode(0);
93  gROOT->ForceStyle();
94 
95  char networkname[256];
96  if(strlen(ifo[0])>0) strcpy(networkname,ifo[0]); else strcpy(networkname,detParms[0].name);
97  for(int i=1;i<nIFO;i++) {
98  if(strlen(ifo[i])>0) sprintf(networkname,"%s-%s",networkname,ifo[i]); // built in detector
99  else sprintf(networkname,"%s-%s",networkname,detParms[i].name); // user define detector
100  }
101 
102  // int gIFACTOR=1;
103  // double FACTORS[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
104 
105  network* net=NULL;
107  strcpy(cfg->tmp_dir,"tmp");
108  // load MDC setup (define MDC set)
109 
110  gROOT->Macro(configPlugin.GetTitle());
111  TString Insp = MDC.GetInspiral();
112 
116 
117 
118  if(MDC.GetInspiralOption("--waveform")!=""){waveform = MDC.GetInspiralOption("--waveform");}
119  if(MDC.GetInspiralOption("--min-mtotal")!=""){minMtot = (float)MDC.GetInspiralOption("--min-mtotal").Atof();bminMtot = 1;}
120  if(MDC.GetInspiralOption("--max-mtotal")!=""){maxMtot = (float)MDC.GetInspiralOption("--max-mtotal").Atof(); bmaxMtot = 1;}
121  if(MDC.GetInspiralOption("--min-distance")!=""){minDistance = (float)MDC.GetInspiralOption("--min-distance").Atof(); bminDistance = 1;}
122  if(MDC.GetInspiralOption("--max-distance")!=""){maxDistance = (float)MDC.GetInspiralOption("--max-distance").Atof(); bmaxDistance = 1;}
123  if(MDC.GetInspiralOption("--min-mratio")!=""){minRatio = (float)MDC.GetInspiralOption("--min-mratio").Atof(); bminRatio = 1;}
124  if(MDC.GetInspiralOption("--max-mratio")!=""){maxRatio = (float)MDC.GetInspiralOption("--max-mratio").Atof(); bmaxRatio = 1;}
125  if(MDC.GetInspiralOption("--d-distr")!=""){
126  if(MDC.GetInspiralOption("--d-distr").CompareTo("uniform")==0){DDistrUniform =1;cout << "Uniform distribution in linear distance"<<endl;}
127  else if(MDC.GetInspiralOption("--d-distr").CompareTo("volume")==0){DDistrVolume =1;cout << "Uniform distribution in volume"<<endl;}
128  else {cout << "No defined distance distribution?????! Or different from uniform and volume???"<<endl;exit(1);}
129  }
130 
131 
132 
133  #ifdef FIXMAXDISTANCE
134  bmaxDistance=1;
135  maxDistance=FIXMAXDISTANCE;
136  #endif
137 
138  #ifdef FIXMINRATIO
139  bminRatio=1;
140  minRatio=FIXMINRATIO;
141  float MINRATIO = minRatio;
142  #else
143  float MINRATIO = 0.0;
144  #endif
145 
146 #ifdef FIXMAXRATIO
147  bmaxRatio=1;
148  maxRatio=FIXMAXRATIO;
149  float MAXRATIO = maxRatio;
150  #else
151  float MAXRATIO = 1.;
152  #endif
153 
154  #ifdef FIXMINTOT
155  bminMtot=1;
156  minMtot=FIXMINTOT;
157  #endif
158 
159 #ifdef FIXMAXTOT
160  bmaxMtot=1;
161  maxMtot=FIXMAXTOT;
162  #endif
163 
164 
165  if(bminMtot){float MINMtot = 0.99*minMtot;}
166  else {float MINMtot = 0.0;}
167 
168  if(bmaxMtot){
169  float MAXMtot = 1.01*maxMtot;
170  int NBINS_MTOT = TMath::FloorNint((maxMtot-minMtot)/MASS_BIN/2.);
171  cout<<"NBINS_MTOT: "<<NBINS_MTOT<<endl;
172 
173  }
174  else {cout <<"Undefined maximal total mass!! Define float MAXMtot"<<endl;exit(1);}
175 
176  if(bminDistance) {float MINDISTANCE = 0.9*minDistance;}
177  else {cout <<"Undefined minimal distance!! Defined float MINDISTANCE"<<endl;exit(1);}
178 
179  if(bmaxDistance) {float MAXDISTANCE = 1.1*maxDistance;}
180  else {cout <<"Undefined maximal distance!! Define float MAXDISTANCE"<<endl;exit(1);}
181 
182 
183 //if((bminRatio)&&(bmaxRatio)){float MINCHIRP = 0.7*minMtot*pow(minRatio,3./5)/pow(1+minRatio,6./5); float MAXCHIRP = 0.6*maxMtot/pow(2.,6./5); }
184 // else {cout <<"Undefined minRatio.. "<<endl;float MINCHIRP = 0.7*minMtot*pow(MINRATIO,3./5)/pow(1+MINRATIO,6./5); float MAXCHIRP = 0.6*maxMtot/pow(2.,6./5);}
185 
186  float MINCHIRP = 0.9*pow(min_mass1*min_mass2,3./5.)/pow(maxMtot,1./5.);
187  float MAXCHIRP = 1.1*pow(max_mass1*max_mass2,3./5.)/pow(minMtot,1./5.);
188 
189  //If not vetoes are defined, then the char string veto_not_vetoed is forced to be equal to ch2
190  if (strlen(veto_not_vetoed) == 0){sprintf(veto_not_vetoed,"%s",ch2);}
191 
192 //###############################################################################################################################
193 // Definitions of Canvas and histograms
194 //###############################################################################################################################
195 
196 
197  TCanvas *c1 = new TCanvas("c1","c1",3,47,1000,802);
198  c1->Range(-1.216392,-477.6306,508.8988,2814.609);
199  c1->SetFillColor(0);
200  c1->SetBorderMode(0);
201  c1->SetBorderSize(2);
202  c1->SetGridx();
203  c1->SetGridy();
204  c1->SetRightMargin(0.1154618);
205  c1->SetTopMargin(0.07642487);
206  c1->SetBottomMargin(0.1450777);
207 
208 
209  TH2F *inj_events = new TH2F("inj_events","D_Minj",NBINS_mass,MIN_MASS,MAX_MASS,NBINS_mass2,min_mass2,max_mass2);
210  inj_events->GetXaxis()->SetRangeUser(MIN_plot_mass1,MAX_plot_mass1);
211  inj_events->GetYaxis()->SetRangeUser(MIN_plot_mass2,MAX_plot_mass2);
212  inj_events->GetXaxis()->SetTitle("Mass 1 (M_{#odot})");
213  inj_events->GetYaxis()->SetTitle("Mass 2 (M_{#odot})");
214  inj_events->GetXaxis()->SetTitleOffset(1.3);
215  inj_events->GetYaxis()->SetTitleOffset(1.25);
216  inj_events->GetXaxis()->CenterTitle(kTRUE);
217  inj_events->GetYaxis()->CenterTitle(kTRUE);
218  inj_events->GetXaxis()->SetNdivisions(410);
219  inj_events->GetYaxis()->SetNdivisions(410);
220  inj_events->GetXaxis()->SetTickLength(0.01);
221  inj_events->GetYaxis()->SetTickLength(0.01);
222  inj_events->GetZaxis()->SetTickLength(0.01);
223  inj_events->GetXaxis()->SetTitleFont(42);
224  inj_events->GetXaxis()->SetLabelFont(42);
225  inj_events->GetYaxis()->SetTitleFont(42);
226  inj_events->GetYaxis()->SetLabelFont(42);
227  inj_events->GetZaxis()->SetLabelFont(42);
228  inj_events->GetZaxis()->SetLabelSize(0.03);
229  //inj_events->GetXaxis()->SetNdivisions(10,kFALSE);
230 //inj_events->GetYaxis()->SetNdivisions(10,kFALSE);
231 
232  inj_events->SetTitle("");
233 
234 
235  inj_events->SetContour(NCont);
236 
237 
238  TH2F *rec_events = new TH2F("rec_events","D_Mrec",NBINS_mass,MIN_MASS,MAX_MASS,NBINS_mass2,min_mass2,max_mass2);
239  rec_events->GetXaxis()->SetRangeUser(MIN_plot_mass1,MAX_plot_mass1);
240  rec_events->GetYaxis()->SetRangeUser(MIN_plot_mass2,MAX_plot_mass2);
241  rec_events->GetXaxis()->SetTitle("Mass 1 (M_{#odot})");
242  rec_events->GetYaxis()->SetTitle("Mass 2 (M_{#odot})");
243  rec_events->GetXaxis()->SetTitleOffset(1.3);
244  rec_events->GetYaxis()->SetTitleOffset(1.25);
245  rec_events->GetXaxis()->CenterTitle(kTRUE);
246  rec_events->GetYaxis()->CenterTitle(kTRUE);
247  rec_events->GetXaxis()->SetNdivisions(410);
248  rec_events->GetYaxis()->SetNdivisions(410);
249  rec_events->GetXaxis()->SetTickLength(0.01);
250  rec_events->GetYaxis()->SetTickLength(0.01);
251  rec_events->GetZaxis()->SetTickLength(0.01);
252  rec_events->GetXaxis()->SetTitleFont(42);
253  rec_events->GetXaxis()->SetLabelFont(42);
254  rec_events->GetYaxis()->SetTitleFont(42);
255  rec_events->GetYaxis()->SetLabelFont(42);
256  rec_events->GetZaxis()->SetLabelFont(42);
257  rec_events->GetZaxis()->SetLabelSize(0.03);
258  rec_events->SetTitle("");
259 
260 
261  rec_events->SetContour(NCont);
263  for(int i=0;i<nfactor;i++){
264  factor_events_inj[i]= new TH2F("factor_events_inj","",NBINS_mass,MIN_MASS,MAX_MASS,NBINS_mass2,min_mass2,max_mass2);
265  }
266  TH2F *factor_events_rec = new TH2F("factor_events_rec","",NBINS_mass,MIN_MASS,MAX_MASS,NBINS_mass2,min_mass2,max_mass2);
267 
268 
269  TH2F *D_Mtot_inj = new TH2F("Distance vs Mtot inj.","",1000,MINMtot,MAXMtot,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000./FACTORS[0]);
270 
271  D_Mtot_inj->GetXaxis()->SetTitle("Total mass (M_{#odot})");
272  D_Mtot_inj->GetYaxis()->SetTitle("Distance (Mpc)");
273  D_Mtot_inj->GetXaxis()->SetTitleOffset(1.3);
274  D_Mtot_inj->GetYaxis()->SetTitleOffset(1.3);
275  D_Mtot_inj->GetXaxis()->SetTickLength(0.01);
276  D_Mtot_inj->GetYaxis()->SetTickLength(0.01);
277  D_Mtot_inj->GetXaxis()->CenterTitle(kTRUE);
278  D_Mtot_inj->GetYaxis()->CenterTitle(kTRUE);
279  D_Mtot_inj->GetXaxis()->SetTitleFont(42);
280  D_Mtot_inj->GetXaxis()->SetLabelFont(42);
281  D_Mtot_inj->GetYaxis()->SetTitleFont(42);
282  D_Mtot_inj->GetYaxis()->SetLabelFont(42);
283  D_Mtot_inj->SetMarkerStyle(20);
284  D_Mtot_inj->SetMarkerSize(0.5);
285  D_Mtot_inj->SetMarkerColor(2);
286  D_Mtot_inj->SetTitle("");
287 // D_Mtot_inj->Draw("ap");
288 D_Mtot_inj->SetName("D_Mtotinj");
289 
290  TH2F *D_Mtot_rec = new TH2F("Distance vs Mtot rec.","",1000,MINMtot,MAXMtot,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000/FACTORS[0]);
291 // D_Mtot_rec->SetName("D_rec");
292  D_Mtot_rec->SetMarkerStyle(20);
293  D_Mtot_rec->SetMarkerSize(0.5);
294  D_Mtot_rec->SetMarkerColor(4);
295 
296  TH2F *D_Mchirp_inj = new TH2F("Distance vs MChirp inj.","",1000,MINCHIRP,MAXCHIRP,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000./FACTORS[0]);
297 
298  D_Mchirp_inj->GetXaxis()->SetTitle("Chirp mass (M_{#odot})");
299  D_Mchirp_inj->GetYaxis()->SetTitle("Distance (Mpc)");
300  D_Mchirp_inj->GetXaxis()->SetTitleOffset(1.3);
301  D_Mchirp_inj->GetYaxis()->SetTitleOffset(1.3);
302  D_Mchirp_inj->GetXaxis()->SetTickLength(0.01);
303  D_Mchirp_inj->GetYaxis()->SetTickLength(0.01);
304  D_Mchirp_inj->GetXaxis()->CenterTitle(kTRUE);
305  D_Mchirp_inj->GetYaxis()->CenterTitle(kTRUE);
306  D_Mchirp_inj->GetXaxis()->SetTitleFont(42);
307  D_Mchirp_inj->GetXaxis()->SetLabelFont(42);
308  D_Mchirp_inj->GetYaxis()->SetTitleFont(42);
309  D_Mchirp_inj->GetYaxis()->SetLabelFont(42);
310  D_Mchirp_inj->SetMarkerStyle(20);
311  D_Mchirp_inj->SetMarkerSize(0.5);
312  D_Mchirp_inj->SetMarkerColor(2);
313  D_Mchirp_inj->SetTitle("");
314 // D_Mchirp_inj->Draw("ap");
315 D_Mchirp_inj->SetName("D_Chirp_inj");
316 
317  TH2F *D_Mchirp_rec = new TH2F("Distance vs MChirp rec.","",1000,MINCHIRP,MAXCHIRP,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000/FACTORS[0]);
318  //D_Mchirp_rec->SetName("D_rec");
319  D_Mchirp_rec->SetMarkerStyle(20);
320  D_Mchirp_rec->SetMarkerSize(0.5);
321  D_Mchirp_rec->SetMarkerColor(4);
322 
323 #ifdef MINCHI
324 
325  int NBINS_SPIN = (int)((MAXCHI-MINCHI)/CHI_BIN); cout<<"NBINS_SPIN: "<<NBINS_SPIN<<endl;
326  TH2F* inj_events_spin_mtot = new TH2F("inj_spin_Mtot","",NBINS_SPIN,MINCHI,MAXCHI,NBINS_MTOT,minMtot,maxMtot);
327  TH2F* rec_events_spin_mtot = new TH2F("rec_spin_Mtot","",NBINS_SPIN,MINCHI,MAXCHI,NBINS_MTOT,minMtot,maxMtot);
329  for(int i=0;i<nfactor;i++){
330  factor_events_spin_mtot_inj[i]= new TH2F("factor_events_spin_mtot_inj","",NBINS_SPIN,MINCHI,MAXCHI,NBINS_MTOT,minMtot,maxMtot);
331  }
332 
333 TH2F *D_Chi_inj = new TH2F("Distance vs Chi inj.","",1000,MINCHI,MAXCHI,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000./FACTORS[0]);
334 
335  D_Chi_inj->GetXaxis()->SetTitle("#chi_{z}");
336  D_Chi_inj->GetYaxis()->SetTitle("Distance (Mpc)");
337  D_Chi_inj->GetXaxis()->SetTitleOffset(1.3);
338  D_Chi_inj->GetYaxis()->SetTitleOffset(1.3);
339  D_Chi_inj->GetXaxis()->SetTickLength(0.01);
340  D_Chi_inj->GetYaxis()->SetTickLength(0.01);
341  D_Chi_inj->GetXaxis()->CenterTitle(kTRUE);
342  D_Chi_inj->GetYaxis()->CenterTitle(kTRUE);
343  D_Chi_inj->GetXaxis()->SetTitleFont(42);
344  D_Chi_inj->GetXaxis()->SetLabelFont(42);
345  D_Chi_inj->GetYaxis()->SetTitleFont(42);
346  D_Chi_inj->GetYaxis()->SetLabelFont(42);
347  D_Chi_inj->SetMarkerStyle(20);
348  D_Chi_inj->SetMarkerSize(0.5);
349  D_Chi_inj->SetMarkerColor(2);
350  D_Chi_inj->SetTitle("");
351 // D_Mchirp_inj->Draw("ap");
352 D_Chi_inj->SetName("D_Chi_inj");
353 
354  TH2F *D_Chi_injx = (TH2F*)D_Chi_inj->Clone("D_Chi_injx");
355  D_Chi_injx->GetXaxis()->SetTitle("#chi_{x}");
356  D_Chi_injx->SetName("D_Chi_injx");
357  TH2F *D_Chi_injy = (TH2F*)D_Chi_inj->Clone("D_Chi_injy");
358  D_Chi_injy->GetXaxis()->SetTitle("#chi_{y}");
359  D_Chi_injy->SetName("D_Chi_injy");
360 
361  TH2F *D_Chi_rec = new TH2F("Distance vs Chi rec.","",1000,MINCHI,MAXCHI,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000/FACTORS[0]);
362  D_Chi_rec->SetMarkerStyle(20);
363  D_Chi_rec->SetMarkerSize(0.5);
364  D_Chi_rec->SetMarkerColor(4);
365 
366  TH2F *D_Chi_recx = (TH2F*)D_Chi_rec->Clone("D_Chi_recx");
367  TH2F *D_Chi_recy = (TH2F*)D_Chi_rec->Clone("D_Chi_recy");
368 
369 #endif
370 
371 
372  TH2F *D_q_inj = new TH2F("Distance vs q inj.","",1000,MINRATIO,MAXRATIO,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000./FACTORS[0]);
373 
374  D_q_inj->GetXaxis()->SetTitle("Mass ratio (M_{#odot})");
375  D_q_inj->GetYaxis()->SetTitle("Distance (Mpc)");
376  D_q_inj->GetXaxis()->SetTitleOffset(1.3);
377  D_q_inj->GetYaxis()->SetTitleOffset(1.3);
378  D_q_inj->GetXaxis()->SetTickLength(0.01);
379  D_q_inj->GetYaxis()->SetTickLength(0.01);
380  D_q_inj->GetXaxis()->CenterTitle(kTRUE);
381  D_q_inj->GetYaxis()->CenterTitle(kTRUE);
382  D_q_inj->GetXaxis()->SetTitleFont(42);
383  D_q_inj->GetXaxis()->SetLabelFont(42);
384  D_q_inj->GetYaxis()->SetTitleFont(42);
385  D_q_inj->GetYaxis()->SetLabelFont(42);
386  D_q_inj->SetMarkerStyle(20);
387  D_q_inj->SetMarkerSize(0.5);
388  D_q_inj->SetMarkerColor(2);
389  D_q_inj->SetTitle("");
390 // D_q_inj->Draw("ap");
391 D_q_inj->SetName("D_q_inj");
392 
393  TH2F *D_q_rec = new TH2F("Distance vs q rec.","",1000,MINRATIO,MAXRATIO,5000,MINDISTANCE/1000./FACTORS[nfactor-1],MAXDISTANCE/1000/FACTORS[0]);
394 // D_q_rec->SetName("D_rec");
395  D_q_rec->SetMarkerStyle(20);
396  D_q_rec->SetMarkerSize(0.5);
397  D_q_rec->SetMarkerColor(4);
398  TH1F *Dt = new TH1F("Dt","",1000,-0.5,0.5);
399  Dt->SetMarkerStyle(20);
400  Dt->SetMarkerSize(0.5);
401  Dt->SetMarkerColor(4);
402 
403 
404  TH2F* rhocc = new TH2F("rhocc","",100,0.,1.,100,pp_rho_min,pp_rho_max);
405  rhocc->SetTitle("0 < cc < 1");
406  rhocc->SetTitleOffset(1.3,"Y");
407  rhocc->GetXaxis()->SetTitle("network correlation");
408  rhocc->GetYaxis()->SetTitle("#rho");
409  rhocc->SetStats(kFALSE);
410  rhocc->SetMarkerStyle(20);
411  rhocc->SetMarkerSize(0.5);
412  rhocc->SetMarkerColor(1);
413 
414  TH2F* rho_pf = new TH2F("rho_pf","",100,-1.,2.,100,pp_rho_min,pp_rho_max);
415  rho_pf->SetTitle("chi2");
416  rho_pf->GetXaxis()->SetTitle("log10(chi2)");
417  rho_pf->SetTitleOffset(1.3,"Y");
418  rho_pf->GetYaxis()->SetTitle("#rho");
419  rho_pf->SetMarkerStyle(20);
420  rho_pf->SetMarkerColor(1);
421  rho_pf->SetMarkerSize(0.5);
422  rho_pf->SetStats(kFALSE);
423 
424 
425  //###############################################################################################################################
426 // Loop over the injected and recovered events
427 //###############################################################################################################################
428 
429  TChain sim("waveburst");
430  TChain mdc("mdc");
431  sim.Add(sim_file_name);
432  mdc.Add(mdc_file_name);
433 
434  int l,m,mt,cz;
435  double time[6];
437  float mass[2];
438  float spin[6];
439  float chi[3];
440  mdc.SetBranchAddress("time",time);
441  mdc.SetBranchAddress("mass",mass);
442  mdc.SetBranchAddress("factor",&factor);
443  mdc.SetBranchAddress("distance",&distance);
444  mdc.SetBranchAddress("mchirp",&mchirp);
445  mdc.SetBranchAddress("spin",spin);
446 
447 bool write_ascii = false;
448 #ifdef WRITE_ASCII
449  write_ascii = true;
450  char fname3[1024];
451  sprintf(fname3,"%s/injected_signals.txt",netdir);
452  FILE* finj = fopen(fname3,"w");
453  fprintf(finj,"#GPS@L1 mass1 mass2 distance spinx1 spiny1 spinz1 spinx2 spiny2 spinz2 \n");
454 #endif
455  for(int g=0;g<(int)mdc.GetEntries();g++){
456  mdc.GetEntry(g);
457 
458 #ifdef MINCHI
459  for(int i=0;i<3;i++){chi[i] = (spin[i]*mass[1]+spin[i+3]*mass[0])/(mass[1]+mass[0]);}
460  D_Chi_inj->Fill(chi[2],distance);
461  D_Chi_injx->Fill(chi[0],distance);
462  D_Chi_injy->Fill(chi[1],distance);
463  inj_events_spin_mtot->Fill(chi[2],mass[1]+mass[0]);
464 #endif
465 
466  for(int i = 0; i<nfactor; i++){
467  if(fabs(factor-FACTORS[i])<TOLERANCE){
468  factor_events_inj[i]->Fill(mass[1],mass[0]);
469 #ifdef MINCHI
470  factor_events_spin_mtot_inj[i]->Fill(chi[2],mass[1]+mass[0]);
471 #endif
472  break;
473  }
474  }
475 
476  inj_events->Fill(mass[1],mass[0]);
477  D_Mtot_inj->Fill(mass[1]+mass[0],distance);
478  D_Mchirp_inj->Fill(mchirp,distance);
479  D_q_inj->Fill(mass[0]/mass[1],distance);
480 #ifdef WRITE_ASCII
481  fprintf(finj,"%10.3f %3.2f %3.2f %3.2f %3.2f %3.2f %3.2f %3.2f %3.2f %3.2f\n",time[0],mass[0],mass[1],distance,spin[0],spin[1],spin[2],spin[3],spin[4],spin[5]);
482 #endif
483 
484  }
485  fclose(finj);
486  cout << "Injected signals: " << mdc.GetEntries() << endl;
487  char cut[512];
488 
489  sprintf(cut,"rho[%d]>%f && %s",pp_irho,T_out,ch2);
490 
491 
492  float ecor,m1,m2,netcc[3],neted,penalty;
493  float rho[2];
494 
495  float chirp[6];
496  float range[2];
497 float iSNR[3],sSNR[3];
498  sim.SetBranchAddress("mass",mass);
499  sim.SetBranchAddress("factor",&factor);
500 #ifdef OLD_TREE
501  sim.SetBranchAddress("distance",&distance);
502  sim.SetBranchAddress("mchirp",&mchirp);
503 #else
504  sim.SetBranchAddress("range",range);
505  sim.SetBranchAddress("chirp",chirp);
506 #endif
507  sim.SetBranchAddress("rho",rho);
508  sim.SetBranchAddress("netcc",netcc);
509  sim.SetBranchAddress("neted",&neted);
510  sim.SetBranchAddress("ecor",&ecor);
511  sim.SetBranchAddress("penalty",&penalty);
512  sim.SetBranchAddress("time",time);
513  sim.SetBranchAddress("iSNR",iSNR);
514  sim.SetBranchAddress("sSNR",sSNR);
515  sim.SetBranchAddress("spin",spin);
516 
517 #ifdef HVETO_VETO
518 cout<<"Adding hveto flags.."<<endl;
519  UChar_t veto_hveto_L1, veto_hveto_H1,veto_hveto_V1;
520  sim.SetBranchAddress("veto_hveto_L1",&veto_hveto_L1);
521  sim.SetBranchAddress("veto_hveto_H1",&veto_hveto_H1);
522  sim.SetBranchAddress("veto_hveto_V1",&veto_hveto_V1);
523 #endif
524 
525 #ifdef CAT3_VETO
526 cout<<"Adding cat3 flags.."<<endl;
527 UChar_t veto_cat3_L1, veto_cat3_H1,veto_cat3_V1;
528  sim.SetBranchAddress("veto_cat3_L1",&veto_cat3_L1);
529  sim.SetBranchAddress("veto_cat3_H1",&veto_cat3_H1);
530  sim.SetBranchAddress("veto_cat3_V1",&veto_cat3_V1);
531 
532 #endif
533 
537  int cnt = 0;
538  int cnt2 = 0;
539 
540 #ifdef MINCHI
541  float spin_mtot_volume[NBINS_MTOT+1][NBINS_SPIN+1], error_spin_mtot_volume[NBINS_MTOT+1][NBINS_SPIN+1];
542  float spin_mtot_radius[NBINS_MTOT+1][NBINS_SPIN+1], error_spin_mtot_radius[NBINS_MTOT+1][NBINS_SPIN+1];
543  for(int i = 0; i <NBINS_MTOT+1; i++){
544  for(int j = 0; j <NBINS_SPIN+1; j++){
545  spin_mtot_volume[i][j] = 0.;
546  error_spin_mtot_volume[i][j] = 0.;
547  // volume_first_shell[i][j] = 0.;
548  //error_volume_first_shell[i][j] = 0.;
549  spin_mtot_radius[i][j] = 0.;
550  error_spin_mtot_radius[i][j] = 0.;
551  }
552  }
553 
554 #endif
555 
556 #ifdef WRITE_ASCII
557  char fname2[1024];
558  sprintf(fname2,"%s/recovered_signals.txt",netdir);
559  FILE* fev = fopen(fname2,"w");
560  fprintf(fev,"#GPS@L1 rho iSNR sSNR \n");
561 #endif
562 
563  for(int i = 0; i <NBINS_mass; i++){
564  for(int j = 0; j <NBINS_mass2; j++){
565  volume[i][j] = 0.;
566  error_volume[i][j] = 0.;
567  volume_first_shell[i][j] = 0.;
568  error_volume_first_shell[i][j] = 0.;
569  radius[i][j] = 0.;
570  error_radius[i][j] = 0.;
571  }
572  }
573  double Vrho[RHO_NBINS],eVrho[RHO_NBINS],Rrho[RHO_NBINS],eRrho[RHO_NBINS],Trho[RHO_NBINS];
574  for(int i=0; i<RHO_NBINS; i++){Vrho[i]=0.;eVrho[i]=0.;Rrho[i]=0.;eRrho[i]=0.;Trho[i]=RHO_MIN+i*RHO_BIN;}
575  double dV,dV1,dV_spin_mtot;
576  int nT;
577  for(int g=0;g<(int)sim.GetEntries();g++){
578  sim.GetEntry(g);
579  if(rho[pp_irho]<=RHO_MIN){continue;}
580  if(netcc[0]<=T_cor){continue;}
581  if((time[0]-time[nIFO])<-T_win || (time[0]-time[nIFO])>2*T_win) {continue;} //NOT checking for detector 1 and 2: very small bias...
582  if (T_vED>0) {if(neted/ecor>=T_vED){continue;}}
583  if (T_pen>0) {if(penalty<=T_vED){continue;}}
584  //if(factor<11.){continue;}
585  if (++cnt%1000==0) {cout << cnt << " - ";}
586  Dt->Fill(time[0]-time[nIFO]);
587 int countv=0;
588 
589 
590 #ifdef OLD_TREE
591 
592  range[1]=distance;
593  chirp[0]=mchirp;
594 
595 #endif
596  if(range[1]==0.0){continue;}
597  D_Mtot_rec->Fill(mass[1]+mass[0],range[1]);
598  D_Mchirp_rec->Fill(chirp[0],range[1]);
599  D_q_rec->Fill(mass[0]/mass[1],range[1]);
600  //D_rec->Fill(range[1]);
601  rhocc->Fill(netcc[0],rho[pp_irho]);
602 
603  float chi2 = penalty>0 ? log10(penalty) : 0;
604  rho_pf->Fill(chi2,rho[pp_irho]);
605 
606  m1=mass[1];
607  m2=mass[0];
608  l = TMath::FloorNint((m2-min_mass2)/MASS_BIN);
609  if(l>NBINS_mass2+1){cout<<"l="<<l<<" is larger than NBINS_mass2="<<NBINS_mass2<< " m2="<<m2<<" min_mass2="<<min_mass2<<endl;exit(1);}
610  m = TMath::FloorNint((m1-MIN_MASS)/MASS_BIN);
611  if(m>NBINS_mass+1){cout<<"m="<<m<<" is larger than NBINS_mass="<<NBINS_mass<< " m1="<<m1<<" MIN_MASS="<<MIN_MASS<<endl;exit(1);}
612  rec_events->Fill(mass[1],mass[0]);
613 #ifdef MINCHI
614  for(int i=0;i<3;i++){chi[i] = (spin[i]*mass[1]+spin[i+3]*mass[0])/(mass[1]+mass[0]);}
615  D_Chi_rec->Fill(chi[2],range[1]);
616  D_Chi_recx->Fill(chi[0],range[1]);
617  D_Chi_recy->Fill(chi[1],range[1]);
618  rec_events_spin_mtot->Fill(chi[2],mass[1]+mass[0]);
619  mt = TMath::FloorNint((m1+m2-minMtot)/MASS_BIN/2.);
620  if(mt>NBINS_MTOT){cout<<"mt="<<mt<<" is larger than NBINS_MTOT="<<NBINS_MTOT<< " Mtot="<<m1+m2<<" minMtot="<<minMtot<<endl;exit(1);}
621  cz = TMath::FloorNint((chi[2]-MINCHI)/CHI_BIN);
622  if(cz>NBINS_SPIN){cout<<"cz="<<cz<<" is larger than NBINS_SPIN="<<NBINS_SPIN<< " chi[2]="<<chi[2]<<" MINCHI="<<MINCHI<<endl;exit(1);}
623 
624 #endif
625 
626 
627  // if(factor==FACTORS[nfactor]){factor_events_rec->Fill(mass[1],mass[0]);dV = 1./(pow(factor,3)*factor_events_inj[i]->GetBinContent(m+1,l+1));}break;} ///BEWARE! Sorted factors!
628  for(int i = 0; i<nfactor; i++){
629  if(fabs(factor-FACTORS[i])<TOLERANCE){
630  if(DDistrVolume){
631  dV = 1./(pow(factor,3));
632  // dV1 = 1./(pow(factor,3)*factor_events_inj[i]->GetEntries());
633  cnt2++;
634  }
635  else if (DDistrUniform){
636 
637  dV = pow(range[1],2)/factor;
638 
639  }
640  // else {cout << "No defined distance distribution?????! Or different from uniform and volume???"<<endl;exit(1);}
641  else {cout << "No defined distance distribution?????! WARNING: Assuming uniform in volume"<<endl;
642  DDistrVolume=1;
643  dV = 1./(pow(factor,3));
644  }
645  dV1 = dV/factor_events_inj[i]->GetEntries();
646 
647 #ifdef MINCHI
648  dV_spin_mtot = dV/factor_events_spin_mtot_inj[i]->GetBinContent(cz+1,mt+1);
649 #endif
650  dV /= factor_events_inj[i]->GetBinContent(m+1,l+1);
651  // if(i==nfactor-1){factor_events_rec->Fill(mass[1],mass[0]);
652  //if(fabs(factor-INTERNAL_FACTOR)<TOLERANCE){factor_events_rec->Fill(mass[1],mass[0]);
653  if(fabs(factor-FACTORS[nfactor-1])<TOLERANCE){factor_events_rec->Fill(mass[1],mass[0]);
654  }
655  break;
656  }
657  }
658 
659  nT=TMath::Min(TMath::Floor((rho[pp_irho]-RHO_MIN)/RHO_BIN),(double)RHO_NBINS)+1;
660  for(int i=0; i<nT; i++){
661 //#ifdef VOL_Mtotq
662  Vrho[i]+=dV1;
663  eVrho[i]+=pow(dV1,2);
664 //#else
665 // Vrho[i]+=dV;
666 // eVrho[i]+=pow(dV,2);
667 //#endif
668  }
669  // cout << "rho[1]="<<rho[1]<<" nT="<<nT<<" Vrho[0]="<<Vrho[0] <<" Vrho[1]="<<Vrho[1]<<" Vrho[2]="<<Vrho[2]<<endl;
670  if(rho[pp_irho]<=T_out){continue;}
671 #ifdef WRITE_ASCII
672 
673  fprintf(fev,"%10.3f %3.2f %3.2f %3.2f\n",time[0],rho[pp_irho],sqrt(iSNR[0]+iSNR[1]+iSNR[2]),sqrt(sSNR[0]+sSNR[1]+sSNR[2]));
674 // if((m !=1)||(l !=1)) cout << "rho[1]="<<rho[1]<<" nT="<<nT<<" range[1]="<<range[1]<<" time[0]="<<time[0] <<" mass[0]="<<mass[0]<<" "<<l<<" mass[1]="<<mass[1]<<" "<<m<<endl;
675 
676 #endif
677 
678  volume[m][l] += dV;
679  error_volume[m][l] += pow(dV,2);
680 #ifdef MINCHI
681 
682  spin_mtot_volume[mt][cz] += dV_spin_mtot;
683  error_spin_mtot_volume[mt][cz] += pow(dV_spin_mtot,2);
684 #endif
685 
686  }
687  cout << endl;
688 fclose(fev);
689 
690  cout << "Recovered entries: " << cnt << endl;
691  cout << "Recovered entries: " << cnt2 << endl;
692  cout<< "Recovered entries vetoed: "<<countv<<endl;
693  cout <<"dV : "<<dV<< " dV1 : "<<dV1<<endl;
694 for(int i=0; i<RHO_NBINS; i++){eVrho[i]=TMath::Sqrt(eVrho[i]);}
695 cout << "Vrho[0] = "<<Vrho[0] <<" +/- "<<eVrho[0]<<endl;
696 cout<< "Vrho[RHO_NBINS-1] = "<<Vrho[RHO_NBINS-1] <<" +/- "<<eVrho[RHO_NBINS-1] <<endl;
697  //for(int i=0;i<mdc_entries;i++){inj_events->Fill(im1[i],im2[i]);iMtot[i]=im1[i]+im2[i];}
698  inj_events->Draw("colz");
699 
700  int MAX_AXIS_Z = inj_events->GetBinContent(inj_events->GetMaximumBin()) + 1;
701  inj_events->GetZaxis()->SetRangeUser(0,MAX_AXIS_Z);
702 
703 
704  char inj_title[256];
705  sprintf(inj_title,"Injected events");
706  //#ifdef MIN_CHI
707  //sprintf(inj_title,"%s (%.1f < #chi < %.1f)",inj_title,MIN_CHI,MAX_CHI);
708  //#endif
709 
710 
711  TPaveText *p_inj = new TPaveText(0.325301,0.926166,0.767068,0.997409,"blNDC");
712  p_inj->SetBorderSize(0);
713  p_inj->SetFillColor(0);
714  p_inj->SetTextColor(1);
715  p_inj->SetTextFont(32);
716  p_inj->SetTextSize(0.045);
717  TText *text = p_inj->AddText(inj_title);
718  p_inj->Draw();
719 
720  char fname[1024];
721  // #ifdef MIN_CHI
722  // sprintf(fname,"%s/Injected_mass1_mass2_chi_%f_%f.eps",netdir,MIN_CHI,MAX_CHI);
723  // #else
724  sprintf(fname,"%s/Injected_mass1_mass2.eps",netdir);
725  // #endif
726  c1->Update();
727  c1->SaveAs(fname);
728 
729 
730  //###############################################################################################################################
731  // Delta time
732  //###############################################################################################################################
733  c1->Clear();
734  c1->SetLogy(1);
735  Dt->Draw();
736  sprintf(fname,"%s/Delta_t.png",netdir);
737  c1->Update();
738  c1->SaveAs(fname);
739 
740  //###############################################################################################################################
741  // RHO CC/chi2 PLOTS
742  //###############################################################################################################################
743 
744 
745  c1->Clear();
746  if(pp_rho_log) c1->SetLogy(kTRUE); else c1->SetLogy(kFALSE);
747  rhocc->Draw("colz");
748  sprintf(fname,"%s/rho_cc.eps",netdir);
749  c1->Update(); c1->SaveAs(fname);
750 
751 
752  c1->Clear();
753  c1->SetLogx(kFALSE);
754  if(pp_rho_log) c1->SetLogy(kTRUE); else c1->SetLogy(kFALSE);
755  rho_pf->Draw("colz");
756  sprintf(fname,"%s/rho_pf.eps",netdir);
757  c1->Update(); c1->SaveAs(fname);
758 
759 
760 
761  //###############################################################################################################################
762  // SNR PLOTS
763  //###############################################################################################################################
764  char sel[1024]="";
765  char newcut[2048]="";
766  char newcut2[2048]="";
767  char title[1024]="";
768  TString ptitle;
769  TString xtitle;
770  TString ytitle;
771  c1->Clear();
772  c1->SetLogx(true);
773  c1->SetLogy(true);
774  ptitle="Number of reconstructed events as a function of injected SNR";
775  gStyle->SetOptStat(0);
776 
777  xtitle = "Injected network snr";
778  //xtitle = "sqrt(snr[0]^2";
779  // for(int i=1;i<nIFO;i++){xtitle += " + snr["+xtitle.Itoa(i,10)+"]^2";}
780  // xtitle +=")";
781  ytitle = "counts";
782 
783  // INJECTED
784  sprintf(newcut,"");
785  sprintf(sel,"sqrt(pow(snr[%d],2)",0);
786  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + pow(snr[%d],2)",sel,i);}
787 
788  sprintf(sel,"%s)>>hist(500)",sel);
789  mdc.Draw(sel,"");
790 
791  int nmdc = mdc.GetSelectedRows();
792  cout << "nmdc : " << nmdc << endl;
793 
794  TH2F *htemp = (TH2F*)gPad->GetPrimitive("hist");
795  htemp->GetXaxis()->SetTitle(xtitle);
796  htemp->GetYaxis()->SetTitle(ytitle);
797  //htemp->GetXaxis()->SetRangeUser(4e-25,1e-21);
798  //htemp->GetXaxis()->SetRangeUser(gTRMDC->GetMinimum("snr"),gTRMDC->GetMaximum("snr"));
799  //htemp->GetXaxis()->SetRangeUser(5, pow(10.,TMath::Ceil(TMath::Log10(htemp->GetMaximum()))));
800  htemp->GetYaxis()->SetRangeUser(1, pow(10.,TMath::Ceil(TMath::Log10(htemp->GetMaximum()))));
801  htemp->SetLineColor(kBlack);
802  htemp->SetLineWidth(3);
803  sprintf(newcut,"(((time[0]-time[%d])>-%g) || (time[0]-time[%d])<%g) && rho[%d]> %g",nIFO,T_win,nIFO,2*T_win,pp_irho,T_out);
804 
805  // DETECTED iSNR
806  sprintf(sel,"sqrt(iSNR[%d]",0);
807  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + iSNR[%d]",sel,i);}
808  sprintf(sel,"%s)>>hist2(500)",sel);
809  cout << "cut " << newcut << endl;
810 
811  // sim.SetFillColor(kBlue);
812  // htemp->SetLineColor(kBlue);
813  sim.Draw(sel,newcut,"same");
814  TH2F *htemp2 = (TH2F*)gPad->GetPrimitive("hist2");
815  htemp2->SetFillColor(kRed);
816  htemp2->SetFillStyle(3017);
817  htemp2->SetLineColor(kRed);
818  htemp2->SetLineWidth(2);
819 
820  int nwave = sim.GetSelectedRows();
821  cout << "nwave : " << nwave << endl;
822  sprintf(title,"%s",newcut);
823 
824  // DETECTED iSNR after final cuts
825  sprintf(newcut2,"%s && %s",newcut,veto_not_vetoed);
826  cout << "final cut " << newcut2 << endl;
827  TString sel_fin = sel;
828  sel_fin.ReplaceAll("hist2","hist3");
829 
830  // sim.SetFillColor(kRed);
831  //htemp->SetLineColor(kRed);
832  sim.Draw(sel_fin,newcut2,"same");
833  TH2F *htemp3 = (TH2F*)gPad->GetPrimitive("hist3");
834  htemp3->SetFillColor(kBlue);
835  htemp3->SetFillStyle(3017);
836  htemp3->SetLineColor(kBlue);
837  htemp3->SetLineWidth(2);
838  int nwave_final = sim.GetSelectedRows();
839  cout << "nwave_final : " << nwave_final << endl;
840  sprintf(title,"%s",newcut);
841 
842  sprintf(title,"%s",ptitle.Data());
843  htemp->SetTitle(title);
844  char lab[256];
845  leg_snr = new TLegend(0.6,0.755,0.885,0.923,"","brNDC");
846  sprintf(lab,"Injections Average SNR: %g",htemp->GetMean());
847  leg_snr->AddEntry("",lab,"a");
848  sprintf(lab,"Injected: %i",nmdc);
849  leg_snr->AddEntry(htemp, lab, "l");
850  sprintf(lab,"Found(minimal cuts): %i",nwave);
851  leg_snr->AddEntry(htemp2, lab, "l");
852  sprintf(lab,"Found(final cuts): %i",nwave_final);
853  leg_snr->AddEntry(htemp3, lab, "l");
854  leg_snr->SetFillColor(0);
855  leg_snr->Draw();
856 
857  sprintf(fname,"%s/Injected_snr_distributions.png",netdir);
858  c1->Update();
859  c1->SaveAs(fname);
860 
861 
862  //Plot estimated network snr vs injected network snr
863  sim.SetMarkerStyle(20);
864  sim.SetMarkerSize(0.5);
865  sim.SetMarkerColor(kRed);
866 
867 
868  sprintf(sel,"sqrt(sSNR[%d]",0);
869  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + sSNR[%d]",sel,i);}
870  sprintf(sel,"%s):sqrt(iSNR[%d]",sel,0);
871  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + iSNR[%d]",sel,i);}
872  sprintf(sel,"%s)>>hist4(500)",sel);
873  sim.Draw(sel,newcut,"");
874  TH2F *htemp4 = (TH2F*)gPad->GetPrimitive("hist4");
875  htemp4->GetXaxis()->SetTitle("Injected network snr");
876  htemp4->GetYaxis()->SetTitle("Estimated network snr");
877  htemp4->GetXaxis()->SetRangeUser(5, pow(10.,TMath::Ceil(TMath::Log10(htemp4->GetMaximum()))));
878  htemp4->GetYaxis()->SetRangeUser(5, pow(10.,TMath::Ceil(TMath::Log10(htemp4->GetMaximum()))));
879 
880  //htemp4->SetMarkerStyle(20);
881  //htemp4->SetMarkerSize(0.5);
882  //htemp4->SetMarkerColor(kRed);
883  sim.SetMarkerColor(kBlue);
884  sim.Draw(sel,newcut2,"same");
885  //TH2F *htemp5 = (TH2F*)gPad->GetPrimitive("hist4");
886  //htemp5->SetMarkerStyle(20);
887  //htemp5->SetMarkerSize(0.5);
888  //htemp5->SetMarkerColor(kBlue);
889 
890  sprintf(fname,"%s/Estimated_snr_vs_Injected_snr.eps",netdir);
891  c1->Update();
892  c1->SaveAs(fname);
893 
894  //Plot relative snr loss
895 
896  sprintf(sel,"(sqrt(sSNR[%d]",0);
897  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + sSNR[%d]",sel,i);}
898  sprintf(sel,"%s)-sqrt(iSNR[%d]",sel,0);
899  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + iSNR[%d]",sel,i);}
900  sprintf(sel,"%s))/sqrt(iSNR[%d]",sel,0);
901  for(int i=1;i<nIFO;i++){sprintf(sel,"%s + iSNR[%d]",sel,i);}
902  sprintf(sel,"%s)>>hist5",sel);
903  cout << "Selection: "<< sel <<endl;
904 
905  gStyle->SetOptStat(1);
906  gStyle->SetOptFit(1);
907  c1->SetLogx(false);
908 
909  sim.Draw(sel,newcut2);
910  TH2F *htemp5 = (TH2F*)gPad->GetPrimitive("hist5");
911  htemp5->GetXaxis()->SetTitle("(Estimated snr -Injected snr)/Injected snr");
912  htemp5->GetYaxis()->SetTitle("Counts");
913  sim.GetHistogram().Fit("gaus");
914 
915  sprintf(fname,"%s/Relative_snr_Loss.png",netdir);
916  c1->Update();
917  c1->SaveAs(fname);
918 
919  gStyle->SetOptStat(0);
920  gStyle->SetOptFit(0);
921  //###############################################################################################################################
922  // DISTANCE PLOTS
923  //###############################################################################################################################
924 
925 
926  c1->Clear();
927 
928 
929 D_Mtot_inj->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
930 D_Mtot_inj->Draw("p");
931 D_Mtot_rec->Draw("p same");
932 
933  leg_D = new TLegend(0.679719,0.156425,0.997992,0.327409,"","brNDC");
934 sprintf(lab,"Injections: %i",(int)mdc.GetEntries());
935  leg_D->AddEntry("",lab,"a");
936  sprintf(lab,"found: %i",cnt);
937  leg_D->AddEntry(D_Mtot_rec, lab, "p");
938  sprintf(lab,"missed: %i",(int)mdc.GetEntries()-cnt);
939 
940  leg_D->AddEntry(D_Mtot_inj,lab, "p");
941 
942  leg_D->SetFillColor(0);
943  leg_D->Draw();
944 
945 
946  sprintf(fname,"%s/Distance_vs_total_mass.eps",netdir);
947  c1->Update();
948  c1->SaveAs(fname);
949 
950 
951  c1->Clear();
952 D_Mchirp_inj->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
953 //D_Mchirp_rec->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
954 D_Mchirp_inj->Draw("p");
955 D_Mchirp_rec->Draw("p same");
956  leg_D->Draw();
957 
958 
959  sprintf(fname,"%s/Distance_vs_chirp_mass.eps",netdir);
960  c1->Update();
961  c1->SaveAs(fname);
962 #ifdef MINCHI
963  c1->Clear();
964 D_Chi_inj->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
965 //D_Chi_rec->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
966 D_Chi_inj->Draw("p");
967 D_Chi_rec->Draw("p same");
968  leg_D->Draw();
969 
970  sprintf(fname,"%s/Distance_vs_Chi.eps",netdir);
971  c1->Update();
972  c1->SaveAs(fname);
973 #endif
974 
975  c1->Clear();
976 D_q_inj->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
977 D_q_inj->Draw("p");
978 D_q_rec->Draw("p same");
979  leg_D->Draw();
980 
981  sprintf(fname,"%s/Distance_vs_mass_ratio.eps",netdir);
982  c1->Update();
983  c1->SaveAs(fname);
984 
985  c1->Clear();
986  leg_D = new TLegend(0.679719,0.826425,0.997992,0.997409,"","brNDC");
987 sprintf(lab,"Injections: %i",(int)mdc.GetEntries());
988  leg_D->AddEntry("",lab,"a");
989  sprintf(lab,"found: %i",cnt);
990  leg_D->AddEntry(D_Mtot_rec, lab, "p");
991  sprintf(lab,"missed: %i",(int)mdc.GetEntries()-cnt);
992 
993  leg_D->AddEntry(D_Mtot_inj,lab, "p");
994 
995 leg_D->SetFillColor(0);
996  leg_D->Draw();
997  c1->SetLogy(1);
998 // D_inj->GetXaxis()->SetRangeUser(10.,2000.);
999 TH1D* D_inj = D_Mtot_inj->ProjectionY();
1000 D_inj->GetXaxis()->SetTitle("Distance (Mpc)");
1001  D_inj->GetYaxis()->SetTitle("Events");
1002  D_inj->GetXaxis()->SetTickLength(0.02);
1003  D_inj->GetYaxis()->SetTickLength(0.01);
1004  D_inj->GetXaxis()->SetTitleOffset(1.3);
1005  D_inj->GetYaxis()->SetTitleOffset(1.3);
1006  D_inj->GetXaxis()->CenterTitle(kTRUE);
1007  D_inj->GetYaxis()->CenterTitle(kTRUE);
1008  D_inj->GetXaxis()->SetTitleFont(42);
1009  D_inj->GetXaxis()->SetLabelFont(42);
1010  D_inj->GetYaxis()->SetTitleFont(42);
1011  D_inj->GetYaxis()->SetLabelFont(42);
1012  D_inj->SetLineWidth(4);
1013  D_inj->SetLineColor(2);
1014  D_inj->SetFillColor(2);
1015  D_inj->SetFillStyle(3017);
1016  D_inj->SetTitle(0);
1017  D_inj->GetXaxis()->SetRangeUser(10.,MAX_EFFECTIVE_RADIUS);
1018  D_inj->GetYaxis()->SetRangeUser(0.1,pow(10.,TMath::Ceil(TMath::Log10(D_inj->GetMaximum()))));
1019 D_inj->Draw("lp");
1020 
1021 TH1D* D_rec=D_Mtot_rec->ProjectionY();
1022 D_rec->SetLineWidth(4);
1023  D_rec->SetLineColor(4);
1024  D_rec->SetFillColor(4);
1025  D_rec->SetFillStyle(3017);
1026  D_rec->SetTitle(0);
1027  D_rec->Draw("lp same");
1028  leg_D->Draw();
1029  TPaveText *p_radius = new TPaveText(0.125301,0.926166,0.567068,0.997409,"blNDC");
1030  p_radius->SetBorderSize(0);
1031  p_radius->SetFillColor(0);
1032  p_radius->SetTextColor(1);
1033  p_radius->SetTextFont(32);
1034  p_radius->SetTextSize(0.045);
1035  TText *text = p_radius->AddText("Found&Lost vs distance");
1036  p_radius->Draw();
1037 
1038  sprintf(fname,"%s/Injected_distances_distribution.png",netdir);
1039  c1->Update();
1040  c1->SaveAs(fname);
1041 
1042  c1->Clear();
1043 
1044  TH1D* Mtot_inj = D_Mtot_inj->ProjectionX();
1045 Mtot_inj->GetXaxis()->SetTitle("Total Mass");
1046  Mtot_inj->GetYaxis()->SetTitle("Events");
1047 Mtot_inj->GetXaxis()->SetTickLength(0.02);
1048  Mtot_inj->GetYaxis()->SetTickLength(0.01);
1049  Mtot_inj->GetXaxis()->SetTitleOffset(1.3);
1050  Mtot_inj->GetYaxis()->SetTitleOffset(1.3);
1051  Mtot_inj->GetXaxis()->CenterTitle(kTRUE);
1052  Mtot_inj->GetYaxis()->CenterTitle(kTRUE);
1053  Mtot_inj->GetXaxis()->SetTitleFont(42);
1054  Mtot_inj->GetXaxis()->SetLabelFont(42);
1055  Mtot_inj->GetYaxis()->SetLabelFont(42);
1056  Mtot_inj->SetLineWidth(4);
1057  Mtot_inj->SetLineColor(2);
1058  Mtot_inj->SetFillColor(2);
1059  Mtot_inj->SetFillStyle(3017);
1060  Mtot_inj->SetTitle(0);
1061 
1062  Mtot_inj->GetXaxis()->SetRangeUser(MINMtot,MAXMtot);
1063  Mtot_inj->GetYaxis()->SetRangeUser(1,pow(10.,TMath::Ceil(TMath::Log10(Mtot_inj->GetMaximum()))));
1064 Mtot_inj->Draw("lp");
1065  TH1D* Mtot_rec = D_Mtot_rec->ProjectionX();
1066  Mtot_rec->SetLineWidth(4);
1067  Mtot_rec->SetLineColor(4);
1068  Mtot_rec->SetFillColor(4);
1069  Mtot_rec->SetFillStyle(3017);
1070  Mtot_rec->SetTitle(0);
1071  Mtot_rec->Draw("lp same");
1072 leg_D->Draw();
1073 sprintf(fname,"%s/Total_mass_distribution.png",netdir);
1074  c1->Update();
1075  c1->SaveAs(fname);
1076 
1077  TH1D* Mchirp_inj = D_Mchirp_inj->ProjectionX();
1078 Mchirp_inj->GetXaxis()->SetTitle("Chirp Mass");
1079  Mchirp_inj->GetYaxis()->SetTitle("Events");
1080 Mchirp_inj->GetXaxis()->SetTickLength(0.02);
1081  Mchirp_inj->GetYaxis()->SetTickLength(0.01);
1082  Mchirp_inj->GetXaxis()->SetTitleOffset(1.3);
1083  Mchirp_inj->GetYaxis()->SetTitleOffset(1.3);
1084  Mchirp_inj->GetXaxis()->CenterTitle(kTRUE);
1085  Mchirp_inj->GetYaxis()->CenterTitle(kTRUE);
1086  Mchirp_inj->GetXaxis()->SetTitleFont(42);
1087  Mchirp_inj->GetXaxis()->SetLabelFont(42);
1088  Mchirp_inj->GetYaxis()->SetTitleFont(42);
1089  Mchirp_inj->GetYaxis()->SetLabelFont(42);
1090  Mchirp_inj->SetLineWidth(4);
1091  Mchirp_inj->SetLineColor(2);
1092  Mchirp_inj->SetFillColor(2);
1093  Mchirp_inj->SetFillStyle(3017);
1094  Mchirp_inj->SetTitle(0);
1095 
1096  Mchirp_inj->GetXaxis()->SetRangeUser(MINCHIRP,MAXCHIRP);
1097  Mchirp_inj->GetYaxis()->SetRangeUser(1,pow(10.,TMath::Ceil(TMath::Log10(Mtot_inj->GetMaximum()))));
1098 Mchirp_inj->Draw("lp");
1099  TH1D* Mchirp_rec = D_Mchirp_rec->ProjectionX();
1100  Mchirp_rec->SetLineWidth(4);
1101  Mchirp_rec->SetLineColor(4);
1102  Mchirp_rec->SetFillColor(4);
1103  Mchirp_rec->SetFillStyle(3017);
1104  Mchirp_rec->SetTitle(0);
1105  Mchirp_rec->Draw("lp same");
1106 leg_D->Draw();
1107 sprintf(fname,"%s/Chirp_mass_distribution.png",netdir);
1108  c1->Update();
1109  c1->SaveAs(fname);
1110 
1111 #ifdef MINCHI
1112  TH1D* Chi_inj = D_Chi_inj->ProjectionX();
1113 Chi_inj->GetXaxis()->SetTitle("#chi_{z}");
1114  Chi_inj->GetYaxis()->SetTitle("Events");
1115  Chi_inj->GetXaxis()->SetTickLength(0.02);
1116  Chi_inj->GetYaxis()->SetTickLength(0.01);
1117  Chi_inj->GetXaxis()->SetTitleOffset(1.3);
1118  Chi_inj->GetYaxis()->SetTitleOffset(1.3);
1119  Chi_inj->GetXaxis()->CenterTitle(kTRUE);
1120  Chi_inj->GetYaxis()->CenterTitle(kTRUE);
1121  Chi_inj->GetXaxis()->SetTitleFont(42);
1122  Chi_inj->GetXaxis()->SetLabelFont(42);
1123  Chi_inj->GetYaxis()->SetTitleFont(42);
1124  Chi_inj->GetYaxis()->SetLabelFont(42);
1125  Chi_inj->SetLineWidth(4);
1126  Chi_inj->SetLineColor(2);
1127  Chi_inj->SetFillColor(2);
1128  Chi_inj->SetFillStyle(3017);
1129  Chi_inj->SetTitle(0);
1130 
1131  Chi_inj->GetXaxis()->SetRangeUser(MINCHI,MAXCHI);
1132  Chi_inj->GetYaxis()->SetRangeUser(1,pow(10.,TMath::Ceil(TMath::Log10(Mtot_inj->GetMaximum()))));
1133  Chi_inj->Draw("lp");
1134  TH1D* Chi_rec = D_Chi_rec->ProjectionX();
1135  Chi_rec->SetLineWidth(4);
1136  Chi_rec->SetLineColor(4);
1137  Chi_rec->SetFillColor(4);
1138  Chi_rec->SetFillStyle(3017);
1139  Chi_rec->SetTitle(0);
1140  Chi_rec->Draw("lp same");
1141 leg_D->Draw();
1142 sprintf(fname,"%s/Chi_distribution.png",netdir);
1143  c1->Update();
1144  c1->SaveAs(fname);
1145 
1146 TH1D* Chi_injx = D_Chi_injx->ProjectionX();
1147 Chi_injx->SetLineWidth(4);
1148  Chi_injx->SetLineColor(2);
1149  Chi_injx->SetFillColor(2);
1150  Chi_injx->SetFillStyle(3017);
1151  Chi_injx->SetTitle(0);
1152 Chi_injx->Draw("lp");
1153  TH1D* Chi_recx = D_Chi_recx->ProjectionX();
1154 Chi_recx->SetLineWidth(4);
1155  Chi_recx->SetLineColor(4);
1156  Chi_recx->SetFillColor(4);
1157  Chi_recx->SetFillStyle(3017);
1158  Chi_recx->SetTitle(0);
1159 Chi_recx->Draw("lp same");
1160 leg_D->Draw();
1161 sprintf(fname,"%s/Chi_x_distribution.png",netdir);
1162  c1->Update();
1163  c1->SaveAs(fname);
1164 
1165 TH1D* Chi_injy = D_Chi_injy->ProjectionX();
1166 Chi_injy->SetLineWidth(4);
1167  Chi_injy->SetLineColor(2);
1168  Chi_injy->SetFillColor(2);
1169  Chi_injy->SetFillStyle(3017);
1170  Chi_injy->SetTitle(0);
1171 Chi_injy->Draw("lp");
1172  TH1D* Chi_recy = D_Chi_recy->ProjectionX();
1173 Chi_recy->SetLineWidth(4);
1174  Chi_recy->SetLineColor(4);
1175  Chi_recy->SetFillColor(4);
1176  Chi_recy->SetFillStyle(3017);
1177  Chi_recy->SetTitle(0);
1178 Chi_recy->Draw("lp same");
1179 leg_D->Draw();
1180 sprintf(fname,"%s/Chi_y_distribution.png",netdir);
1181  c1->Update();
1182  c1->SaveAs(fname);
1183 
1184 #endif
1185 
1186  TH1D* q_inj = D_q_inj->ProjectionX();
1187 q_inj->GetXaxis()->SetTitle("Mass Ratio");
1188  q_inj->GetYaxis()->SetTitle("Events");
1189 q_inj->GetXaxis()->SetTickLength(0.02);
1190  q_inj->GetYaxis()->SetTickLength(0.01);
1191  q_inj->GetXaxis()->SetTitleOffset(1.3);
1192  q_inj->GetYaxis()->SetTitleOffset(1.3);
1193  q_inj->GetXaxis()->CenterTitle(kTRUE);
1194  q_inj->GetYaxis()->CenterTitle(kTRUE);
1195  q_inj->GetXaxis()->SetTitleFont(42);
1196  q_inj->GetXaxis()->SetLabelFont(42);
1197  q_inj->GetYaxis()->SetTitleFont(42);
1198  q_inj->GetYaxis()->SetLabelFont(42);
1199  q_inj->SetLineWidth(4);
1200  q_inj->SetLineColor(2);
1201  q_inj->SetFillColor(2);
1202  q_inj->SetFillStyle(3017);
1203  q_inj->SetTitle(0);
1204 
1205  q_inj->GetXaxis()->SetRangeUser(MINRATIO,MAXRATIO);
1206  q_inj->GetYaxis()->SetRangeUser(1,pow(10.,TMath::Ceil(TMath::Log10(Mtot_inj->GetMaximum()))));
1207 q_inj->Draw("lp");
1208  TH1D* q_rec = D_q_rec->ProjectionX();
1209  q_rec->SetLineWidth(4);
1210  q_rec->SetLineColor(4);
1211  q_rec->SetFillColor(4);
1212  q_rec->SetFillStyle(3017);
1213  q_rec->SetTitle(0);
1214  q_rec->Draw("lp same");
1215 leg_D->Draw();
1216 sprintf(fname,"%s/Mass_ratio_distribution.png",netdir);
1217  c1->Update();
1218  c1->SaveAs(fname);
1219 
1220 
1221 
1222  c1->SetLogy(0);
1223 
1224  //###############################################################################################################################
1225  // Efficiency
1226  //###############################################################################################################################
1227 
1228 
1229  c1->Clear();
1230  TH2F *efficiency = (TH2F*) rec_events->Clone();
1231  efficiency->SetName("efficiency");
1232  efficiency->Divide(inj_events);
1233  efficiency->GetZaxis()->SetRangeUser(0,1.0);
1234  efficiency->SetTitle("");
1235 
1236 
1237  efficiency->Draw("colz text");
1238 
1239 
1240  TExec *ex1 = new TExec("ex1","gStyle->SetPaintTextFormat(\".2f\");");
1241  ex1->Draw();
1242 
1243 
1244  char eff_title[256];
1245  sprintf(eff_title,"Efficiency");
1246  #ifdef MIN_CHI
1247  sprintf(eff_title,"%s (%.1f < #chi < %.1f)",eff_title,MIN_CHI,MAX_CHI);
1248  #endif
1249 
1250 
1251  TPaveText *p_eff = new TPaveText(0.325301,0.926166,0.767068,0.997409,"blNDC");
1252  p_eff->SetBorderSize(0);
1253  p_eff->SetFillColor(0);
1254  p_eff->SetTextColor(1);
1255  p_eff->SetTextFont(32);
1256  p_eff->SetTextSize(0.045);
1257  TText *text = p_eff->AddText(eff_title);
1258  p_eff->Draw();
1259 
1260 
1261 
1262  #ifdef MIN_CHI
1263  sprintf(fname,"%s/Efficiency_mass1_mass2_chi_%f_%f.png",netdir,MIN_CHI,MAX_CHI);
1264  #else
1265  sprintf(fname,"%s/Efficiency_mass1_mass2.png",netdir);
1266  #endif
1267  c1->Update();
1268  c1->SaveAs(fname);
1269 
1270 
1271 
1272 
1273  //###############################################################################################################################
1274  // Entries in the shell corresponding to the largest factor
1275  //###############################################################################################################################
1276 
1277  if(DDistrVolume){
1278  float shell_volume = 4.*TMath::Pi()*(pow(maxDistance/1000.,3) - pow(minDistance/1000.,3))/3;
1279  float volume_internal_sphere = 4.*TMath::Pi()*pow(minDistance/1000.,3)/(3*pow(FACTORS[nfactor-1],3));
1280  }
1281  else if (DDistrUniform){
1282  float shell_volume = 4.*TMath::Pi()*(maxDistance/1000. - minDistance/1000.);
1283  float volume_internal_sphere = 4.*TMath::Pi()*minDistance/1000./FACTORS[nfactor-1];
1284  }
1285  else {cout << "No defined distance distribution?????!"<<endl;exit(1);}
1286 
1287 
1288  cout << "Shell volume: " << shell_volume << ", internal sphere volume: " << volume_internal_sphere << endl;
1289 
1290  TH2F *efficiency_first_shell = (TH2F*) factor_events_rec->Clone();
1291  efficiency_first_shell->Divide(factor_events_inj[nfactor-1]);
1292 
1293 
1294  //###############################################################################################################################
1295  // Effective radius calculation
1296  //###############################################################################################################################
1297 
1298 
1299  int massbins = 0;
1300  double V0 = 0.0;
1301  for(int j=1; j<NBINS_mass; j++){
1302  for(int k=0; k<NBINS_mass2; k++){
1303 
1304  volume_first_shell[j][k] = efficiency_first_shell->GetBinContent(j+1,k+1);
1305  if(factor_events_rec->GetBinContent(j+1,k+1) != 0.) {
1306  error_volume_first_shell[j][k] = 1./TMath::Sqrt(factor_events_rec->GetBinContent(j+1,k+1));
1307  massbins++;
1308  }
1309 
1310  volume[j][k] = shell_volume*volume[j][k] + volume_internal_sphere*volume_first_shell[j][k];
1311  V0 += volume[j][k];
1312  error_volume[j][k] = shell_volume*TMath::Sqrt(error_volume[j][k]) + volume_internal_sphere*volume_first_shell[j][k]*error_volume_first_shell[j][k];
1313 
1314  radius[j][k] = pow(3.*volume[j][k]/(4*TMath::Pi()),1./3);
1315  if(volume[j][k] != 0.) error_radius[j][k] = (1./3)*pow(3./(4*TMath::Pi()),1./3)*pow(1./pow(volume[j][k],2),1./3)*error_volume[j][k];
1316 
1317  }
1318  }
1319  cout << "Average Volume at threshold V0 = "<<V0/massbins <<endl;
1320  c1->Clear();
1321 
1322  TH2F* h_radius = new TH2F("h_radius","",NBINS_mass,MIN_MASS,MAX_MASS,NBINS_mass2,min_mass2,max_mass2);
1323  h_radius->GetXaxis()->SetRangeUser(MIN_plot_mass1,MAX_plot_mass1);
1324  h_radius->GetYaxis()->SetRangeUser(MIN_plot_mass2,MAX_plot_mass2);
1325  h_radius->GetXaxis()->SetTitle("Mass 1 (M_{#odot})");
1326  h_radius->GetYaxis()->SetTitle("Mass 2 (M_{#odot})");
1327  h_radius->GetXaxis()->SetTitleOffset(1.3);
1328  h_radius->GetYaxis()->SetTitleOffset(1.25);
1329  h_radius->GetXaxis()->CenterTitle(kTRUE);
1330  h_radius->GetYaxis()->CenterTitle(kTRUE);
1331  h_radius->GetXaxis()->SetNdivisions(410);
1332  h_radius->GetYaxis()->SetNdivisions(410);
1333  h_radius->GetXaxis()->SetTickLength(0.01);
1334  h_radius->GetYaxis()->SetTickLength(0.01);
1335  h_radius->GetZaxis()->SetTickLength(0.01);
1336  h_radius->GetXaxis()->SetTitleFont(42);
1337  h_radius->GetXaxis()->SetLabelFont(42);
1338  h_radius->GetYaxis()->SetTitleFont(42);
1339  h_radius->GetYaxis()->SetLabelFont(42);
1340  h_radius->GetZaxis()->SetLabelFont(42);
1341  h_radius->GetZaxis()->SetLabelSize(0.03);
1342  h_radius->SetTitle("");
1343 
1344 
1345  for(int i=1; i<=NBINS_mass; i++){
1346  for(int j=1; j<=NBINS_mass2; j++){
1347  h_radius->SetBinContent(i,j,radius[i-1][j-1]);
1348  h_radius->SetBinError(i,j,error_radius[i-1][j-1]);
1349  }
1350  }
1351 
1352 
1353  h_radius->SetContour(NCont);
1354  h_radius->SetEntries(1); // This option needs to be enabled when filling 2D histogram with SetBinContent
1355  h_radius->Draw("colz text colsize=2"); // Option to write error associated to the bin content
1356  // h_radius->Draw("colz text");
1357  h_radius->GetZaxis()->SetRangeUser(0,MAX_EFFECTIVE_RADIUS/2.);
1358 
1359  TExec *ex2 = new TExec("ex2","gStyle->SetPaintTextFormat(\".0f\");");
1360  ex2->Draw();
1361 
1362 
1363  char radius_title[256];
1364 
1365  sprintf(radius_title,"%s : Effective radius (Mpc)",networkname);
1366 
1367 
1368  TPaveText *p_radius = new TPaveText(0.325301,0.926166,0.767068,0.997409,"blNDC");
1369  p_radius->SetBorderSize(0);
1370  p_radius->SetFillColor(0);
1371  p_radius->SetTextColor(1);
1372  p_radius->SetTextFont(32);
1373  p_radius->SetTextSize(0.045);
1374  TText *text = p_radius->AddText(radius_title);
1375  p_radius->Draw();
1376 
1377 
1378  #ifdef PLOT_CHIRP
1379 
1380  TH2F *chirp_mass = new TH2F("Chirp_Mass","",NBINS_mass*10,MIN_plot_mass1,MAX_plot_mass1*1.1,NBINS_mass2*10,MIN_plot_mass2,MAX_plot_mass2*1.1);
1381  float M1,M2;
1382  for (Int_t i = 0; i < NBINS_mass*10; i++) {
1383  for(Int_t j = 0; j < NBINS_mass2*10; j++){
1384  M1=MIN_plot_mass1+i*(MAX_plot_mass1-MIN_plot_mass1*1.1)/NBINS_mass/10.;
1385  M2=MIN_plot_mass2+j*(MAX_plot_mass2-MIN_plot_mass2*1.1)/NBINS_mass/10.;
1386  chirp_mass->SetBinContent(i,j,(float) TMath::Power(pow(M1*M2,3.)/(M1+M2),1./5));
1387  }
1388  }
1389  Double_t contours[CONTOURS];
1390  //contours[0] = 35.;
1391  for(int i=0;i<CONTOURS;i++){contours[i] = TMath::Nint((i+1)*MAXCHIRP/CONTOURS);}
1392  chirp_mass->GetZaxis()->SetRangeUser(0,MAXCHIRP);
1393  chirp_mass->SetContour(CONTOURS, contours);
1394 
1395 
1396  TCanvas* c1t = new TCanvas("c1t","Contour List",610,0,600,600);
1397  c1t->SetTopMargin(0.15);
1398 
1399  // Draw contours as filled regions, and Save points
1400  chirp_mass->Draw("CONT Z LIST");
1401  c1t->Update();
1402 
1403  // Get Contours
1404  TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
1405  TList* contLevel = NULL;
1406  TGraph* curv = NULL;
1407  TGraph* gc = NULL;
1408 
1409  Int_t nGraphs = 0;
1410  Int_t TotalConts = 0;
1411 
1412  if (conts == NULL){
1413  printf("*** No Contours Were Extracted!\n");
1414  TotalConts = 0;
1415  return;
1416  } else {
1417  TotalConts = conts->GetSize();
1418  }
1419 
1420  printf("TotalConts = %d\n", TotalConts);
1421 
1422  for(i = 0; i < TotalConts; i++){
1423  contLevel = (TList*)conts->At(i);
1424  printf("Contour %d has %d Graphs\n", i, contLevel->GetSize());
1425  nGraphs += contLevel->GetSize();
1426  }
1427 
1428  nGraphs = 0;
1429 
1430  Double_t x0, y0, z0;
1431  TLatex Tl;
1432  Tl.SetTextSize(0.02);
1433  char val[20];
1434  //break;
1435  for(i = 0; i < TotalConts-1; i++){
1436  contLevel = (TList*)conts->At(i);
1437  z0 = contours[i+1]; //Temporary Patch
1438 
1439  printf("Z-Level Passed in as: Z = %f\n", z0);
1440 
1441  // Get first graph from list on curves on this level
1442  curv = (TGraph*)contLevel->First();
1443  int point;
1444  for(j = 0; j < contLevel->GetSize(); j++){
1445  point = curv->GetN()-1;
1446  curv->GetPoint(point, x0, y0);
1447  point--;
1448  //printf("\tCoordinate Point # %d : x0 = %f y0 = %f \n",point,x0,y0);
1449 
1450  while (((x0 < MIN_plot_mass1)||(x0 > MAX_plot_mass1)||(y0 < MIN_plot_mass2)||(y0 > MAX_plot_mass2))&&(point>0)){
1451 
1452  curv->GetPoint(point, x0, y0);
1453  //printf("\tCoordinate Point # %d : x0 = %f y0 = %f \n",point,x0,y0);
1454  point--;
1455 
1456 
1457  }
1458  curv->SetLineWidth(1);
1459  curv->SetLineStyle(3);
1460  //curv->SetLineColor(2);
1461  nGraphs ++;
1462  printf("\tGraph: %d -- %d Elements\n", nGraphs,curv->GetN());
1463  printf("\tCoordinate Point # %d : x0 = %f y0 = %f \n",point,x0,y0);
1464  // Draw clones of the graphs to avoid deletions in case the 1st
1465  // pad is redrawn.
1466  c1->cd();
1467  gc = (TGraph*)curv->Clone();
1468  gc->Draw("C");
1469  //sprintf(val,"#color[2]{#scale[0.9]{#it{M_{chirp}}=%g}}",z0);
1470  //sprintf(val,"#scale[0.9]{#it{M_{chirp}}=%g}",z0);
1471  sprintf(val,"#it{M_{c}}=%g",z0);
1472  // if(i==0){x0 = 110.;y0=10v.;}
1473  Tl.DrawLatex(x0-1.,y0+0.5,val);
1474 
1475  //x0=0.;
1476  //y0=0.;
1477  //curv = (TGraph*)contLevel->After(curv); // Get Next graph
1478  }
1479  }
1480  c1->Update();
1481 
1482 #endif
1483 
1484 
1485 #ifdef PLOT_MASS_RATIO
1486  TH2F *mass_ratio = new TH2F("Mass_Ratio","",NBINS_mass*10,MIN_plot_mass1,MAX_plot_mass1*1.1,NBINS_mass2*10,MIN_plot_mass2,MAX_plot_mass2*1.1);
1487 
1488  float M1,M2;
1489 
1490  for(int i=0; i<NBINS_mass*10; i++) {
1491  for(int j=0; j<NBINS_mass*10; j++){
1492 
1493  M1 = MIN_MASS + i*(MAX_plot_mass1-MIN_plot_mass1*1.1)/NBINS_mass/10.;
1494  M2 = MIN_MASS + j*(MAX_plot_mass2-MIN_plot_mass2*1.1)/NBINS_mass/10.;
1495 
1496  mass_ratio->SetBinContent(i,j,(float) M1/M2);
1497 
1498  }
1499  }
1500 
1501  Double_t contours2[5];
1502  contours2[0] = 1.;
1503  contours2[1] = 1.5;
1504  contours2[2] = 2.;
1505  contours2[3] = 3.;
1506  contours2[4] = 4.;
1507 
1508 
1509  mass_ratio->SetContour(5, contours2);
1510 
1511  TCanvas* c1t2 = new TCanvas("c1t2","Contour List2",610,0,600,600);
1512  c1t2->SetTopMargin(0.15);
1513 
1514  // Draw contours as filled regions, and Save points
1515  mass_ratio->Draw("CONT Z LIST");
1516  c1t2->Update();
1517 
1518  // Get Contours
1519  TObjArray *conts2 = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
1520  TList* contLevel2 = NULL;
1521  TGraph* curv2 = NULL;
1522  TGraph* gc2 = NULL;
1523 
1524  int nGraphs = 0;
1525  int TotalConts = 0;
1526 
1527  if(conts2 == NULL){ printf("*** No Contours Were Extracted!\n"); TotalConts = 0; return; }
1528  else{ TotalConts = conts2->GetSize(); }
1529  printf("TotalConts = %d\n", TotalConts);
1530 
1531  for(i=0; i<TotalConts; i++){
1532  contLevel2 = (TList*)conts2->At(i);
1533  printf("Contour %d has %d Graphs\n", i, contLevel2->GetSize());
1534  nGraphs += contLevel2->GetSize();
1535  }
1536 
1537  nGraphs = 0;
1538  Double_t x0, y0, z0;
1539  TLatex Tl2;
1540  Tl2.SetTextSize(0.02);
1541  char val[20];
1542 
1543  for(i=0; i<TotalConts; i++){
1544 
1545  contLevel2 = (TList*)conts2->At(i);
1546  z0 = contours2[i];
1547 
1548  printf("Z-Level Passed in as: Z = %f\n", z0);
1549 
1550  // Get first graph from list on curves on this level
1551  curv2 = (TGraph*)contLevel2->First();
1552  int point;
1553 
1554  for(j=0; j<contLevel2->GetSize(); j++){
1555 
1556  point = curv2->GetN()-1;
1557  curv2->GetPoint(point, x0, y0);
1558  point--;
1559  printf("\tCoordinate Point # %d : x0 = %f y0 = %f \n",point,x0,y0);
1560 
1561  while(((x0 < MIN_plot_mass1)||(x0 > MAX_plot_mass1-4.)||(y0 < MIN_plot_mass2)||(y0 > MAX_plot_mass2 -2))&&(point>0)){
1562  curv2->GetPoint(point, x0, y0);
1563  point--;
1564 
1565  }
1566 
1567  curv2->SetLineWidth(1);
1568  curv2->SetLineStyle(3);
1569 
1570  nGraphs ++;
1571  printf("\tGraph: %d -- %d Elements\n", nGraphs,curv2->GetN());
1572  printf("\tCoordinate Point # %d : x0 = %f y0 = %f \n",point,x0,y0);
1573 
1574  //Draw clones of the graphs to avoid deletions in case the 1st pad is redrawn.
1575  c1->cd();
1576  gc2 = (TGraph*)curv2->Clone();
1577  gc2->Draw("C");
1578 
1579  sprintf(val,"#it{q}=%.2g",1./z0);
1580  Tl2.DrawLatex(x0,y0,val);
1581 
1582  // curv2 = (TGraph*)contLevel2->After(curv2); // Get Next graph
1583 
1584  }
1585  }
1586  c1->Update();
1587  #endif
1588 
1589 
1590 
1591  // #ifdef MIN_CHI
1592  // sprintf(fname,"%s/Effective_radius_chi_%f_%f.png",netdir,MIN_CHI,MAX_CHI);
1593  // #else
1594  sprintf(fname,"%s/Effective_radius.png",netdir);
1595  //sprintf(fname,"%s/Effective_radius.png",netdir);
1596  #endif
1597  c1->Update();
1598  c1->SaveAs(fname);
1599 
1600 
1601 
1602 //###############################################################################################################################
1603  // Effective radius spin-mtot calculation
1604  //###############################################################################################################################
1605 
1606 #ifdef MINCHI
1607  int spin_mtot_bins = 0;
1608  double V0_spin_mtot = 0.0;
1609  for(int j=0; j<NBINS_MTOT+1; j++){
1610  for(int k=0; k<NBINS_SPIN+1; k++){
1611 
1612  // volume_first_shell[j][k] = efficiency_first_shell->GetBinContent(j+1,k+1);
1613  // if(factor_events_rec->GetBinContent(j+1,k+1) != 0.) {
1614  // error_volume_first_shell[j][k] = 1./TMath::Sqrt(factor_events_rec->GetBinContent(j+1,k+1));
1615  // massbins++;
1616  //}
1617  if(spin_mtot_volume[j][k] != 0.) {
1618  spin_mtot_volume[j][k] = shell_volume*spin_mtot_volume[j][k]; /// Warning: neglecting the internal volume... + volume_internal_sphere*volume_first_shell[j][k];
1619  V0_spin_mtot += spin_mtot_volume[j][k];
1620  error_spin_mtot_volume[j][k] = shell_volume*TMath::Sqrt(error_spin_mtot_volume[j][k]); ///Warning: neglecting the internal volume...+ volume_internal_sphere*volume_first_shell[j][k]*error_volume_first_shell[j][k];
1621 
1622  spin_mtot_radius[j][k] = pow(3.*spin_mtot_volume[j][k]/(4*TMath::Pi()),1./3);
1623 
1624  error_spin_mtot_radius[j][k] = (1./3)*pow(3./(4*TMath::Pi()),1./3)*pow(1./pow(spin_mtot_volume[j][k],2),1./3)*error_spin_mtot_volume[j][k];
1625  spin_mtot_bins++;
1626  }
1627  cout<<j<< " "<<k<< " "<<spin_mtot_volume[j][k]<<" "<<spin_mtot_radius[j][k]<<endl;
1628  }
1629  }
1630  cout << "Average Spin-Mtot Volume at threshold V0 = "<<V0_spin_mtot/spin_mtot_bins <<endl;
1631  c1->Clear();
1632 
1633  TH2F* h_spin_mtot_radius = new TH2F("h_spin_mtot_radius","",NBINS_SPIN,MINCHI,MAXCHI,NBINS_MTOT,minMtot,maxMtot);
1634  //h_spin_mtot_radius->GetXaxis()->SetRangeUser(MIN_plot_mass1,MAX_plot_mass1);
1635  //h_spin_mtot_radius->GetYaxis()->SetRangeUser(MIN_plot_mass2,MAX_plot_mass2);
1636  h_spin_mtot_radius->GetXaxis()->SetTitle("#chi_{z}");
1637  h_spin_mtot_radius->GetYaxis()->SetTitle("Total Mass (M_{#odot})");
1638  h_spin_mtot_radius->GetXaxis()->SetTitleOffset(1.3);
1639  h_spin_mtot_radius->GetYaxis()->SetTitleOffset(1.25);
1640  h_spin_mtot_radius->GetXaxis()->CenterTitle(kTRUE);
1641  h_spin_mtot_radius->GetYaxis()->CenterTitle(kTRUE);
1642  h_spin_mtot_radius->GetXaxis()->SetNdivisions(410);
1643  h_spin_mtot_radius->GetYaxis()->SetNdivisions(410);
1644  h_spin_mtot_radius->GetXaxis()->SetTickLength(0.01);
1645  h_spin_mtot_radius->GetYaxis()->SetTickLength(0.01);
1646  h_spin_mtot_radius->GetZaxis()->SetTickLength(0.01);
1647  h_spin_mtot_radius->GetXaxis()->SetTitleFont(42);
1648  h_spin_mtot_radius->GetXaxis()->SetLabelFont(42);
1649  h_spin_mtot_radius->GetYaxis()->SetTitleFont(42);
1650  h_spin_mtot_radius->GetYaxis()->SetLabelFont(42);
1651  h_spin_mtot_radius->GetZaxis()->SetLabelFont(42);
1652  h_spin_mtot_radius->GetZaxis()->SetLabelSize(0.03);
1653  h_spin_mtot_radius->SetTitle("");
1654 
1655 
1656  for(int i=1; i<=NBINS_MTOT+1; i++){
1657  for(int j=1; j<=NBINS_SPIN+1; j++){
1658  h_spin_mtot_radius->SetBinContent(j,i,spin_mtot_radius[i-1][j-1]);
1659  h_spin_mtot_radius->SetBinError(j,i,error_spin_mtot_radius[i-1][j-1]);
1660  cout<<j<< " "<<i<< " "<<h_spin_mtot_radius->GetBinContent(j,i)<<endl;
1661  }
1662  }
1663 
1664 
1665  h_spin_mtot_radius->SetContour(NCont);
1666  h_spin_mtot_radius->SetEntries(1); // This option needs to be enabled when filling 2D histogram with SetBinContent
1667  h_spin_mtot_radius->Draw("colz text colsize=2"); // Option to write error associated to the bin content
1668  //h_spin_mtot_radius->Draw("colz text");
1669  h_spin_mtot_radius->GetZaxis()->SetRangeUser(0,MAX_EFFECTIVE_RADIUS/2.);
1670 
1671  TExec *ex2 = new TExec("ex2","gStyle->SetPaintTextFormat(\".0f\");");
1672  ex2->Draw();
1673 
1674 
1675  //char radius_title[256];
1676  sprintf(radius_title,"Effective radius Mtot vs #chi_{z} (Mpc, %.1f < #chi_{z} < %.1f)",MINCHI,MAXCHI);
1677 
1678 
1679 
1680  TPaveText *p_radius = new TPaveText(0.325301,0.926166,0.767068,0.997409,"blNDC");
1681  p_radius->SetBorderSize(0);
1682  p_radius->SetFillColor(0);
1683  p_radius->SetTextColor(1);
1684  p_radius->SetTextFont(32);
1685  p_radius->SetTextSize(0.045);
1686  TText *text = p_radius->AddText(radius_title);
1687  p_radius->Draw();
1688 
1689  sprintf(fname,"%s/Effective_radius_chi.png",netdir);
1690 
1691  c1->Update();
1692  c1->SaveAs(fname);
1693 
1694 
1695 #endif
1696  //###############################################################################################################################
1697 // Calculating Radius vs rho
1698 //###############################################################################################################################
1699 
1700 
1701  for(int i=0; i<RHO_NBINS; i++){
1702  Vrho[i]*=shell_volume;
1703  eVrho[i]*=shell_volume;
1704 //#ifndef VOL_Mtotq
1705 // Vrho[i]/=massbins;
1706 //#endif
1707  Rrho[i] = pow(3.*Vrho[i]/(4*TMath::Pi()),1./3);
1708  eRrho[i] = pow(3./(4*TMath::Pi()),1./3)*1./3*pow(Vrho[i],-2./3)*eVrho[i];
1709  }
1710 TF1* f2 = cbcTool.doRangePlot(RHO_NBINS, Trho, Rrho, eRrho, RHO_MIN, T_out, c1, networkname, netdir, write_ascii);
1711 
1712 
1713  //###############################################################################################################################
1714 // Calculating background outlier
1715 //###############################################################################################################################
1716 
1717  #ifdef BKG_NTUPLE
1718 
1719 
1720  char net_file_name[1024];
1721  char liv_file_name[1024];
1722 
1723  if(BKG_NTUPLE=="True"){
1725 
1726  if(gSystem->Getenv("CWB_BKG_TREE")==NULL) {
1727  cout << "Error : environment CWB_BKG_TREE is not defined!!!" << endl;exit(1);
1728  } else {
1729  Compare_tree=TString(gSystem->Getenv("CWB_BKG_TREE"));
1730 
1731 
1732  }
1733 
1734  sprintf(net_file_name,"%s",Compare_tree.Data());
1735  TString Liv_file_name = Compare_tree;
1736  }
1737  else {
1738  sprintf(net_file_name,"%s",BKG_NTUPLE);
1739  TString Liv_file_name = BKG_NTUPLE;
1740 
1741  }
1742 
1743  Liv_file_name.ReplaceAll("wave_","live_");
1744  sprintf(liv_file_name,"%s",Liv_file_name.Data());
1745 
1746  cout << net_file_name << endl;
1747  cout << liv_file_name << endl;
1748 
1749  // Check if files exist
1750  TB.checkFile(net_file_name);
1751  TB.checkFile(liv_file_name);
1752  TChain wave("waveburst");
1753  TChain liv("liveTime");
1754  wave.Add(net_file_name);
1755  liv.Add(liv_file_name);
1756 
1757 
1758 
1759  Int_t bkg_size = (Int_t)wave.GetEntries();
1760  cout << "bkg size : " << bkg_size << endl;
1761 
1762 //###############################################################################################################################
1763 // Live time calculation (after CAT2)
1764 //###############################################################################################################################
1765 
1766  #ifdef SKIP_GETLIVETIME
1767  cout <<"WARNING!!! Skipping actual calculation of BKG and zero lag live time"<<endl;
1768  #else
1769  double liveTot = 0.;
1770  double liveZero = 0.;
1771 
1772  TStopwatch w;
1773  w.Start();
1774  cout << "Start CWB::CBCTool::getLiveTime of current ntuple : " << endl;
1775  double liveTot=cbcTool.getLiveTime2(liv);
1776  cout<<"Elapsed time: "<<w.RealTime()<<" s"<<endl;
1777  cout << "Start CWB::CBCTool::getZero(lag)LiveTime of current ntuple : " << endl;
1778  w.Start();
1779  double zerolagtime=cbcTool.getZeroLiveTime2(nIFO,liv);
1780  cout<<"Elapsed time: "<<w.RealTime()<<" s"<<endl;
1781 
1782  liveTot-=zerolagtime; //The live time calculated by CWB::CBCTool::getLiveTime2 is inclusive of zero-lag
1783 
1784  #endif
1785 
1786  double OLIVETIME_yr = liveZero/86400./365.;
1787  double BKG_LIVETIME_yr = liveTot/86400./365.;
1788  double OLIVETIME_Myr = OLIVETIME_yr/(1.e6);
1789  double BKG_LIVETIME_Myr = BKG_LIVETIME_yr/(1.e6);
1790 
1791  cout.precision(14);
1792  cout << "Total live time ---> zero lag: " << liveZero << " s, background: " << liveTot << " s" << endl;
1793  cout << "Total live time ---> zero lag: " << OLIVETIME_yr << " yr, background: " << BKG_LIVETIME_yr << " yr" << endl;
1794  cout << "Total live time ---> zero lag: " << OLIVETIME_Myr << " Myr, background: " << BKG_LIVETIME_Myr << " Myr" << endl;
1795 
1796 
1797 
1798 
1799 
1800  char bkg_threshold[256];
1801  char s1[64];
1802  char SEL[64];
1803  sprintf(s1,"!(lag[%i]==0&&slag[%i]==0)", nIFO, nIFO);
1804  #if defined(HVETO_VETO) || defined(CAT3_VETO)
1805  sprintf(bkg_threshold,"rho[%d] > %f && %s && %s",pp_irho,RHO_MIN,veto_not_vetoed,s1);
1806 #else
1807  sprintf(bkg_threshold,"rho[%d] > %f && %s && %s",pp_irho,RHO_MIN,ch2,s1);
1808 #endif
1809 
1810  sprintf(SEL,"rho[%d]",pp_irho);
1811  wave.Draw(SEL,bkg_threshold,"goff");
1812  int sel_events = wave.GetSelectedRows();
1813  const int bkg_entries = sel_events;
1814 
1815 
1816  double* rho_bkg=wave.GetV1();
1817 
1818 
1819  cout << "Threshold on background events:" << endl;
1820  cout << bkg_threshold << endl;
1821  cout << "Background events above thresholds: " << bkg_entries << endl;
1822 
1823  Int_t *index = new Int_t[bkg_entries];
1824 
1825  TMath::Sort(bkg_entries,rho_bkg,index);
1826  double RHO_MAX = ceil(rho_bkg[index[0]]*1.05);
1827  double VRHO[bkg_entries];
1828 for(int i = 0 ; i < bkg_entries ; i++) {VRHO[i] = rho_bkg[index[i]];}
1829 
1830  cout << "Background loudest event: " << VRHO[0] << ", largest considered rho: " << RHO_MAX << endl;
1831 
1832 
1833  #ifdef OPEN_LAG
1834  cout << "Opening lag: " << OPEN_LAG << endl;
1835  char s2[64];
1836 
1837  sprintf(s2,"(lag[%i]==%i&&slag[%i]==0)", nIFO, OPEN_LAG, nIFO);
1838 
1839  TString lag_cut = bkg_threshold;
1840  lag_cut.ReplaceAll(s1,s2); cout<<lag_cut.Data()<<endl;
1841 
1842  wave.Draw(SEL,lag_cut,"goff");
1843  int lag_events = wave.GetSelectedRows();
1844  const int fl_entries = lag_events;
1845 
1846  double* rho_lag = wave.GetV1();
1847 
1848  Int_t *indexl = new Int_t[fl_entries];
1849  TMath::Sort(fl_entries,rho_lag,indexl);
1850 
1851  cout << "Total entries in lag["<<OPEN_LAG <<"]: " << lag_events << endl;
1852 // cout << "Loudest event: " << rho_lag[indexl[0]] << endl
1853 
1854  double LFAD[fl_entries];
1855  double LRHO[fl_entries];
1856  double LMU[fl_entries];
1857 
1858  for(int i = 0 ; i < fl_entries ; i++) {LRHO[i] = rho_lag[indexl[i]]; LFAD[i] = 0.0;LMU[i] = 0.0;}
1859 #endif
1860 
1861 
1862 
1863 // double C = BKG_LIVETIME_Myr*shell_volume; //BEWARE!!! no volume_internal_sphere
1864  double C = OLIVETIME_Myr*shell_volume; //BEWARE!!! no volume_internal_sphere
1865 for(int i=0; i<RHO_NBINS; i++){
1866 
1867  Vrho[i]*=C;
1868  }
1869  cout<<"Mass averaged Visible volume @rho="<<RHO_MIN<<" : "<<Vrho[0]<<" on "<<massbins<< " mass bins" <<endl;
1870  cout<<"OLIVETIME_Myr : "<<OLIVETIME_Myr<<" shell_volume : "<<shell_volume<<endl;
1871 
1872 
1873 ///ROC Curve
1874 cbcTool.doROCPlot(bkg_entries, rho_bkg, index, RHO_BIN, liveTot, f2, c1, netdir, write_ascii);
1875 
1876 
1877 #ifdef FAD
1878 
1879 TGraph* grtmp = new TGraph(RHO_NBINS,Trho,Vrho);
1880  // f1->SetParameters(500.,-2.3);
1881  f2->SetParameters(500.,-2.5,0.0);
1882  grtmp->Fit("f2","R");
1883 
1884  grtmp->SetMarkerStyle(20);
1885  grtmp->SetMarkerSize(1.0);
1886  //grtmp->Draw("alp");
1887  TMultiGraph *multigraph = new TMultiGraph();
1888  multigraph->Add(grtmp);
1889  multigraph->Paint("AP");
1890  multigraph->GetHistogram()->GetXaxis()->SetTitle("#rho");
1891  multigraph->GetHistogram()->GetXaxis()->SetRangeUser(RHO_MIN,RHO_MAX);
1892  multigraph->GetHistogram()->GetYaxis()->SetRangeUser(f2->Eval(RHO_MAX),f2->Eval(RHO_MIN));
1893  multigraph->GetHistogram()->GetYaxis()->SetTitle("Productivity(Mpc^{3} Myr)");
1894  multigraph->GetXaxis()->SetTitleOffset(1.3);
1895  multigraph->GetYaxis()->SetTitleOffset(1.25);
1896  multigraph->GetXaxis()->SetTickLength(0.01);
1897  multigraph->GetYaxis()->SetTickLength(0.01);
1898  multigraph->GetXaxis()->CenterTitle(kTRUE);
1899  multigraph->GetYaxis()->CenterTitle(kTRUE);
1900  multigraph->GetXaxis()->SetTitleFont(42);
1901  multigraph->GetXaxis()->SetLabelFont(42);
1902  multigraph->GetYaxis()->SetTitleFont(42);
1903  multigraph->GetYaxis()->SetLabelFont(42);
1904  multigraph->GetYaxis()->SetMoreLogLabels(kTRUE);
1905  multigraph->GetYaxis()->SetNoExponent(kTRUE);
1906 
1907  c1->Clear();
1908  c1->SetLogy(1);
1909  gStyle->SetOptFit(kTRUE);
1910 
1911  multigraph->Draw("ALP");
1912  sprintf(radius_title,"%s : Productivity (%s) ", networkname, waveform.Data());
1913  p_radius->Clear();
1914  TText *text = p_radius->AddText(radius_title);
1915  p_radius->Draw();
1916 
1917 sprintf(fname,"%s/Productivity.png",netdir);
1918  c1->Update();
1919  c1->SaveAs(fname);
1920 
1921 
1922 
1923  double VFAD[bkg_entries];
1924  //double VRHO[bkg_entries];
1925  double MU[bkg_entries];
1926  double FAD = 0.0;
1927  double dFAD = 0.0;
1928  int cnt2 = 0;
1929  int lag_cnt = 0;
1930 
1931 double K = OLIVETIME_Myr/BKG_LIVETIME_Myr;
1932 for (int i = 0;i<bkg_entries;i++){
1933 
1934 //#ifdef OLD_FAD
1935  dFAD = 1.0/f2->Eval(VRHO[i]);
1936  FAD += dFAD;
1937  //FAD = (i+1)*dFAD;
1938 
1939  VFAD[i] = FAD*K;
1940  //VRHO[i] = rho_bkg[index[i]];
1941  MU[i] = VFAD[i]/dFAD ;
1942  // if ((++cnt2%500==0)||fabs(FAD-0.001) < 0.001)){cout << i << " Rho : " << rho_bkg[index[i]] << " Productivity : " << MU[i]/FAD <<" FAD : "<<FAD<<" MU : "<<MU[i]<< endl;}
1943 #ifdef OPEN_LAG
1944 
1945  while((VRHO[i]<=LRHO[lag_cnt])&&lag_cnt < fl_entries){
1946 
1947  // for(int j = -1;j<1;j++) { cout << i+j << " Rho : " << VRHO[i+j] << " Productivity : " << MU[i+j]/VFAD[i+j] <<" FAD : "<<VFAD[i+j]<<" MU : "<<MU[i+j]<< endl;}
1948  cout << i+j << " Rho : " << VRHO[i] << " Productivity : " << MU[i]/VFAD[i] <<" FAD : "<<VFAD[i]<<" MU : "<<MU[i]<< endl;
1949  if (i==0){cout << "BEWARE!!! Loudest lag event larger than loudest bkg event: "<< LRHO[lag_cnt] <<" "<<VRHO[0] <<endl;LFAD[lag_cnt]=FAD;LMU[lag_cnt] = MU[i];lag_cnt++;}
1950  else {
1951  #ifdef OLD_FAD
1952  LFAD[lag_cnt] = (VFAD[i]-VFAD[i-1])/(VRHO[i]-VRHO[i-1])*(LRHO[lag_cnt]-VRHO[i-1])+VFAD[i-1];
1953  LMU[lag_cnt] = (MU[i]-MU[i-1])/(VRHO[i]-VRHO[i-1])*(LRHO[lag_cnt]-VRHO[i-1])+MU[i-1];
1954  cout << "Lag["<< OPEN_LAG << "]: #"<< lag_cnt <<" rho=" << LRHO[lag_cnt] << " FAD="<<LFAD[lag_cnt] << " MU="<<LMU[lag_cnt]<<endl;
1955  lag_cnt++;
1956  #endif
1957  }
1958  }
1959 
1960 #endif
1961 
1962  //
1963  //RhoH->Fill(rho_bkg[index[i]],FAD);
1964 
1965 }
1966 
1967 cout << "Loop ends..."<<endl;
1968 
1969  c1->Clear();
1970  c1->SetLogy(1);
1971  gStyle->SetOptFit(kFALSE);
1972 TGraph* grFAD = new TGraph(bkg_entries,VRHO,VFAD);
1973 
1974 grFAD->GetHistogram()->GetXaxis()->SetRangeUser(RHO_MIN,RHO_MAX);
1975  // multigraph->GetHistogram()->GetYaxis()->SetRangeUser(0.01,10.0);
1976  grFAD->GetXaxis()->SetTitle("#rho");
1977  grFAD->GetYaxis()->SetTitle("FAD (Mpc^{-3} Myr^{-1})");
1978  grFAD->GetXaxis()->SetTitleOffset(1.3);
1979  grFAD->GetYaxis()->SetTitleOffset(1.25);
1980  grFAD->GetXaxis()->SetTickLength(0.01);
1981  grFAD->GetYaxis()->SetTickLength(0.01);
1982  grFAD->GetXaxis()->CenterTitle(kTRUE);
1983  grFAD->GetYaxis()->CenterTitle(kTRUE);
1984  grFAD->GetXaxis()->SetTitleFont(42);
1985  grFAD->GetXaxis()->SetLabelFont(42);
1986  grFAD->GetYaxis()->SetTitleFont(42);
1987  grFAD->GetYaxis()->SetLabelFont(42);
1988  grFAD->SetMarkerStyle(20);
1989  grFAD->SetMarkerSize(1.0);
1990  grFAD->SetMarkerColor(1);
1991  grFAD->SetLineColor(kBlue);
1992  grFAD->SetTitle("");
1993  grFAD->Draw("alp");
1994 
1995  #ifdef OPEN_LAG
1996  TGraph* gr_LAG = new TGraph(fl_entries,LRHO,LFAD);
1997  gr_LAG->SetMarkerStyle(20);
1998  gr_LAG->SetMarkerSize(1.0);
1999  gr_LAG->SetMarkerColor(2);
2000  gr_LAG->Draw("p,SAME");
2001  #endif
2002 
2003  char leg[128];
2004  sprintf(leg,"Background");
2005  leg_FAD = new TLegend(0.707831,0.735751,0.940763,0.897668,"","brNDC");
2006  leg_FAD->AddEntry(grFAD,leg,"p");
2007 
2008  #ifdef OPEN_LAG
2009 
2010  sprintf(leg,"Events from lag[%i]",OPEN_LAG);
2011  leg_FAD->AddEntry(gr_LAG,leg,"p");
2012  #endif
2013  leg_FAD->SetFillColor(0);
2014  leg_FAD->Draw();
2015 
2016 
2017 
2018  char FAD_title[256];
2019  sprintf(FAD_title,"%s : FAD distribution (%s) ", networkname, waveform.Data());
2020 
2021 
2022 
2023  TPaveText *title1 = new TPaveText(0.2941767,0.9274611,0.7359438,0.9974093,"blNDC");
2024  title1->SetBorderSize(0);
2025  title1->SetFillColor(0);
2026  title1->SetTextColor(1);
2027  title1->SetTextFont(32);
2028  title1->SetTextSize(0.045);
2029  TText *text = title1->AddText(FAD_title);
2030  title1->Draw();
2031 
2032 
2033 
2034  //RhoH->Draw("LP");
2035 sprintf(fname,"%s/FAD.eps",netdir);
2036  c1->Update();
2037  c1->SaveAs(fname);
2038 
2039 
2040  c1->Clear();
2041  c1->SetLogy(1);
2042  gStyle->SetOptFit(kFALSE);
2043 TGraph* grMU = new TGraph(bkg_entries,VRHO,MU);
2044 grMU->GetHistogram()->GetYaxis()->SetTitle("#mu");
2045 grMU->GetHistogram()->GetXaxis()->SetTitle("#rho");
2046 grMU->GetHistogram()->GetXaxis()->SetRangeUser(RHO_MIN,RHO_MAX);
2047 
2048 grMU->GetXaxis()->SetTitleOffset(1.3);
2049 grMU->GetYaxis()->SetTitleOffset(1.25);
2050 grMU->GetXaxis()->SetTickLength(0.01);
2051 grMU->GetYaxis()->SetTickLength(0.01);
2052 grMU->GetXaxis()->CenterTitle(kTRUE);
2053 grMU->GetYaxis()->CenterTitle(kTRUE);
2054 grMU->GetXaxis()->SetTitleFont(42);
2055 grMU->GetXaxis()->SetLabelFont(42);
2056 grMU->GetYaxis()->SetTitleFont(42);
2057 grMU->GetYaxis()->SetLabelFont(42);
2058 grMU->SetMarkerStyle(20);
2059 grMU->SetMarkerSize(1.0);
2060 grMU->SetMarkerColor(1);
2061 grMU->SetLineColor(kBlue);
2062 grMU->SetTitle("");
2063 grMU->Draw("alp");
2064 
2065  #ifdef OPEN_LAG
2066  TGraph* gr_LAG_MU = new TGraph(fl_entries,LRHO,LMU);
2067  gr_LAG_MU->SetMarkerStyle(20);
2068  gr_LAG_MU->SetMarkerSize(1.0);
2069  gr_LAG_MU->SetMarkerColor(2);
2070  gr_LAG_MU->Draw("p,SAME");
2071  #endif
2072 
2073 
2074  sprintf(FAD_title,"%s : Expected events per Observation time ", networkname);
2075 
2076 
2077  title1->Clear();
2078  TText *text = title1->AddText(FAD_title);
2079  title1->Draw();
2080  leg_FAD->Draw();
2081 sprintf(fname,"%s/MU.eps",netdir);
2082  c1->Update();
2083  c1->SaveAs(fname);
2084 
2085 #ifdef OPEN_LAG
2086  wave.Scan("rho[1]:time[0]:netcc[0]:run",lag_cut,"colsize=12");
2087  #endif
2088 
2089  #endif
2090 #endif
2091 
2092  exit(0);
2093 }
float spin[6]
Definition: cbc_plots.C:438
char ch2[2048]
#define NRGBs
TText * text
Definition: cbc_plots.C:717
double rho
char cut[512]
exit(1)
static const double C
Definition: GNGen.cc:10
char val[20]
TCut chirp
char xtitle[1024]
static double g(double e)
Definition: GNGen.cc:98
double T_pen
CWB::config * cfg
Definition: cbc_plots.C:106
char eff_title[256]
delete[] spin_mtot_volume
printf("total live time: non-zero lags = %10.1f \n", liveTot)
float min_mass1
Definition: cbc_plots.C:45
int nwave
Definition: cbc_plots.C:820
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
TString waveform
Definition: cbc_plots.C:115
char networkname[256]
Definition: cbc_plots.C:95
int point
float maxRatio
Definition: cbc_plots.C:113
bool bminDistance
Definition: cbc_plots.C:114
sprintf(inj_title,"Injected events")
TH2F * efficiency
TH2F * factor_events_rec
Definition: cbc_plots.C:266
par[0] name
TH1F * ecor
static double getZeroLiveTime2(int nIFO, TChain &liv)
Definition: CBCTool.cc:563
netevent W & wave
TString("c")
delete[] error_radius
TH2F * rec_events
Definition: cbc_plots.C:238
int NBINS_mass1
Definition: cbc_plots.C:62
#define RHO_NBINS
Definition: cbc_plots.C:4
double red[NRGBs]
Definition: cbc_plots.C:24
double T_cor
TExec * ex1
TLegend * leg
Definition: compare_bkg.C:246
Double_t z0
TPaveText * p_radius
Definition: cbc_plots.C:1368
delete[] radius
int mt
Definition: cbc_plots.C:434
ofstream finj
Int_t TotalConts
TH2F * D_Chi_injy
float mchirp
Definition: cbc_plots.C:436
CWB::mdc * MDC
TLatex Tl2
TH2F * efficiency_first_shell
TH1D * D_rec
TH1F * Dt
Definition: cbc_plots.C:398
static TF1 * doRangePlot(int RHO_NBINS, double *Trho, double *Rrho, double *eRrho, float RHO_MIN, float T_out, TCanvas *c1, TString networkname, TString odir, bool write_ascii)
Definition: CBCTool.cc:362
TH1D * Mchirp_rec
int j
Definition: cwb_net.C:10
static double getLiveTime2(TChain &liv)
Definition: CBCTool.cc:519
i drho i
network * net
Definition: cbc_plots.C:105
float min_mass2
Definition: cbc_plots.C:57
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:3956
int nwave_final
Definition: cbc_plots.C:838
CWB::Toolbox TB
Definition: ComputeSNR.C:5
delete[] volume_first_shell
char ifo[NIFO_MAX][8]
double green[NRGBs]
Definition: cbc_plots.C:25
TCut netcc("netcc","netcc[0]>0.7")
nT
Definition: cbc_plots.C:659
bool bminRatio
Definition: cbc_plots.C:114
char inj_title[256]
Definition: cbc_plots.C:704
float distance
Definition: cbc_plots.C:436
bool DDistrUniform
Definition: cbc_plots.C:114
#define NCont
TChain mdc("mdc")
dV
Temporay patch for redshifted mass distributions, i.e. point-like in source frame and spread over mul...
TChain sim("waveburst")
wavearray< double > w
Definition: Test1.C:27
#define nIFO
int countv
Definition: cbc_plots.C:587
bool bminMtot
Definition: cbc_plots.C:114
TH1D * Mtot_rec
float M1
int nmdc
Definition: cbc_plots.C:791
delete[] factor_events_spin_mtot_inj
TH1D * Mtot_inj
TH2F * rho_pf
Definition: cbc_plots.C:414
TH2F * D_Chi_recx
char fname[1024]
Definition: cbc_plots.C:720
double time[6]
Definition: cbc_plots.C:435
float shell_volume[nfactor]
int NBINS_mass
Definition: cbc_plots.C:33
TLatex Tl
double blue[NRGBs]
Definition: cbc_plots.C:26
TH2F * D_Mchirp_rec
Definition: cbc_plots.C:317
float MINCHIRP
Definition: cbc_plots.C:186
double pp_rho_min
float mass[2]
Definition: cbc_plots.C:437
bool DDistrVolume
Definition: cbc_plots.C:114
bool bmaxRatio
Definition: cbc_plots.C:114
Int_t nGraphs
i() int(T_cor *100))
double Pi
TF1 * f2
Definition: cbc_plots.C:1710
float max_mass2
Definition: cbc_plots.C:39
char netdir[1024]
float maxMtot
Definition: cbc_plots.C:113
char tmp_dir[1024]
Definition: config.hh:307
delete[] error_spin_mtot_radius
TH2F * D_Mchirp_inj
Definition: cbc_plots.C:296
float MAXCHIRP
Definition: cbc_plots.C:187
bool pp_rho_log
float maxDistance
Definition: cbc_plots.C:113
delete[] error_volume_first_shell
float minMtot
Definition: cbc_plots.C:113
std::vector< double > iSNR[NIFO_MAX]
TH2F * D_q_rec
Definition: cbc_plots.C:393
TExec * ex2
Definition: cbc_plots.C:1359
double liveTot
int k
m1
Definition: cbc_plots.C:606
leg_snr
Definition: cbc_plots.C:845
int m
Definition: cbc_plots.C:434
TH2F * h_radius
Definition: cbc_plots.C:1322
float chi2
Definition: cbc_plots.C:603
TH2F * htemp5
Definition: cbc_plots.C:910
float FACTORS[nfactor]
TH2F * D_Chi_rec
double pp_rho_max
int NBINS_mass2
Definition: cbc_plots.C:63
delete[] error_spin_mtot_volume
char net_file_name[256]
char fname3[2048]
char liv_file_name[256]
char sim_file_name[1024]
float max_mass1
Definition: cbc_plots.C:51
TH2F * inj_events
Definition: cbc_plots.C:209
delete[] spin_mtot_radius
Double_t x0
double stops[NRGBs]
Definition: cbc_plots.C:23
error_volume[m][l]
Definition: cbc_plots.C:679
double liveZero
TH1D * Mchirp_inj
static void doROCPlot(int bkg_entries, double *rho_bkg, int *index, float RHO_BIN, double liveTot, TF1 *AverageRad, TCanvas *c1, TString odir, bool write_ascii)
Definition: CBCTool.cc:88
char lab[256]
Definition: cbc_plots.C:844
char title[256]
Definition: SSeriesExample.C:1
float MAXRATIO
Definition: cbc_plots.C:151
TH2F * factor_events_inj[nfactor]
Definition: cbc_plots.C:262
float factor
Definition: cbc_plots.C:436
float MINRATIO
Definition: cbc_plots.C:143
TH1D * q_inj
int NBINS_SPIN
TH2F * D_Mtot_rec
Definition: cbc_plots.C:290
wavearray< int > index
double T_win
TH1D * q_rec
TH1D * D_inj
TH2F * rhocc
Definition: cbc_plots.C:404
double fabs(const Complex &x)
Definition: numpy.cc:37
TString sel("slag[1]:rho[1]")
TH2F * D_Mtot_inj
Definition: cbc_plots.C:269
static void mkDir(TString dir, bool question=false, bool remove=true)
Definition: Toolbox.cc:4000
int cz
Definition: cbc_plots.C:434
char radius_title[256]
Definition: cbc_plots.C:1363
int cnt
strcpy(cfg->tmp_dir,"tmp")
int l
Definition: cbc_plots.C:434
nfactor[0]
Definition: cwb_eced.C:10
int MAX_AXIS_Z
Definition: cbc_plots.C:700
char veto_not_vetoed[1024]
TCanvas * c1
Definition: cbc_plots.C:197
float minDistance
Definition: cbc_plots.C:113
float M2
TH2F * D_Chi_injx
float minRatio
Definition: cbc_plots.C:113
TMacro configPlugin
double T_out
fclose(finj)
bool write_ascii
Definition: cbc_plots.C:447
bool bmaxDistance
Definition: cbc_plots.C:114
int massbins
Definition: cbc_plots.C:1299
#define RHO_BIN
Definition: cbc_plots.C:3
Int_t GetEntries()
Definition: netevent.cc:381
double T_vED
TString Insp
Definition: cbc_plots.C:111
char mdc_file_name[1024]
#define CONTOURS
Definition: cbc_plots.C:5
TString config
detectorParams detParms[4]
double Start
Definition: wavearray.hh:307
float chi[3]
Definition: cbc_plots.C:439
volume[m][l]
Definition: cbc_plots.C:678
#define RHO_MIN
Definition: cbc_plots.C:2
TH2F * D_q_inj
Definition: cbc_plots.C:372
TPaveText * p_inj
Definition: cbc_plots.C:711
TH2F * D_Chi_recy
i drho pp_irho
bool bmaxMtot
Definition: cbc_plots.C:114
TPaveText * p_eff
TH2F * htemp4
Definition: cbc_plots.C:874
TString Compare_tree
Definition: compare_bkg.C:99
m2
Definition: cbc_plots.C:607
double V0
Definition: cbc_plots.C:1300