Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Namespaces Files Functions Variables Macros Pages
Make_PP_IFAR.C
Go to the documentation of this file.
1 
2 // ---------------------------------------------------------------------------------
3 // INFOS
4 // ---------------------------------------------------------------------------------
5 
6 // this macro produce the plot number of events in the observation period vs IFAR
7 
8 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9 // The macro Make_PP_IFAR_GW151226_RATE.C is the standard macro in cWB config master.2.8 + add GW151226 event in O1 data
10 // set a superior limit @ GW170104 IFAR -> set a limit for BBH IFAR, ecluded BBH IFAR estimated as lower limits (IFAR_TAG)
11 // compute the index which contains the maximum common IFAR for rates -> limits the sim green plot (RATE_TAG)
12 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13 
14 // input list chunk root files format
15 //
16 // wave_merged_root_file_name_with_ifar.root lag slag chunk chunkLabel
17 // Ex: K02,K03 open box
18 // wave_O2_K02_merged_root_file_name_with_ifar.root 0 0 2 K02
19 // wave_O2_K03_merged_root_file_name_with_ifar.root 0 0 3 K03
20 
21 // ---------------------------------------------------------------------------------
22 // HOW RUN THIS MACRO (EXAMPLES)
23 // ---------------------------------------------------------------------------------
24 
25 /*
26 
27 // uses a list of input root merged files generated by different chunks
28 root 'Make_PP_IFAR.C("pp_lag0_slag0_M2.lst", "--search=BurstLF:bin2 \\
29  --exit=1 --pfname=pp_plot_bin2_kall_lag0_slag0_M2.png")'
30 
31 // uses a single chunk root merged file
32 root 'Make_PP_IFAR.C("wave_merged_output_root_file.root", "--search=BurstLF:bin2 --lag=0 --slag=0 \\
33  --chunk=0 --exit=1 --pfname=pp_plot_bin2_kall_lag0_slag0_M2.png")'
34 
35 */
36 
37 // ---------------------------------------------------------------------------------
38 // USER DEFINES : Default values used to initialize the user options input parameters
39 // ---------------------------------------------------------------------------------`
40 
41 #define PP_RUN "O2"
42 #define PP_SEARCH ""
43 
44 #define PP_NIFO 2
45 #define PP_LAG -1
46 #define PP_SLAG -1
47 #define PP_CHUNK -1
48 
49 #define PP_TRIALS 1
50 
51 #define PP_PFNAME ""
52 #define PP_ZL_STYLE "step"
53 #define PP_BELT_STYLE "continuos"
54 
55 #define PP_IFAR_MAX -1
56 #define PP_RHO_MIN 4.5
57 
58 #define PP_SFACTOR 1
59 
60 #define PP_BBH false
61 #define PP_EXIT true
62 
63 #define PP_GW151226_IFAR 0
64 #define PP_ZL_MAX_IFAR 0
65 
66 // ---------------------------------------------------------------------------------`
67 // DEFINES
68 // ---------------------------------------------------------------------------------`
69 
70 #define MAX_RUNS 2 // O1,O2
71 
72 #define MAX_LIST_SIZE 100 // maximum number of input root files in the input file list
73 #define IFAR_NBIN 2000 // number of bins used for the computation of belts and background
74  // NOTE: if such value is low it introduces some discretization artifacts for low IFAR values
75 
76 //#define FAP0 0.1 // FAP used to define the first belt : SIGMA = 1.644
77 //#define FAP1 0.01 // FAP used to define the second belt : SIGMA = 2.575
78 //#define FAP2 0.001 // FAP used to define the third belt : SIGMA = 3.290
79 
80 //
81 // Probability = TMath::Erf(Sigma/sqrt(2))
82 // Sigma = TMath::ErfInverse(Probability)*sqrt(2)
83 
84 #define FAP0 1-0.682689 // FAP used to define the first belt : SIGMA = 1.
85 #define FAP1 1-0.954499 // FAP used to define the second belt : SIGMA = 2.
86 #define FAP2 1-0.997300 // FAP used to define the third belt : SIGMA = 3.
87 
88 #define DAY (24.*3600.)
89 #define YEAR (24.*3600.*365.)
90 
91 #define CHUNK_FILE_LIST "Chunk_List.txt"
92 #define MACRO_READ_CHUNK "ReadChunkList.C"
93 #define CHUNK_MAX_SIZE 100
94 
95 //#define APPLY_MAX_SIM_IFAR // if enabled the detected & sim events are cut to the common superior limit
96 
97 //#define PLOT_FOREGROUND_NOBBH // plot ZL BBH && ZL NOBBH
98 
99 // ---------------------------------------------------------------------------------
100 // INPUT USER OPTIONS
101 // ---------------------------------------------------------------------------------
102 
103 struct uoptions {
104  TString run; // run type: allowed types are O1/O2/O1O2
105  TString search; // search type: allowed types are BBH,BurstLF,BurstHF,BurstLD,IMBHB + bin1/bin2/...
106 
107  int nifo; // number of detectors
108  int lag; // lag -> used only when the input file name is a root file
109  int slag; // slag -> used only when the input file name is a root file
110  int chunk; // chunk -> used only when the input file name is a root file
111 
112  int trials; // trials factor
113 
114  TString pfname; // output plot file name, included the directory
115  TString zl_style; // zero lag plot style: step/marker -> zero lag is displayed with step/marker
116  TString belt_style; // belt plot style: step/continuos -> belts are displayed with step/continuos poisson distribution
117 
118  double ifar_max; // max ifar (year) used for plot, if =-1 then it is the maximum ifar from root files
119  double rho_min; // min rho[0/1] read from output root files
120 
121  double sfactor; // factor used to normalize the simulation events for the computation of astrophysical rates
122  // SIM_NEVT/=gOPT.sfactor
123  bool bbh; // if false the known bbh events are excluded from analysis
124  bool exit; // if true program exit ayt the end of the execution
125 
126  double gw151226_ifar; // GW151226 IFAR(YEARS): if > 0 the GW151226 event is added to the foreground (def=0, not added)
127  // This is a special option introduced for the GW151226 event which has been found only into the extended segments
128 
129  double zl_max_ifar; // It is used to set the maximum IFAR of foreground events (YEARS).
130  // The value to be used is the minimum common foreground IFAR for which the IFAR is well estimated
131  // The foreground with IFAR>zl_max_ifar is set to zl_max_ifar
132 };
133 
134 
135 // ---------------------------------------------------------------------------------
136 // Global Variables
137 // ---------------------------------------------------------------------------------
138 
139 // global User Options
140 
141 uoptions gOPT;
142 
143 // plot objects
144 
145 TCanvas* canvas;
146 TGraph* gr;
147 TGraph* grsim;
148 TGraph* grsimbkg;
149 TGraphAsymmErrors* egrsimbkg0;
150 TGraph* gr0;
151 TLegend* leg0;
152 TLegend* leg;
153 TMarker* marker;
154 TArrow* arrow;
155 TLatex *latex;
156 TGraphAsymmErrors* egr0;
157 TGraphAsymmErrors* egr1;
158 TGraphAsymmErrors* egr2;
159 
160 // wavearray used for plots
161 
162 wavearray<double> xIFAR;
163 wavearray<double> xFAR;
164 wavearray<double> bRATE;
165 wavearray<double> yBKG;
166 
167 wavearray<double> SIM_IFAR; // sorted ifar simulated events
168 wavearray<double> SIM_NEVT; // number of simulated events with ifar > SIM_IFAR[k]
169 wavearray<double> SIM_BKG_NEVT; // number of simulated events + expected from BKG with ifar > SIM_IFAR[k]
170 
171 // SIM+BKG FAP0 belts
172 wavearray<double> ySIM_BKG;
173 wavearray<double> exSIM_BKG_IFAR;
174 wavearray<double> xSIM_BKG_IFAR;
175 wavearray<double> eySIM_BKG_inf0;
176 wavearray<double> eySIM_BKG_sup0;
177 
178 wavearray<double> ZL_IFAR; // sorted ifar events
179 wavearray<double> ZL_NEVT; // number of events with ifar > ZL_IFAR[k]
180 wavearray<double> ZL_GPS; // event GPS time
181 
182 // FAP0/FAP1/FAP2 error belts
183 wavearray<double> exFAR;
184 wavearray<double> eyRATEinf0;
185 wavearray<double> eyRATEsup0;
186 wavearray<double> eyRATEinf1;
187 wavearray<double> eyRATEsup1;
188 wavearray<double> eyRATEinf2;
189 wavearray<double> eyRATEsup2;
190 
191 // O1 known BBH events
192 vector<TString> O1_BBH_NAME;
193 vector<double> O1_BBH_GPS;
194 
195 // O2 known BBH events
196 vector<TString> O2_BBH_NAME;
197 vector<double> O2_BBH_GPS;
198 
199 // chunks id,start,stop
200 int chunk_size[MAX_RUNS];
201 int chunk_id[MAX_RUNS][CHUNK_MAX_SIZE];
202 double chunk_start[MAX_RUNS][CHUNK_MAX_SIZE];
203 double chunk_stop[MAX_RUNS][CHUNK_MAX_SIZE];
204 
205 // ---------------------------------------------------------------------------------
206 // FUNCTIONS
207 // ---------------------------------------------------------------------------------
208 
209 void ResetUserOptions();
210 void ReadUserOptions(TString options);
211 void PrintUserOptions();
212 void DumpUserOptions();
213 
214 void InitBBH(); // init array Ox_BBH_NAME,Ox_BBH_GPS with known BBH events
215 TString GetNameBBH(double gps); // return the BBH name
216 double GetGPSBBH(TString name); // return the BBH GPS
217 TString GetCutBBH(); // return the cut string used to exclude BBH events from tree selection
218 
219 int GetChunkID(double gps, TString& run); // return chunk ID for a gine gps time
220 int GetPercentiles(double mu, double* q); // return percentiles in q for FAP0,FAP1,FAP2 (continuos)
221 int GetPoissonNsup(double mu, double fap); // return sup number expeted events of for a gine fap
222 int GetPoissonNinf(double mu, double fap); // return inf number expeted events of for a gine fap
223 int ReadFileList(TString ifName, TChain** live, TChain** wave, TChain** mdc, int* lag, int* slag, int* chunk, int* bin, TString* run); // read input root file list
224 void GetLiveTime(int nwave, TChain** live, int* lag, int* slag, int* bin, TString* run, double* liveZero, double* liveTot); // retun live time
225 
226 void MakePlot(TString title, double obsTime, double ifar_cmax); // make final IFAR plot
227 void DumpLoudest(double obsTime); // dump to ascii file the loudest zero lag event list
228 void DumpPeriod(double* obsTime);
229 void DumpBackground(); // Dumps Background & Belts in ASCII format
230 void DumpForeground(); // Dumps Foreground in ASCII format
231 void DumpSimulation(double ifar_cmax); // Dumps Simulation & Belts in ASCII format
232 
233 
234 void Make_PP_IFAR(TString ifName, TString options) {
235 
236  // -----------------------------------------------------
237  // Read Input Options
238  // -----------------------------------------------------
239 
240  ResetUserOptions();
241  ReadUserOptions(options);
242  PrintUserOptions();
243 
244  if((gOPT.run!="O1")&&(gOPT.run!="O2")&&(gOPT.run!="O1O2")) {
245  cout << "Make_PP_IFAR - Error : run option not allowed, must be O1/O2/O1O2" << endl;
246  gSystem->Exit(1);
247  }
248 
249  TString SEARCH = gOPT.search;
250  if(gOPT.search.Contains(":")) SEARCH = gOPT.search(0,gOPT.search.Index(":",0));
251  if((SEARCH!="BBH")&&(SEARCH!="IMBHB")&&(SEARCH!="BurstLF")&&(SEARCH!="BurstHF")&&(SEARCH!="BurstLD")) {
252  cout << "Make_PP_IFAR - Error : search option not allowed, must be BBH/IMBHB/BurstLF/BurstHF/BurstLD + ':bin1/bin2/...'" << endl;
253  gSystem->Exit(1);
254  }
255 
256  if(gOPT.trials<=0) {
257  cout << "Make_PP_IFAR - Error : trials value not allowed, must be >0" << endl;
258  gSystem->Exit(1);
259  }
260 
261  if(gOPT.sfactor<=0) {
262  cout << "Make_PP_SFACTOR - Error : sfactor value not allowed, must be >0" << endl;
263  gSystem->Exit(1);
264  }
265 
266  if(gOPT.nifo!=2 && gOPT.nifo!=3) {
267  cout << "Make_PP_IFAR - Error : nifo value not allowed, must 2/3" << endl;
268  gSystem->Exit(1);
269  }
270 
271  if((gOPT.zl_style!="step")&&(gOPT.zl_style!="marker")) {
272  cout << "Make_PP_IFAR - Error : zl_style value not allowed, must be step/marker" << endl;
273  gSystem->Exit(1);
274  }
275 
276  if((gOPT.belt_style!="step")&&(gOPT.belt_style!="continuos")) {
277  cout << "Make_PP_IFAR - Error : belt_style value not allowed, must be step/continuos" << endl;
278  gSystem->Exit(1);
279  }
280 
281  // check if is a root file
282  TString fext = ifName(ifName.Last('.')+1,4);
283  if(fext=="root") {
284  if(gOPT.lag<0 || gOPT.slag<0 || gOPT.chunk<0) {
285  cout << "Make_PP_IFAR - Error : when input file ext=root lag,slag,chunk must be >=0" << endl;
286  gSystem->Exit(1);
287  }
288  }
289 
290  DumpUserOptions();
291 
292  if(gOPT.exit) gROOT->SetBatch(kTRUE); // NOTE: batch mode save plot with high quality !!!
293 
294  // initialize wavearray
295  xIFAR.resize(IFAR_NBIN);
296  xFAR.resize(IFAR_NBIN);
297  bRATE.resize(IFAR_NBIN);
298  yBKG.resize(IFAR_NBIN);
299 
300  exFAR.resize(IFAR_NBIN); exFAR=0;
301  eyRATEinf0.resize(IFAR_NBIN);
302  eyRATEsup0.resize(IFAR_NBIN);
303  eyRATEinf1.resize(IFAR_NBIN);
304  eyRATEsup1.resize(IFAR_NBIN);
305  eyRATEinf2.resize(IFAR_NBIN);
306  eyRATEsup2.resize(IFAR_NBIN);
307 
308  ySIM_BKG.resize(IFAR_NBIN);
309  exSIM_BKG_IFAR.resize(IFAR_NBIN); exSIM_BKG_IFAR=0;
310  xSIM_BKG_IFAR.resize(IFAR_NBIN);
311  eySIM_BKG_inf0.resize(IFAR_NBIN);
312  eySIM_BKG_sup0.resize(IFAR_NBIN);
313 
314  InitBBH(); // init known BBH
315 
316 
317  // -----------------------------------------------------
318  // Read Chunk start stop GPS times
319  // -----------------------------------------------------
320 
321  // get CWB_CONFIG
322  char cwb_config_env[1024] = "";
323  if(gSystem->Getenv("CWB_CONFIG")!=NULL) {
324  strcpy(cwb_config_env,TString(gSystem->Getenv("CWB_CONFIG")).Data());
325  }
326 
327  for(int n=0;n<MAX_RUNS;n++) {
328 
329  TString RUN;
330 
331  if(n==0 && gOPT.run.Contains("O1")) RUN="O1";
332  else
333  if(n==1 && gOPT.run.Contains("O2")) RUN="O2";
334  else continue;
335 
336  char chunk_file_list[1024];
337  sprintf(chunk_file_list,"%s/%s/CHUNKS/%s",cwb_config_env,RUN.Data(),CHUNK_FILE_LIST);
338  cout << chunk_file_list << endl;
339 
340  char macro_read_chunk_path[1024];
341  sprintf(macro_read_chunk_path,"%s/MACROS/%s",cwb_config_env,MACRO_READ_CHUNK);
342 
343  CWB::Toolbox::checkFile(macro_read_chunk_path);
344 
345  // load macro
346  gROOT->LoadMacro(gSystem->ExpandPathName(macro_read_chunk_path));
347 
348  cout << endl;
349  cout << "----------------------------------------------------------------------------------------------------" << endl;
350  cout << "chunk" << "\t\t" << "start" << "\t\t" << "stop" << "\t\t"
351  << " days\t\t" << "beg-date" << " - " << "end-date" << "\t" << "run" << endl;
352  cout << "----------------------------------------------------------------------------------------------------" << endl;
353  cout << endl;
354 
355  chunk_size[n] = ReadChunkList(chunk_file_list, chunk_id[n], chunk_start[n], chunk_stop[n]);
356 
357  cout << "----------------------------------------------------------------------------------------------------" << endl;
358  cout << endl;
359  }
360 
361  // -----------------------------------------------------
362  // open wave root trees
363  // -----------------------------------------------------
364 
365  int rho_id;
366  if(SEARCH=="BBH" || SEARCH=="IMBHB") rho_id=1;
367  else rho_id=0;
368 
369  TString run[MAX_LIST_SIZE];
370  int bin[MAX_LIST_SIZE];
371  int chunk[MAX_LIST_SIZE];
372  int lag[MAX_LIST_SIZE];
373  int slag[MAX_LIST_SIZE];
374 
375  TChain* live[MAX_LIST_SIZE];
376  for(int i=0;i<MAX_LIST_SIZE;i++) live[i] = new TChain("liveTime");
377 
378  TChain* wave[MAX_LIST_SIZE];
379  for(int i=0;i<MAX_LIST_SIZE;i++) wave[i] = new TChain("waveburst");
380 
381  TChain* mdc[MAX_LIST_SIZE];
382  for(int i=0;i<MAX_LIST_SIZE;i++) mdc[i] = new TChain("mdc");
383 
384  int nwave = 1;
385  if(fext=="root") {
386  chunk[0] = gOPT.chunk;
387  lag[0] = gOPT.lag;
388  slag[0] = gOPT.slag;
389  run[0] = gOPT.run;
390  bin[0] = 0;
391  TString fwave = ifName;
392  TString flive = ifName;
393  flive.ReplaceAll("wave_","live_");
394  wave[0]->Add(fwave);
395  live[0]->Add(flive);
396  wave[0]->SetTitle(fwave);
397  live[0]->SetTitle(flive);
398  } else {
399  // Read File List & Fill TChain
400  nwave = ReadFileList(ifName, live, wave, mdc, lag, slag, chunk, bin, run);
401  gOPT.lag = lag[0];
402  gOPT.slag = slag[0];
403  gOPT.chunk = chunk[0];
404  }
405  cout << endl;
406  if(nwave<=0) {
407  cout << "Make_PP_IFAR - Error : no files read from input list, check input list format" << endl;
408  gSystem->Exit(1);
409  }
410 
411  // check if slag and lag are consistent
412  int xlag=lag[0];
413  int xslag=slag[0];
414  for(int i=1;i<nwave;i++) {
415  if(bin[i]<0) continue; // simulation
416  if(lag[i]!=xlag || slag[i]!=xslag) {
417  cout << "Make_PP_IFAR - Error : lag,slag defined in input list root files are not consistent" << endl;
418  gSystem->Exit(1);
419  }
420  }
421 
422  // check duplicated entries
423  for(int i=0;i<nwave;i++) {
424  if(bin[i]<0) continue; // simulation
425  for(int j=0;j<nwave;j++) {
426  if(bin[j]<0) continue; // simulation
427  if((i!=j) && (lag[i]==lag[j] && slag[i]==slag[j] && chunk[i]==chunk[j] && bin[i]==bin[j] && TString(wave[i]->GetTitle())==TString(wave[j]->GetTitle()))) {
428  cout << "Make_PP_IFAR - Error : input file list contains duplicated entries for lag,slag,chunk,bin = "
429  << lag[i] << "," << slag[i] << "," << chunk[i] << "," << bin[i] << endl;
430  cout << " : "
431  << wave[i]->GetTitle() << endl;
432  gSystem->Exit(1);
433  }
434  }
435  }
436 
437  // check consistency bin1 bin2 entries
438  vector<int> vbin; // contains the list of unique bins
439  for(int i=0;i<nwave;i++) {
440  if(bin[i]<0) continue; // simulation
441  bool exist=false;
442  for(int j=0;j<vbin.size();j++) if(bin[i]==vbin[j]) exist=true;
443  if(!exist) vbin.push_back(bin[i]);
444  }
445  for(int i=0;i<nwave;i++) {
446  if(bin[i]<0) continue; // simulation
447  int nbins=1;
448  for(int j=0;j<nwave;j++) {
449  if(bin[j]<0) continue; // simulation
450  if((i!=j) && (lag[i]==lag[j] && slag[i]==slag[j] && chunk[i]==chunk[j])) {
451  for(int k=0;k<vbin.size();k++) if(bin[j]!=vbin[k]) nbins++;
452  }
453  }
454  if(nbins!=vbin.size()) {
455  cout << "Make_PP_IFAR - Error : input file list must contains all bins for each chunk: "
456  << wave[i]->GetTitle() << endl;
457  cout << " lag,slag,chunk,bin = "
458  << lag[i] << "," << slag[i] << "," << chunk[i] << "," << bin[i] << endl;
459  gSystem->Exit(1);
460  }
461  }
462 
463  // -----------------------------------------------------
464  // compute live time
465  // -----------------------------------------------------
466 
467  double liveTot[MAX_RUNS];
468  double liveZero[MAX_RUNS];
469  for(int n=0;n<MAX_RUNS;n++) {liveTot[n]=0;liveZero[n]=0;}
470 
471  GetLiveTime(nwave, live, lag, slag, bin, run, liveZero, liveTot);
472 
473  double obsTime=0;
474  for(int n=0;n<MAX_RUNS;n++) {
475 
476  obsTime+=liveTot[n]/gOPT.trials; // apply trials factor
477 
478  TString RUN;
479  if(n==0) RUN="O1";
480  if(n==1) RUN="O2";
481  printf("%s : total live time: non-zero lags = %10.1f sec, %10.2f days, %10.4f years \n\n",RUN.Data(),liveTot[n],liveTot[n]/DAY,liveTot[n]/YEAR);
482  }
483 
484  printf("%s : total observational time: non-zero lags = %10.1f sec, %10.2f days, %10.4f years \n\n",gOPT.run.Data(),obsTime,obsTime/DAY,obsTime/YEAR);
485 
486 
487  // =======================================================================
488  // DETECTED EVENTS STUFF
489  // =======================================================================
490 
491  // -----------------------------------------------------
492  // Get BBH cut string
493  // -----------------------------------------------------
494 
495  TString CutBBH = GetCutBBH();
496 
497  // -----------------------------------------------------
498  // compute min.max ifar in wave files
499  // -----------------------------------------------------
500 
501  double ifar_cmin = 0; // minimum common IFAR
502  double ifar_cmax = 1e100; // maximum common IFAR (used for rate plots) (RATE_TAG)
503  double ifar_max = 0;
504  double value;
505  for(int m=0;m<nwave;m++) {
506  if(bin[m]<0) continue; // simulation
507  int isize = (int)wave[m]->GetEntries();
508  wave[m]->SetEstimate(isize);
509  char cut[64]; sprintf(cut,"rho[%d]>%f",rho_id,gOPT.rho_min);
510  wave[m]->Draw("ifar",cut,"goff",isize);
511  double* ifar = wave[m]->GetV1();
512  int sel_entries = wave[m]->GetSelectedRows();
513  value = log10(TMath::MinElement(sel_entries,ifar));
514  if(fext=="root")
515  cout << "tMax events per year : " << (1./pow(10,value))*YEAR << endl;
516  else
517  cout << "chunk\t" << chunk[m] << "\tbin\t" << bin[m] << "\tMax events per year : " << (1./pow(10,value))*YEAR << endl;
518  if(value>ifar_cmin) ifar_cmin=value; // this condition is used to select the common IFAR inferior limit
519  value = log10(TMath::MaxElement(sel_entries,ifar));
520  if(value>ifar_max) ifar_max=value; // this condition is used to select the maximum IFAR limit
521  if(value<ifar_cmax) ifar_cmax=value; // this condition is used to select the common IFAR superior limit (RATE_TAG)
522  }
523 
524  cout << endl;
525  cout << "ifar_cmin : " << pow(10,ifar_cmin)/YEAR << " years" << endl;
526  cout << "ifar_cmax : " << pow(10,ifar_cmax)/YEAR << " years" << endl; // (RATE_TAG)
527  cout << "ifar_max : " << pow(10,ifar_max)/YEAR << " years" << endl;
528 
529  double inf = ifar_cmin;
530  double sup = TMath::Nint(ifar_max+0.5);
531 
532  TString Cuts="";
533 
534  // -----------------------------------------------------
535  // read zero ifar lag events
536  // -----------------------------------------------------
537 
538  char sel[1024];
539  char tmp[1024];
540  cout << endl;
541  vector<double> vIFAR;
542  vector<double> vGPS;
543  int nsel=0;
544  for(int m=0;m<nwave;m++) {
545  if(bin[m]<0) continue; // simulation
546  Cuts = CutBBH;
547  // ------------------------> APPLY TRIALS FACTOR !!!
548  sprintf(sel,"ifar/%f:time[0]",gOPT.trials);
549  sprintf(tmp,"ifar>%f",pow(10,ifar_cmin)*gOPT.trials);
550  if(lag[m]==-1 && slag[m]==-1) sprintf(tmp,"ifar>%f && lag[%d]!=0 && rho[%d]>%f",pow(10,ifar_cmin)*gOPT.trials,gOPT.nifo,rho_id,gOPT.rho_min);
551  if(lag[m]>= 0 && slag[m]==-1) sprintf(tmp,"ifar>%f && lag[%d]==%d && rho[%d]>%f",pow(10,ifar_cmin)*gOPT.trials,gOPT.nifo,lag[m],rho_id,gOPT.rho_min);
552  if(lag[m]>= 0 && slag[m]>= 0) sprintf(tmp,"ifar>%f && lag[%d]==%d && slag[%d]==%d && rho[%d]>%f",pow(10,ifar_cmin)*gOPT.trials,gOPT.nifo,lag[m],gOPT.nifo,slag[m],rho_id,gOPT.rho_min);
553  if(gOPT.bbh) {
554  Cuts = TString::Format("%s ",tmp);
555  } else {
556  Cuts += TString::Format(" && %s ",tmp);
557  }
558  cout << "CutBBH : " << Cuts << endl;
559  wave[m]->SetEstimate(wave[m]->GetEntries());
560  wave[m]->Draw(sel,Cuts.Data(),"goff");
561  nsel+= wave[m]->GetSelectedRows();
562  double* V1 = wave[m]->GetV1();
563  double* V2 = wave[m]->GetV2();
564  for(int k=0;k<wave[m]->GetSelectedRows();k++) {vIFAR.push_back(V1[k]);vGPS.push_back(V2[k]);}
565  }
566 
567  // IF gOPT.gw151226_ifar>0 -> ADD GW151226 Event !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
568  // This is a special option introduced for the GW151226 event which has been found only into the extended segments
569  if(gOPT.bbh && gOPT.gw151226_ifar>0) {
570  vIFAR.push_back(gOPT.gw151226_ifar*YEAR);vGPS.push_back(GetGPSBBH("GW151226"));
571  }
572 
573  // -----------------------------------------------------
574  // sort ifar zero lag events
575  // -----------------------------------------------------
576 
577  if(vIFAR.size()>0) {
578  int size = vIFAR.size();
579  int *index = new int[size];
580  double *IFAR = new double[size];
581  for(int k=0;k<size;k++) IFAR[k]=vIFAR[k];
582  TMath::Sort(size,IFAR,index,true);
583  if(gOPT.zl_style=="marker") {
584  ZL_GPS.resize(size);
585  ZL_IFAR.resize(size);
586  ZL_NEVT.resize(size);
587  for(int k=0;k<size;k++) { // invert the ifar order
588  ZL_GPS[size-k-1] = vGPS[index[k]]; // event GPS time
589  ZL_IFAR[size-k-1] = IFAR[index[k]]/YEAR; // sorted ifar events
590  ZL_NEVT[size-k-1] = k+1; // number of events with ifar > ZL_IFAR[k]
591  }
592  } else { // step
593  ZL_GPS.resize(2*size-1);
594  ZL_IFAR.resize(2*size-1);
595  ZL_NEVT.resize(2*size-1);
596  int n=0;
597  for(int k=size-1;k>=0;k--) { // invert the ifar order
598  ZL_GPS[n] = vGPS[index[k]]; // event GPS
599  ZL_IFAR[n] = IFAR[index[k]]/YEAR; // sorted ifar events
600  ZL_NEVT[n] = k+1; // number of events with ifar > ZL_IFAR[k]
601  n++;
602  if(n>=ZL_NEVT.size()-1) break;
603  // add point to make ZL STYLE STEPWISE
604  ZL_GPS[n] = ZL_GPS[n-1];
605  ZL_IFAR[n] = ZL_IFAR[n-1];
606  ZL_NEVT[n] = ZL_NEVT[n-1]-1;
607  n++;
608  }
609  }
610  delete [] index;
611  delete [] IFAR;
612  // The foreground with IFAR>zl_max_ifar is set to zl_max_ifar
613  if(gOPT.zl_max_ifar>0) {
614  for(int k=0;k<ZL_NEVT.size();k++) if(ZL_IFAR[k]>gOPT.zl_max_ifar) ZL_IFAR[k]=gOPT.zl_max_ifar;
615  }
616  }
617 
618  // adjust sup limit according to the loudest zero lag event
619  double loudest_zl_ifar_sec = ZL_NEVT.size()>0 ? ZL_IFAR[ZL_NEVT.size()-1]*YEAR : 0;
620  if(log10(loudest_zl_ifar_sec)>sup) sup=log10(loudest_zl_ifar_sec);
621  if(inf>0 && sup/inf<2) sup=2*inf;
622  else sup*=1.01;
623 
624  if(pow(10, inf)/YEAR>0.01) inf=log10(0.01*YEAR);
625 
626  if(gOPT.ifar_max>0) sup = log10(gOPT.ifar_max*YEAR);
627 
628  cout << endl;
629  cout << "ifar inf : " << pow(10, inf)/YEAR << " years" << endl;
630  cout << "ifar sup : " << pow(10, sup)/YEAR << " years" << endl;
631 
632  double difar_bin = (sup-inf)/IFAR_NBIN;
633 
634 
635  // =======================================================================
636  // SIMULATED EVENTS STUFF USED TO COMPUTE ASTROPHYSICAL RATES
637  // =======================================================================
638 
639  // -----------------------------------------------------
640  // read ifar simulated events
641  // -----------------------------------------------------
642 
643  double ifar_sim_min = 0;
644  double ifar_sim_max = 1.e20;
645 
646  while(!vIFAR.empty()) vIFAR.pop_back();
647  vIFAR.clear(); std::vector<double>().swap(vIFAR);
648  vector<int> vNINJ; // store injections per chunk in each event (used for normalization)
649 
650  cout << endl;
651  for(int m=0;m<nwave;m++) {
652  if(bin[m]>=0) continue; // select simulation root files
653 
654  int ninj = mdc[m]->GetEntries();
655  cout << run[m] << "\t" << chunk[m] << "\tmdc injected entries : " << ninj << endl;
656 
657  wave[m]->SetEstimate(wave[m]->GetEntries());
658  // ------------------------> APPLY TRIALS FACTOR !!!
659  sprintf(sel,"ifar/%f",gOPT.trials);
660  sprintf(tmp,"ifar>%f",pow(10,ifar_cmin)*gOPT.trials);
661  sprintf(tmp,"ifar>%f && rho[%d]>%f",pow(10,ifar_cmin)*gOPT.trials,rho_id,gOPT.rho_min);
662  Cuts = TString::Format("%s ",tmp);
663  wave[m]->Draw(sel,Cuts.Data(),"goff");
664  double* V1 = wave[m]->GetV1();
665  for(int k=0;k<wave[m]->GetSelectedRows();k++) {vIFAR.push_back(V1[k]);vNINJ.push_back(ninj);}
666  double efficiency = ninj>0 ? (double)wave[m]->GetSelectedRows()/(double)ninj : 0;
667  cout << run[m] << "\t" << chunk[m] << "\tselected wave sim detected entries : "
668  << wave[m]->GetSelectedRows() << "\tefficiency(%): " << int(100*efficiency) << "\tcumulative: " << vIFAR.size() << endl;
669  if(vIFAR.size()==0) continue;
670  int sel_entries = wave[m]->GetSelectedRows();
671  value = TMath::MinElement(sel_entries,V1);
672  //cout << "min sim IFAR per year : " << value/YEAR << endl;
673  if(value>ifar_sim_min) ifar_sim_min=value; // this condition is used to select the common sim IFAR inferior limit
674  value = TMath::MaxElement(sel_entries,V1);
675  //cout << "max sim IFAR per year : " << value/YEAR << endl;
676  if(value<ifar_sim_max) ifar_sim_max=value; // this condition is used to select the common sim IFAR superior limit
677  }
678 
679  if(vIFAR.size()>0) {
680  ifar_sim_min/=YEAR;
681  ifar_sim_max/=YEAR;
682 
683  cout << endl;
684  cout << "ifar_sim_min : " << ifar_sim_min << " years" << endl;
685  cout << "ifar_sim_max : " << ifar_sim_max << " years" << endl;
686  }
687 
688  // -----------------------------------------------------
689  // sort ifar simulated events
690  // -----------------------------------------------------
691 
692  if(vIFAR.size()>0) {
693  int size = vIFAR.size();
694  int *index = new int[size];
695  double *IFAR = new double[size];
696  for(int k=0;k<size;k++) IFAR[k]=vIFAR[k];
697  TMath::Sort(size,IFAR,index,true);
698  SIM_IFAR.resize(2*size-1);
699  SIM_NEVT.resize(2*size-1);
700  SIM_BKG_NEVT.resize(2*size-1);
701  int n=0;
702 // vector<double> vWEIGHT(size); // store weigths used for normalization
703 // vWEIGHT[0]=1./vNINJ[0];
704 // for(int k=1;k<size;k++) vWEIGHT[k]=vWEIGHT[k-1]+1./vNINJ[k-1]; // init weigths
705  for(int k=size-1;k>=0;k--) { // invert the ifar order
706  SIM_IFAR[n] = IFAR[index[k]]/YEAR; // sorted ifar simulated events
707 #ifdef APPLY_MAX_SIM_IFAR
708  if(SIM_IFAR[n]<ifar_sim_min) SIM_IFAR[n]=ifar_sim_min;
709  if(SIM_IFAR[n]>ifar_sim_max) SIM_IFAR[n]=ifar_sim_max;
710 #endif
711 // SIM_NEVT[n] = vWEIGHT[k]; // number of simulated events with ifar > SIM_IFAR[k]
712  SIM_NEVT[n] = k+1; // number of simulated events with ifar > SIM_IFAR[k]
713  n++;
714  if(n>=SIM_NEVT.size()-1) break;
715  // add point to make ZL STYLE STEPWISE
716  SIM_IFAR[n] = SIM_IFAR[n-1];
717  SIM_NEVT[n] = SIM_NEVT[n-1]-1;
718 // SIM_NEVT[n] = vWEIGHT[k-1];
719  n++;
720  }
721  delete [] index;
722  delete [] IFAR;
723  }
724  SIM_NEVT*=1./gOPT.sfactor; // expected events from astrophysical reates
725 // SIM_NEVT*=gOPT.sfactor; // expected events from astrophysical reates
726  for(int k=0;k<SIM_NEVT.size();k++) SIM_BKG_NEVT[k]=SIM_NEVT[k]+obsTime/(SIM_IFAR[k]*YEAR); // add expected from background
727 
728 #ifdef APPLY_MAX_SIM_IFAR
729  // set cut off of zero lag events according to the common sim IFAR inferior limit
730  if(SIM_NEVT.size()>0) {
731  for(int i=0;i<ZL_NEVT.size();i++) if(ZL_IFAR[i]>ifar_sim_max) ZL_IFAR[i]=ifar_sim_max;
732  }
733 #endif
734 
735  // -----------------------------------------------------
736  // make sim plot
737  // -----------------------------------------------------
738 
739  grsim=NULL;
740  grsimbkg=NULL;
741  if(SIM_NEVT.size()>0) {
742 
743 // (RATE_TAG)
744  // compute the index which contains the maximum common IFAR
745  int ifar_cmax_index=0;
746  for(int i=0;i<SIM_IFAR.size();i++) if(SIM_IFAR.data[i]<pow(10,ifar_cmax)/YEAR) ifar_cmax_index=i;
747 
748 // (RATE_TAG)
749  grsim = new TGraph(ifar_cmax_index,SIM_IFAR.data,SIM_NEVT.data); // expected events from astrophysical reates
750  grsimbkg = new TGraph(ifar_cmax_index,SIM_IFAR.data,SIM_BKG_NEVT.data); // expected events from astrophysical reates + expected from background
751  }
752 
753  // -----------------------------------------------------
754  // make sim belts at FAP0
755  // -----------------------------------------------------
756 
757  if(SIM_NEVT.size()>0) {
758  int M=0;
759  for(int n=0;n<IFAR_NBIN;n++) {
760  double x = inf+n*difar_bin;
761  if(x<(ifar_cmin-difar_bin) && x<inf) continue; // RATE_TAG
762  if(x>ifar_cmax) continue; // RATE_TAG
763  xSIM_BKG_IFAR[M] = pow(10,x)/YEAR;
764  if((xSIM_BKG_IFAR[M]<SIM_IFAR[0]) || (xSIM_BKG_IFAR[M]>SIM_IFAR[SIM_IFAR.size()-1])) continue;
765  double mu = grsimbkg->Eval(xSIM_BKG_IFAR[M]); // expected events in ZERO_LAG_TIME
766  ySIM_BKG[M]=mu;
767 
768  double q[6];
769  if(gOPT.belt_style=="continuos") GetPercentiles(mu, q);
770 
771  double nevt=0;
772  nevt = gOPT.belt_style=="continuos" ? q[3] : GetPoissonNsup(mu,FAP0/2.);
773  eySIM_BKG_sup0[M] = (nevt-mu);
774  nevt = gOPT.belt_style=="continuos" ? q[2] : GetPoissonNinf(mu,FAP0/2.);
775  eySIM_BKG_inf0[M] = (mu-nevt);
776  M++;
777  }
778  xSIM_BKG_IFAR.resize(M);
779  ySIM_BKG.resize(M);
780  eySIM_BKG_inf0.resize(M);
781  eySIM_BKG_sup0.resize(M);
782  }
783 
784 
785  // =======================================================================
786  // MAKE BELTS FOR EXPECTED BACKGROUND
787  // =======================================================================
788 
789  // -----------------------------------------------------
790  // make belts at FAP0,FAP1,FAP2
791  // -----------------------------------------------------
792 
793  if(gOPT.belt_style=="step") {
794  // add points to make BELT STYLE STEPWISE
795 
796  xIFAR.resize(2*IFAR_NBIN);
797  xFAR.resize(2*IFAR_NBIN);
798  bRATE.resize(2*IFAR_NBIN);
799  yBKG.resize(2*IFAR_NBIN);
800 
801  exFAR.resize(2*IFAR_NBIN); exFAR=0;
802  eyRATEinf0.resize(2*IFAR_NBIN);
803  eyRATEsup0.resize(2*IFAR_NBIN);
804  eyRATEinf1.resize(2*IFAR_NBIN);
805  eyRATEsup1.resize(2*IFAR_NBIN);
806  eyRATEinf2.resize(2*IFAR_NBIN);
807  eyRATEsup2.resize(2*IFAR_NBIN);
808  }
809 
810  int N=0;
811  for(int n=0;n<IFAR_NBIN;n++) {
812  double x = inf+n*difar_bin;
813  //if(x<(ifar_cmin-difar_bin) || x>(ifar_max+2*difar_bin)) continue; // RATE_TAG
814  if(x<(ifar_cmin-difar_bin) && x<inf) continue; // RATE_TAG
815  xIFAR[N] = pow(10,x);
816  xFAR[N] = 1./xIFAR[N];
817  bRATE[N] = xFAR[N];
818  double mu = xFAR[N]*obsTime; // expected bakground events in ZERO_LAG_TIME
819 
820  double q[6];
821  if(gOPT.belt_style=="continuos") GetPercentiles(mu, q);
822 
823  yBKG[N] = gOPT.belt_style=="continuos" ? q[3] : GetPoissonNsup(mu,FAP0/2.);
824  eyRATEsup0[N] = (yBKG[N]-mu)/obsTime;
825  yBKG[N] = gOPT.belt_style=="continuos" ? q[2] : GetPoissonNinf(mu,FAP0/2.);
826  eyRATEinf0[N] = (mu-yBKG[N])/obsTime;
827  yBKG[N] = gOPT.belt_style=="continuos" ? q[4] : GetPoissonNsup(mu,FAP1/2.);
828  eyRATEsup1[N] = (yBKG[N]-mu)/obsTime;
829  yBKG[N] = gOPT.belt_style=="continuos" ? q[1] : GetPoissonNinf(mu,FAP1/2.);
830  eyRATEinf1[N] = (mu-yBKG[N])/obsTime;
831  yBKG[N] = gOPT.belt_style=="continuos" ? q[5] : GetPoissonNsup(mu,FAP2/2.);
832  eyRATEsup2[N] = (yBKG[N]-mu)/obsTime;
833  yBKG[N] = gOPT.belt_style=="continuos" ? q[0] : GetPoissonNinf(mu,FAP2/2.);
834  eyRATEinf2[N] = (mu-yBKG[N])/obsTime;
835  N++;
836 
837  if(gOPT.belt_style=="step") { // add points to make BELT STYLE STEPWISE
838 
839  if(N>1) {
840 
841  xIFAR[N] = xIFAR[N-1];
842  xFAR[N] = xFAR[N-1];
843  bRATE[N] = bRATE[N-1];
844 
845  xIFAR[N-1] = xIFAR[N-2];
846  xFAR[N-1] = xFAR[N-2];
847  bRATE[N-1] = bRATE[N-2];
848 
849  yBKG[N-1] = mu;
850  yBKG[N] = mu;
851 
852  eyRATEsup0[N] = eyRATEsup0[N-1];
853  eyRATEinf0[N] = eyRATEinf0[N-1];
854  eyRATEsup1[N] = eyRATEsup1[N-1];
855  eyRATEinf1[N] = eyRATEinf1[N-1];
856  eyRATEsup2[N] = eyRATEsup2[N-1];
857  eyRATEinf2[N] = eyRATEinf2[N-1];
858 
859  N++;
860 
861  } else {
862  yBKG[N-1]=mu;
863  }
864  } else {
865  yBKG[N-1]=bRATE[N-1]*obsTime;
866  }
867  }
868  xIFAR.resize(N);
869  xFAR.resize(N);
870 
871  // normalize -> per SEC -> per YEAR
872  xFAR *=YEAR;
873  xIFAR*=1./YEAR;
874 
875  // compute number of events in the obsTime period
876  bRATE*=obsTime;
877  eyRATEsup0*=obsTime;
878  eyRATEinf0*=obsTime;
879  eyRATEsup1*=obsTime;
880  eyRATEinf1*=obsTime;
881  eyRATEsup2*=obsTime;
882  eyRATEinf2*=obsTime;
883 
884 
885  // -----------------------------------------------------
886  // print sigma belts
887  // -----------------------------------------------------
888 
889  double sigma[3];
890  sigma[0] = sqrt(2)*TMath::ErfcInverse(FAP0);
891  sigma[1] = sqrt(2)*TMath::ErfcInverse(FAP1);
892  sigma[2] = sqrt(2)*TMath::ErfcInverse(FAP2);
893 
894  cout << endl << "--------------------------------------------------------" << endl << endl;;
895  cout << "FAP = " << FAP0 << "\t -> SIGMA = " << sigma[0] << endl;
896  cout << "FAP = " << FAP1 << "\t -> SIGMA = " << sigma[1] << endl;
897  cout << "FAP = " << FAP2 << "\t -> SIGMA = " << sigma[2] << endl;
898  cout << endl << "--------------------------------------------------------" << endl << endl;;
899 
900 
901  // =======================================================================
902  // DUMP PLOTS AND ASCII FILES
903  // =======================================================================
904 
905  // -----------------------------------------------------
906  // Make Plot
907  // -----------------------------------------------------
908 
909  char title[1024];
910  if(gOPT.bbh) {
911  sprintf(title,"cWB %s (run=%s : lag=%d : slag=%d : included BBH)",gOPT.search.Data(),gOPT.run.Data(),gOPT.lag,gOPT.slag);
912  } else {
913  sprintf(title,"cWB %s (run=%s : lag=%d : slag=%d : excluded BBH)",gOPT.search.Data(),gOPT.run.Data(),gOPT.lag,gOPT.slag);
914  }
915 
916  sprintf(title,"");
917 
918  MakePlot(title, obsTime, ifar_cmax);
919 
920  // -----------------------------------------------------
921  // Dump Zero Lag Events to ascii file
922  // -----------------------------------------------------
923 
924  DumpPeriod(liveTot);
925  DumpLoudest(obsTime);
926  DumpBackground();
927  DumpForeground();
928  DumpSimulation(ifar_cmax);
929 
930  if(gOPT.exit && gOPT.pfname!="") exit(0);
931 }
932 
933 void MakePlot(TString title, double obsTime, double ifar_cmax) {
934 
935  canvas = new TCanvas("FAR", "FAR", 300,40, 800, 600);
936  canvas->SetGridx();
937  canvas->SetGridy();
938  canvas->SetLogx();
939  canvas->SetLogy();
940 // canvas->SetFrameBorderMode(0);
941 // canvas->SetFrameLineColor(17);
942 // canvas->SetFrameBorderMode(0);
943 
944  gr0 = new TGraph(xIFAR.size(),xIFAR.data,bRATE.data);
945  gr0->SetTitle(title);
946  gr0->Draw("LA3");
947  gr0->SetLineWidth(2);
948  gr0->SetLineStyle(1);
949  gr0->SetLineColor(kGray+2);
950  gr0->GetHistogram()->GetXaxis()->CenterTitle();
951  gr0->GetHistogram()->GetYaxis()->CenterTitle();
952  gr0->GetHistogram()->GetXaxis()->SetLabelOffset(0.0001);
953  gr0->GetHistogram()->GetXaxis()->SetTickLength(0.02);
954  gr0->GetHistogram()->GetXaxis()->SetLabelSize(0.035);
955  gr0->GetHistogram()->GetYaxis()->SetLabelSize(0.035);
956  gr0->GetHistogram()->GetXaxis()->SetTitleSize(0.035);
957  gr0->GetHistogram()->GetYaxis()->SetTitleSize(0.035);
958  gr0->GetHistogram()->GetXaxis()->SetTitleOffset(1.7);
959  gr0->GetHistogram()->GetYaxis()->SetTitleOffset(1.5);
960  gr0->GetHistogram()->GetXaxis()->SetTitleOffset(1.4);
961  gr0->GetHistogram()->GetYaxis()->SetTitleOffset(1.2);
962  gr0->GetHistogram()->GetXaxis()->SetAxisColor(17);
963  gr0->GetHistogram()->GetYaxis()->SetAxisColor(17);
964  canvas->RedrawAxis();
965  canvas->Update();
966 
967  // set axis ranges
968  if(xIFAR.size()>0) {
969  double ymin = 0.7;
970  double ymax = gr0->GetHistogram()->GetMaximum();
971  if(ZL_NEVT.size()>0) ymax=2*ZL_NEVT[0];
972  if(ymax/ymin<10) ymax=10*ymin;
973  if(ymax>gr0->GetHistogram()->GetMaximum()) ymax=gr0->GetHistogram()->GetMaximum();
974  gr0->GetHistogram()->GetYaxis()->SetRangeUser(ymin,ymax);
975  gr0->GetHistogram()->GetXaxis()->SetRangeUser(0.01,1e5);
976  if((xIFAR[0]>0) && (xIFAR[xIFAR.size()-1]/xIFAR[0])<10) {
977  canvas->SetLogx(false);
978  canvas->SetLogy(false);
979  }
980  } else {
981  canvas->SetLogx(false);
982  canvas->SetLogy(false);
983  }
984  gr0->GetHistogram()->GetXaxis()->SetTitle("Inverse False Alarm Rate (yr)");
985  gr0->GetHistogram()->GetYaxis()->SetTitle("Cumulative Number of Events");
986  gr0->GetHistogram()->GetXaxis()->SetLabelFont(132);
987  gr0->GetHistogram()->GetYaxis()->SetLabelFont(132);
988  gr0->GetHistogram()->GetXaxis()->SetTitleFont(132);
989  gr0->GetHistogram()->GetYaxis()->SetTitleFont(132);
990  gr0->GetHistogram()->GetXaxis()->SetTickLength(0);
991 
992  // Plot FAP2
993  egr2 = new TGraphAsymmErrors(xIFAR.size(),xIFAR.data,yBKG.data,exFAR.data,exFAR.data,eyRATEinf2.data,eyRATEsup2.data);
994  egr2->SetMarkerStyle(20);
995  egr2->SetMarkerSize(0);
996  egr2->SetLineColor(15);
997  egr2->SetFillStyle(1001);
998  egr2->SetFillColor(18); // light gray
999  egr2->Draw("3SAME");
1000 
1001  // Plot FAP1
1002  egr1 = new TGraphAsymmErrors(xIFAR.size(),xIFAR.data,yBKG.data,exFAR.data,exFAR.data,eyRATEinf1.data,eyRATEsup1.data);
1003  egr1->SetMarkerStyle(20);
1004  egr1->SetMarkerSize(0);
1005  egr1->SetLineColor(15);
1006  egr1->SetFillStyle(1001);
1007  egr1->SetFillColor(17); // gray
1008  egr1->Draw("3SAME");
1009 
1010  // Plot FAP0
1011  egr0 = new TGraphAsymmErrors(xIFAR.size(),xIFAR.data,yBKG.data,exFAR.data,exFAR.data,eyRATEinf0.data,eyRATEsup0.data);
1012  egr0->SetMarkerStyle(20);
1013  egr0->SetMarkerSize(0);
1014  egr0->SetLineColor(15);
1015  egr0->SetFillStyle(1001);
1016  egr0->SetFillColor(16); // dark gray
1017  egr0->Draw("3SAME");
1018 
1019 // BEG TEST P-ASTRO
1020  // Plot Simulated
1021  if(grsimbkg!=NULL) {
1022 
1023  // compute the index which contains the maximum common IFAR (RATE_TAG)
1024  int ifar_cmax_index=0;
1025  for(int i=0;i<xSIM_BKG_IFAR.size();i++) if(xSIM_BKG_IFAR.data[i]<pow(10,ifar_cmax)/YEAR) ifar_cmax_index=i;
1026 
1027  // Plot sim FAP0
1028  egrsimbkg0 = new TGraphAsymmErrors(ifar_cmax_index,xSIM_BKG_IFAR.data,ySIM_BKG.data,exSIM_BKG_IFAR.data,exSIM_BKG_IFAR.data,eySIM_BKG_inf0.data,eySIM_BKG_sup0.data); // (RATE_TAG)
1029  egrsimbkg0->SetMarkerSize(0);
1030  egrsimbkg0->SetLineColor(15);
1031  egrsimbkg0->SetFillColorAlpha(kGreen, 0.1);
1032  egrsimbkg0->Draw("3SAME"); // belt
1033 
1034  grsim->SetLineColor(kGreen);
1035  grsim->SetLineWidth(1);
1036  grsim->SetMarkerSize(0);
1037  grsim->Draw("LSAME"); // expected signal from astrophysical rates
1038 
1039  grsimbkg->SetLineColor(kGreen+2);
1040  grsimbkg->SetLineWidth(1);
1041  grsimbkg->SetMarkerSize(0);
1042  grsimbkg->Draw("LSAME"); // expected signal from astrophysical rates + expected from backgrouns
1043  }
1044 // END TEST P-ASTRO
1045 
1046  // Plot Foreground
1047  if(ZL_NEVT.size()>0) {
1048  gr = new TGraph(ZL_NEVT.size(),ZL_IFAR.data,ZL_NEVT.data);
1049 //gr->SetLineStyle(2);
1050 //gr->SetLineWidth(2);
1051  gr->SetLineColor(kRed);
1052  gr->SetMarkerColor(kRed);
1053  gr->SetMarkerStyle(23);
1054  if(gOPT.zl_style=="marker") {
1055  gr->SetMarkerSize(0.4);
1056  gr->Draw("PSAME"); // marker
1057  } else {
1058  gr->SetMarkerSize(0);
1059  gr->Draw("SAME"); // stepwise
1060  }
1061  } else {
1062  gr = NULL;
1063  //canvas->SetLogy(false);
1064  }
1065 
1066 #ifdef PLOT_FOREGROUND_NOBBH
1067 
1068  wavearray<double> ZL_IFAR_NOBBH(ZL_NEVT.size()); // sorted ifar events
1069 
1070  // select event without known BBH
1071  int nobbh_size=0;
1072  for(int i=0;i<ZL_NEVT.size();i++) {
1073  if((GetNameBBH(ZL_GPS[i])=="")) ZL_IFAR_NOBBH[nobbh_size++]=ZL_IFAR[i];
1074  }
1075  if(ZL_IFAR_NOBBH[nobbh_size-1]==ZL_IFAR_NOBBH[nobbh_size-2]) nobbh_size-=1;
1076  ZL_IFAR_NOBBH.resize(nobbh_size);
1077  wavearray<double> ZL_NEVT_NOBBH(nobbh_size); // number of nobbh events with ifar > ZL_IFAR[k]
1078  for(int i=0;i<nobbh_size;i+=2) {
1079  ZL_NEVT_NOBBH[i]=(nobbh_size-i)/2+1;
1080  ZL_NEVT_NOBBH[i+1]=(nobbh_size-i)/2;
1081  }
1082  ZL_NEVT_NOBBH[0]=(nobbh_size)/2+1;
1083 
1084  // Plot Foreground no bbh
1085  TGraph* gr_nobbh;
1086  if(ZL_NEVT_NOBBH.size()>0) {
1087  gr_nobbh = new TGraph(ZL_NEVT_NOBBH.size(),ZL_IFAR_NOBBH.data,ZL_NEVT_NOBBH.data);
1088  //gr_nobbh->SetLineColor(kRed);
1089  //gr_nobbh->SetLineColor(kGreen+2);
1090  //gr_nobbh->SetLineColor(kOrange+7);
1091  gr_nobbh->SetLineColor(kBlack);
1092  gr_nobbh->SetLineWidth(1);
1093  gr_nobbh->SetLineStyle(3);
1094  gr_nobbh->SetMarkerColor(kRed);
1095  gr_nobbh->SetMarkerStyle(23);
1096  if(gOPT.zl_style=="marker") {
1097  gr_nobbh->SetMarkerSize(0.4);
1098  gr_nobbh->Draw("PSAME"); // marker
1099  } else {
1100  gr_nobbh->SetMarkerSize(0);
1101  gr_nobbh->Draw("SAME"); // stepwise
1102  }
1103  }
1104 
1105 #endif
1106 
1107  gr0->Draw("LSAME");
1108 
1109  // Plot blue marker & name for known BBH
1110  bool anyBBH=false;
1111  for(int i=0;i<ZL_NEVT.size();i++) {
1112  if((GetNameBBH(ZL_GPS[i])!="")&&(i%2==0)) {
1113  marker = new TMarker(ZL_IFAR[i],ZL_NEVT[i], 20);
1114  bool selected=false;
1115  if((GetNameBBH(ZL_GPS[i])=="GW150914")) selected=true;
1116  if((GetNameBBH(ZL_GPS[i])=="GW170814")) selected=true;
1117  if((GetNameBBH(ZL_GPS[i])=="GW170104")) selected=true;
1118  if((GetNameBBH(ZL_GPS[i])=="GW170608")) selected=true;
1119  if(selected) {
1120  arrow = new TArrow(ZL_IFAR[i],ZL_NEVT[i],ZL_IFAR[i]*1.5,ZL_NEVT[i],0.05,"|>");
1121  arrow->SetAngle(60);
1122  arrow->SetLineWidth(2);
1123  arrow->SetLineColor(kBlue);
1124  arrow->SetFillColor(kBlue);
1125  arrow->SetArrowSize(0.012);
1126  arrow->Draw();
1127  }
1128  marker->SetMarkerSize(1.2);
1129  marker->SetMarkerColor(kBlue);
1130  marker->Draw();
1131 
1132  latex = new TLatex(ZL_IFAR[i]*(4./3.),ZL_NEVT[i]*1.1, GetNameBBH(ZL_GPS[i]).Data());
1133  latex->SetTextFont(132);
1134  latex->SetTextSize(0.028);
1135  latex->SetTextColor(kBlack);
1136  latex->Draw();
1137 
1138  anyBBH=true;
1139  }
1140  }
1141  // if there no BBH events then a red marker is draw for loudest event
1142  if((!anyBBH)&&(ZL_NEVT.size()>0)&&(GetNameBBH(ZL_GPS[ZL_NEVT.size()-1])=="")) {
1143  marker = new TMarker(ZL_IFAR[ZL_NEVT.size()-1],ZL_NEVT[ZL_NEVT.size()-1], 23);
1144  marker->SetMarkerSize(1.2);
1145  marker->SetMarkerColor(kRed);
1146  marker->Draw();
1147  }
1148 
1149  // draw legend
1150 #ifdef PLOT_FOREGROUND_NOBBH
1151  leg = new TLegend(0.6,0.6,0.90,0.90,NULL,"brNDC");
1152 #else
1153  if(grsimbkg!=NULL) {
1154  leg = new TLegend(0.6,0.630,0.90,0.90,NULL,"brNDC");
1155  } else {
1156  leg = new TLegend(0.6,0.6,0.90,0.90,NULL,"brNDC");
1157  }
1158 #endif
1159 
1160  leg->SetBorderSize(1);
1161  leg->SetTextAlign(22);
1162  leg->SetTextFont(132);
1163  leg->SetLineColor(1);
1164  leg->SetLineStyle(1);
1165  leg->SetLineWidth(1);
1166  leg->SetFillColor(0);
1167  leg->SetFillStyle(1001);
1168  if(grsimbkg!=NULL) leg->SetTextSize(0.035); else leg->SetTextSize(0.035);
1169  leg->SetLineColor(kBlack);
1170  leg->SetFillColor(kWhite);
1171 
1172  char legLabel[64];
1173  sprintf(legLabel,"Foreground"); if(gr!=NULL) leg->AddEntry(gr,legLabel,"lp");
1174 #ifdef PLOT_FOREGROUND_NOBBH
1175  sprintf(legLabel,"Residual Foreground"); if(gr_nobbh!=NULL) leg->AddEntry(gr_nobbh,legLabel,"lp");
1176 #endif
1177  sprintf(legLabel,"Expected Background"); leg->AddEntry(gr0,legLabel,"lp");
1178  if(gOPT.bbh && anyBBH) leg->AddEntry(marker,"BBH Events","p");
1179  else if(ZL_NEVT.size()>0) leg->AddEntry(marker,"Loudest","p");
1180  if(grsimbkg!=NULL) {
1181  leg->AddEntry(grsim,"Signal Model","lp");;
1182  leg->AddEntry(grsimbkg,"Noise+Signal Model","lp");;
1183  } else {
1184  // add background belts labels
1185  leg->AddEntry(egr2,"< 3 #sigma","f");
1186  leg->AddEntry(egr1,"< 2 #sigma","f");
1187  leg->AddEntry(egr0,"< 1 #sigma","f");
1188  }
1189  leg->Draw();
1190 
1191  // fix frame lines color
1192  double xmax = gr0->GetHistogram()->GetXaxis()->GetXmax();
1193  double xmin = gr0->GetHistogram()->GetXaxis()->GetXmin();
1194  //double ymin = pow(10,gPad->GetUymin());
1195  double ymin = 0.0;
1196  double ymax = pow(10,gPad->GetUymax());
1197  cout << "xmin=" << xmin << " xmax=" << xmax << endl;
1198  cout << "ymin=" << ymin << " ymax=" << ymax << endl;
1199 
1200  TLine* frameLine[4];
1201  frameLine[0] = new TLine(xmin,ymin,xmax,ymin);
1202  frameLine[1] = new TLine(xmin,ymin,xmin,ymax);
1203  frameLine[2] = new TLine(xmax,ymin,xmax,ymax);
1204  frameLine[3] = new TLine(0.,ymax,xmax,ymax);
1205  for(int i=0;i<4;i++) frameLine[i]->Draw();
1206 
1207 
1208  // dump plot to disk
1209  if(gOPT.pfname) {
1210  char ofname[1024];
1211  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1212  if(gOPT.bbh) sprintf(ofname,"%s_bbh_plot.png",PFNAME.Data());
1213  else sprintf(ofname,"%s_nobbh_plot.png",PFNAME.Data());
1214  if(gOPT.bbh) sprintf(ofname,"%s_bbh_plot.root",PFNAME.Data());
1215  else sprintf(ofname,"%s_nobbh_plot.root",PFNAME.Data());
1216  if(gOPT.bbh) sprintf(ofname,"%s_bbh_plot.pdf",PFNAME.Data());
1217  else sprintf(ofname,"%s_nobbh_plot.pdf",PFNAME.Data());
1218  canvas->Print(ofname);
1219  }
1220 }
1221 
1222 int GetPercentiles(double mu, double* q) {
1223 
1224  double max = mu>10 ? 10*mu : 100;
1225 
1226  TF1 f("poisson","TMath::Poisson(x,[1])",0,max);
1227  f.SetParameter(1,mu);
1228  f.SetNpx(10000);
1229 
1230  double probSum[6] = {FAP2/2,FAP1/2,FAP0/2,1-FAP0/2,1-FAP1/2,1-FAP2/2};
1231  return f.GetQuantiles(6,q,probSum);
1232 }
1233 
1234 int GetPoissonNsup(double mu, double fap) {
1235 
1236  double FAP = 1.;
1237  for(int j=0;j<1000000;j++) {
1238  //FAP -= pow(mu,j)*exp(-mu)/TMath::Factorial(j);
1239  FAP -= TMath::PoissonI(j,mu);
1240  //cout << " N : " << j << " MU : " << mu << " FAP " << FAP << " fap : " << fap << endl;
1241  if(FAP<fap) return j==0 ? 0 : j+1;
1242  }
1243  return 0;
1244 }
1245 
1246 int GetPoissonNinf(double mu, double fap) {
1247 
1248  double FAP = 0.;
1249  for(int j=0;j<1000000;j++) {
1250  //FAP -= pow(mu,j)*exp(-mu)/TMath::Factorial(j);
1251  FAP += TMath::PoissonI(j,mu);
1252  //cout << " N : " << j << " MU : " << mu << " FAP " << FAP << " fap : " << fap << endl;
1253  if(FAP>fap) return j==0 ? 0 : j-1;
1254  }
1255  return 0;
1256 }
1257 
1258 int ReadFileList(TString ifName, TChain** live, TChain** wave, TChain** mdc, int* lag, int* slag, int* chunk, int* bin, TString* run) {
1259 
1260  // Read WAVE/LIVE file list
1261 
1262  ifstream in;
1263  in.open(ifName.Data(),ios::in);
1264  if (!in.good()) {cout << "ReadFileList Error Opening File : " << ifName.Data() << endl;exit(1);}
1265 
1266  int size=0;
1267  char str[1024];
1268  int fpos=0;
1269  while(true) {
1270  in.getline(str,1024);
1271  if (!in.good()) break;
1272  if(str[0] != '#') size++;
1273  }
1274  in.clear(ios::goodbit);
1275  in.seekg(0, ios::beg);
1276  if (size==0) {cout << "ReadFileList Error : File " << ifName.Data() << " is empty" << endl;exit(1);}
1277  if (size>MAX_LIST_SIZE) {cout << "ReadFileList Error : Files in " << ifName.Data() << " > " << MAX_LIST_SIZE << endl;exit(1);}
1278 
1279  char sfile[1024];
1280  int xlag,xslag,ichunk,ibin;
1281  char srun[256];
1282  int nwave=0;
1283  while(true) {
1284  fpos=in.tellg();
1285  in.getline(str,1024);
1286  if (!in.good()) break;
1287  if(str[0] == '#' || str[0]=='\0') continue;
1288  in.seekg(fpos, ios::beg);
1289  in >> sfile >> xlag >> xslag >> ichunk >> ibin >> srun;
1290  if(!in.good()) break;
1291  lag[nwave]=xlag;
1292  slag[nwave]=xslag;
1293  TString fwave = sfile;
1294  if(!CWB::Toolbox::isFileExisting(fwave)) {cout << "ReadFileList Error : File not exist : " << fwave << endl;exit(1);}
1295  wave[nwave]->Add(fwave);
1296  wave[nwave]->SetTitle(fwave);
1297  if(ibin>=0) { // ibin<=0/ibin>=0 -> simulation/background
1298  TString flive = sfile;
1299  flive.ReplaceAll("wave_","live_");
1300  if(!CWB::Toolbox::isFileExisting(flive)) {cout << "ReadFileList Error : File not exist : " << flive << endl;exit(1);}
1301  live[nwave]->Add(flive);
1302  live[nwave]->SetTitle(flive);
1303  } else { // simulation -> read mdc
1304  TString fmdc = sfile;
1305  fmdc.ReplaceAll("wave_","mdc_");
1306  if(!CWB::Toolbox::isFileExisting(fmdc)) {cout << "ReadFileList Error : File not exist : " << fmdc << endl;exit(1);}
1307  mdc[nwave]->Add(fmdc);
1308  mdc[nwave]->SetTitle(fmdc);
1309  }
1310  chunk[nwave]=ichunk;
1311  bin[nwave]=ibin;
1312  run[nwave]=srun;
1313  cout << nwave+1 << "\t" << fwave << "\t" << xlag << "\t" << xslag << "\t" << ichunk << "\t" << ibin << "\t" << srun << endl;
1314  nwave++;
1315  }
1316  in.close();
1317 
1318  return nwave;
1319 }
1320 
1321 void GetLiveTime(int nwave, TChain** live, int* lag, int* slag, int* bin, TString* run, double* liveZero, double* liveTot) {
1322 
1323  CWB::Toolbox TB;
1324 
1325  wavearray<double> Trun(500000); Trun = 0.;
1326  wavearray<double> Wlag[NIFO_MAX+1];
1327  wavearray<double> Wslag[NIFO_MAX+1];
1328  wavearray<double> Tdlag;
1329  wavearray<double> Tlag;
1330  wavearray<double> Rlag;
1331  wavearray<double> Elag;
1332  wavearray<double> Olag;
1333 
1334  cout << "Start CWB::Toolbox::getLiveTime : be patient, it takes a while ..." << endl;
1335  for(int m=0;m<nwave;m++) {
1336  if(bin[m]<0) continue; // simulation
1337  cout << "Compute liveTot : Read : " << live[m]->GetTitle() << endl;
1338  if(run[m]=="O1") liveTot[0]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,lag[m],slag[m]);
1339  if(run[m]=="O2") liveTot[1]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,lag[m],slag[m]);
1340  }
1341 
1342  // get live time zero (when is not included into Tlag array)
1343  for(int m=0;m<nwave;m++) {
1344  if(bin[m]<0) continue; // simulation
1345  if((lag[m]!=0)&&(slag[m]!=0)) {
1346  cout << "Compute liveZero : Read : " << live[m]->GetTitle() << endl;
1347  if(run[m]=="O1") liveZero[0]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,0,0);
1348  if(run[m]=="O2") liveZero[1]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,0,0);
1349  //liveZero+=TB.getZeroLiveTime(gOPT.nifo,*live[m]); // crash in ATLAS ???
1350  if(run[m]=="O1") printf("live time: zero lags = %10.1f, %10.1f days \n",liveZero[0],liveZero[0]/DAY);
1351  if(run[m]=="O2") printf("live time: zero lags = %10.1f, %10.1f days \n",liveZero[1],liveZero[1]/DAY);
1352  }
1353  }
1354 }
1355 
1356 void InitBBH() {
1357 
1358 // ---------------------------------------------------------------------------------
1359 // INIT KNOWN O1,O2 BBH
1360 // ---------------------------------------------------------------------------------
1361 
1362  // O1 known BBH events
1363  O1_BBH_NAME.resize(3);
1364  O1_BBH_GPS.resize(3);
1365 
1366  O1_BBH_NAME[0]="GW150914"; O1_BBH_GPS[0]=1126259462.40;
1367  O1_BBH_NAME[1]="LVT151012"; O1_BBH_GPS[1]=1128678900.55;
1368  O1_BBH_NAME[2]="GW151226"; O1_BBH_GPS[2]=1135136350.56;
1369 
1370  // O2 known BBH events
1371  O2_BBH_NAME.resize(7);
1372  O2_BBH_GPS.resize(7);
1373 
1374  O2_BBH_NAME[0]="GW170104"; O2_BBH_GPS[0]=1167559936.59;
1375  O2_BBH_NAME[1]="GW170608"; O2_BBH_GPS[1]=1180922494.38;
1376  O2_BBH_NAME[2]="GW170729"; O2_BBH_GPS[2]=1185389807.33;
1377  O2_BBH_NAME[3]="GW170809"; O2_BBH_GPS[3]=1186302519.73;
1378  O2_BBH_NAME[4]="GW170814"; O2_BBH_GPS[4]=1186741861.50;
1379  O2_BBH_NAME[5]="GW170818"; O2_BBH_GPS[5]=1187058327.07;
1380  O2_BBH_NAME[6]="GW170823"; O2_BBH_GPS[6]=1187529256.51;
1381 
1382 }
1383 
1384 TString GetCutBBH() {
1385 
1386  TString CutBBH="";
1387 
1388  if(gOPT.run.Contains("O1")) {
1389  CutBBH += TString::Format(" !(abs(time[0]-%f)<0.2) ",O1_BBH_GPS[0]);
1390  }
1391  if(gOPT.run.Contains("O1")) {
1392  for(int i=1;i<O1_BBH_GPS.size();i++) CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.2) ",O1_BBH_GPS[i]);
1393  }
1394 
1395  if(gOPT.run=="O1O2") {
1396  CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.2) ",O2_BBH_GPS[0]);
1397  } else if(gOPT.run=="O2") {
1398  CutBBH += TString::Format(" !(abs(time[0]-%f)<0.2) ",O2_BBH_GPS[0]);
1399  }
1400  if(gOPT.run.Contains("O2")) {
1401  for(int i=1;i<O2_BBH_GPS.size();i++) CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.2) ",O2_BBH_GPS[i]);
1402  }
1403 
1404  return CutBBH;
1405 }
1406 
1407 TString GetNameBBH(double gps) {
1408 
1409  // known O1 BBH
1410  for(int i=0;i<O1_BBH_GPS.size();i++) if(fabs(gps-O1_BBH_GPS[i])<0.2) return O1_BBH_NAME[i];
1411 
1412  // known O2 BBH
1413  for(int i=0;i<O2_BBH_GPS.size();i++) if(fabs(gps-O2_BBH_GPS[i])<0.2) return O2_BBH_NAME[i];
1414 
1415  return "";
1416 }
1417 
1418 double GetGPSBBH(TString name) {
1419 
1420  // known O1 BBH
1421  for(int i=0;i<O1_BBH_NAME.size();i++) if(O1_BBH_NAME[i]==name) return O1_BBH_GPS[i];
1422 
1423  // known O2 BBH
1424  for(int i=0;i<O2_BBH_NAME.size();i++) if(O2_BBH_NAME[i]==name) return O2_BBH_GPS[i];
1425 
1426  return 0;
1427 }
1428 
1429 void ResetUserOptions() {
1430 
1431  gOPT.run = PP_RUN;
1432  gOPT.search = PP_SEARCH;
1433 
1434  gOPT.nifo = PP_NIFO;
1435  gOPT.lag = PP_LAG;
1436  gOPT.slag = PP_SLAG;
1437  gOPT.chunk = PP_CHUNK;
1438 
1439  gOPT.trials = PP_TRIALS;
1440 
1441  gOPT.pfname = PP_PFNAME;
1442  gOPT.zl_style = PP_ZL_STYLE;
1443  gOPT.belt_style = PP_BELT_STYLE;
1444 
1445  gOPT.ifar_max = PP_IFAR_MAX;
1446  gOPT.rho_min = PP_RHO_MIN;
1447 
1448  gOPT.sfactor = PP_SFACTOR;
1449 
1450  gOPT.bbh = PP_BBH;
1451  gOPT.exit = PP_EXIT;
1452 
1453  gOPT.gw151226_ifar = PP_GW151226_IFAR;
1454  gOPT.zl_max_ifar = PP_ZL_MAX_IFAR;
1455 }
1456 
1457 void PrintUserOptions() {
1458 
1459  cout.precision(14);
1460 
1461  TString sexit = gOPT.exit?"true":"false";
1462  TString bbh = gOPT.bbh?"true":"false";
1463 
1464  cout << endl;
1465 
1466  cout << "-----------------------------------------" << endl;
1467  cout << "PP IFAR config options " << endl;
1468  cout << "-----------------------------------------" << endl << endl;
1469 
1470  cout << "PP_RUN " << gOPT.run << endl;
1471  cout << "PP_SEARCH " << gOPT.search << endl;
1472 
1473  cout << "PP_NIFO " << gOPT.nifo << endl;
1474  cout << "PP_LAG " << gOPT.lag << endl;
1475  cout << "PP_SLAG " << gOPT.slag << endl;
1476  cout << "PP_CHUNK " << gOPT.chunk << endl;
1477 
1478  cout << "PP_TRIALS " << gOPT.trials << endl;
1479 
1480  cout << "PP_PFNAME " << gOPT.pfname << endl;
1481  cout << "PP_ZL_STYLE " << gOPT.zl_style << endl;
1482  cout << "PP_BELT_STYLE " << gOPT.belt_style << endl;
1483 
1484  cout << "PP_IFAR_MAX " << gOPT.ifar_max << " yr"<< endl;
1485  cout << "PP_RHO_MIN " << gOPT.rho_min << endl;
1486 
1487  cout << "PP_SFACTOR " << gOPT.sfactor << endl;
1488 
1489  cout << "PP_BBH " << bbh << endl;
1490  cout << "PP_EXIT " << sexit << endl;
1491 
1492  cout << "PP_GW151226_IFAR " << gOPT.gw151226_ifar << endl;
1493  cout << "PP_ZL_MAX_IFAR " << gOPT.zl_max_ifar << endl;
1494 
1495  cout << endl;
1496 }
1497 
1498 void ReadUserOptions(TString options) {
1499 
1500  // get plugin options
1501 
1502  if(TString(options)!="") {
1503 
1504  //cout << "PP_IFAR options : " << options << endl;
1505  TObjArray* token = TString(options).Tokenize(TString(' '));
1506  for(int j=0;j<token->GetEntries();j++) {
1507 
1508  TObjString* tok = (TObjString*)token->At(j);
1509  TString stok = tok->GetString();
1510 
1511  if(stok.Contains("run=")) {
1512  TString pp_run=stok;
1513  pp_run.Remove(0,pp_run.Last('=')+1);
1514  gOPT.run=pp_run;
1515  }
1516 
1517  if(stok.Contains("search=")) {
1518  TString pp_search=stok;
1519  pp_search.Remove(0,pp_search.Last('=')+1);
1520  gOPT.search=pp_search;
1521  }
1522 
1523  if(stok.Contains("nifo=")) {
1524  TString pp_nifo=stok;
1525  pp_nifo.Remove(0,pp_nifo.Last('=')+1);
1526  if(pp_nifo.IsDigit()) gOPT.nifo=pp_nifo.Atoi();
1527  }
1528 
1529  if(stok.Contains("lag=") && !stok.Contains("slag=")) {
1530  TString pp_lag=stok;
1531  pp_lag.Remove(0,pp_lag.Last('=')+1);
1532  if(pp_lag.IsDigit()) gOPT.lag=pp_lag.Atoi();
1533  }
1534 
1535  if(stok.Contains("slag=")) {
1536  TString pp_slag=stok;
1537  pp_slag.Remove(0,pp_slag.Last('=')+1);
1538  if(pp_slag.IsDigit()) gOPT.slag=pp_slag.Atoi();
1539  }
1540 
1541  if(stok.Contains("chunk=")) {
1542  TString pp_chunk=stok;
1543  pp_chunk.Remove(0,pp_chunk.Last('=')+1);
1544  if(pp_chunk.IsDigit()) gOPT.chunk=pp_chunk.Atoi();
1545  }
1546 
1547  if(stok.Contains("trials=")) {
1548  TString pp_trials=stok;
1549  pp_trials.Remove(0,pp_trials.Last('=')+1);
1550  if(pp_trials.IsDigit()) gOPT.trials=pp_trials.Atoi();
1551  }
1552 
1553  if(stok.Contains("pfname=")) {
1554  TString pp_pfname=stok;
1555  pp_pfname.Remove(0,pp_pfname.Last('=')+1);
1556  gOPT.pfname=pp_pfname;
1557  }
1558 
1559  if(stok.Contains("zl_style=")) {
1560  TString pp_zl_style=stok;
1561  pp_zl_style.Remove(0,pp_zl_style.Last('=')+1);
1562  gOPT.zl_style=pp_zl_style;
1563  }
1564 
1565  if(stok.Contains("belt_style=")) {
1566  TString pp_belt_style=stok;
1567  pp_belt_style.Remove(0,pp_belt_style.Last('=')+1);
1568  gOPT.belt_style=pp_belt_style;
1569  }
1570 
1571  if(stok.Contains("ifar_max=")) {
1572  TString pp_ifar_max=stok;
1573  pp_ifar_max.Remove(0,pp_ifar_max.Last('=')+1);
1574  if(pp_ifar_max.IsFloat()) gOPT.ifar_max=pp_ifar_max.Atof();
1575  }
1576 
1577  if(stok.Contains("rho_min=")) {
1578  TString pp_rho_min=stok;
1579  pp_rho_min.Remove(0,pp_rho_min.Last('=')+1);
1580  if(pp_rho_min.IsFloat()) gOPT.rho_min=pp_rho_min.Atof();
1581  }
1582 
1583  if(stok.Contains("sfactor=")) {
1584  TString pp_sfactor=stok;
1585  pp_sfactor.Remove(0,pp_sfactor.Last('=')+1);
1586  if(pp_sfactor.IsFloat()) gOPT.sfactor=pp_sfactor.Atof();
1587  }
1588 
1589  if(stok.Contains("bbh=")) {
1590  TString pp_bbh=stok;
1591  pp_bbh.Remove(0,pp_bbh.Last('=')+1);
1592  if(pp_bbh=="true") gOPT.bbh=true;
1593  if(pp_bbh=="false") gOPT.bbh=false;
1594  }
1595 
1596  if(stok.Contains("exit=")) {
1597  TString pp_exit=stok;
1598  pp_exit.Remove(0,pp_exit.Last('=')+1);
1599  if(pp_exit=="true") gOPT.exit=true;
1600  if(pp_exit=="false") gOPT.exit=false;
1601  }
1602 
1603  if(stok.Contains("gw151226_ifar=")) {
1604  TString pp_gw151226_ifar=stok;
1605  pp_gw151226_ifar.Remove(0,pp_gw151226_ifar.Last('=')+1);
1606  if(pp_gw151226_ifar.IsFloat()) gOPT.gw151226_ifar=pp_gw151226_ifar.Atof();
1607  }
1608 
1609  if(stok.Contains("zl_max_ifar=")) {
1610  TString pp_zl_max_ifar=stok;
1611  pp_zl_max_ifar.Remove(0,pp_zl_max_ifar.Last('=')+1);
1612  if(pp_zl_max_ifar.IsFloat()) gOPT.zl_max_ifar=pp_zl_max_ifar.Atof();
1613  }
1614  }
1615  }
1616 }
1617 
1618 void DumpLoudest(double obsTime) {
1619 
1620  if(gOPT.pfname=="") return;
1621 
1622  char ofname[1024];
1623  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1624  if(gOPT.bbh) sprintf(ofname,"%s_bbh_loudest.txt",PFNAME.Data());
1625  else sprintf(ofname,"%s_nobbh_loudest.txt",PFNAME.Data());
1626 
1627  char sout[1024];
1628 
1629  ofstream out;
1630  out.open(ofname,ios::out);
1631  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
1632  cout << "Create Output File : " << ofname << endl;
1633  out.precision(19);
1634 
1635  // write header
1636  sprintf(sout,"%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s",
1637  4, "#run", 6, "chunk", 18, "gps_time", 10, "known.BBH", 16, "ifar(sec)", 16, "ifar(year)", 16, "obs_time(sec)",
1638  16, "obs_time(days)", 16, "expected", 10, "observed", 11, "cumul.FAP", 8, "sigma");
1639  out << sout << endl;
1640 
1641  // write data
1642  int observed=1;
1643  for (int i=ZL_IFAR.size()-1; i>=0; i--) {
1644 
1645  double expected = obsTime/(ZL_IFAR[i]*YEAR);
1646 
1647  // compute FAP
1648  double FAP = 1.;
1649  for(int j=0;j<observed;j++) {
1650  FAP -= TMath::PoissonI(j,expected);
1651  //cout << " N : " << j << " MU : " << expected << " fap : " << fap << endl;
1652  }
1653 
1654  // compute sigma
1655  double sigma = sqrt(2)*TMath::ErfcInverse(FAP);
1656 
1657  TString bbh_name = GetNameBBH(ZL_GPS[i]);
1658  if(bbh_name=="") bbh_name=".";
1659 
1660  TString run;
1661  int chunkID = GetChunkID(ZL_GPS[i],run);
1662 
1663  sprintf(sout,"%*s %*d %*.6f %*s %*.1f %*.4f %*d %*.1f %*.9f %*d %*.3e %*.2f",
1664  4, run.Data(), 6, chunkID, 18, ZL_GPS[i], 10, bbh_name.Data(), 16,
1665  ZL_IFAR[i]*YEAR, 16, ZL_IFAR[i], 16, obsTime,
1666  16, obsTime/DAY, 16, expected, 10, observed, 11, FAP, 8, sigma);
1667  if(gOPT.zl_style=="marker") {
1668  out << sout << endl;
1669  observed++;
1670  } else {
1671  if(i%2==0) {out << sout << endl; observed++;}
1672  }
1673  }
1674 
1675  out.close();
1676 }
1677 
1678 int GetChunkID(double gps, TString& run) {
1679 
1680  for(int n=0;n<MAX_RUNS;n++) {
1681  for(int i=0;i<chunk_size[n];i++) {
1682  if(gps>chunk_start[n][i] && gps<chunk_stop[n][i]) {
1683  if(n==0) run="O1";
1684  if(n==1) run="O2";
1685  return chunk_id[n][i];
1686  }
1687  }
1688  }
1689 
1690  return -1;
1691 }
1692 
1693 void DumpPeriod(double* obsTime) {
1694 
1695  if(gOPT.pfname=="") return;
1696 
1697  char ofname[1024];
1698  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1699  if(gOPT.bbh) sprintf(ofname,"%s_bbh_period.txt",PFNAME.Data());
1700  else sprintf(ofname,"%s_nobbh_period.txt",PFNAME.Data());
1701 
1702  char sout[1024];
1703 
1704  ofstream out;
1705  out.open(ofname,ios::out);
1706  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
1707  cout << "Create Output File : " << ofname << endl;
1708  out.precision(19);
1709 
1710  // write header
1711  sprintf(sout,"%*s %*s %*s %*s %*s %*s %*s %*s",
1712  4, "#run", 20, "gps_start", 20, "date_start", 20, "gps_stop", 20, "date_stop", 20, "interval(days)", 20, "obs_time(sec)", 20, "obs_time(days)");
1713  out << sout << endl;
1714 
1715  for(int n=0;n<MAX_RUNS;n++) {
1716 
1717  TString RUN;
1718 
1719  if(n==0 && gOPT.run.Contains("O1")) RUN="O1";
1720  else
1721  if(n==1 && gOPT.run.Contains("O2")) RUN="O2";
1722  else continue;
1723 
1724  wat::Time beg_date(chunk_start[n][0]);
1725  wat::Time end_date(chunk_start[n][chunk_size[n]-1]);
1726 
1727  double interval = (end_date.GetGPS()-beg_date.GetGPS())/DAY;
1728 
1729  TString sbeg_date = beg_date.GetDateString();sbeg_date.Resize(19);
1730  TString send_date = end_date.GetDateString();send_date.Resize(19);
1731 
1732  sbeg_date.ReplaceAll(" ","-");
1733  send_date.ReplaceAll(" ","-");
1734 
1735  sprintf(sout,"%*s %*.0f %*s %*.0f %*s %*.2f %*d %*.1f", 4, RUN.Data(), 20, beg_date.GetGPS(), 20, sbeg_date.Data(),
1736  20, end_date.GetGPS(), 20, send_date.Data(), 20, interval, 20, obsTime[n]/gOPT.trials, 20, obsTime[n]/DAY/gOPT.trials);
1737  out << sout << endl;
1738  }
1739 
1740  out.close();
1741 }
1742 
1743 void DumpUserOptions() {
1744 
1745  if(gOPT.pfname=="") return;
1746 
1747  char ofname[1024];
1748  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1749  if(gOPT.bbh) sprintf(ofname,"%s_bbh_config.txt",PFNAME.Data());
1750  else sprintf(ofname,"%s_nobbh_config.txt",PFNAME.Data());
1751 
1752  ofstream out;
1753  out.open(ofname,ios::out);
1754  if(!out.good()) {cout << "DumpUserOptions : Error Opening File : " << ofname << endl;exit(1);}
1755  cout << "dump: " << ofname << endl;
1756 
1757  TString info="";
1758  TString tabs="\t\t\t\t";
1759 
1760  // cWB library
1761  char version[128];
1762  sprintf(version,"WAT Version : %s - SVN Revision : %s - Tag/Branch : %s",watversion('f'),watversion('r'),watversion('b'));
1763 
1764  // cWB config
1765 
1766  // get CWB_CONFIG
1767  char cwb_config_env[1024] = "";
1768  if(gSystem->Getenv("CWB_CONFIG")!=NULL) {
1769  strcpy(cwb_config_env,TString(gSystem->Getenv("CWB_CONFIG")).Data());
1770  }
1771  char cfg_version[256];
1772  TString cfg_branch = GetGitInfos("branch",cwb_config_env);
1773  TString cfg_tag = GetGitInfos("tag",cwb_config_env);
1774  TString cfg_diff = GetGitInfos("diff",cwb_config_env);
1775  if(cfg_branch!="") {
1776  if(cfg_diff!="") cfg_branch=cfg_branch+"/M";
1777  sprintf(cfg_version,"cWB Config Branch\t%s",cfg_branch.Data());
1778  } else if(cfg_tag!="") {
1779  if(cfg_diff!="") cfg_branch=cfg_branch+"/M";
1780  sprintf(cfg_version,"cWB Config Tag\t%s",cfg_tag.Data());
1781  } else {
1782  sprintf(cfg_version,"cWB Config\t%s","Undefined");
1783  }
1784 
1785  out << endl;
1786  out << "--------------------------------" << endl;
1787  out << "PP IFAR version " << endl;
1788  out << "--------------------------------" << endl;
1789  out << endl;
1790  out << version << endl;
1791  out << endl;
1792  out<< cfg_version <<endl;
1793  out<< "cWB Config Hash\t" << GetGitInfos("hash",cwb_config_env) <<endl;
1794  out << endl;
1795 
1796  out << "--------------------------------" << endl;
1797  out << "PP IFAR config options " << endl;
1798  out << "--------------------------------" << endl;
1799  out << endl;
1800 
1801  out << "pp_run " << gOPT.run << endl;
1802  info = "// run type: allowed types are O1/O2/O1O2 (Default: O2)";
1803  out << tabs << info << endl << endl;
1804 
1805  out << "pp_search " << gOPT.search << endl;
1806  info = "// search type: allowed types are BBH,BurstLF,BurstHF,BurstLD,IMBHB + bin1/bin2/... (Default: none)";
1807  out << tabs << info << endl << endl;
1808 
1809  out << "pp_nifo " << gOPT.nifo << endl;
1810  info = "// number of detectors (Default: 2)";
1811  out << tabs << info << endl << endl;
1812 
1813  out << "pp_lag " << gOPT.lag << endl;
1814  info = "// lag -> used only when the input file name is a root file (Default: -1)";
1815  out << tabs << info << endl << endl;
1816 
1817  out << "pp_slag " << gOPT.slag << endl;
1818  info = "// slag -> used only when the input file name is a root file (Default: -1)";
1819  out << tabs << info << endl << endl;
1820 
1821  out << "pp_chunk " << gOPT.chunk << endl;
1822  info = "// chunk -> used only when the input file name is a root file (Default: -1)";
1823  out << tabs << info << endl << endl;
1824 
1825  out << "pp_trials " << gOPT.trials << endl;
1826  info = "// trials factor (Default: 1)";
1827  out << tabs << info << endl << endl;
1828 
1829  out << "pp_pfname " << gOPT.pfname << endl;
1830  info = "// output plot file name, included the directory (Default: none)";
1831  out << tabs << info << endl << endl;
1832 
1833  out << "pp_zl_style " << gOPT.zl_style << endl;
1834  info = "// zero lag plot style: step/marker -> zero lag is displayed with step/marker (Default: step)";
1835  out << tabs << info << endl << endl;
1836 
1837  out << "pp_belt_style " << gOPT.belt_style << endl;
1838  info = "// belt plot style: step/continous -> belts are displayed with step/continuos poisson distribution (Default: continuos)";
1839  out << tabs << info << endl << endl;
1840 
1841  out << "pp_ifar_max " << gOPT.ifar_max << endl;
1842  info = "// max ifar (year) used for plot, if =-1 then it is the maximum ifar from root files (Default: -1)";
1843  out << tabs << info << endl << endl;
1844 
1845  out << "pp_rho_min " << gOPT.rho_min << endl;
1846  info = "// min rho[0/1] read from output root files (Default: 4.5)";
1847  out << tabs << info << endl << endl;
1848 
1849  out << "pp_sfactor " << gOPT.sfactor << endl;
1850  info = "// factor used to normalize the simulation events for the computation of astrophysical rates (Default: 1)";
1851  out << tabs << info << endl << endl;
1852 
1853  out << "bbh " << gOPT.bbh << endl;
1854  info = "// false/true -> if false the known bbh events are excluded from analysis (Default: false)";
1855  out << tabs << info << endl << endl;
1856 
1857  out << "exit " << gOPT.exit << endl;
1858  info = "// false/true -> if true program exit ayt the end of the execution (Default: true)";
1859  out << tabs << info << endl << endl;
1860 
1861  out << "pp_gw151226_ifar " << gOPT.gw151226_ifar << endl;
1862  info = "// GW151226 IFAR(YEARS): if > 0 the GW151226 event is added to the foreground (Default: 0, not added)";
1863  out << tabs << info << endl << endl;
1864 
1865  out << "pp_zl_max_ifar " << gOPT.zl_max_ifar << endl;
1866  info = "// ZL_MAX_IFAR(YEARS): if > 0 all foreground with IFAR>zl_max_ifar is set to zl_max_ifar (Default: 0, not used)";
1867  out << tabs << info << endl << endl;
1868 
1869 }
1870 
1871 // Dumps Background & Belts in ASCII
1872 void DumpBackground() {
1873 
1874  if(gOPT.pfname=="") return;
1875 
1876  char ofname[1024];
1877  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1878  if(gOPT.bbh) sprintf(ofname,"%s_bbh_background.txt",PFNAME.Data());
1879  else sprintf(ofname,"%s_nobbh_background.txt",PFNAME.Data());
1880 
1881  ofstream out;
1882  out.open(ofname,ios::out);
1883  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
1884  cout << "Create Output File : " << ofname << endl;
1885  out.precision(19);
1886 
1887  // write header
1888  out << "#background data : ifar(years), expected, lower_sigma1, lower_sigma2, lower_sigma3, upper_sigma1, upper_sigma2, upper_sigma3" << endl;
1889 
1890  for(int i=0;i<xIFAR.size();i++) {
1891 
1892  double med = bRATE[i];
1893 
1894  double l900 = yBKG[i]-eyRATEinf0[i];
1895  double u900 = yBKG[i]-eyRATEsup0[i];
1896 
1897  double l990 = yBKG[i]-eyRATEinf1[i];
1898  double u990 = yBKG[i]+eyRATEsup1[i];
1899 
1900  double l999 = yBKG[i]-eyRATEinf2[i];
1901  double u999 = yBKG[i]+eyRATEsup2[i];
1902 
1903  out << xIFAR[i] << " " << med
1904  << " " << l900 << " " << l990 << " " << l999
1905  << " " << u900 << " " << u990 << " " << u999
1906  << endl;
1907  }
1908 
1909  out.close();
1910 }
1911 
1912 // Dumps Foreground in ASCII format
1913 void DumpForeground() {
1914 
1915  if(gOPT.pfname=="") return;
1916 
1917  char ofname[1024];
1918  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1919  if(gOPT.bbh) sprintf(ofname,"%s_bbh_foreground.txt",PFNAME.Data());
1920  else sprintf(ofname,"%s_nobbh_foreground.txt",PFNAME.Data());
1921 
1922  ofstream out;
1923  out.open(ofname,ios::out);
1924  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
1925  cout << "Create Output File : " << ofname << endl;
1926  out.precision(19);
1927 
1928  // write header
1929  out << "#foreground data : ifar(years), observed" << endl;
1930 
1931  for(int i=0;i<ZL_NEVT.size();i++) {
1932 
1933  out << ZL_IFAR[i] << " " << ZL_NEVT[i] << endl;
1934  }
1935 
1936  out.close();
1937 }
1938 
1939 // Dumps Simulation & Belts in ASCII format
1940 void DumpSimulation(double ifar_cmax) {
1941 
1942  if(SIM_NEVT.size()<=0) return;
1943  if(gOPT.pfname=="") return;
1944  if(xSIM_BKG_IFAR.size()<=0) return;
1945 
1946  char ofname[1024];
1947  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
1948  if(gOPT.bbh) sprintf(ofname,"%s_bbh_simulation.txt",PFNAME.Data());
1949  else sprintf(ofname,"%s_nobbh_simulation.txt",PFNAME.Data());
1950 
1951  ofstream out;
1952  out.open(ofname,ios::out);
1953  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
1954  cout << "Create Output File : " << ofname << endl;
1955  out.precision(19);
1956 
1957  // write header
1958  out << "#simulation data : ifar(years), estimated+background, lower_90.0, upper_90.0" << endl;
1959 
1960  // compute the index which contains the maximum common IFAR (RATE_TAG)
1961  int ifar_cmax_index=0;
1962  for(int i=0;i<xSIM_BKG_IFAR.size();i++) if(xSIM_BKG_IFAR.data[i]<pow(10,ifar_cmax)/YEAR) ifar_cmax_index=i;
1963 
1964  for(int i=0;i<ifar_cmax_index;i++) {
1965 
1966  double med = xSIM_BKG_IFAR[i];
1967 
1968  double l90 = med-eySIM_BKG_inf0[i];
1969  double u90 = med+eySIM_BKG_sup0[i];
1970 
1971  out << SIM_IFAR[i] << " " << med << " " << l90 << " " << u90 << endl;
1972  }
1973 
1974  out.close();
1975 }
1976 
egr1
TGraphAsymmErrors * egr1
Definition: Make_PP_IFAR.C:157
strcpy
strcpy(analysis,"2G")
PP_NIFO
#define PP_NIFO
Definition: Make_PP_IFAR.C:44
GetPercentiles
int GetPercentiles(double mu, double *q)
Definition: Make_PP_IFAR.C:1222
chunk_size
int chunk_size[MAX_RUNS]
Definition: Make_PP_IFAR.C:200
CHUNK_FILE_LIST
#define CHUNK_FILE_LIST
Definition: Make_PP_IFAR.C:91
InitBBH
void InitBBH()
Definition: Make_PP_IFAR.C:1356
latex
TLatex * latex
Definition: Make_PP_IFAR.C:155
ZL_GPS
wavearray< double > ZL_GPS
Definition: Make_PP_IFAR.C:180
DumpLoudest
void DumpLoudest(double obsTime)
Definition: Make_PP_IFAR.C:1618
chunk_id
int chunk_id[MAX_RUNS][CHUNK_MAX_SIZE]
Definition: Make_PP_IFAR.C:201
PP_BELT_STYLE
#define PP_BELT_STYLE
Definition: Make_PP_IFAR.C:53
O1_BBH_NAME
vector< TString > O1_BBH_NAME
Definition: Make_PP_IFAR.C:192
DumpPeriod
void DumpPeriod(double *obsTime)
Definition: Make_PP_IFAR.C:1693
eySIM_BKG_sup0
wavearray< double > eySIM_BKG_sup0
Definition: Make_PP_IFAR.C:176
CHUNK_MAX_SIZE
#define CHUNK_MAX_SIZE
Definition: Make_PP_IFAR.C:93
ReadChunkList
int ReadChunkList(TString ifile, int *chunk=NULL, double *start=NULL, double *stop=NULL)
Definition: ReadChunkList.C:4
n
shift breaksw case n
Definition: cwb_clchunk.csh:70
PP_ZL_STYLE
#define PP_ZL_STYLE
Definition: Make_PP_IFAR.C:52
MAX_RUNS
#define MAX_RUNS
Definition: Make_PP_IFAR.C:70
PP_ZL_MAX_IFAR
#define PP_ZL_MAX_IFAR
Definition: Make_PP_IFAR.C:64
ySIM_BKG
wavearray< double > ySIM_BKG
Definition: Make_PP_IFAR.C:172
O1_BBH_GPS
vector< double > O1_BBH_GPS
Definition: Make_PP_IFAR.C:193
yBKG
wavearray< double > yBKG
Definition: Make_PP_IFAR.C:165
PP_PFNAME
#define PP_PFNAME
Definition: Make_PP_IFAR.C:51
arrow
TArrow * arrow
Definition: Make_PP_IFAR.C:154
leg0
TLegend * leg0
Definition: Make_PP_IFAR.C:151
GetCutBBH
TString GetCutBBH()
Definition: Make_PP_IFAR.C:1384
grsim
TGraph * grsim
Definition: Make_PP_IFAR.C:147
O2_BBH_NAME
vector< TString > O2_BBH_NAME
Definition: Make_PP_IFAR.C:196
MACRO_READ_CHUNK
#define MACRO_READ_CHUNK
Definition: Make_PP_IFAR.C:92
chunk
int chunk[CHUNK_MAX_SIZE]
Definition: CreateChunkCuts.C:14
PP_IFAR_MAX
#define PP_IFAR_MAX
Definition: Make_PP_IFAR.C:55
sprintf
sprintf(tag,"wave_%s", data_label)
eyRATEsup0
wavearray< double > eyRATEsup0
Definition: Make_PP_IFAR.C:185
ResetUserOptions
void ResetUserOptions()
Definition: Make_PP_IFAR.C:1429
PP_GW151226_IFAR
#define PP_GW151226_IFAR
Definition: Make_PP_IFAR.C:63
grsimbkg
TGraph * grsimbkg
Definition: Make_PP_IFAR.C:148
GetNameBBH
TString GetNameBBH(double gps)
Definition: Make_PP_IFAR.C:1407
gr
TGraph * gr
Definition: Make_PP_IFAR.C:146
DumpSimulation
void DumpSimulation(double ifar_cmax)
Definition: Make_PP_IFAR.C:1940
GetChunkID
int GetChunkID(double gps, TString &run)
Definition: Make_PP_IFAR.C:1678
Make_PP_IFAR
void Make_PP_IFAR(TString ifName, TString options)
Definition: Make_PP_IFAR.C:234
DAY
#define DAY
Definition: Make_PP_IFAR.C:88
size
Long_t size
Definition: cwb_mklinks.C:8
gOPT
uoptions gOPT
Definition: Make_PP_IFAR.C:141
PP_SLAG
#define PP_SLAG
Definition: Make_PP_IFAR.C:46
PP_SFACTOR
#define PP_SFACTOR
Definition: Make_PP_IFAR.C:58
IFAR_NBIN
#define IFAR_NBIN
Definition: Make_PP_IFAR.C:73
eyRATEinf1
wavearray< double > eyRATEinf1
Definition: Make_PP_IFAR.C:186
egr2
TGraphAsymmErrors * egr2
Definition: Make_PP_IFAR.C:158
FAP2
#define FAP2
Definition: Make_PP_IFAR.C:86
GetPoissonNsup
int GetPoissonNsup(double mu, double fap)
Definition: Make_PP_IFAR.C:1234
egrsimbkg0
TGraphAsymmErrors * egrsimbkg0
Definition: Make_PP_IFAR.C:149
leg
TLegend * leg
Definition: Make_PP_IFAR.C:152
O2_BBH_GPS
vector< double > O2_BBH_GPS
Definition: Make_PP_IFAR.C:197
xFAR
wavearray< double > xFAR
Definition: Make_PP_IFAR.C:163
eyRATEinf0
wavearray< double > eyRATEinf0
Definition: Make_PP_IFAR.C:184
PP_RHO_MIN
#define PP_RHO_MIN
Definition: Make_PP_IFAR.C:56
egr0
TGraphAsymmErrors * egr0
Definition: Make_PP_IFAR.C:156
eySIM_BKG_inf0
wavearray< double > eySIM_BKG_inf0
Definition: Make_PP_IFAR.C:175
chunk_start
double chunk_start[MAX_RUNS][CHUNK_MAX_SIZE]
Definition: Make_PP_IFAR.C:202
name
par[0] name
Definition: CWB_Plugin_EBBH_Config.C:14
GetLiveTime
void GetLiveTime(int nwave, TChain **live, int *lag, int *slag, int *bin, TString *run, double *liveZero, double *liveTot)
Definition: Make_PP_IFAR.C:1321
token
TObjArray * token
Definition: cwb_mklinks.C:77
tmp
char tmp[256]
Definition: cwb_mklinks.C:73
pp_rho_min
pp_rho_min
Definition: user_pparameters.C:15
exit
exit(0)
SIM_BKG_NEVT
wavearray< double > SIM_BKG_NEVT
Definition: Make_PP_IFAR.C:169
q
shift breaksw case q
Definition: cwb_obchunk.csh:120
DumpBackground
void DumpBackground()
Definition: Make_PP_IFAR.C:1872
canvas
TCanvas * canvas
Definition: Make_PP_IFAR.C:145
GetPoissonNinf
int GetPoissonNinf(double mu, double fap)
Definition: Make_PP_IFAR.C:1246
PP_RUN
#define PP_RUN
Definition: Make_PP_IFAR.C:41
YEAR
#define YEAR
Definition: Make_PP_IFAR.C:89
eyRATEsup1
wavearray< double > eyRATEsup1
Definition: Make_PP_IFAR.C:187
marker
TMarker * marker
Definition: Make_PP_IFAR.C:153
exSIM_BKG_IFAR
wavearray< double > exSIM_BKG_IFAR
Definition: Make_PP_IFAR.C:173
FAP0
#define FAP0
Definition: Make_PP_IFAR.C:84
PP_EXIT
#define PP_EXIT
Definition: Make_PP_IFAR.C:61
PP_LAG
#define PP_LAG
Definition: Make_PP_IFAR.C:45
m
shift breaksw case m
Definition: cwb_obchunk.csh:110
PP_BBH
#define PP_BBH
Definition: Make_PP_IFAR.C:60
eyRATEinf2
wavearray< double > eyRATEinf2
Definition: Make_PP_IFAR.C:188
SIM_NEVT
wavearray< double > SIM_NEVT
Definition: Make_PP_IFAR.C:168
SIM_IFAR
wavearray< double > SIM_IFAR
Definition: Make_PP_IFAR.C:167
ZL_IFAR
wavearray< double > ZL_IFAR
Definition: Make_PP_IFAR.C:178
exFAR
wavearray< double > exFAR
Definition: Make_PP_IFAR.C:183
PrintUserOptions
void PrintUserOptions()
Definition: Make_PP_IFAR.C:1457
PP_CHUNK
#define PP_CHUNK
Definition: Make_PP_IFAR.C:47
GetGPSBBH
double GetGPSBBH(TString name)
Definition: Make_PP_IFAR.C:1418
value
par[0] value
Definition: CWB_Plugin_EBBH_Config.C:35
MAX_LIST_SIZE
#define MAX_LIST_SIZE
Definition: Make_PP_IFAR.C:72
DumpUserOptions
void DumpUserOptions()
Definition: Make_PP_IFAR.C:1743
ReadUserOptions
void ReadUserOptions(TString options)
Definition: Make_PP_IFAR.C:1498
xIFAR
wavearray< double > xIFAR
Definition: Make_PP_IFAR.C:162
xSIM_BKG_IFAR
wavearray< double > xSIM_BKG_IFAR
Definition: Make_PP_IFAR.C:174
ZL_NEVT
wavearray< double > ZL_NEVT
Definition: Make_PP_IFAR.C:179
PP_SEARCH
#define PP_SEARCH
Definition: Make_PP_IFAR.C:42
ReadFileList
int ReadFileList(TString ifName, TChain **live, TChain **wave, TChain **mdc, int *lag, int *slag, int *chunk, int *bin, TString *run)
Definition: Make_PP_IFAR.C:1258
eyRATEsup2
wavearray< double > eyRATEsup2
Definition: Make_PP_IFAR.C:189
gr0
TGraph * gr0
Definition: Make_PP_IFAR.C:150
IFAR
#define IFAR
Definition: MakePlotEgw.C:13
chunk_stop
double chunk_stop[MAX_RUNS][CHUNK_MAX_SIZE]
Definition: Make_PP_IFAR.C:203
FAP1
#define FAP1
Definition: Make_PP_IFAR.C:85
MakePlot
void MakePlot(TString title, double obsTime, double ifar_cmax)
Definition: Make_PP_IFAR.C:933
bRATE
wavearray< double > bRATE
Definition: Make_PP_IFAR.C:164
PP_TRIALS
#define PP_TRIALS
Definition: Make_PP_IFAR.C:49
search
search
Definition: user_parameters.C:10
DumpForeground
void DumpForeground()
Definition: Make_PP_IFAR.C:1913