Logo coherent WaveBurst  
Library Reference Guide
Logo
DrawSearchAreaPRC.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
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 // Draw PRC Search Area
21 // Note : this macro is used to generate the PRC report (cwb_report merge_label prc)
22 // Author : Gabriele Vedovato
23 
24 #include "TROOT.h"
25 #include "TSystem.h"
26 #include "TFile.h"
27 #include "TTree.h"
28 #include <fstream>
29 #include <iostream>
30 #include "TGraph.h"
31 #include "TMultiGraph.h"
32 #include "TCanvas.h"
33 #include "TH1F.h"
34 #include "TMath.h"
35 #include "TH1F.h"
36 #include "TH2F.h"
37 #include <TComplex.h>
38 #include <TStyle.h>
39 #include <TRandom.h>
40 #include "TVector3.h"
41 #include "TRotation.h"
42 #include "Math/Vector3Dfwd.h"
43 #include "Math/Rotation3D.h"
44 
45 TH1F* sa_hist;
46 TCanvas* sa_canvas;
47 
48 void
50  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {
51 
52  using namespace ROOT::Math;
53 
54  char wave_file_name[1024];
55  sprintf(wave_file_name,"merge/wave_%s.%s.root",data_label.Data(),merge_label.Data());
56 
57  TFile* ifile = TFile::Open(wave_file_name);
58  TTree* itree = (TTree *) gROOT->FindObject("waveburst");
59 
60  char cut[1024];
61  char tmp[1024];
62  sprintf(cut,"abs(time[0]-time[%d])<%f && netcc[%d]>%f && rho[%d]>%f",
63  nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
64  if(T_vED>0) {strcpy(tmp,cut);sprintf(cut,"%s && neted[0]/ecor<%f",tmp,T_vED);}
65  if(T_pen>0) {strcpy(tmp,cut);sprintf(cut,"%s && penalty>%f",tmp,T_pen);}
66  if(T_ifar>0) {strcpy(tmp,cut);sprintf(cut,"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}
67 
68  itree->Draw("erA[0]",cut,"goff");
69  int size = itree->GetSelectedRows();
70  cout << "detected events : " << size << endl;
71 
72  double* erA0 = itree->GetV1();
73  Int_t *index = new Int_t[size];
74 
75  TMath::Sort(size,itree->GetV1(),index,false);
76 
77  double area_min = pow(erA0[index[0]],2);
78  double area_max = pow(erA0[index[size-1]],2);
79  area_min = area_min ? area_min : 0.1;
80  int nbin = TMath::Ceil(area_max/area_min);
81 
82  sa_hist = new TH1F("sa_hist","sa_hist",nbin,area_min,area_max);
83 
84  for(int i=0;i<size;i++) {
85  double sarea = erA0[index[i]]*erA0[index[i]];
86  sa_hist->Fill(sarea);
87  }
88 
89  // Normalization
90  double integral = 0;
91  for (int i=0;i<=sa_hist->GetNbinsX();i++) integral+=sa_hist->GetBinContent(i);
92  for (int i=0;i<=sa_hist->GetNbinsX();i++) sa_hist->SetBinContent(i,sa_hist->GetBinContent(i)/integral);
93  for (int i=2;i<=sa_hist->GetNbinsX();i++) sa_hist->SetBinContent(i,sa_hist->GetBinContent(i)+sa_hist->GetBinContent(i-1));
94  // find area at 50%
95  double search_area50=0;
96  for (int i=0;i<=sa_hist->GetNbinsX();i++) if(sa_hist->GetBinContent(i)>0.5) {search_area50=sa_hist->GetBinCenter(i);break;}
97  cout << "search_area50 : " << search_area50 << endl;
98  // find area at 90%
99  double search_area90=0;
100  for (int i=0;i<=sa_hist->GetNbinsX();i++) if(sa_hist->GetBinContent(i)>0.9) {search_area90=sa_hist->GetBinCenter(i);break;}
101  cout << "search_area90 : " << search_area90 << endl;
102 
103 
104  sa_canvas = new TCanvas("fom", "PRC", 300,40, 600, 600);
105  sa_canvas->SetGridx();
106  sa_canvas->SetGridy();
107  sa_canvas->SetLogx();
108 
109  // dump search area
110 
111  char ofname[1024];
112  ofstream out;
113  sprintf(ofname,"%s/search_area.txt",odir.Data());
114  out.open(ofname,ios::out);
115  if (!out.good()) {cout << "DrawSearchAreaPRC.C: Error Opening Output File : " << ofname << endl;exit(1);}
116  cout << "Create Dump Search Area Output File : " << ofname << endl;
117  int hsize = sa_hist->GetNbinsX()+1;
118  double* fraction = new double[hsize];
119  double* perc = new double[hsize];
120  double* eperc = new double[hsize];
121  char ostring[256];
122  for (int i=0;i<=sa_hist->GetNbinsX();i++) {
123  // binomial error
124  double N = integral;
125  double k = integral*sa_hist->GetBinContent(i);
126  if(k>N) k=N;
127  fraction[i] = sa_hist->GetBinCenter(i);
128  perc[i] = sa_hist->GetBinContent(i);
129  eperc[i] = (1/N)*sqrt(k*(1-k/N));
130  sprintf(ostring,"%.6f\t%.6f\t%.6f\t%d\n",sa_hist->GetBinCenter(i), perc[i], eperc[i], int(N));
131  out << ostring;
132  }
133  out.close();
134 
135  TGraphErrors* sa_graph = new TGraphErrors(hsize, fraction, perc, 0, eperc);
136  sa_graph->GetXaxis()->SetTitleOffset(1.3);
137  sa_graph->GetYaxis()->SetTitleOffset(1.25);
138  sa_graph->GetXaxis()->SetTickLength(0.01);
139  sa_graph->GetYaxis()->SetTickLength(0.01);
140  sa_graph->GetXaxis()->CenterTitle(kTRUE);
141  sa_graph->GetYaxis()->CenterTitle(kTRUE);
142  sa_graph->GetXaxis()->SetTitleFont(42);
143  sa_graph->GetXaxis()->SetLabelFont(42);
144  sa_graph->GetYaxis()->SetTitleFont(42);
145  sa_graph->GetYaxis()->SetLabelFont(42);
146  sa_graph->SetMarkerStyle(20);
147  sa_graph->SetMarkerSize(0.2);
148  sa_graph->SetMarkerColor(1);
149  sa_graph->SetLineColor(kOrange);
150  sa_graph->SetLineWidth(3);
151  sa_graph->SetTitle("");
152  sa_graph->SetFillColor(kOrange);
153  sa_graph->SetFillStyle(3001);
154 
155  char title[256];
156  sprintf(title,"Search Area - Ndet = %d",size);
157  sa_graph->SetTitle(title);
158  sa_graph->Paint("alp");
159  sa_graph->GetXaxis()->SetTitle("Searched area deg^{2}");
160  sa_graph->GetYaxis()->SetTitle("Cumulative fraction of events");
161  sa_graph->GetXaxis()->SetLimits(1, 10000);
162  sa_graph->GetYaxis()->SetRangeUser(0, 1);
163  sa_graph->GetXaxis()->SetTitleOffset(1.3);
164  sa_graph->GetYaxis()->SetTitleOffset(1.25);
165  sa_graph->GetXaxis()->SetTickLength(0.01);
166  sa_graph->GetYaxis()->SetTickLength(0.01);
167  sa_graph->GetXaxis()->CenterTitle(kTRUE);
168  sa_graph->GetYaxis()->CenterTitle(kTRUE);
169  sa_graph->GetXaxis()->SetTitleFont(132);
170  sa_graph->GetXaxis()->SetLabelFont(132);
171  sa_graph->GetYaxis()->SetTitleFont(132);
172  sa_graph->GetYaxis()->SetLabelFont(132);
173 
174  sa_graph->Draw("aple3");
175 
176  // draw the legend
177  TLegend *legend = new TLegend(0.1, 0.75, 0.5, 0.9,NULL,"brNDC");
178  legend->SetBorderSize(1);
179  legend->SetTextFont(132);
180  legend->SetTextSize(0.03);
181  legend->SetLineColor(kBlack);
182  legend->SetFillColor(kWhite);
183  char label[256];
184  sprintf(label,"50%% : %0.2f (deg^{2})",search_area50);
185  legend->AddEntry(sa_graph,label,"");
186  sprintf(label,"90%% : %0.2f (deg^{2})",search_area90);
187  legend->AddEntry(sa_graph,label,"");
188  legend->Draw();
189 
190  sprintf(ofname,"%s/search_area.gif",odir.Data());
191  sa_canvas->Update();
192  sa_canvas->SaveAs(ofname);
193  char cmd[1024];
194  TString pfname(ofname);
195  pfname.ReplaceAll(".gif",".png");
196  sprintf(cmd,"convert %s %s",ofname,pfname.Data());
197  cout << cmd << endl;
198  gSystem->Exec(cmd);
199  sprintf(cmd,"rm %s",ofname);
200  cout << cmd << endl;
201  gSystem->Exec(cmd);
202 
203  delete sa_graph;
204  delete [] fraction;
205  delete [] perc;
206  delete [] eperc;
207 }
208 
double T_ifar
double T_pen
TString("c")
double T_cor
ofstream out
Definition: cwb_merge.C:214
char odir[1024]
i pp_inetcc
Long_t size
i drho i
#define N
#define nIFO
char data_label[512]
Definition: test_config1.C:160
i() int(T_cor *100))
TString label
Definition: MergeTrees.C:21
double * tmp
Definition: testWDM_5.C:31
char cut[512]
void DrawSearchAreaPRC(TString data_label, TString odir, TString merge_label, int nIFO, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar)
TString merge_label
int k
TFile * ifile
TCanvas * sa_canvas
char title[256]
Definition: SSeriesExample.C:1
wavearray< int > index
double T_win
char cmd[1024]
strcpy(RunLabel, RUN_LABEL)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double T_cut
TTree * itree
TH1F * sa_hist
double T_vED
i drho pp_irho
exit(0)