Logo coherent WaveBurst  
Library Reference Guide
Logo
CreateDistanceParplots.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 // Creates various distance vs pars plots
21 // Note : this macro is used to generate the CBC report
22 // Author : Francesco Salemi
23 
24 #include <TComplex.h>
25 #include <TRandom.h>
26 #include <TStyle.h>
27 #include <fstream>
28 #include <iostream>
29 #include "Math/Rotation3D.h"
30 #include "Math/Vector3Dfwd.h"
31 #include "TCanvas.h"
32 #include "TFile.h"
33 #include "TGraphErrors.h"
34 #include "TH1F.h"
35 #include "TH2F.h"
36 #include "TLegend.h"
37 #include "TMath.h"
38 #include "TMultiGraph.h"
39 #include "TROOT.h"
40 #include "TRatioPlot.h"
41 #include "TRotation.h"
42 #include "TSystem.h"
43 #include "TTree.h"
44 #include "TTreeIndex.h"
45 #include "TVector3.h"
46 
47 #define MAXY 50000.0
48 #define MAXSNR 150.0
49 
50 // TGraphErrors* CreateGraphRadiusIFAR (char* sim_file_name, char*
51 // mdc_file_name, float shell_volume, double liveTot=1e6, float T_ifar=0.0, float
52 // T_win=0.2, int TRIALS = 1, TString SEL="", Color_t color=kBlue, TString
53 // opt="default"){
55  char* mdc_file_name,
56  char* netdir,
57  TString opt = "",
58  double MINX = 0.0,
59  double MAXX = 1.0,
60  double MAXDISTANCE = 5000.,
61  int NBIN_DIST = 10,
62  float T_ifar = 0.0,
63  float T_win = 0.2,
64  int nIFO = 2) {
65  TCanvas* co_canvas3 = new TCanvas("sd3", "SD3", 3, 47, 1000, 802);
66  co_canvas3->SetGridx();
67  co_canvas3->SetGridy();
68  co_canvas3->SetLogy();
69 
70  float myifar, netcc[3];
71  float rho[2];
72  double mytime[6];
73  float factor, mydistance, mchirp;
74  float mass[2];
75  float spin[6];
76  float chip;
77  float iSNR[nIFO];
78  float snr[nIFO];
79  float chirp[6];
80  float range[2];
81  float iota[2];
82 
83  CWB::CBCTool cbcTool;
84 
85  TFile* filein = new TFile(sim_file_name);
86  TTree* sim = nullptr;
87  filein->GetObject("waveburst", sim);
88  if (!sim->GetListOfBranches()->FindObject("chip")) {
89  cout << "Adding Chi_p branch to wave tree" << endl;
90  cbcTool.AddChip(sim_file_name, "waveburst");
91  }
92  sim->SetBranchAddress("mass", mass);
93  sim->SetBranchAddress("factor", &factor);
94  sim->SetBranchAddress("range", range);
95  sim->SetBranchAddress("chirp", chirp);
96  sim->SetBranchAddress("rho", rho);
97  sim->SetBranchAddress("netcc", netcc);
98  sim->SetBranchAddress("ifar", &myifar);
99  sim->SetBranchAddress("time", mytime);
100  sim->SetBranchAddress("spin", spin);
101  sim->SetBranchAddress("chip", &chip);
102  sim->SetBranchAddress("iSNR", iSNR);
103  sim->SetBranchAddress("iota", iota);
104 
105  TFile* filein2 = new TFile(mdc_file_name);
106  TTree* mdc = nullptr;
107  filein2->GetObject("mdc", mdc);
108  gROOT->cd();
109  if (!mdc->GetListOfBranches()->FindObject("chip")) {
110  cout << "Adding Chi_p branch to mdc tree" << endl;
111  cbcTool.AddChip(mdc_file_name, "mdc");
112  }
113 
114  mdc->SetBranchAddress("time", mytime);
115  mdc->SetBranchAddress("mass", mass);
116  mdc->SetBranchAddress("factor", &factor);
117  mdc->SetBranchAddress("distance", &mydistance);
118  mdc->SetBranchAddress("mchirp", &mchirp);
119  mdc->SetBranchAddress("spin", spin);
120  mdc->SetBranchAddress("chip", &chip);
121  mdc->SetBranchAddress("snr", snr);
122  mdc->SetBranchAddress("iota", iota);
123 
124  int nevts = (int)mdc->GetEntries();
125  float CYS = 86400. * 365.25;
126 
127  cout << nevts << " injected signals " << sim->GetEntries("ifar>0")
128  << " recovered signals" << endl;
129  int countv = 0;
130  int countvifar = 0;
131  int cnt = 0;
132  std::vector<double> xi, xr, yi, yr;
133 
134  auto inj = new TH2F("Injected snr vs stat inj", "", NBIN_DIST, MINX,
135  MAXX, 100, 0.0, MAXSNR);
136  // auto SNR= new TH2F("iSNR vs stat", "", 10, MINX, MAXX, 100,0.0,
137  // 100.0);
138  TMultiGraph* mg = new TMultiGraph();
139 
140  if (opt.Contains("chieff")) {
141  inj->GetXaxis()->SetTitle("#chi_{eff}");
142  mg->GetXaxis()->SetTitle("#chi_{eff}");
143  } else if (opt.Contains("chip")) {
144  inj->GetXaxis()->SetTitle("#chi_{p}");
145  mg->GetXaxis()->SetTitle("#chi_{p}");
146  } else if (opt.Contains("chirp")) {
147  inj->GetXaxis()->SetTitle("Chirp Mass (M_{#odot})");
148  mg->GetXaxis()->SetTitle("Chirp Mass (M_{#odot})");
149  } else if (opt.Contains("mtot")) {
150  inj->GetXaxis()->SetTitle("Total Mass (M_{#odot})");
151  mg->GetXaxis()->SetTitle("Total Mass (M_{#odot})");
152  } else if (opt.Contains("eta")) {
153  inj->GetXaxis()->SetTitle("#eta, Symmetric Mass Ratio");
154  mg->GetXaxis()->SetTitle("#eta, Symmetric Mass Ratio");
155  } else if (opt.Contains("iota")) {
156  inj->GetXaxis()->SetTitle("Cos(#iota), Inclination");
157  mg->GetXaxis()->SetTitle("Cos(#iota), Inclination");
158  } else if (opt.Contains("distance")) {
159  inj->GetXaxis()->SetTitle("Sensitive Distance (Mpc)");
160  mg->GetXaxis()->SetTitle("Injected snr");
161  // cout <<"Distance plots" << endl;
162  } else {
163  cout << "Not a valid option! "
164  "opt={\"chip\",\"chirp\",\"eta\",\"mtot\", \"chieff\", \"iota\"}"
165  << endl;
166  exit(1);
167  }
168 
169  // inj->GetYaxis()->SetRangeUser(10., MAXDISTANCE);
170  inj->GetYaxis()->SetTitle("Injected snr");
171  inj->GetXaxis()->SetTitleOffset(1.3);
172  inj->GetYaxis()->SetTitleOffset(1.3);
173  inj->GetXaxis()->SetTickLength(0.01);
174  inj->GetYaxis()->SetTickLength(0.01);
175  inj->GetXaxis()->CenterTitle(kTRUE);
176  inj->GetYaxis()->CenterTitle(kTRUE);
177  inj->GetXaxis()->SetTitleFont(42);
178  inj->GetXaxis()->SetLabelFont(42);
179  inj->GetYaxis()->SetTitleFont(42);
180  inj->GetYaxis()->SetLabelFont(42);
181  inj->SetMarkerStyle(20);
182  inj->SetMarkerSize(0.5);
183  inj->SetMarkerColor(2);
184  inj->SetTitle("");
185  TH2F* rec = (TH2F*)inj->Clone("Injected snr vs stat rec");
186  // rec->GetYaxis()->SetRangeUser(10., MAXDISTANCE/1000);
187  // TH2F *SNRrec = (TH2F *)SNR->Clone("iSNR vs stat rec");
188  rec->SetMarkerColor(4);
189 
190  mg->GetYaxis()->SetTitle("Sensitive Distance (Mpc)");
191  // mg->GetXaxis()->SetTitle("Inverse False Alarm Rate [yr]");
192  mg->GetXaxis()->SetTitleOffset(1.3);
193  mg->GetYaxis()->SetTitleOffset(1.3);
194  mg->GetXaxis()->SetTickLength(0.01);
195  mg->GetYaxis()->SetTickLength(0.01);
196  mg->GetXaxis()->CenterTitle(kTRUE);
197  mg->GetYaxis()->CenterTitle(kTRUE);
198  mg->GetXaxis()->SetTitleFont(42);
199  mg->GetXaxis()->SetLabelFont(42);
200  mg->GetYaxis()->SetTitleFont(42);
201  mg->GetYaxis()->SetLabelFont(42);
202  mg->SetTitle("");
203 
204  //Loop over mdc TTree
205  float SNR2, mSNR;
206  for (int g = 0; g < (int)mdc->GetEntries(); g++) {
207  mdc->GetEntry(g);
208  SNR2 = pow(snr[0], 2.0) + pow(snr[1], 2.0);
209  for (int i = 2; i < nIFO; i++) {
210  SNR2 += pow(snr[i], 2.0);
211  }
212  mSNR = TMath::Sqrt(SNR2);
213  yi.push_back(mydistance);
214  if (opt.Contains("chieff")) {
215  xi.push_back((spin[2] * mass[0] + spin[5] * mass[1]) /
216  (mass[1] + mass[0]));
217  } else if (opt.Contains("chip")) {
218  xi.push_back(chip);
219  } else if (opt.Contains("chirp")) {
220  xi.push_back(mchirp);
221  } else if (opt.Contains("mtot")) {
222  xi.push_back(mass[0] + mass[1]);
223  } else if (opt.Contains("eta")) {
224  xi.push_back(mass[0] * mass[1] /
225  pow(mass[0] + mass[1], 2.0));
226  } else if (opt.Contains("iota")) {
227  xi.push_back(iota[1]);
228  } else if (opt.Contains("distance")) {
229  xi.push_back(mydistance);
230  }
231 
232  inj->Fill(xi[xi.size() - 1], mSNR);
233  // SNR->Fill(xi[xi.size()-1], mSNR);
234  }
235  // cout << MINX << " " << MAXX << " " << NBIN_DIST << " " << MAXDISTANCE
236  // << endl;
237  inj->Draw("p");
238  // Loop over sim TTree
239  for (int g = 0; g < (int)sim->GetEntries(); g++) {
240  sim->GetEntry(g);
241  // ifactor = (int)factor - 1;
242  // cout << "g=" << g << " trueindex=" <<index1[g] <<" IFAR=" <<
243  // myifar/CYS << " RHO1=" << rho[1] <<endl;
244  if (myifar <= T_ifar * CYS) {
245  countvifar++;
246  // cout << g << " " <<index1[g] <<" " << myifar/CYS << "
247  // ";
248  continue;
249  }
250 
251  if ((mytime[0] - mytime[nIFO]) < -T_win ||
252  (mytime[0] - mytime[nIFO]) > 2 * T_win) {
253  countv++;
254  continue;
255  } // NOT checking for detector 1 and 2: very small bias...
256  SNR2 = 0.0;
257  for (int i = 0; i < nIFO; i++) {
258  SNR2 += iSNR[i];
259  }
260  mSNR = TMath::Sqrt(SNR2);
261  yr.push_back(range[1]);
262  if (opt.Contains("chieff")) {
263  xr.push_back((spin[2] * mass[0] + spin[5] * mass[1]) /
264  (mass[1] + mass[0]));
265  cnt++;
266  } else if (opt.Contains("chip")) {
267  xr.push_back(chip);
268  cnt++;
269  } else if (opt.Contains("chirp")) {
270  xr.push_back(chirp[0]);
271  cnt++;
272  } else if (opt.Contains("mtot")) {
273  xr.push_back(mass[1] + mass[0]);
274  cnt++;
275  } else if (opt.Contains("eta")) {
276  xr.push_back(mass[0] * mass[1] /
277  pow(mass[0] + mass[1], 2.0));
278  cnt++;
279  } else if (opt.Contains("iota")) {
280  xr.push_back(iota[1]);
281  cnt++;
282  } else if (opt.Contains("distance")) {
283  xr.push_back(range[1]);
284  cnt++;
285  }
286  rec->Fill(xr[xr.size() - 1], mSNR);
287  }
288 
289  cout << endl;
290  // cout << countvifar << " events vetoed by T_ifar : " << T_ifar << endl;
291  //cout << countv << " events vetoed by T_win" << endl;
292  //cout << rec->GetEntries() << " events selected" << endl;
293 
294  char lab[1024];
295  char fname[1024];
296  char fname2[1024];
297  char fname3[1024];
298 
299  sprintf(fname, "%s/iSNR_vs_%s.eps", netdir, opt.Data());
300  sprintf(fname3, "%s/Distance_vs_%s.eps", netdir, opt.Data());
301  sprintf(fname2, "%s/%s_distribution.png", netdir, opt.Data());
302 
303  // D_Chi_rec->GetYaxis()->SetRangeUser(10.,3*MAXDISTANCE);
304  // co_canvas3->SetLogx();
305  inj->Draw("p");
306  rec->Draw("p same");
307  auto leg_D = new TLegend(0.6, 0.1, 0.9, 0.25, "", "brNDC");
308  sprintf(lab, "Injections: %i", (int)mdc->GetEntries());
309  leg_D->AddEntry("", lab, "a");
310  sprintf(lab, "found: %i", cnt);
311  leg_D->AddEntry(rec, lab, "p");
312  sprintf(lab, "missed: %i", (int)mdc->GetEntries() - cnt);
313  leg_D->AddEntry(inj, lab, "p");
314  leg_D->SetFillColor(0);
315  leg_D->SetFillColorAlpha(0, 0.9);
316  leg_D->Draw();
317  co_canvas3->SaveAs(fname);
318  co_canvas3->SetLogx(0);
319  // Test new scatter plot
320  TGraph* rec_gr = new TGraph(xr.size(), &xr[0], &yr[0]);
321  TGraph* inj_gr = new TGraph(xi.size(), &xi[0], &yi[0]);
322  TH1D* injx;
323  TH1D* recx;
324  TH1D* snrx;
325  if (!opt.Contains("distance")) {
326 
327  inj_gr->SetMarkerColor(2);
328  inj_gr->SetMarkerStyle(20);
329  inj_gr->SetMarkerSize(0.5);
330 
331  rec_gr->SetMarkerColor(4);
332  rec_gr->SetMarkerStyle(20);
333  rec_gr->SetMarkerSize(0.5);
334 
335  mg->Add(inj_gr);
336  mg->Add(rec_gr);
337  // mg->Draw("aple3");
338  co_canvas3->Clear();
339  mg->GetXaxis()->SetLimits(MINX, MAXX);
340  mg->GetYaxis()->SetRangeUser(10., MAXDISTANCE);
341  mg->Draw("ap");
342  // rec_gr->Draw("ap");
343  leg_D->Draw();
344  co_canvas3->SaveAs(fname3);
345  injx = (TH1D*)inj->ProjectionX();
346  recx = (TH1D*)rec->ProjectionX();
347  snrx = (TH1D*)inj->ProfileX();
348  } else {
349  injx = (TH1D*)inj->ProjectionX();
350  recx = (TH1D*)rec->ProjectionX();
351  snrx = (TH1D*)inj->ProfileX();
352  }
353 
354  injx->Sumw2();
355  // injx->SetFillColor(kRed);
356  injx->SetFillColorAlpha(kRed, 0.3);
357  injx->SetFillStyle(3004);
358  recx->Sumw2();
359  recx->SetFillColorAlpha(kBlue, 0.99);
360  recx->SetFillStyle(3001);
361 
362  co_canvas3->Clear();
363  co_canvas3->SetLogx(0);
364 
365  //Ratio plot
366  auto rp1 = new TRatioPlot(recx, injx, "divsym");
367  rp1->SetH1DrawOpt("PEBAR");
368  rp1->SetH2DrawOpt("PEBAR");
369  rp1->SetGraphDrawOpt("PB");
370  rp1->SetSeparationMargin(0.01);
371 
372  // rp1->SetConfidenceIntervalColors();
373  rp1->Draw("PEBAR");
374  rp1->SetSplitFraction(0.5);
375  // if (!opt.Contains("distance")){
376  rp1->GetUpperRefYaxis()->SetRangeUser(1., MAXY);
377  //rp1->GetLowerRefGraph()->SetMinimum(0.0);
378  //rp1->GetLowerRefGraph()->SetMaximum(1.0);
379  /* } else {
380  //co_canvas3->SetLog();
381  rp1->GetUpperRefYaxis()->SetRangeUser(1.,MAXY);
382  //rp1->GetUpperRefXaxis()->SetRangeUser(1.,2000.);
383  rp1->GetLowerRefGraph()->SetMinimum(0.0);
384  rp1->GetLowerRefGraph()->SetMaximum(1.0);
385  //rp1->GetLowYaxis()->SetNdivisions(510);
386  }*/
387 
388  rp1->GetLowerRefYaxis()->SetTitle("ratio");
389  rp1->GetLowerRefYaxis()->CenterTitle(kTRUE);
390  rp1->GetUpperRefYaxis()->SetTitle("entries");
391  rp1->GetUpperRefYaxis()->CenterTitle(kTRUE);
392  // rp1->GetUpperRefYaxis()->SetLineColor(kBlue);
393  rp1->GetLowerRefYaxis()->SetLabelColor(kBlue);
394  rp1->GetLowerRefXaxis()->CenterTitle(kTRUE);
395  rp1->GetLowerRefXaxis()->SetTitleOffset(1.3);
396  rp1->GetLowerRefGraph()->SetFillColor(kBlue);
397  rp1->GetLowerRefGraph()->SetFillStyle(3001);
398  // leg_D2->Draw();
399 
400  auto leg_D3 = new TLegend(0.6, 0.8455, 0.8965, 0.94555, "", "brNDC");
401  // sprintf(lab, "Injections: %i", (int)mdc->GetEntries());
402  // leg_D3->AddEntry("", lab, "a");
403  sprintf(lab, "Injected: %i", (int)mdc->GetEntries());
404  leg_D3->AddEntry(inj, lab, "p");
405  sprintf(lab, "Found: %i", cnt);
406  leg_D3->AddEntry(rec, lab, "p");
407  // leg_D3->SetFillColor(0);
408  leg_D3->SetFillColorAlpha(0, 0.9);
409  leg_D3->Draw();
410  // TGraph *g = rp1->GetLowerRefGraph();
411  // rp11 = TRatioPlot(h1,h0)
412  // rp1->GetUpperPad()->cd()
413  // h2->Draw("same")
414  rp1->GetLowerPad()->cd();
415  // TPad* mp = (TPad*)rp1->GetLowerPad();
416 
417  Float_t rightmax = 1.2 * snrx->GetMaximum();
418  Float_t rightmin = 0.8 * snrx->GetMinimum();
419  double ymax = rp1->GetLowerRefYaxis()->GetXmax();
420  double ymin = rp1->GetLowerRefYaxis()->GetXmin();
421  double xmax = rp1->GetLowerRefXaxis()->GetXmax();
422  Float_t scale = ( ymax ) / (rightmax );
423  // cout << "rightmax: " << rightmax << " Uymax: "<<gPad->GetUymax() << "
424  // scale: " << scale << endl; g->Draw(); rp1->GetLowerPad()->cd();
425  snrx->Scale(scale);
426  snrx->SetLineColor(kGreen + 2);
427  snrx->SetMarkerColor(kGreen + 2);
428  // snrx->SetFillColor(kGreen+2);
429  snrx->SetFillColorAlpha(kGreen + 2, 0.4);
430  snrx->SetFillStyle(3001);
431  snrx->SetMarkerStyle(20);
432  snrx->SetMarkerSize(0.5);
433 
434  snrx->Draw("SAME PEBAR");
435  // TGaxis*myaxis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
436  //gPad->GetUxmax(),gPad->GetUymax(), 0,rightmax,510,"+L");
437  TGaxis* myaxis =
438  new TGaxis(xmax, ymin, xmax, ymax, rightmin, rightmax, 510, "+SL");
439  myaxis->SetLineColor(kGreen + 2);
440  myaxis->SetLabelColor(kGreen + 2);
441  myaxis->SetTitle("Average Injected SNR");
442  myaxis->SetTitleOffset(0.7);
443  myaxis->SetTickLength(0.01);
444  myaxis->CenterTitle(kTRUE);
445  myaxis->SetTitleFont(42);
446  myaxis->SetLabelFont(42);
447  myaxis->SetLabelSize(0.06);
448  myaxis->SetTitleSize(0.06);
449 
450  myaxis->Draw();
451 
452  auto leg_D4 = new TLegend(0.6, 0.75, 0.9, 0.975, "", "brNDC");
453  // sprintf(lab, "Injections: %i", (int)mdc->GetEntries());
454  // leg_D3->AddEntry("", lab, "a");
455  sprintf(lab, "Ratio (i.e. Found / Injected)");
456  leg_D4->AddEntry(recx, lab, "p");
457  sprintf(lab, "Average Injected SNR");
458  leg_D4->AddEntry(snrx, lab, "p");
459  // leg_D3->SetFillColor(0);
460  leg_D4->SetFillColorAlpha(0, 0.9);
461  // leg_D4->Draw();
462  // rp1->GetLowYaxis()->SetNdivisions(505);
463  co_canvas3->Update();
464 
465  if (!opt.Contains("I")) {
466  rp1->GetLowerPad()->SetEditable(kFALSE);
467  co_canvas3->SaveAs(fname2);
468  filein->Close();
469  filein2->Close();
470  delete filein, filein2;
471  delete co_canvas3;
472  delete inj, rec, injx, recx, snrx, leg_D, leg_D3, leg_D4, rp1,
473  myaxis;
474  xi.clear(), xr.clear(), yi.clear(), yr.clear();
475  delete inj_gr, rec_gr, mg;
476  }
477 
478  // exit(0);
479 }
double rho
double T_ifar
TChain sim("waveburst")
static double g(double e)
Definition: GNGen.cc:116
float factor
TString opt
TString("c")
float CYS
TString mdc[4]
float SNR2
i drho i
void CreateDistanceParplots(char *sim_file_name, char *mdc_file_name, char *netdir, TString opt="", double MINX=0.0, double MAXX=1.0, double MAXDISTANCE=5000., int NBIN_DIST=10, float T_ifar=0.0, float T_win=0.2, int nIFO=2)
char fname3[2048]
char lab[256]
#define nIFO
i() int(T_cor *100))
char netdir[1024]
char fname[1024]
float mchirp
std::vector< double > iSNR[NIFO_MAX]
static void AddChip(TString filein, TString treename)
Definition: CBCTool.cc:1129
float spin[6]
char sim_file_name[1024]
double T_win
double mytime[6]
int cnt
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
#define MAXY
float mass[2]
char mdc_file_name[1024]
TMultiGraph * mg
#define MAXSNR
exit(0)