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