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