Logo coherent WaveBurst  
Library Reference Guide
Logo
DrawCoverageVsPercentagePRC.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 Error Regions Percentage vs Coverage
21 // Note : this macro is used to generate the PRC report (cwb_report merge_label prc)
22 // Author : Gabriele Vedovato
23 
24 
25 #define NMDC 1
26 
27 //#define DUMP_ASCII
28 
29 double GetPercentage(TString gtype, int iderA, TString fname, int nIFO, float T_win, int pp_inetcc,
30  float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar);
31 
32 void
34  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {
35 
36  if(gtype!="nre" && gtype!="prc" && gtype!="fre") {
37  cout << "DrawCoverageVsPercentagePRC : Error - wrong input gtype (prc/nre/fre)" << endl;
38  gSystem->Exit(1);
39  }
40 
41 
42  char wave_file_name[1024];
43  sprintf(wave_file_name,"merge/wave_%s.%s.root",data_label.Data(),merge_label.Data());
44 
45  TObjArray* token = TString(wave_file_name).Tokenize(TString('/'));
46  TObjString* sfile = (TObjString*)token->At(token->GetEntries()-1);
47  TString ofile = sfile->GetString();
48  ofile.ReplaceAll(".root",".txt");
49  //TString TITLE = sfile->GetString();
50  //TITLE.ReplaceAll(".root","");
51  TString TITLE = "";
52  if(gtype=="prc") TITLE = "pp-plot-prc";
53  if(gtype=="nre") TITLE = "pp-plot-nre";
54  if(gtype=="fre") TITLE = "pp-plot-fre";
55 
56 #ifdef DUMP_ASCII
57  char ofname[1024];
58  if(gtype=="prc") sprintf(ofname,"%s/pp_plot_prc.txt",odir.Data());
59  if(gtype=="nre") sprintf(ofname,"%s/pp_plot_nre.txt",odir.Data());
60  if(gtype=="fre") sprintf(ofname,"%s/pp_plot_fre.txt",odir.Data());
61  ofstream out;
62  out.open(ofname,ios::out);
63  if (!out.good()) {cout << "Error Opening File : " << ofname << endl;exit(1);}
64  cout << "Create file : " << ofname << endl;
65 #endif
66 
67  // if coverage > percentage then erA is too large, the true erA is smaller (erA is conservative)
68  // if coverage < percentage then erA is too small, the true erA is greater (erA is optimistic)
69 
70  double coverage[NMDC][11];
71  for (int i=0;i<NMDC;i++) {
72  coverage[i][0]=0.;
73  coverage[i][10]=100.;
74  for (int j=1;j<10;j++) {
75  coverage[i][j]=GetPercentage(gtype, j,wave_file_name, nIFO, T_win, pp_inetcc,
76  T_cor, pp_irho, T_cut, T_vED, T_pen, T_ifar);
77  if(coverage[i][j]<0) return; // tree do not contains the branch (for example the nre erR array)
78 #ifdef DUMP_ASCII
79  out << j*10 << " " << coverage[i][j] << endl;
80 #endif
81  cout << j*10 << " " << coverage[i][j] << endl;
82  }
83  }
84 #ifdef DUMP_ASCII
85  out.close();
86 #endif
87 
88  TCanvas* canvas = new TCanvas("fom", "PRC", 300,40, 600, 600);
89  canvas->Clear();
90  canvas->ToggleEventStatus();
91  canvas->SetFillColor(0);
92  canvas->SetGridx();
93  canvas->SetGridy();
94  canvas->SetLogy(false);
95  canvas->SetLogx(false);
96 
97  gStyle->SetTitleH(0.062);
98  gStyle->SetTitleW(0.98);
99  gStyle->SetTitleY(0.98);
100  gStyle->SetTitleFont(72);
101  gStyle->SetLineColor(kWhite);
102  gStyle->SetPalette(1,0);
103  gStyle->SetNumberContours(256);
104 
105  Style_t markers[32]= {20,21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,
106  21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,20 };
107 
108  Color_t colors[32] = { 6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7,
109  6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7 };
110 
111  double perc[11]={0,10,20,30,40,50,60,70,80,90,100};
112  TGraph* gr[NMDC+11];
113  for(int i=0;i<NMDC;i++) {
114  gr[i] = new TGraph(11,perc,coverage[i]);
115  gr[i]->SetLineColor(colors[i]);
116  gr[i]->SetLineWidth(1);
117  gr[i]->SetMarkerSize(1);
118  gr[i]->SetMarkerColor(colors[i]);
119  gr[i]->SetMarkerStyle(markers[i]);
120  }
121  double x[2]={0,100};
122  double y[2]={0,100};
123  gr[NMDC+1] = new TGraph(2,x,y);
124  gr[NMDC+1]->SetLineColor(1);
125  gr[NMDC+1]->SetLineWidth(2);
126  gr[NMDC+1]->SetLineStyle(2);
127 
128  TMultiGraph* mg = new TMultiGraph();
129  for(int i=0;i<NMDC;i++) mg->Add(gr[i],"lp");
130  mg->Add(gr[NMDC+1],"lp");
131  mg->Paint("ap");
132  char title[256];
133  sprintf(title,"%s",TITLE.Data());
134  mg->GetHistogram()->SetTitle(title);
135  mg->GetHistogram()->GetXaxis()->SetTitle("Percentage");
136  mg->GetHistogram()->GetYaxis()->SetTitle("Coverage");
137  mg->GetHistogram()->GetXaxis()->SetRangeUser(0,100);
138  mg->GetHistogram()->GetYaxis()->SetRangeUser(0,100);
139  mg->GetHistogram()->GetXaxis()->SetTitleOffset(1.3);
140  mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.3);
141  mg->GetHistogram()->GetXaxis()->CenterTitle(true);
142  mg->GetHistogram()->GetYaxis()->CenterTitle(true);
143  mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
144  mg->GetHistogram()->GetXaxis()->SetTitleFont(42);
145  mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
146  mg->GetHistogram()->GetYaxis()->SetTitleFont(42);
147  mg->GetHistogram()->GetXaxis()->SetNdivisions(511);
148  mg->Draw("a");
149 
150  TLegend* leg;
151  double hleg = 0.8-NMDC*0.05;
152  leg = new TLegend(0.1291513,hleg,0.6244966,0.8738739,NULL,"brNDC");
153 
154  leg->SetBorderSize(1);
155  leg->SetTextAlign(22);
156  leg->SetTextFont(12);
157  leg->SetLineColor(1);
158  leg->SetLineStyle(1);
159  leg->SetLineWidth(1);
160  leg->SetFillColor(0);
161  leg->SetFillStyle(1001);
162  leg->SetTextSize(0.04);
163  leg->SetLineColor(kBlack);
164  leg->SetFillColor(kWhite);
165 
166  char label[256];
167  for(int i=0;i<NMDC;i++) {
168  sprintf(label,"%s ",data_label.Data());
169  leg->AddEntry(gr[i],label,"lp");
170  }
171 // leg->Draw();
172 
173  char ofname[1024];
174  if(gtype=="prc") sprintf(ofname,"%s/pp_plot_prc.gif",odir.Data());
175  if(gtype=="nre") sprintf(ofname,"%s/pp_plot_nre.gif",odir.Data());
176  if(gtype=="fre") sprintf(ofname,"%s/pp_plot_fre.gif",odir.Data());
177  TString gfileName=ofname;
178  canvas->Print(gfileName);
179  TString pfileName=gfileName;
180  pfileName.ReplaceAll(".gif",".png");
181  char cmd[1024];
182  sprintf(cmd,"convert %s %s",gfileName.Data(),pfileName.Data());
183  cout << cmd << endl;
184  gSystem->Exec(cmd);
185  sprintf(cmd,"rm %s",gfileName.Data());
186  cout << cmd << endl;
187  gSystem->Exec(cmd);
188  //exit(0);
189 }
190 
191 double GetPercentage(TString gtype, int iderA, TString fname, int nIFO, float T_win, int pp_inetcc,
192  float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {
193 
194  TFile *ifile = TFile::Open(fname.Data());
195  TTree* itree = (TTree *) gROOT->FindObject("waveburst");
196  itree->SetEstimate(itree->GetEntries());
197 
198  // check if erR exist
199  if(gtype=="nre") {
200  TString nre_name="erR";
201  TBranch* branch;
202  bool nre_exist=false;
203  TIter next(itree->GetListOfBranches());
204  while ((branch=(TBranch*)next())) {
205  if (nre_name.CompareTo(branch->GetName())==0) nre_exist=true;
206  }
207  if(!nre_exist) return -1;
208  }
209 
210  // check if erF exist
211  if(gtype=="fre") {
212  TString fre_name="erF";
213  TBranch* branch;
214  bool fre_exist=false;
215  TIter next(itree->GetListOfBranches());
216  while ((branch=(TBranch*)next())) {
217  if (fre_name.CompareTo(branch->GetName())==0) fre_exist=true;
218  }
219  if(!fre_exist) return -1;
220  }
221 
222  char selection[1024];
223  char cut[1024];
224  char tmp[1024];
225  sprintf(cut,"abs(time[0]-time[%d])<%f && netcc[%d]>%f && rho[%d]>%f",
226  nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
227  if(T_vED>0) {strcpy(tmp,cut);sprintf(cut,"%s && neted[0]/ecor<%f",tmp,T_vED);}
228  if(T_pen>0) {strcpy(tmp,cut);sprintf(cut,"%s && penalty>%f",tmp,T_pen);}
229  if(T_ifar>0) {strcpy(tmp,cut);sprintf(cut,"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}
230  //cout << cut << endl;
231 
232  if(gtype=="prc") sprintf(selection,"erA[0]-erA[%d]>>hist_cumulative_erA1(2000,-100,100)",iderA);
233  if(gtype=="nre") sprintf(selection,"erR[0]-erR[%d]>>hist_cumulative_erA1(2000,-1,1)",iderA);
234  if(gtype=="fre") sprintf(selection,"erF[0]-erF[%d]>>hist_cumulative_erA1(2000,-1,1)",iderA);
235 
236  itree->Draw(selection,cut,"goff");
237  TH1F *hist_cumulative_erA1 = (TH1F*)gDirectory->Get("hist_cumulative_erA1");
238  int size = itree->GetSelectedRows();
239  //cout << "size : " << size << endl;
240  if(size==0) {
241  delete hist_cumulative_erA1;
242  delete itree;
243  delete ifile;
244  return 0;
245  }
246 
247  double integral = hist_cumulative_erA1->ComputeIntegral();
248  if (integral==0) {cout << "Empty histogram !!!" << endl;exit(0);}
249  double* cumulative = hist_cumulative_erA1->GetIntegral();
250  for (int i=0;i<hist_cumulative_erA1->GetNbinsX();i++) hist_cumulative_erA1->SetBinContent(i,cumulative[i]/integral);
251  //cout << "integral " << integral << endl;
252 
253 
254  double perc = 100.*hist_cumulative_erA1->GetBinContent(1001);
255  delete hist_cumulative_erA1;
256  delete itree;
257  delete ifile;
258  //ifile->Close();
259 
260  return perc;
261 }
262 
char ofile[1024]
double T_ifar
double T_pen
TString("c")
#define TITLE
Definition: TestSiteJ1.C:44
double T_cor
ofstream out
Definition: cwb_merge.C:214
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ 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
char odir[1024]
return wmap canvas
i pp_inetcc
Long_t size
Color_t colors[16]
int j
Definition: cwb_net.C:28
i drho i
char data_label[512]
Definition: test_config1.C:160
TGraph * gr
void DrawCoverageVsPercentagePRC(TString gtype, 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 label
Definition: MergeTrees.C:21
double * tmp
Definition: testWDM_5.C:31
char cut[512]
TIter next(twave->GetListOfBranches())
char fname[1024]
TString merge_label
TObjArray * token
#define nIFO
#define NMDC
TFile * ifile
char title[256]
Definition: SSeriesExample.C:1
double T_win
TBranch * branch
double GetPercentage(TString gtype, int iderA, TString fname, 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)
char cmd[1024]
strcpy(RunLabel, RUN_LABEL)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double T_cut
TTree * itree
double T_vED
wavearray< double > y
Definition: Test10.C:31
TMultiGraph * mg
i drho pp_irho
exit(0)