Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CBCTool.cc
Go to the documentation of this file.
1 #include "CBCTool.hh"
2 #include "watversion.hh"
3 #include "TTreeFormula.h"
4 #include "TStopwatch.h"
5 //#include "GCBCTool.hh"
6 #include "TThread.h"
7 #include <Riostream.h>
8 #include "TF1.h"
9 #include "Math/WrappedTF1.h"
10 #include "Math/BrentRootFinder.h"
11 
12 #define CCAST(PAR) const_cast<char*>(PAR) \
13 
14 // Without this macro the THtml doc for CWB::CBCTool can not be generated
15 ClassImp(CWB::CBCTool)
16 
17 //int compareSegments(waveSegment a, waveSegment b) {return a.start < b.start;}
18 
19 // definitions used by CWB::CBCTool::mergeCWBTrees with threads
20 
21 #define MAX_THREADS 8
22 
23 ////______________________________________________________________________________
24 
25 Long64_t CWB::CBCTool::GetTreeIndex(TTree *tree, const char* selection){
26 
27  tree->Draw("Entry$>>hist(Entries$,0,Entries$)",selection,"goff");
28  TH1I *hist = (TH1I*)gDirectory->Get("hist");
29  Long64_t iEntry = hist->GetBinLowEdge(hist->FindFirstBinAbove(0));
30  delete hist;
31  return iEntry;
32 }
33 
34 //______________________________________________________________________________
35 /*
36 TH1* CWB::CBCTool::GetCumulative(TH1* in, Bool_t forward)
37 {
38  // Return a pointer to an histogram containing the cumulative The
39  // cumulative can be computed both in the forward (default) or backward
40  // direction; the name of the new histogram is constructed from
41  // the name of this histogram with the suffix suffix appended.
42  //
43  // The cumulative distribution is formed by filling each bin of the
44  // resulting histogram with the sum of that bin and all previous
45  // (forward == kTRUE) or following (forward = kFALSE) bins.
46  //
47  // note: while cumulative distributions make sense in one dimension, you
48  // may not be getting what you expect in more than 1D because the concept
49  // of a cumulative distribution is much trickier to define; make sure you
50  // understand the order of summation before you use this method with
51  // histograms of dimension >= 2.
52 
53  const Int_t nbinsx = in->GetNbinsX();
54  const Int_t nbinsy = in->GetNbinsY();
55  const Int_t nbinsz = in->GetNbinsZ();
56  TH1* hintegrated = (TH1*) in->Clone();
57  hintegrated->Reset();
58  if (forward) { // Forward computation
59  Double_t sum = 0.;
60  for (Int_t binz = 1; binz <= nbinsz; ++binz) {
61  for (Int_t biny = 1; biny <= nbinsy; ++biny) {
62  for (Int_t binx = 1; binx <= nbinsx; ++binx) {
63  const Int_t bin = hintegrated->GetBin(binx, biny, binz);
64  sum += in->GetBinContent(bin);
65  hintegrated->SetBinContent(bin, sum);
66  }
67  }
68  }
69  } else { // Backward computation
70  Double_t sum = 0.;
71  for (Int_t binz = nbinsz; binz >= 1; --binz) {
72  for (Int_t biny = nbinsy; biny >= 1; --biny) {
73  for (Int_t binx = nbinsx; binx >= 1; --binx) {
74  const Int_t bin = hintegrated->GetBin(binx, biny, binz);
75  sum += in->GetBinContent(bin);
76  hintegrated->SetBinContent(bin, sum);
77  }
78  }
79  }
80  }
81  return hintegrated;
82 }
83 
84 */
85 
86 //______________________________________________________________________________
87 
88 void CWB::CBCTool::doROCPlot(int bkg_entries, double* rho_bkg, int* index, float RHO_BIN, double liveTot, TF1* AverageRad, TCanvas* c1, TString odir, bool write_ascii) {
89 
90 //
91 // It produces ROC Plot from background and simulation ntuples
92 //
93 // Input
94 // bkg_entries : number of background events
95 // rho_bkg : array of doubles containing the events rho
96 // index : array of integers containing the sorted entries in ascending order // MY: too much useless stuff, need some clening and rearrangement!
97 // RHO_BIN : float containing the rho bin //MY: remove it? it should be defined in user-pp conf file
98 // liveTot : double containing the total background live time in seconds as calculated from the liveTime ntuple
99 // AverageRad : fitting function of the measured average radius as a function of rho (a power-law with some exponential, A*rho^B*exp(-C))
100 // calculated by the following doRangePlot method
101 // BEWARE, the fit does a decent job on all sofar tested waveforms, but no check is done! Look at the Range.png.
102 // The reasoning in favour of using the radius from fitting in place of the measured one (which produces some overestimate bias at low rho)
103 // is to have a better estimate at high rho, where the statistics is poor.
104 // c1 : the formatted canvas used for all plots by cbc_plots.C
105 // odir : output dir
106 // write_ascii : boolean to toggle on/off the writing of ascii output
107 
108 
109  int nbins = -TMath::Nint(rho_bkg[index[bkg_entries-1]]-rho_bkg[index[0]])/RHO_BIN;
110  cout << "Rho max : "<< rho_bkg[index[0]] << " Rho min : "<<rho_bkg[index[bkg_entries-1]] << " nbins : "<<nbins<<endl;
111 
112  //definition and filling of the events rho histogram
113  TH1D* h = new TH1D("h", "h", nbins,rho_bkg[index[bkg_entries-1]],rho_bkg[index[0]]*1.05);
114  for (int i = 0;i<bkg_entries;i++){h->Fill(rho_bkg[i]);}
115 //hc is the cumulative histogram of h
116 // TH1* hc = h->GetCumulative(kFALSE);
117  const Int_t nbinsx = h->GetNbinsX();
118  const Int_t nbinsy = h->GetNbinsY();
119  const Int_t nbinsz = h->GetNbinsZ();
120  TH1* hc = (TH1*) h->Clone();
121  hc->Reset();
122  Double_t sum = 0.;
123  for (Int_t binz = nbinsz; binz >= 1; --binz) {
124  for (Int_t biny = nbinsy; biny >= 1; --biny) {
125  for (Int_t binx = nbinsx; binx >= 1; --binx) {
126  const Int_t bin = hc->GetBin(binx, biny, binz);
127  sum += h->GetBinContent(bin);
128  hc->SetBinContent(bin, sum);
129  }
130  }
131  }
132 
133 
134  double far[nbins],efar[nbins],rad[nbins],rhobin[nbins],erad[nbins];
135  for (int i = 0; i<nbins; i++) {
136  far[i] = (double)hc->GetBinContent(i)/liveTot;
137  efar[i] = (double)hc->GetBinError(i)/liveTot;
138  rhobin[i] = hc->GetBinCenter(i); //MY : should I consider the left bound of the bin? this choice should be more conservative
139  rad[i] = AverageRad->Eval(rhobin[i] ); //BEWARE: no guarantee that the fitting function returns correct estimates
140  erad[i] = 0.0;
141 // cout << "(far, rho, radius) - "<<far[i]<<","<<hc->GetBinCenter(i)<<","<<rad[i]<<" ";
142  }
143  //cout<<endl;
144 // cout << "FAR max : "<< far[0] << " FAR min : "<<far[nbins-1] <<endl;
145  c1->Clear();
146  c1->SetLogy(1);
147  c1->SetLogx(1);
148 
149  hc->Scale(1./liveTot);
150  hc->GetYaxis()->SetTitle("FAR [Hz]");
151  hc->GetXaxis()->SetTitle("#rho");
152  hc->Draw("LP");
153  char fname[1024];
154  sprintf(fname,"%s/FAR.png",odir.Data());
155  c1->Update();
156  c1->SaveAs(fname); //A dump of the FAR used for the ROC
157 
158  c1->Clear();
159  c1->SetLogx(1);
160  c1->SetLogy(0); //MY: log?
161  c1->SetTheta(89.999);// the setting of the theta and phi is for the viewing angle: there is no TGraphError with "pcol" option, so the trick is to use the 2D version and use a viewing angle the closest as possible to the vertical (Drawbacks: some artefacts due to that can be noticed; no gridlines)
162  c1->SetPhi(0.0001);
163 
164  TGraph2DErrors* gr2=new TGraph2DErrors(nbins,far,rad,rhobin,efar,erad,erad);
165  gr2->SetMarkerStyle(20);
166  gr2->SetLineColor(20);
167  gr2->SetMinimum(0.0);
168  double RHO_MAX = ceil(rho_bkg[index[0]]*1.05);
169  gr2->SetMaximum(RHO_MAX);
170  // gr1->Draw("ALP");
171  //TMultiGraph *multi = new TMultiGraph();
172  // multi->Add(gr2);
173  // multi->Paint("AP");
174  gr2->SetTitle("ROC Curve: Average Range vs FAR @rho threshold");
175  gr2->GetHistogram()->GetXaxis()->SetTitle("FAR [Hz]");
176  gr2->GetHistogram()->GetZaxis()->SetTitle("#rho");
177 // multi->GetHistogram()->GetXaxis()->SetRangeUser(RHO_MIN,10.);
178 // multi->GetHistogram()->GetXaxis()->SetRangeUser(1e-12, 1e-7);
179  gr2->GetHistogram()->GetXaxis()->SetLimits(1e-12, 1e-7); //MY: Add similar line for the y axis? declare the X and Y bounds in the definition
180  gr2->GetHistogram()->GetYaxis()->SetTitle("Average Range [Mpc]");
181  gr2->GetXaxis()->SetTitleOffset(1.3);
182  gr2->GetYaxis()->SetTitleOffset(1.25);
183  gr2->GetXaxis()->SetTickLength(0.02);
184  gr2->GetYaxis()->SetTickLength(0.01);
185  gr2->GetXaxis()->CenterTitle(kTRUE);
186  gr2->GetYaxis()->CenterTitle(kTRUE);
187  gr2->GetZaxis()->CenterTitle(kTRUE);
188  gr2->GetXaxis()->SetTitleFont(42);
189  gr2->GetXaxis()->SetLabelFont(42);
190  gr2->GetYaxis()->SetTitleFont(42);
191  gr2->GetYaxis()->SetLabelFont(42);
192 // gr2->GetYaxis()->SetMoreLogLabels(kTRUE);
193  gr2->GetYaxis()->SetNoExponent(kTRUE);
194 
195  gr2->Draw("err zcolpcol");
196  sprintf(fname,"%s/ROC.png",odir.Data());
197  c1->Update();
198  c1->SaveAs(fname);
199 
200  if(write_ascii){
201  sprintf(fname,"%s/ROC.txt",odir.Data());
202  FILE* frange = fopen(fname,"w");
203  fprintf(frange,"#rho FAR[Hz] eFAR range[Mpc] erange\n");
204  for(int i=0; i<nbins;i++){fprintf(frange,"%3.3g %4.3g %4.3g %4.4g %4.4g\n",hc->GetBinCenter(i),far[i],efar[i],rad[i],erad[i]);}
205  fclose(frange);
206 
207  }
208 
209  delete h, hc,gr2;
210 
211  return;
212 }
213 
214 
215 //______________________________________________________________________________
216 
217 void CWB::CBCTool::doROCPlot(int bkg_entries, double* rho_bkg, int* index, float RHO_BIN, float RHO_NBINS, float RHO_MIN, double liveTot, double* Rrho, double* eRrho, double* Trho, TCanvas* c1, TString odir, bool write_ascii) {
218 
219 //
220 // It produces ROC Plot from background and simulation ntuples
221 //
222 // Input
223 // bkg_entries : number of background events
224 // rho_bkg : array of doubles containing the events rho
225 // index : array of integers containing the sorted entries in ascending order // MY: too much useless stuff, need some clening and rearrangement!
226 // RHO_BIN : float containing the rho bin //MY: remove it? it should be defined in user-pp conf file
227 // liveTot : double containing the total background live time in seconds as calculated from the liveTime ntuple
228 // Rrho and eRrho : arrays of doubles containing the Radius and error radius estimates at fixed rho thresholds Trho
229 //
230 // c1 : the formatted canvas used for all plots by cbc_plots.C
231 // odir : output dir
232 // write_ascii : boolean to toggle on/off the writing of ascii output
233 
234 
235  //int nbins = -TMath::Nint(rho_bkg[index[bkg_entries-1]]-rho_bkg[index[0]])/RHO_BIN;
236  int nbins = (int)RHO_NBINS;
237  cout << "Rho max : "<< rho_bkg[index[0]] << " Rho min : "<<rho_bkg[index[bkg_entries-1]] << " nbins : "<<nbins<<endl;
238 
239  //definition and filling of the events rho histogram
240  TH1D* h = new TH1D("h", "h", nbins,RHO_MIN,RHO_NBINS*RHO_BIN);
241  for (int i = 0;i<bkg_entries;i++){h->Fill(rho_bkg[i]);}
242 //hc is the cumulative histogram of h
243 // TH1* hc = h->GetCumulative(kFALSE);
244 
245  const Int_t nbinsx = h->GetNbinsX();
246  const Int_t nbinsy = h->GetNbinsY();
247  const Int_t nbinsz = h->GetNbinsZ();
248  TH1* hc = (TH1*) h->Clone();
249  hc->Reset();
250  Double_t sum = 0.;
251  for (Int_t binz = nbinsz; binz >= 1; --binz) {
252  for (Int_t biny = nbinsy; biny >= 1; --biny) {
253  for (Int_t binx = nbinsx; binx >= 1; --binx) {
254  const Int_t bin = hc->GetBin(binx, biny, binz);
255  sum += h->GetBinContent(bin);
256  hc->SetBinContent(bin, sum);
257  }
258  }
259  }
260 
261  double far[nbins],efar[nbins],rad[nbins],rhobin[nbins],erad[nbins], vol[nbins], evol[nbins], erhobin[nbins];
262  int j=0;
263  for (int i = 0; i<nbins; i++) {
264  far[i] = (double)hc->GetBinContent(i+1)/liveTot*365.*86400.;
265  efar[i] = (double)hc->GetBinError(i+1)/liveTot*365.*86400.;
266  rhobin[i] = hc->GetBinCenter(i+1); //MY : should I consider the left bound of the bin? this choice should be more conservative
267  while(rhobin[i]>Trho[j]) {j++;}
268  erhobin[i] = RHO_BIN;
269  rad[i] = Rrho[j];
270  erad[i] = eRrho[j];
271  vol[i] = 4./3.*TMath::Pi()*pow(Rrho[j],3.);
272  evol[i] = 4.*TMath::Pi()*pow(Rrho[j],2.)*eRrho[j];
273 // cout << "(far, rho, radius) - "<<far[i]<<","<<hc->GetBinCenter(i)<<","<<rad[i]<<" ";
274  }
275  //cout<<endl;
276 // cout << "FAR max : "<< far[0] << " FAR min : "<<far[nbins-1] <<endl;
277  c1->Clear();
278  c1->SetLogy(1);
279  c1->SetLogx(1);
280 
281  hc->Scale(365.*86400./liveTot);
282  hc->GetYaxis()->SetTitle("FAR [year^{-1}]");
283  hc->GetXaxis()->SetTitle("#rho");
284  hc->Draw("LP");
285  char fname[1024];
286  sprintf(fname,"%s/FAR.png",odir.Data());
287  c1->Update();
288  c1->SaveAs(fname); //A dump of the FAR used for the ROC
289 
290  c1->Clear();
291  c1->SetLogx(1);
292  c1->SetLogy(0); //MY: log?
293  c1->SetTheta(89.999);// the setting of the theta and phi is for the viewing angle: there is no TGraphError with "pcol" option, so the trick is to use the 2D version and use a viewing angle the closest as possible to the vertical (Drawbacks: some artefacts due to that can be noticed; no gridlines)
294  c1->SetPhi(0.0001);
295 
296  TGraph2DErrors* gr2=new TGraph2DErrors(nbins,far,rad,rhobin,efar,erad,erhobin);
297  gr2->SetMarkerStyle(20);
298  gr2->SetLineColor(20);
299  gr2->SetMinimum(0.0);
300  double RHO_MAX = ceil(rho_bkg[index[0]]*1.05);
301  gr2->SetMaximum(RHO_MAX);
302  // gr1->Draw("ALP");
303  //TMultiGraph *multi = new TMultiGraph();
304  // multi->Add(gr2);
305  // multi->Paint("AP");
306  gr2->SetTitle("ROC Curve: Average Range vs FAR @rho threshold");
307  gr2->GetHistogram()->GetXaxis()->SetTitle("FAR [year^{-1}]");
308  gr2->GetHistogram()->GetZaxis()->SetTitle("#rho");
309 // multi->GetHistogram()->GetXaxis()->SetRangeUser(RHO_MIN,10.);
310 // multi->GetHistogram()->GetXaxis()->SetRangeUser(1e-12, 1e-7);
311  gr2->GetHistogram()->GetXaxis()->SetLimits(1e-5, 1e1); //MY: Add similar line for the y axis? declare the X and Y bounds in the definition
312  gr2->GetHistogram()->GetYaxis()->SetTitle("Average Range [Mpc]");
313  gr2->GetXaxis()->SetTitleOffset(1.3);
314  gr2->GetYaxis()->SetTitleOffset(1.25);
315  gr2->GetXaxis()->SetTickLength(0.02);
316  gr2->GetYaxis()->SetTickLength(0.01);
317  gr2->GetXaxis()->CenterTitle(kTRUE);
318  gr2->GetYaxis()->CenterTitle(kTRUE);
319  gr2->GetZaxis()->CenterTitle(kTRUE);
320  gr2->GetXaxis()->SetTitleFont(42);
321  gr2->GetXaxis()->SetLabelFont(42);
322  gr2->GetYaxis()->SetTitleFont(42);
323  gr2->GetYaxis()->SetLabelFont(42);
324 // gr2->GetYaxis()->SetMoreLogLabels(kTRUE);
325  gr2->GetYaxis()->SetNoExponent(kTRUE);
326 
327  gr2->Draw("err zcolpcol");
328  sprintf(fname,"%s/ROC.png",odir.Data());
329  c1->Update();
330  c1->SaveAs(fname);
331 
332  if(write_ascii){
333  sprintf(fname,"%s/ROC.txt",odir.Data());
334  FILE* frange = fopen(fname,"w");
335  fprintf(frange,"#rho FAR[year^{-1}] eFAR range[Mpc] erange\n");
336  for(int i=0; i<nbins;i++){fprintf(frange,"%3.3g %4.3g %4.3g %4.4g %4.4g\n",hc->GetBinCenter(i),far[i],efar[i],rad[i],erad[i]);}
337  fclose(frange);
338 
339  }
340 
341  TGraph2DErrors* gr3=new TGraph2DErrors(nbins,far,vol,rhobin,efar,erad,erhobin);
342  gr3->SetTitle("ROC Curve: Average Volume vs FAR @rho threshold");
343  gr3->GetHistogram()->GetYaxis()->SetTitle("Average volume [Mpc^3]");
344  gr3->GetHistogram()->GetXaxis()->SetTitle("FAR [year^{-1}]");
345  gr3->GetHistogram()->GetZaxis()->SetTitle("#rho");
346  gr3->GetHistogram()->GetXaxis()->SetLimits(1e-5, 1e2); //MY: Add similar line for the y axis? declare the X and Y bounds in the definition
347 
348  gr3->Draw("err zcolpcol");
349  sprintf(fname,"%s/ROCV.png",odir.Data());
350  c1->Update();
351  c1->SaveAs(fname);
352 
353  delete h, hc,gr2;
354 
355  return;
356 }
357 
358 
359 
360 //______________________________________________________________________________
361 
362 TF1* CWB::CBCTool::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){
363 
364 //
365 // It produces Range Plot from simulation ntuple
366 // AverageRad : fitting function of the measured average radius as a function of rho (a power-law with some exponential, A*rho^B*exp(-C))
367 // BEWARE, the fit does a decent job on all sofar tested waveforms, but no check is done! Look at the Range.png.
368 // The reasoning in favour of using the radius from fitting in place of the measured one (which produces some overestimate bias at low rho)
369 // is to have a better estimate at high rho, where the statistics is poor.
370 //
371 // Input
372 // RHO_NBINS : number of rho bins // MY: remove this in favour of the standard user_pp variables?
373 // Trho : array of doubles containg the rho thresholds at the left margin of the bins
374 // Rrho : array of doubles containing the average radius as calculated at the left margin of the bin
375 // eRrho : array of doubles containing the error on the average radius
376 // RHO_MIN : float with minimal rho // MY: do I really need this?
377 // T_out : float with output threashold
378 // c1 : the formatted canvas used for all plots by cbc_plots.C
379 // networkname : TString with the network name //MY: read it from the ntuple?
380 // odir : output dir
381 // write_ascii : boolean to toggle on/off the writing of ascii output
382 
383 
384  char fname[1024];
385 
386  if(write_ascii){
387  sprintf(fname,"%s/range_threshold.txt",odir.Data());
388  FILE* frange = fopen(fname,"w");
389  fprintf(frange,"#rho range[Mpc] \n");
390  for(int i=0; i<RHO_NBINS; i++){fprintf(frange,"%3.2f %4.3e\n",Trho[i],Rrho[i]);}
391  fclose(frange);
392  }
393 
394 
395 //TF1 *f1 = new TF1("f1","[0]*pow(x,[1])",T_out,10.);
396 TF1 *f2 = new TF1("f2","[0]*pow(x,[1])*TMath::Exp(-x*[2])",T_out,Trho[RHO_NBINS-1]); //Empirical function to fit the radius vs rho
397 //TF1 *f3 = new TF1("f3","[0]+[1]/pow(x,1)",T_out,10.);
398 f2->SetParameters(500.,-1.,0.0);
399 
400  TGraphErrors* grtmp = new TGraphErrors(RHO_NBINS,Trho,Rrho,NULL,eRrho);
401  grtmp->SetMarkerStyle(20);
402  grtmp->SetMarkerSize(1.0);
403  grtmp->GetHistogram()->GetXaxis()->SetTitle("#rho");
404  grtmp->GetHistogram()->GetXaxis()->SetRangeUser(RHO_MIN,Trho[RHO_NBINS-1]);
405  grtmp->GetHistogram()->GetYaxis()->SetRangeUser(f2->Eval(Trho[RHO_NBINS-1]),f2->Eval(T_out));
406  grtmp->GetHistogram()->GetYaxis()->SetTitle("Average Range [Mpc]");
407  grtmp->GetXaxis()->SetTitleOffset(1.3);
408  grtmp->GetYaxis()->SetTitleOffset(1.25);
409  grtmp->GetXaxis()->SetTickLength(0.01);
410  grtmp->GetYaxis()->SetTickLength(0.01);
411  grtmp->GetXaxis()->CenterTitle(kTRUE);
412  grtmp->GetYaxis()->CenterTitle(kTRUE);
413  grtmp->GetXaxis()->SetTitleFont(42);
414  grtmp->GetXaxis()->SetLabelFont(42);
415  grtmp->GetYaxis()->SetTitleFont(42);
416  grtmp->GetYaxis()->SetLabelFont(42);
417  grtmp->GetYaxis()->SetMoreLogLabels(kTRUE);
418  grtmp->GetYaxis()->SetNoExponent(kTRUE);
419 
420 
421  grtmp->Fit("f2","R");
422 
423 
424  c1->Clear();
425  c1->SetLogy();
426  c1->SetLogx();
427  gStyle->SetOptFit(kTRUE);
428  sprintf(fname,"%s : Range (%s) ", networkname.Data(), f2->GetExpFormula().Data());
429  grtmp->SetTitle(fname);
430  grtmp->Draw("ALP");
431 
432 
433  sprintf(fname,"%s/Range.png",odir.Data());
434  c1->Update();
435  c1->SaveAs(fname);
436 
437  delete grtmp;
438 
439  return f2;
440 }
441 
442 //______________________________________________________________________________
443 
444 double CWB::CBCTool::getLiveTime(int nIFO, TChain& liv,int lag_number, int slag_number, int dummy) {
445 //
446 // It calculates the total live time of all time shifts.
447 // If defining a lag_number/slag_number, it calculates live time for these only.
448 // This is a semplified version of the method (same name) that can be found in the standard toolbox: it's much faster as it does less things...
449 // Input: nIFO : detector number
450 // liv : tree containing live time information
451 // lag_number : specific lag number
452 // slag_number : specific slag number
453 // dummy : used to increase calculation speed
454 //
455 // Examples : (lag_number==-1) && (slag_number==-1) -> all non zero lags : excluded (lag==0&&slag==0)
456 // (lag_number>=0) && (slag_number>=0) -> select (lag==lag_number) && (slag==slag_number)
457 // (lag_number>=0) && (slag_number==-1) -> select lag=lag_number : excluded (lag==0&&slag==0)
458 // (lag_number==-1) && (slag_number>=0) -> select slag=slag_number : excluded (lag==0&&slag==0)
459 //
460 // Note: the variable dummy has been introduced to optimize
461 // speed of the method (the reason is unknown!!!)
462 
463  float xlag[NIFO_MAX+1];
464  float xslag[NIFO_MAX+1];
465  double xlive;
466  double liveTot = 0.;
467  double Live=0.;
468 
469  // check if slag is presents in liv tree
470  TBranch* branch;
471  bool slagFound=false;
472  TIter next(liv.GetListOfBranches());
473  while ((branch=(TBranch*)next())) {
474  if (TString("slag").CompareTo(branch->GetName())==0) slagFound=true;
475  }
476  next.Reset();
477  if(!slagFound) {
478  cout << "CWB::Toolbox::getLiveTime : Error - live tree do not contains slag leaf" << endl;
479  gSystem->Exit(1);
480  }
481 
482 // liv.SetCacheSize(10000000);
483 // liv.AddBranchToCache("slag");
484  // liv.AddBranchToCache("lag");
485  // liv.AddBranchToCache("live");
486  liv.SetBranchAddress("slag",xslag);
487  liv.SetBranchAddress("lag",xlag);
488 
489  liv.SetBranchAddress("live",&xlive);
490  liv.SetBranchStatus("gps",false);
491 
492 
493 Long64_t ntrg = liv.GetEntries();
494  for(Long64_t i=0; i<ntrg; i++) {
495  liv.GetEntry(i);
496 
497  if((lag_number>=0)&&(slag_number>=0)) {
498  Live = (xlag[nIFO]==lag_number&&xslag[nIFO]==slag_number) ? xlive : 0.; // lag/slag live time
499  }
500  if((lag_number>=0)&&(slag_number==-1)) {
501  Live = ((xlag[nIFO]==lag_number)&&!((xlag[nIFO]==0)&&(xslag[nIFO]==0))) ? xlive : 0.; // lag live time
502  }
503  if((lag_number==-1)&&(slag_number>=0)) {
504  Live = ((xslag[nIFO]==slag_number)&&!((xlag[nIFO]==0)&&(xslag[nIFO]==0))) ? xlive : 0.; // slag live time
505  }
506  if((lag_number==-1)&&(slag_number==-1)) {
507  Live = !((xlag[nIFO]==0)&&(xslag[nIFO]==0)) ? xlive : 0.; // non-zero live time
508  }
509 
510  liveTot += Live;
511  }
512  return liveTot;
513 
514 
515 }
516 
517 //______________________________________________________________________________
518 
519 double CWB::CBCTool::getLiveTime2(TChain& liv){
520 // It calculates the total live time of all time shifts (both zero lag and non-zero lags) : remember to subtract the zero-lag livetime if you need precise times!
521 // This is a semplified version of the previous method: it only reads&sums the live branch...faster than the speed of light...:) (very useful for large liveTime ntuples)
522 // Input:
523 // liv : tree containing live time information
524 
525  double xlive;
526  double liveTot = 0.;
527  TBranch* LIVE = liv.GetBranch("live");
528  LIVE->SetAddress(&xlive);
529  Long64_t nentries = LIVE->GetEntries();
530  for (Long64_t i=0;i<nentries;i++) {
531  LIVE->GetEntry(i);
532  liveTot += xlive;
533  }
534  return liveTot;
535 
536 }
537 //______________________________________________________________________________
538 
539 double CWB::CBCTool::getZeroLiveTime(int nIFO, TChain& liv){
540 // It calculates the total live time of zero lag
541 // This is mostly needed to correct the background estimate from the previous method
542 
543 // Input: nIFO : detector number
544 // liv : tree containing live time information
545  double liveTot = 0.;
546  char cut[256];
547 
548  sprintf(cut,"(lag[%d]==0&&slag[%d]==0)",nIFO,nIFO);
549  liv.SetEstimate(-1);
550  Long64_t nentries = liv.Draw("live",cut,"goff");
551  double* xlive = liv.GetV1();
552  for (Long64_t i=0;i<nentries;i++) {
553  liveTot += xlive[i];
554  }
555  return liveTot;
556 
557 }
558 
559 //______________________________________________________________________________
560 
561 //______________________________________________________________________________
562 
563 double CWB::CBCTool::getZeroLiveTime2(int nIFO, TChain& liv){
564 // It calculates the total live time of zero lag
565 // This is mostly needed to correct the background estimate from the previous method. Yet another emplementation to improve speed wrt the previous method (best so far) : it reads from branches just when it is really needed.
566 //
567 
568 // Input: nIFO : detector number
569 // liv : tree containing live time information
570 
571  float xlag[NIFO_MAX+1];
572  float xslag[NIFO_MAX+1];
573  double xlive;
574  double liveTot = 0.;
575 
576  Long64_t nentries = liv.GetEntries();
577  TBranch *lag = liv.GetBranch("lag");
578  lag->SetAddress(xlag);
579  TBranch *slag = liv.GetBranch("slag");
580  slag->SetAddress(xslag);
581  TBranch *live = liv.GetBranch("live");
582  live->SetAddress(&xlive);
583 
584  for (Long64_t i=0;i<nentries;i++) {
585  // T->LoadTree(i);
586  slag->GetEntry(i);
587  if (xslag[nIFO]!=0) continue;
588  lag->GetEntry(i);
589  if (xlag[nIFO]!=0) continue;
590  live->GetEntry(i);
591  liveTot += xlive;
592  }
593  return liveTot;
594 
595 
596 }
597 
598 //______________________________________________________________________________
599 
600 //______________________________________________________________________________
601 
602  double CWB::CBCTool::getFAR(float rho, TH1* hc, double liveTot){
603  int w=0;
604  while((rho>(float)hc->GetBinCenter(w))&&(w<hc->GetNbinsX()-1)){w++;}
605  double far = (double)hc->GetBinContent(w+1)/liveTot;
606  double efar = (double)hc->GetBinError(w+1)/liveTot;
607 
608  return far,efar;
609 }
610 
611 //______________________________________________________________________________
612 
613  TH2F* CWB::CBCTool::FillSLagHist(int NIFO_MAX, TChain& live, int NSlag, double SlagMin, double SlagMax){
614 
615 
616  TH2F* lSlag = new TH2F("LSLAG","Live time distribution over slags",NSlag,SlagMin/86400.,SlagMax/86400.,NSlag,SlagMin/86400.,SlagMax/86400.);
617  cout << "Start filling lSlag histogram : " << endl;
618  // float xlag[NIFO_MAX+1];
619  float xslag[NIFO_MAX+1];
620  double xlive;
621  live.SetBranchStatus("*",0); //disable all branches
622  live.SetBranchStatus("slag",1);
623  live.SetBranchStatus("live",1);
624  live.SetBranchAddress("slag",xslag);
625  // live.SetBranchAddress("lag",xlag);
626  live.SetBranchAddress("live",&xlive);
627  // live.SetBranchStatus("gps",false);
628 
629  long int ntrg = live.GetEntries();
630  for(long int l=0; l<ntrg; l++) {
631  live.GetEntry(l);
632  lSlag->Fill(xslag[0]/86400.,xslag[1]/86400.,xlive);
633  if(l%10000000==0) {cout << scientific << (double)l/ntrg <<"% - ";}
634 
635  }
636 
637  return lSlag;
638 }
639 //______________________________________________________________________________
640 
641 
642 
643 void CWB::CBCTool::doChirpFARPlot(int sel_events, double* recMchirp, double* injMchirp, double* far, TCanvas* c1, TString odir) {
644 
645 
646  c1->Clear();
647  c1->SetLogx(0);
648  c1->SetLogy(0);
649  c1->SetLogz(1);
650  c1->SetTheta(89.999);// the setting of the theta and phi is for the viewing angle: there is no TGraphError with "pcol" option, so the trick is to use the 2D version and use a viewing angle the closest as possible to the vertical (Drawbacks: some artefacts due to that can be noticed; no gridlines)
651  c1->SetPhi(0.0001);
652  double efar[sel_events];
653  for(int i=0;i<sel_events;i++) {efar[i]=0.0;}
654 
655 
656 
657  TGraph2DErrors* gr2=new TGraph2DErrors(sel_events,recMchirp,injMchirp,far,efar,efar,efar);
658  gr2->SetMarkerStyle(20);
659  gr2->SetLineColor(20);
660  //gr2->SetMinimum(1e-4);
661  //gr2->SetMaximum(1e1);
662  gr2->SetTitle("");
663  gr2->GetHistogram()->GetXaxis()->SetTitle("Injected mchirp");
664  gr2->GetHistogram()->GetYaxis()->SetTitle("Recovered mchirp");
665  gr2->GetHistogram()->GetXaxis()->SetLimits(0.,70.); //MY: Add similar line for the y axis? declare the X and Y bounds in the definition
666  gr2->GetHistogram()->GetYaxis()->SetLimits(0.,70.); //MY: Add similar line for the y axis? declare the X and Y bounds in the definition
667  gr2->GetHistogram()->GetZaxis()->SetTitle("False Alarm Rate [years^{-1}]");
668  gr2->GetXaxis()->SetTitleOffset(1.3);
669  gr2->GetYaxis()->SetTitleOffset(1.25);
670  gr2->GetXaxis()->SetTickLength(0.02);
671  gr2->GetYaxis()->SetTickLength(0.01);
672  gr2->GetXaxis()->CenterTitle(kTRUE);
673  gr2->GetYaxis()->CenterTitle(kTRUE);
674  gr2->GetZaxis()->CenterTitle(kTRUE);
675  gr2->GetXaxis()->SetTitleFont(42);
676  gr2->GetXaxis()->SetLabelFont(42);
677  gr2->GetYaxis()->SetTitleFont(42);
678  gr2->GetYaxis()->SetLabelFont(42);
679 // gr2->GetYaxis()->SetMoreLogLabels(kTRUE);
680  gr2->GetYaxis()->SetNoExponent(kTRUE);
681 
682  gr2->Draw("err zcolpcol");
683  char fname[1024];
684  sprintf(fname,"%s/ChirpFAR.png",odir.Data());
685  c1->Update();
686  c1->SaveAs(fname);
687 
688  return;
689 
690 }
691 
692 //______________________________________________________________________________
693 double CWB::CBCTool::calc_isco_radius(double a) {
694 
695  //Calculate the ISCO radius of a Kerr BH as a function of the Kerr parameter
696  //a : Kerr parameter
697  // Ref. Eq. (2.5) of Ori, Thorne Phys Rev D 62 124022 (2000)
698 
699  double z1 = 1.+pow((1.-pow(a,2.)),1./3)*(pow(1.+a,1./3) + pow(1.-a,1./3));
700  double z2 = pow(3.*pow(a,2.) + pow(z1,2.),1./2);
701  double a_sign = TMath::Sign(1.,a);
702  return 3+z2 - pow((3.-z1)*(3.+z1+2.*z2),1./2)*a_sign;
703 
704 }
705 
706 
707 //______________________________________________________________________________
708 
709 double CWB::CBCTool::calc_isco_freq(double a) {
710  /*
711  Calculate the ISCO frequency of a Kerr BH as a function of the Kerr parameter
712 
713  a : Kerr parameter (numpy array)
714  */
715 
716  double r_isco = CWB::CBCTool::calc_isco_radius(a);
717  double u_isco = pow(r_isco,-0.5);
718  double v_isco = u_isco*pow(1.-a*pow(u_isco,3.)+pow(a,2.)*pow(u_isco,6.),1./3.);
719  return pow(v_isco,3.)/TMath::Pi();
720 
721 }
722 //______________________________________________________________________________
723 
724 
725 double CWB::CBCTool::_final_spin_diff(double a_f, double eta, double delta_m, double S, double Delta){
726  // Internal function: the final spin is determined by minimizing this function
727 
728  // calculate ISCO radius
729  double r_isco = CWB::CBCTool::calc_isco_radius(a_f);
730 
731  // angular momentum at ISCO -- Eq.(2.8) of Ori, Thorne Phys Rev D 62 124022 (2000)
732  double J_isco = (3*pow(r_isco,1./2)-2*a_f)*2./pow(3*r_isco,1./2);
733 
734  // fitting coefficients - Table X1 of Healy et al Phys Rev D 90, 104004 (2014)
735  // [forth order fits]
736  double L0 = 0.686710;
737  double L1 = 0.613247;
738  double L2a = -0.145427;
739  double L2b = -0.115689;
740  double L2c = -0.005254;
741  double L2d = 0.801838;
742  double L3a = -0.073839;
743  double L3b = 0.004759;
744  double L3c = -0.078377;
745  double L3d = 1.585809;
746  double L4a = -0.003050;
747  double L4b = -0.002968;
748  double L4c = 0.004364;
749  double L4d = -0.047204;
750  double L4e = -0.053099;
751  double L4f = 0.953458;
752  double L4g = -0.067998;
753  double L4h = 0.001629;
754  double L4i = -0.066693;
755 
756  double a_f_new = pow((4.*eta),2.)*(L0 + L1*S + L2a*Delta*delta_m + L2b*pow(S,2.) + L2c*pow(Delta,2.) \
757  + L2d*pow(delta_m,2.) + L3a*Delta*S*delta_m + L3b*S*pow(Delta,2.) + L3c*pow(S,3.) \
758  + L3d*S*pow(delta_m,2.) + L4a*Delta*pow(S,2)*delta_m + L4b*pow(Delta,3.)*delta_m \
759  + L4c*pow(Delta,4.) + L4d*pow(S,4.) + L4e*pow(Delta,2.)*pow(S,2.) + L4f*pow(delta_m,4.) + L4g*Delta*pow(delta_m,3.) \
760  + L4h*pow(Delta,2.)*pow(delta_m,2.) + L4i*pow(S,2.)*pow(delta_m,2.)) \
761  + S*(1. + 8.*eta)*pow(delta_m,4.) + eta*J_isco*pow(delta_m,6.);
762 
763  return TMath::Abs(a_f-a_f_new);
764 
765 
766 
767 
768 }
769 
770 //______________________________________________________________________________
771 
772 double CWB::CBCTool::_final_spin_diff(Double_t *x, Double_t *par){
773 
774  double a_f = x[0];
775  double eta = par[0];
776  double delta_m = par[1];
777  double S = par[2];
778  double Delta = par[3];
779  // Internal function: the final spin is determined by minimizing this function
780 
781  // calculate ISCO radius
782  double r_isco = CWB::CBCTool::calc_isco_radius(a_f);
783 
784  // angular momentum at ISCO -- Eq.(2.8) of Ori, Thorne Phys Rev D 62 124022 (2000)
785  double J_isco = (3*pow(r_isco,1./2)-2*a_f)*2./pow(3*r_isco,1./2);
786 
787  // fitting coefficients - Table X1 of Healy et al Phys Rev D 90, 104004 (2014)
788  // [forth order fits]
789  double L0 = 0.686710;
790  double L1 = 0.613247;
791  double L2a = -0.145427;
792  double L2b = -0.115689;
793  double L2c = -0.005254;
794  double L2d = 0.801838;
795  double L3a = -0.073839;
796  double L3b = 0.004759;
797  double L3c = -0.078377;
798  double L3d = 1.585809;
799  double L4a = -0.003050;
800  double L4b = -0.002968;
801  double L4c = 0.004364;
802  double L4d = -0.047204;
803  double L4e = -0.053099;
804  double L4f = 0.953458;
805  double L4g = -0.067998;
806  double L4h = 0.001629;
807  double L4i = -0.066693;
808 
809  double a_f_new = pow((4.*eta),2.)*(L0 + L1*S + L2a*Delta*delta_m + L2b*pow(S,2.) + L2c*pow(Delta,2.) \
810  + L2d*pow(delta_m,2.) + L3a*Delta*S*delta_m + L3b*S*pow(Delta,2.) + L3c*pow(S,3.) \
811  + L3d*S*pow(delta_m,2.) + L4a*Delta*pow(S,2)*delta_m + L4b*pow(Delta,3.)*delta_m \
812  + L4c*pow(Delta,4.) + L4d*pow(S,4.) + L4e*pow(Delta,2.)*pow(S,2.) + L4f*pow(delta_m,4.) + L4g*Delta*pow(delta_m,3.) \
813  + L4h*pow(Delta,2.)*pow(delta_m,2.) + L4i*pow(S,2.)*pow(delta_m,2.)) \
814  + S*(1. + 8.*eta)*pow(delta_m,4.) + eta*J_isco*pow(delta_m,6.);
815 
816  return TMath::Abs(a_f-a_f_new);
817 
818 
819 
820 
821 }
822 
823 //______________________________________________________________________________
824 
825 double CWB::CBCTool::bbh_final_mass_and_spin_non_precessing(double m1, double m2, double chi1, double chi2){
826  /*
827  Calculate the mass and spin of the final BH resulting from the
828  merger of two black holes with non-precessing spins
829 
830  m1, m2: component masses
831  chi1, chi2: dimensionless spins of two BHs
832  */
833 
834  // binary parameters
835  double m = m1+m2;
836  double q = m1/m2;
837  double eta = q/pow((1.+q),2.);
838  double delta_m = (m1-m2)/m;
839 
840  double S1 = chi1*pow(m1,2.); // spin angular momentum 1
841  double S2 = chi2*pow(m2,2.); // spin angular momentum 2
842  double S = (S1+S2)/pow(m,2.); // symmetric spin (dimensionless -- called \tilde{S} in the paper)
843  double Delta = (S2/m2-S1/m1)/m; // antisymmetric spin (dimensionless -- called tilde{Delta} in the paper
844 
845  //# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
846  //# compute the final spin
847  //# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
848  //# res = so.minimize_scalar(_final_spin_diff, bounds=(-0.99999, 0.99999), args=(eta, delta_m, S, Delta), method='Bounded', tol=1e-6, options={'maxiter':100, 'disp':False})
849  // a_f = res.x
850 
851  //x, cov_x = so.leastsq(_final_spin_diff, 0., args=(eta, delta_m, S, Delta))
852  TF1 f("spin_diff",CWB::CBCTool::_final_spin_diff,-1.0,1.0,4);
853  f.SetParameters(eta,delta_m, S, Delta);
854  ROOT::Math::WrappedTF1 wf1(f);
855 
856  // Create the Integrator
857  ROOT::Math::BrentRootFinder brf;
858 
859  // Set parameters of the method
860  brf.SetFunction( wf1, -1.0, 1.0 );
861  brf.Solve();
862 
863  //cout << brf.Root() << endl;
864  double a_f = brf.Root();
865 
866  //# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
867  //# now compute the final mass
868  //# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
869  double r_isco = CWB::CBCTool::calc_isco_radius(a_f);
870 
871  //# fitting coefficients - Table X1 of Healy et al Phys Rev D 90, 104004 (2014)
872  //# [forth order fits]
873  double M0 = 0.951507;
874  double K1 = -0.051379;
875  double K2a = -0.004804;
876  double K2b = -0.054522;
877  double K2c = -0.000022;
878  double K2d = 1.995246;
879  double K3a = 0.007064;
880  double K3b = -0.017599;
881  double K3c = -0.119175;
882  double K3d = 0.025000;
883  double K4a = -0.068981;
884  double K4b = -0.011383;
885  double K4c = -0.002284;
886  double K4d = -0.165658;
887  double K4e = 0.019403;
888  double K4f = 2.980990;
889  double K4g = 0.020250;
890  double K4h = -0.004091;
891  double K4i = 0.078441;
892 
893  //# binding energy at ISCO -- Eq.(2.7) of Ori, Thorne Phys Rev D 62 124022 (2000)
894  double E_isco = (1. - 2./r_isco + a_f/pow(r_isco,1.5))/pow(1. - 3./r_isco + 2.*a_f/pow(r_isco,1.5),1./2.);
895 
896  //# final mass -- Eq. (14) of Healy et al Phys Rev D 90, 104004 (2014)
897  double mf = pow(4.*eta,2.)*(M0 + K1*S + K2a*Delta*delta_m + K2b*pow(S,2.) + K2c*pow(Delta,2.) + K2d*pow(delta_m,2.) \
898  + K3a*Delta*S*delta_m + K3b*S*pow(Delta,2.) + K3c*pow(S,3.) + K3d*S*pow(delta_m,2.) \
899  + K4a*Delta*pow(S,2.)*delta_m + K4b*pow(Delta,3.)*delta_m + K4c*pow(Delta,4.) + K4d*pow(S,4.) \
900  + K4e*pow(Delta,2.)*pow(S,2.) + K4f*pow(delta_m,4.) + K4g*Delta*pow(delta_m,3.) + K4h*pow(Delta,2.)*pow(delta_m,2.) \
901  + K4i*pow(S,2.)*pow(delta_m,2.)) + (1+eta*(E_isco+11.))*pow(delta_m,6.);
902 
903  return mf*m;
904 
905 }
906 //______________________________________________________________________________
907 
908 /*
909 void CWB::CBCTool::doEffectiveRadiusChiPlot(){
910 
911 
912  int spin_mtot_bins = 0;
913  double V0_spin_mtot = 0.0;
914  for(int j=0; j<NBINS_MTOT+1; j++){
915  for(int k=0; k<NBINS_SPIN+1; k++){
916 
917  // volume_first_shell[j][k] = efficiency_first_shell->GetBinContent(j+1,k+1);
918  // if(factor_events_rec->GetBinContent(j+1,k+1) != 0.) {
919  // error_volume_first_shell[j][k] = 1./TMath::Sqrt(factor_events_rec->GetBinContent(j+1,k+1));
920  // massbins++;
921  //}
922  if(spin_mtot_volume[j][k] != 0.) {
923  spin_mtot_volume[j][k] = shell_volume*spin_mtot_volume[j][k]; /// Warning: neglecting the internal volume... + volume_internal_sphere*volume_first_shell[j][k];
924  V0_spin_mtot += spin_mtot_volume[j][k];
925  error_spin_mtot_volume[j][k] = shell_volume*TMath::Sqrt(error_spin_mtot_volume[j][k]); ///Warning: neglecting the internal volume...+ volume_internal_sphere*volume_first_shell[j][k]*error_volume_first_shell[j][k];
926 
927  spin_mtot_radius[j][k] = pow(3.*spin_mtot_volume[j][k]/(4*TMath::Pi()),1./3);
928 
929  error_spin_mtot_radius[j][k] = (1./3)*pow(3./(4*TMath::Pi()),1./3)*pow(1./pow(spin_mtot_volume[j][k],2),1./3)*error_spin_mtot_volume[j][k];
930  spin_mtot_bins++;
931  }
932  cout<<j<< " "<<k<< " "<<spin_mtot_volume[j][k]<<" "<<spin_mtot_radius[j][k]<<endl;
933  }
934  }
935  cout << "Average Spin-Mtot Volume at threshold V0 = "<<V0_spin_mtot/spin_mtot_bins <<endl;
936  c1->Clear();
937 
938  TH2F* h_spin_mtot_radius = new TH2F("h_spin_mtot_radius","",NBINS_SPIN,MINCHI,MAXCHI,NBINS_MTOT,minMtot,maxMtot);
939  //h_spin_mtot_radius->GetXaxis()->SetRangeUser(MIN_plot_mass1,MAX_plot_mass1);
940  //h_spin_mtot_radius->GetYaxis()->SetRangeUser(MIN_plot_mass2,MAX_plot_mass2);
941  h_spin_mtot_radius->GetXaxis()->SetTitle("#chi_{z}");
942  h_spin_mtot_radius->GetYaxis()->SetTitle("Total Mass (M_{#odot})");
943  h_spin_mtot_radius->GetXaxis()->SetTitleOffset(1.3);
944  h_spin_mtot_radius->GetYaxis()->SetTitleOffset(1.25);
945  h_spin_mtot_radius->GetXaxis()->CenterTitle(kTRUE);
946  h_spin_mtot_radius->GetYaxis()->CenterTitle(kTRUE);
947  h_spin_mtot_radius->GetXaxis()->SetNdivisions(410);
948  h_spin_mtot_radius->GetYaxis()->SetNdivisions(410);
949  h_spin_mtot_radius->GetXaxis()->SetTickLength(0.01);
950  h_spin_mtot_radius->GetYaxis()->SetTickLength(0.01);
951  h_spin_mtot_radius->GetZaxis()->SetTickLength(0.01);
952  h_spin_mtot_radius->GetXaxis()->SetTitleFont(42);
953  h_spin_mtot_radius->GetXaxis()->SetLabelFont(42);
954  h_spin_mtot_radius->GetYaxis()->SetTitleFont(42);
955  h_spin_mtot_radius->GetYaxis()->SetLabelFont(42);
956  h_spin_mtot_radius->GetZaxis()->SetLabelFont(42);
957  h_spin_mtot_radius->GetZaxis()->SetLabelSize(0.03);
958  h_spin_mtot_radius->SetTitle("");
959 
960 
961  for(int i=1; i<=NBINS_MTOT+1; i++){
962  for(int j=1; j<=NBINS_SPIN+1; j++){
963  h_spin_mtot_radius->SetBinContent(j,i,spin_mtot_radius[i-1][j-1]);
964  h_spin_mtot_radius->SetBinError(j,i,error_spin_mtot_radius[i-1][j-1]);
965  cout<<j<< " "<<i<< " "<<h_spin_mtot_radius->GetBinContent(j,i)<<endl;
966  }
967  }
968 
969 
970  h_spin_mtot_radius->SetContour(NCont);
971  h_spin_mtot_radius->SetEntries(1); // This option needs to be enabled when filling 2D histogram with SetBinContent
972  h_spin_mtot_radius->Draw("colz text colsize=2"); // Option to write error associated to the bin content
973  //h_spin_mtot_radius->Draw("colz text");
974  h_spin_mtot_radius->GetZaxis()->SetRangeUser(0,MAX_EFFECTIVE_RADIUS/2.);
975 
976  TExec *ex2 = new TExec("ex2","gStyle->SetPaintTextFormat(\".0f\");");
977  ex2->Draw();
978 
979 
980  //char radius_title[256];
981  sprintf(radius_title,"Effective radius Mtot vs #chi_{z} (Mpc, %.1f < #chi_{z} < %.1f)",MINCHI,MAXCHI);
982 
983 
984 
985  TPaveText *p_radius = new TPaveText(0.325301,0.926166,0.767068,0.997409,"blNDC");
986  p_radius->SetBorderSize(0);
987  p_radius->SetFillColor(0);
988  p_radius->SetTextColor(1);
989  p_radius->SetTextFont(32);
990  p_radius->SetTextSize(0.045);
991  TText *text = p_radius->AddText(radius_title);
992  p_radius->Draw();
993 
994  sprintf(fname,"%s/Effective_radius_chi.png",netdir);
995 
996  c1->Update();
997  c1->SaveAs(fname);
998 
999 
1000 }
1001 
1002 */
1003 
1004 
1005 
tree
TTree * tree
Definition: TimeSortTree.C:20
rho
double rho
Definition: cwb_report_prod_2.C:149
cut
char cut[512]
Definition: check_slag_bkg.C:62
cwb_online.f
tuple f
Definition: cwb_online.py:91
live
TString live
Definition: HistoryParamFileTest.C:18
fprintf
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
Live
double Live
Definition: cwb_report_prod_2.C:197
networkname
char networkname[256]
Definition: cbc_plots.C:95
CWB::CBCTool::bbh_final_mass_and_spin_non_precessing
static double bbh_final_mass_and_spin_non_precessing(double m1, double m2, double chi1, double chi2)
Definition: CBCTool.cc:825
nentries
Int_t nentries
Definition: TimeSortTree.C:24
a
wavearray< double > a(hp.size())
CWB::CBCTool::getZeroLiveTime2
static double getZeroLiveTime2(int nIFO, TChain &liv)
Definition: CBCTool.cc:563
TString
TString("c")
Definition: cwb_report_skymap.C:111
CWB::CBCTool::getZeroLiveTime
static double getZeroLiveTime(int nIFO, TChain &liv)
Definition: CBCTool.cc:539
RHO_NBINS
#define RHO_NBINS
Definition: cbc_plots.C:4
x
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
Definition: BaseFunctionWDM.C:51
odir
char odir[1024]
Definition: cwb_mkhtml_header.C:148
xslag
float xslag[6]
Definition: cwb_fix_live_slag_missed.C:114
CWB::CBCTool::doChirpFARPlot
static void doChirpFARPlot(int sel_events, double *recMchirp, double *injMchirp, double *far, TCanvas *c1, TString odir)
Definition: CBCTool.cc:643
xlag
float xlag
Definition: cwb_report_prod_2.C:192
CWB::CBCTool::doRangePlot
static TF1 * doRangePlot(int RHO_NBINS, double *Trho, double *Rrho, double *eRrho, float RHO_MIN, float T_out, TCanvas *c1, TString networkname, TString odir, bool write_ascii)
Definition: CBCTool.cc:362
m
int m
Definition: cwb_net.C:10
j
int j
Definition: cwb_net.C:10
CWB::CBCTool::getLiveTime2
static double getLiveTime2(TChain &liv)
Definition: CBCTool.cc:519
i
i drho i
Definition: cwb_epparameters.C:85
CWB::CBCTool::GetTreeIndex
static Long64_t GetTreeIndex(TTree *tree, const char *selection)
Definition: CBCTool.cc:25
CWB::CBCTool::getLiveTime
static double getLiveTime(int nIFO, TChain &liv, int lag_number, int slag_number, int dummy)
Definition: CBCTool.cc:444
w
wavearray< double > w
Definition: Test1.C:27
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
S2
Meyer< double > * S2
Definition: TestTObjectMeyer.C:26
CWB::CBCTool::_final_spin_diff
static double _final_spin_diff(double a_f, double eta, double delta_m, double S, double Delta)
Definition: CBCTool.cc:725
h
wavearray< double > h
Definition: Regression_H1.C:25
c1
TCanvas * c1
Definition: cwb_report_prod_2.C:69
int
i() int(T_cor *100))
Pi
double Pi
Definition: DrawPhaseShift.C:12
f2
TF1 * f2
Definition: cbc_plots.C:1710
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:4
L1
gNET add & L1
Definition: DrawGNET.C:9
next
TIter next(twave->GetListOfBranches())
fname
char fname[1024]
Definition: cwb_report_prod_1.C:110
CWB::CBCTool::getFAR
static double getFAR(float rho, TH1 *hc, double liveTot)
Definition: CBCTool.cc:602
xlive
double xlive
Definition: cwb_report_prod_2.C:194
liveTot
double liveTot
Definition: cwb_report_prod_2.C:195
m1
m1
Definition: cbc_plots.C:606
lSlag
TH1F * lSlag
Definition: check_slag_bkg.C:107
chi2
float chi2
Definition: cbc_plots.C:603
CWB::CBCTool::calc_isco_radius
static double calc_isco_radius(double a)
Definition: CBCTool.cc:693
par
vector< mdcpar > par
Definition: CWB_Plugin_HEN_SIM_Config.C:83
e
double e
Definition: CreateListEBBH.C:35
slagFound
bool slagFound
Definition: cwb_report_prod_2.C:108
gr2
TGraphErrors * gr2
Definition: cwb_report_prod_2.C:1023
CWB::CBCTool::calc_isco_freq
static double calc_isco_freq(double a)
Definition: CBCTool.cc:709
CWB::CBCTool::FillSLagHist
static TH2F * FillSLagHist(int NIFO_MAX, TChain &live, int NSlag, double SlagMin, double SlagMax)
Definition: CBCTool.cc:613
CWB::CBCTool::doROCPlot
static void doROCPlot(int bkg_entries, double *rho_bkg, int *index, float RHO_BIN, double liveTot, TF1 *AverageRad, TCanvas *c1, TString odir, bool write_ascii)
Definition: CBCTool.cc:88
index
wavearray< int > index
Definition: TestEulerCharacteristic.C:5
branch
TBranch * branch
Definition: cwb_fix_live_slag_missed.C:42
slag
Definition: Toolbox.hh:81
S
Meyer< double > S(1024, 2)
l
int l
Definition: cbc_plots.C:434
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
ntrg
int ntrg
Definition: cwb_fix_live_slag_missed.C:124
T_out
double T_out
Definition: cwb_pparameters.C:144
write_ascii
bool write_ascii
Definition: cbc_plots.C:447
RHO_BIN
#define RHO_BIN
Definition: cbc_plots.C:3
fclose
fclose(ftrig)
RHO_MIN
#define RHO_MIN
Definition: cbc_plots.C:2
Delta
static double Delta
Definition: geodesics.cc:8
m2
m2
Definition: cbc_plots.C:607