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