29 #include "Math/Rotation3D.h" 30 #include "Math/Vector3Dfwd.h" 33 #include "TGraphErrors.h" 38 #include "TMultiGraph.h" 40 #include "TRatioPlot.h" 41 #include "TRotation.h" 44 #include "TTreeIndex.h" 60 double MAXDISTANCE = 5000.,
65 TCanvas* co_canvas3 =
new TCanvas(
"sd3",
"SD3", 3, 47, 1000, 802);
66 co_canvas3->SetGridx();
67 co_canvas3->SetGridy();
68 co_canvas3->SetLogy();
70 float myifar, netcc[3];
85 TFile* filein =
new TFile(sim_file_name);
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");
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);
105 TFile* filein2 =
new TFile(mdc_file_name);
106 TTree*
mdc =
nullptr;
107 filein2->GetObject(
"mdc", mdc);
109 if (!mdc->GetListOfBranches()->FindObject(
"chip")) {
110 cout <<
"Adding Chi_p branch to mdc tree" << endl;
111 cbcTool.
AddChip(mdc_file_name,
"mdc");
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);
124 int nevts = (
int)mdc->GetEntries();
125 float CYS = 86400. * 365.25;
127 cout << nevts <<
" injected signals " << sim->GetEntries(
"ifar>0")
128 <<
" recovered signals" << endl;
132 std::vector<double> xi, xr, yi, yr;
134 auto inj =
new TH2F(
"Injected snr vs stat inj",
"", NBIN_DIST, MINX,
138 TMultiGraph*
mg =
new TMultiGraph();
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");
163 cout <<
"Not a valid option! " 164 "opt={\"chip\",\"chirp\",\"eta\",\"mtot\", \"chieff\", \"iota\"}" 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);
185 TH2F* rec = (TH2F*)inj->Clone(
"Injected snr vs stat rec");
188 rec->SetMarkerColor(4);
190 mg->GetYaxis()->SetTitle(
"Sensitive Distance (Mpc)");
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);
206 for (
int g = 0;
g < (
int)mdc->GetEntries();
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);
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")) {
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);
232 inj->Fill(xi[xi.size() - 1], mSNR);
239 for (
int g = 0;
g < (
int)sim->GetEntries();
g++) {
244 if (myifar <=
T_ifar * CYS) {
251 if ((mytime[0] - mytime[
nIFO]) < -
T_win ||
252 (mytime[0] - mytime[nIFO]) > 2 *
T_win) {
257 for (
int i = 0;
i <
nIFO;
i++) {
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]));
266 }
else if (
opt.Contains(
"chip")) {
269 }
else if (
opt.Contains(
"chirp")) {
270 xr.push_back(chirp[0]);
272 }
else if (
opt.Contains(
"mtot")) {
273 xr.push_back(mass[1] + mass[0]);
275 }
else if (
opt.Contains(
"eta")) {
276 xr.push_back(mass[0] * mass[1] /
277 pow(mass[0] + mass[1], 2.0));
279 }
else if (
opt.Contains(
"iota")) {
280 xr.push_back(iota[1]);
282 }
else if (
opt.Contains(
"distance")) {
283 xr.push_back(range[1]);
286 rec->Fill(xr[xr.size() - 1], mSNR);
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());
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);
317 co_canvas3->SaveAs(fname);
318 co_canvas3->SetLogx(0);
320 TGraph* rec_gr =
new TGraph(xr.size(), &xr[0], &yr[0]);
321 TGraph* inj_gr =
new TGraph(xi.size(), &xi[0], &yi[0]);
325 if (!
opt.Contains(
"distance")) {
327 inj_gr->SetMarkerColor(2);
328 inj_gr->SetMarkerStyle(20);
329 inj_gr->SetMarkerSize(0.5);
331 rec_gr->SetMarkerColor(4);
332 rec_gr->SetMarkerStyle(20);
333 rec_gr->SetMarkerSize(0.5);
339 mg->GetXaxis()->SetLimits(MINX, MAXX);
340 mg->GetYaxis()->SetRangeUser(10., MAXDISTANCE);
344 co_canvas3->SaveAs(fname3);
345 injx = (TH1D*)inj->ProjectionX();
346 recx = (TH1D*)rec->ProjectionX();
347 snrx = (TH1D*)inj->ProfileX();
349 injx = (TH1D*)inj->ProjectionX();
350 recx = (TH1D*)rec->ProjectionX();
351 snrx = (TH1D*)inj->ProfileX();
356 injx->SetFillColorAlpha(kRed, 0.3);
357 injx->SetFillStyle(3004);
359 recx->SetFillColorAlpha(kBlue, 0.99);
360 recx->SetFillStyle(3001);
363 co_canvas3->SetLogx(0);
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);
374 rp1->SetSplitFraction(0.5);
376 rp1->GetUpperRefYaxis()->SetRangeUser(1.,
MAXY);
388 rp1->GetLowerRefYaxis()->SetTitle(
"ratio");
389 rp1->GetLowerRefYaxis()->CenterTitle(kTRUE);
390 rp1->GetUpperRefYaxis()->SetTitle(
"entries");
391 rp1->GetUpperRefYaxis()->CenterTitle(kTRUE);
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);
400 auto leg_D3 =
new TLegend(0.6, 0.8455, 0.8965, 0.94555,
"",
"brNDC");
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");
408 leg_D3->SetFillColorAlpha(0, 0.9);
414 rp1->GetLowerPad()->cd();
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 );
426 snrx->SetLineColor(kGreen + 2);
427 snrx->SetMarkerColor(kGreen + 2);
429 snrx->SetFillColorAlpha(kGreen + 2, 0.4);
430 snrx->SetFillStyle(3001);
431 snrx->SetMarkerStyle(20);
432 snrx->SetMarkerSize(0.5);
434 snrx->Draw(
"SAME PEBAR");
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);
452 auto leg_D4 =
new TLegend(0.6, 0.75, 0.9, 0.975,
"",
"brNDC");
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");
460 leg_D4->SetFillColorAlpha(0, 0.9);
463 co_canvas3->Update();
465 if (!
opt.Contains(
"I")) {
466 rp1->GetLowerPad()->SetEditable(kFALSE);
467 co_canvas3->SaveAs(fname2);
470 delete filein, filein2;
472 delete inj, rec, injx, recx, snrx, leg_D, leg_D3, leg_D4, rp1,
474 xi.clear(), xr.clear(), yi.clear(), yr.clear();
475 delete inj_gr, rec_gr,
mg;
static double g(double e)
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)
std::vector< double > iSNR[NIFO_MAX]
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)