Make_PP_IFAR.C

Go to the documentation of this file.

/*
# Copyright (C) 2019 Gabriele Vedovato
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/
 
 
 
// ---------------------------------------------------------------------------------
// INFOS
// ---------------------------------------------------------------------------------
 
// this macro produce the plot number of events in the observation period vs IFAR
 
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// The macro Make_PP_IFAR_GW151226_RATE.C is the standard macro in cWB config master.2.8 + add GW151226 event in O1 data
// set a superior limit @ GW170104 IFAR -> set a limit for BBH IFAR, ecluded BBH IFAR estimated as lower limits (IFAR_TAG) 
// compute the index which contains the maximum common IFAR for rates -> limits the sim green plot (RATE_TAG)
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
// input list chunk root files format
// 
// wave_merged_root_file_name_with_ifar.root lag slag chunk chunkLabel
// Ex: K02,K03 open box 
//     wave_O2_K02_merged_root_file_name_with_ifar.root 0 0 2 K02
//     wave_O2_K03_merged_root_file_name_with_ifar.root 0 0 3 K03
 
// ---------------------------------------------------------------------------------
// HOW RUN THIS MACRO (EXAMPLES)
// ---------------------------------------------------------------------------------
 
/*
 
// uses a list of input root merged files generated by different chunks
root 'Make_PP_IFAR.C("pp_lag0_slag0_M2.lst", "--search=BurstLF:bin2   \\
                      --exit=1 --pfname=pp_plot_bin2_kall_lag0_slag0_M2.png")'
 
// uses a single chunk root merged file
root 'Make_PP_IFAR.C("wave_merged_output_root_file.root", "--search=BurstLF:bin2  --lag=0 --slag=0 \\
                      --chunk=0 --exit=1 --pfname=pp_plot_bin2_kall_lag0_slag0_M2.png")'
 
*/
 
#include "ReadChunkList.C"
 
// ---------------------------------------------------------------------------------
// USER DEFINES : Default values used to initialize the user options input parameters
// ---------------------------------------------------------------------------------`
 
#define PP_RUN      "O3"
#define PP_SEARCH   ""
 
#define PP_NIFO           2
#define PP_LAG      -1
#define PP_SLAG        -1
#define PP_CHUNK    -1
 
#define PP_TRIALS         1
 
#define PP_PFNAME   ""
#define PP_ZL_STYLE    "step"
#define PP_BELT_STYLE     "continuos"
 
#define PP_IFAR_MAX    -1
#define PP_RHO_MIN     4.5
 
#define PP_SFACTOR     1
 
#define PP_BBH      false
#define PP_EXIT        true
 
#define PP_GW151226_IFAR  0
#define PP_ZL_MAX_IFAR    0
 
// ---------------------------------------------------------------------------------`
// DEFINES
// ---------------------------------------------------------------------------------`
 
#define MAX_RUNS  3  // O1,O2,O3
 
#define MAX_LIST_SIZE   100   // maximum number of input root files in the input file list
#define IFAR_NBIN       2000    // number of bins used for the computation of belts and background
                                // NOTE: if such value is low it introduces some discretization artifacts for low IFAR values
 
//#define FAP0 0.1      // FAP used to define the first belt  : SIGMA = 1.644
//#define FAP1 0.01     // FAP used to define the second belt : SIGMA = 2.575
//#define FAP2 0.001    // FAP used to define the third belt  : SIGMA = 3.290
 
// 
// Probability = TMath::Erf(Sigma/sqrt(2))
// Sigma = TMath::ErfInverse(Probability)*sqrt(2)
 
#define FAP0   (1.-0.682689)     // FAP used to define the first belt  : SIGMA = 1.
#define FAP1   (1.-0.954499)     // FAP used to define the second belt : SIGMA = 2.
#define FAP2   (1.-0.997300)     // FAP used to define the third belt  : SIGMA = 3.
 
#define DAY    (24.*3600.)
#define YEAR   (24.*3600.*365.)
 
#define CHUNK_FILE_LIST         "Chunk_List.txt"
#define CHUNK_MAX_SIZE          100
 
//#define APPLY_MAX_SIM_IFAR  // if enabled the detected & sim events are cut to the common superior limit
 
//#define PLOT_FOREGROUND_NOBBH  // plot ZL BBH && ZL NOBBH
 
#define PA_FILE_LIST            "Public_Alerts.lst"
 
// ---------------------------------------------------------------------------------
// INPUT USER OPTIONS
// ---------------------------------------------------------------------------------
 
struct uoptions {
  TString   run;              // run type: allowed types are O1/O2/O1O2/O3
  TString   search;           // search type: allowed types are BBH,BurstLF,BurstHF,BurstLD,IMBHB + bin1/bin2/...
 
  int       nifo;             // number of detectors 
  int       lag;              // lag   -> used only when the input file name is a root file 
  int       slag;             // slag  -> used only when the input file name is a root file 
  int       chunk;            // chunk -> used only when the input file name is a root file 
 
  int       trials;           // trials factor
 
  TString   pfname;           // output plot file name, included the directory 
  TString   zl_style;         // zero lag plot style: step/marker -> zero lag is displayed with step/marker
  TString   belt_style;       // belt plot style: step/continuos -> belts are displayed with step/continuos poisson distribution
 
  double    ifar_max;         // max ifar (year) used for plot, if =-1 then it is the maximum ifar from root files
  double    rho_min;       // min rho[0/1] read from output root files
 
  double    sfactor;       // factor used to normalize the simulation events for the computation of astrophysical rates
                              // SIM_NEVT/=gOPT.sfactor
  bool      bbh;        // if false the known bbh events are excluded from analysis
  bool      exit;       // if true program exit ayt the end of the execution
 
  double    gw151226_ifar;    // GW151226 IFAR(YEARS): if > 0 the GW151226 event is added to the foreground (def=0, not added)
               // This is a special option introduced for the GW151226 event which has been found only into the extended segments
 
  double    zl_max_ifar;      // It is used to set the maximum IFAR of foreground events (YEARS). 
                              // The value to be used is the minimum common foreground IFAR for which the IFAR is well estimated    
                              // The foreground with IFAR>zl_max_ifar is set to zl_max_ifar 
};
 
 
// ---------------------------------------------------------------------------------
// Global Variables
// ---------------------------------------------------------------------------------
 
// global User Options
 
uoptions   gOPT;        
 
// plot objects
 
TCanvas* canvas;
TGraph*  gr;
TGraph*  grsim;
TGraph*  grsimbkg;
TGraphAsymmErrors* egrsimbkg0; 
TGraph*  gr0;
TLegend* leg0;
TLegend* leg;
TMarker* marker;
TArrow* arrow;
TLatex *latex;
TGraphAsymmErrors* egr0; 
TGraphAsymmErrors* egr1; 
TGraphAsymmErrors* egr2; 
 
// wavearray used for plots
 
wavearray<double> xIFAR;
wavearray<double> xFAR;
wavearray<double> bRATE;
wavearray<double> yBKG;
 
wavearray<double> SIM_IFAR;   // sorted ifar simulated events
wavearray<double> SIM_NEVT;   // number of simulated events with ifar > SIM_IFAR[k]
wavearray<double> SIM_BKG_NEVT; // number of simulated events + expected from BKG with ifar > SIM_IFAR[k]
 
// SIM+BKG FAP0 belts
wavearray<double> ySIM_BKG;
wavearray<double> exSIM_BKG_IFAR; 
wavearray<double> xSIM_BKG_IFAR;
wavearray<double> eySIM_BKG_inf0;
wavearray<double> eySIM_BKG_sup0;
 
wavearray<double> ZL_IFAR;    // sorted ifar events
wavearray<double> ZL_NEVT;    // number of events with ifar > ZL_IFAR[k]
wavearray<double> ZL_GPS;     // event GPS time
 
// FAP0/FAP1/FAP2 error belts
wavearray<double> exFAR; 
wavearray<double> eyRATEinf0;
wavearray<double> eyRATEsup0;
wavearray<double> eyRATEinf1;
wavearray<double> eyRATEsup1;
wavearray<double> eyRATEinf2;
wavearray<double> eyRATEsup2;
 
// O1 known BBH events
vector<TString> O1_BBH_NAME;
vector<double>  O1_BBH_GPS;
 
// O2 known BBH events
vector<TString> O2_BBH_NAME;
vector<double>  O2_BBH_GPS;
 
// O3 known BBH events
vector<TString> O3_BBH_NAME;
vector<double>  O3_BBH_GPS;
 
// chunks id,start,stop
int    chunk_size[MAX_RUNS];  
int    chunk_id[MAX_RUNS][CHUNK_MAX_SIZE];
double chunk_start[MAX_RUNS][CHUNK_MAX_SIZE];
double chunk_stop[MAX_RUNS][CHUNK_MAX_SIZE];
 
// ---------------------------------------------------------------------------------
// FUNCTIONS
// ---------------------------------------------------------------------------------
 
void ResetUserOptions();
void ReadUserOptions(TString options);
void PrintUserOptions();
void DumpUserOptions();
 
void    InitBBH();            // init array Ox_BBH_NAME,Ox_BBH_GPS with known BBH events
TString GetNameBBH(double gps);        // return the BBH name
double  GetGPSBBH(TString name);    // return the BBH GPS
TString GetCutBBH();          // return the cut string used to exclude BBH events from tree selection
int     ReadPublicAlerts();         // read public alerts, implemented only for O3  
 
int  GetChunkID(double gps, TString& run);   // return chunk ID for a gine gps time
int  GetPercentiles(double mu, double* q);   // return percentiles in q for FAP0,FAP1,FAP2 (continuos)
int  GetPoissonNsup(double mu, double fap);  // return sup number expeted events of for a gine fap
int  GetPoissonNinf(double mu, double fap);  // return inf number expeted events of for a gine fap
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
void GetLiveTime(int nwave, TChain** live, int* lag, int* slag, int* bin, TString* run, double* liveZero, double* liveTot);   // retun live time 
 
void MakePlot(TString title, double obsTime, double ifar_cmax);   // make final IFAR plot
void DumpLoudest(double obsTime);      // dump to ascii file the loudest zero lag event list
void DumpPeriod(double* obsTime);
void DumpBackground();           // Dumps Background & Belts in ASCII format
void DumpForeground();           // Dumps Foreground in ASCII format
void DumpSimulation(double ifar_cmax);       // Dumps Simulation & Belts in ASCII format 
 
 
void Make_PP_IFAR(TString ifName, TString options) {
 
  // ----------------------------------------------------- 
  // Read Input Options
  // ----------------------------------------------------- 
 
  ResetUserOptions();
  ReadUserOptions(options);
  PrintUserOptions();
 
  if((gOPT.run!="O1")&&(gOPT.run!="O2")&&(gOPT.run!="O1O2")&&(gOPT.run!="O3")) {
    cout << "Make_PP_IFAR - Error : run option not allowed, must be O1/O2/O1O2/O3" << endl;
    gSystem->Exit(1);
  } 
 
  TString SEARCH = gOPT.search; 
  if(gOPT.search.Contains(":")) SEARCH = gOPT.search(0,gOPT.search.Index(":",0));
  if((SEARCH!="BBH")&&(SEARCH!="IMBHB")&&(SEARCH!="BurstLF")&&(SEARCH!="BurstHF")&&(SEARCH!="BurstLD")) {
    cout << "Make_PP_IFAR - Error : search option not allowed, must be BBH/IMBHB/BurstLF/BurstHF/BurstLD + ':bin1/bin2/...'" << endl;
    gSystem->Exit(1);
  } 
 
  if(gOPT.trials<=0) {
    cout << "Make_PP_IFAR - Error : trials value not allowed, must be >0" << endl;
    gSystem->Exit(1);
  } 
 
  if(gOPT.sfactor<=0) {
    cout << "Make_PP_SFACTOR - Error : sfactor value not allowed, must be >0" << endl;
    gSystem->Exit(1);
  } 
 
  if(gOPT.nifo!=2 && gOPT.nifo!=3) {
    cout << "Make_PP_IFAR - Error : nifo value not allowed, must 2/3" << endl;
    gSystem->Exit(1);
  } 
 
  if((gOPT.zl_style!="step")&&(gOPT.zl_style!="marker")) {
    cout << "Make_PP_IFAR - Error : zl_style value not allowed, must be step/marker" << endl;
    gSystem->Exit(1);
  } 
 
  if((gOPT.belt_style!="step")&&(gOPT.belt_style!="continuos")) {
    cout << "Make_PP_IFAR - Error : belt_style value not allowed, must be step/continuos" << endl;
    gSystem->Exit(1);
  } 
 
  // check if is a root file
  TString fext = ifName(ifName.Last('.')+1,4);   
  if(fext=="root") {
    if(gOPT.lag<0 || gOPT.slag<0 || gOPT.chunk<0) {
      cout << "Make_PP_IFAR - Error : when input file ext=root lag,slag,chunk must be >=0" << endl;
      gSystem->Exit(1);
    }
  }
 
  DumpUserOptions();
 
  if(gOPT.exit) gROOT->SetBatch(kTRUE);   // NOTE: batch mode save plot with high quality !!!
 
  // initialize wavearray
  xIFAR.resize(IFAR_NBIN);
  xFAR.resize(IFAR_NBIN);
  bRATE.resize(IFAR_NBIN);
  yBKG.resize(IFAR_NBIN);
 
  exFAR.resize(IFAR_NBIN); exFAR=0;
  eyRATEinf0.resize(IFAR_NBIN);
  eyRATEsup0.resize(IFAR_NBIN);
  eyRATEinf1.resize(IFAR_NBIN);
  eyRATEsup1.resize(IFAR_NBIN);
  eyRATEinf2.resize(IFAR_NBIN);
  eyRATEsup2.resize(IFAR_NBIN);
 
  ySIM_BKG.resize(IFAR_NBIN); 
  exSIM_BKG_IFAR.resize(IFAR_NBIN); exSIM_BKG_IFAR=0;
  xSIM_BKG_IFAR.resize(IFAR_NBIN);
  eySIM_BKG_inf0.resize(IFAR_NBIN);
  eySIM_BKG_sup0.resize(IFAR_NBIN);
 
  InitBBH();   // init known BBH 
 
 
  // ----------------------------------------------------- 
  // Read Chunk start stop GPS times
  // ----------------------------------------------------- 
 
  // get CWB_CONFIG
  char cwb_config_env[1024] = "";
  if(gSystem->Getenv("CWB_CONFIG")!=NULL) {
    strcpy(cwb_config_env,TString(gSystem->Getenv("CWB_CONFIG")).Data());
  } 
 
  for(int n=0;n<MAX_RUNS;n++) {  
 
    TString RUN;
 
    if(n==0 && gOPT.run.Contains("O1")) RUN="O1"; 
    else 
    if(n==1 && gOPT.run.Contains("O2")) RUN="O2"; 
    else 
    if(n==2 && gOPT.run.Contains("O3")) RUN="O3"; 
    else continue;
 
    char chunk_file_list[1024];
    sprintf(chunk_file_list,"%s/%s/CHUNKS/%s/%s",cwb_config_env,RUN.Data(),SEARCH.Data(),CHUNK_FILE_LIST);
    cout << chunk_file_list << endl;
  
    cout << endl;
    cout << "----------------------------------------------------------------------------------------------------" << endl;
    cout << "chunk" << "\t\t" << "start" << "\t\t" << "stop" << "\t\t" 
       << " days\t\t" << "beg-date" << " - " << "end-date" << "\t" << "run" << endl;
    cout << "----------------------------------------------------------------------------------------------------" << endl;
    cout << endl;
  
    chunk_size[n] = ReadChunkList(chunk_file_list, chunk_id[n], chunk_start[n], chunk_stop[n]);
  
    cout << "----------------------------------------------------------------------------------------------------" << endl;
    cout << endl;
  }
 
  // ----------------------------------------------------- 
  // open wave root trees
  // ----------------------------------------------------- 
 
  int rho_id;
  if(SEARCH=="BBH" || SEARCH=="IMBHB")    rho_id=1;
  else               rho_id=0;
 
  TString run[MAX_LIST_SIZE];
  int bin[MAX_LIST_SIZE];
  int chunk[MAX_LIST_SIZE];
  int lag[MAX_LIST_SIZE];
  int slag[MAX_LIST_SIZE];
 
  TChain* live[MAX_LIST_SIZE];
  for(int i=0;i<MAX_LIST_SIZE;i++) live[i] = new TChain("liveTime");
 
  TChain* wave[MAX_LIST_SIZE];
  for(int i=0;i<MAX_LIST_SIZE;i++) wave[i] = new TChain("waveburst");
 
  TChain* mdc[MAX_LIST_SIZE];
  for(int i=0;i<MAX_LIST_SIZE;i++) mdc[i] = new TChain("mdc");
 
  int nwave = 1; 
  if(fext=="root") {
    chunk[0]      = gOPT.chunk;
    lag[0]        = gOPT.lag;
    slag[0]       = gOPT.slag;
    run[0]        = gOPT.run; 
    bin[0]        = 0; 
    TString fwave = ifName;
    TString flive = ifName;
    flive.ReplaceAll("wave_","live_");
    wave[0]->Add(fwave);
    live[0]->Add(flive);
    wave[0]->SetTitle(fwave);
    live[0]->SetTitle(flive);
  } else {
    // Read File List & Fill TChain
    nwave = ReadFileList(ifName, live, wave, mdc, lag, slag, chunk, bin, run);
    gOPT.lag   = lag[0];
    gOPT.slag  = slag[0];
    gOPT.chunk = chunk[0];
  }
  cout << endl;
  if(nwave<=0) {
    cout << "Make_PP_IFAR - Error : no files read from input list, check input list format" << endl;
    gSystem->Exit(1);
  }
 
  // check if slag and lag are consistent
  int xlag=lag[0];
  int xslag=slag[0];
  for(int i=1;i<nwave;i++) {
    if(bin[i]<0) continue;    // simulation
    if(lag[i]!=xlag || slag[i]!=xslag) {
      cout << "Make_PP_IFAR - Error : lag,slag defined in input list root files are not consistent" << endl;
      gSystem->Exit(1);
    }
  }
 
  // check duplicated entries
  for(int i=0;i<nwave;i++) {
    if(bin[i]<0) continue;    // simulation
    for(int j=0;j<nwave;j++) {
      if(bin[j]<0) continue;  // simulation
      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()))) {
        cout << "Make_PP_IFAR - Error : input file list contains duplicated entries for lag,slag,chunk,bin = " 
             << lag[i] << "," << slag[i] << "," << chunk[i] << "," << bin[i] << endl;
        cout << "                     : " 
             << wave[i]->GetTitle() << endl;
        gSystem->Exit(1);
      }
    }
  }
 
  // check consistency bin1 bin2 entries
  vector<int> vbin;  // contains the list of unique bins
  for(int i=0;i<nwave;i++) {
    if(bin[i]<0) continue;    // simulation
    bool exist=false;
    for(int j=0;j<vbin.size();j++) if(bin[i]==vbin[j]) exist=true;
    if(!exist) vbin.push_back(bin[i]);
  }
  for(int i=0;i<nwave;i++) {
    if(bin[i]<0) continue;    // simulation
    int nbins=1;
    for(int j=0;j<nwave;j++) {
      if(bin[j]<0) continue;  // simulation
      if((i!=j) && (lag[i]==lag[j] && slag[i]==slag[j] && chunk[i]==chunk[j])) {
        for(int k=0;k<vbin.size();k++) if(bin[j]!=vbin[k]) nbins++;
      }
    }
    if(nbins!=vbin.size()) {
      cout << "Make_PP_IFAR - Error : input file list must contains all bins for each chunk: "  
           << wave[i]->GetTitle() << endl;
      cout << "                       lag,slag,chunk,bin = " 
           << lag[i] << "," << slag[i] << "," << chunk[i] << "," << bin[i] << endl;
      gSystem->Exit(1);
    }
  }
 
  // ----------------------------------------------------- 
  // compute live time
  // ----------------------------------------------------- 
 
  double liveTot[MAX_RUNS];
  double liveZero[MAX_RUNS];
  for(int n=0;n<MAX_RUNS;n++) {liveTot[n]=0;liveZero[n]=0;}  
 
  GetLiveTime(nwave, live, lag, slag, bin, run, liveZero, liveTot);
 
  double obsTime=0;
  for(int n=0;n<MAX_RUNS;n++) {  
 
    obsTime+=liveTot[n]/gOPT.trials;   // apply trials factor
 
    TString RUN;
    if(n==0) RUN="O1";
    if(n==1) RUN="O2";
    if(n==2) RUN="O3";
    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);
  }
 
  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);
 
 
  // ======================================================================= 
  // DETECTED EVENTS STUFF
  // ======================================================================= 
 
  // ----------------------------------------------------- 
  // Get BBH cut string
  // ----------------------------------------------------- 
 
  TString CutBBH = GetCutBBH();
 
  // ----------------------------------------------------- 
  // compute min.max ifar in wave files
  // ----------------------------------------------------- 
 
  double ifar_cmin = 0;    // minimum common IFAR
  double ifar_cmax = 1e100;     // maximum common IFAR (used for rate plots) (RATE_TAG)
  double ifar_max  = 0;
  double value;
  for(int m=0;m<nwave;m++) {
    if(bin[m]<0) continue;    // simulation
    int isize = (int)wave[m]->GetEntries();
    wave[m]->SetEstimate(isize);
    char cut[64]; sprintf(cut,"rho[%d]>%f",rho_id,gOPT.rho_min);
    wave[m]->Draw("ifar",cut,"goff",isize);
    double* ifar = wave[m]->GetV1();
    int sel_entries = wave[m]->GetSelectedRows();
    value = log10(TMath::MinElement(sel_entries,ifar));
    if(fext=="root")
      cout << "tMax events per year : " << (1./pow(10,value))*YEAR << endl;
    else
      cout << "chunk\t" << chunk[m] << "\tbin\t" << bin[m] << "\tMax events per year : " << (1./pow(10,value))*YEAR << endl;
    if(value>ifar_cmin) ifar_cmin=value;        // this condition is used to select the common IFAR inferior limit
    value = log10(TMath::MaxElement(sel_entries,ifar));
    if(value>ifar_max) ifar_max=value;          // this condition is used to select the maximum IFAR limit
    if(value<ifar_cmax) ifar_cmax=value;        // this condition is used to select the common IFAR superior limit (RATE_TAG)
  }
 
  cout << endl;
  cout << "ifar_cmin : " << pow(10,ifar_cmin)/YEAR << " years" << endl;
  cout << "ifar_cmax : " << pow(10,ifar_cmax)/YEAR << " years" << endl; // (RATE_TAG)
  cout << "ifar_max  : " << pow(10,ifar_max)/YEAR  << " years" << endl;
 
  double inf = ifar_cmin;
  double sup = TMath::Nint(ifar_max+0.5);
 
  TString Cuts="";
 
  // ----------------------------------------------------- 
  // read zero ifar lag events
  // ----------------------------------------------------- 
 
  char sel[1024];
  char tmp[1024];
  cout << endl;
  vector<double> vIFAR;
  vector<double> vGPS;
  int nsel=0;
  for(int m=0;m<nwave;m++) {
    if(bin[m]<0) continue;    // simulation
    Cuts = CutBBH;
    // ------------------------> APPLY TRIALS FACTOR !!!
    sprintf(sel,"ifar/%d:time[0]",gOPT.trials);
    sprintf(tmp,"ifar>%f",pow(10,ifar_cmin)*gOPT.trials);
    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);
    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);
    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);
    if(gOPT.bbh) {
      Cuts  = TString::Format("%s ",tmp);
    } else {
      Cuts += TString::Format(" && %s ",tmp);
    }
    cout << "CutBBH : " << Cuts << endl;
    wave[m]->SetEstimate(wave[m]->GetEntries());
    wave[m]->Draw(sel,Cuts.Data(),"goff");
    nsel+= wave[m]->GetSelectedRows();
    double* V1 = wave[m]->GetV1();
    double* V2 = wave[m]->GetV2();
    for(int k=0;k<wave[m]->GetSelectedRows();k++) {vIFAR.push_back(V1[k]);vGPS.push_back(V2[k]);}
  }
 
  // IF gOPT.gw151226_ifar>0 -> ADD GW151226 Event !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // This is a special option introduced for the GW151226 event which has been found only into the extended segments
  if(gOPT.bbh && gOPT.gw151226_ifar>0) {
    vIFAR.push_back(gOPT.gw151226_ifar*YEAR);vGPS.push_back(GetGPSBBH("GW151226"));
  }
 
  // ----------------------------------------------------- 
  // sort ifar zero lag events
  // ----------------------------------------------------- 
 
  if(vIFAR.size()>0) {
    int size = vIFAR.size();
    int *index = new int[size];
    double *IFAR = new double[size];
    for(int k=0;k<size;k++) IFAR[k]=vIFAR[k];
    TMath::Sort(size,IFAR,index,true);
    if(gOPT.zl_style=="marker") {
      ZL_GPS.resize(size);
      ZL_IFAR.resize(size);
      ZL_NEVT.resize(size);
      for(int k=0;k<size;k++) {        // invert the ifar order
        ZL_GPS[size-k-1]  = vGPS[index[k]];  // event GPS time
        ZL_IFAR[size-k-1] = IFAR[index[k]]/YEAR;   // sorted ifar events
        ZL_NEVT[size-k-1] = k+1;    // number of events with ifar > ZL_IFAR[k]
      }
    } else {   // step
      ZL_GPS.resize(2*size-1);
      ZL_IFAR.resize(2*size-1);
      ZL_NEVT.resize(2*size-1);
      int n=0;
      for(int k=size-1;k>=0;k--) {     // invert the ifar order
        ZL_GPS[n]  = vGPS[index[k]];      // event GPS
        ZL_IFAR[n] = IFAR[index[k]]/YEAR; // sorted ifar events
        ZL_NEVT[n] = k+1;        // number of events with ifar > ZL_IFAR[k]
        n++;
        if(n>=ZL_NEVT.size()-1) break;
        // add point to make ZL STYLE STEPWISE
        ZL_GPS[n]  = ZL_GPS[n-1];      
        ZL_IFAR[n] = ZL_IFAR[n-1];     
        ZL_NEVT[n] = ZL_NEVT[n-1]-1;            
        n++;
      }
    }
    delete [] index; 
    delete [] IFAR; 
    // The foreground with IFAR>zl_max_ifar is set to zl_max_ifar 
    if(gOPT.zl_max_ifar>0) {
      for(int k=0;k<ZL_NEVT.size();k++) if(ZL_IFAR[k]>gOPT.zl_max_ifar) ZL_IFAR[k]=gOPT.zl_max_ifar;
    }
  }
 
  // adjust sup limit according to the loudest zero lag event
  double loudest_zl_ifar_sec = ZL_NEVT.size()>0 ? ZL_IFAR[ZL_NEVT.size()-1]*YEAR : 0;
  if(log10(loudest_zl_ifar_sec)>sup) sup=log10(loudest_zl_ifar_sec);
  if(inf>0 && sup/inf<2) sup=2*inf; 
  else          sup*=1.01;
 
  if(pow(10, inf)/YEAR>0.01) inf=log10(0.01*YEAR);
 
  if(gOPT.ifar_max>0) sup = log10(gOPT.ifar_max*YEAR);
 
  cout << endl;
  cout << "ifar inf : " << pow(10, inf)/YEAR << " years" << endl;
  cout << "ifar sup : " << pow(10, sup)/YEAR << " years" << endl;
 
  double difar_bin = (sup-inf)/IFAR_NBIN;
 
 
  // ======================================================================= 
  // SIMULATED EVENTS STUFF USED TO COMPUTE ASTROPHYSICAL RATES
  // ======================================================================= 
 
  // ----------------------------------------------------- 
  // read ifar simulated events
  // ----------------------------------------------------- 
 
  double ifar_sim_min = 0;
  double ifar_sim_max = 1.e20;
 
  while(!vIFAR.empty()) vIFAR.pop_back();
  vIFAR.clear(); std::vector<double>().swap(vIFAR);
  vector<int> vNINJ;    // store injections per chunk in each event (used for normalization)
 
  cout << endl;
  for(int m=0;m<nwave;m++) {
    if(bin[m]>=0) continue;   // select simulation root files
 
    int ninj = mdc[m]->GetEntries();
    cout << run[m] << "\t" << chunk[m] << "\tmdc injected entries : " << ninj << endl;
 
    wave[m]->SetEstimate(wave[m]->GetEntries());
    // ------------------------> APPLY TRIALS FACTOR !!!
    sprintf(sel,"ifar/%d",gOPT.trials);
    sprintf(tmp,"ifar>%f",pow(10,ifar_cmin)*gOPT.trials);
    sprintf(tmp,"ifar>%f && rho[%d]>%f",pow(10,ifar_cmin)*gOPT.trials,rho_id,gOPT.rho_min);
    Cuts  = TString::Format("%s ",tmp);
    wave[m]->Draw(sel,Cuts.Data(),"goff");
    double* V1 = wave[m]->GetV1();
    for(int k=0;k<wave[m]->GetSelectedRows();k++) {vIFAR.push_back(V1[k]);vNINJ.push_back(ninj);}
    double efficiency = ninj>0 ? (double)wave[m]->GetSelectedRows()/(double)ninj : 0;
    cout << run[m] << "\t" << chunk[m] << "\tselected wave sim detected entries : " 
         << wave[m]->GetSelectedRows() << "\tefficiency(%): " << int(100*efficiency) << "\tcumulative: " << vIFAR.size() << endl;
    if(vIFAR.size()==0) continue;
    int sel_entries = wave[m]->GetSelectedRows();
    value = TMath::MinElement(sel_entries,V1);
    //cout << "min sim IFAR per year : " << value/YEAR << endl;
    if(value>ifar_sim_min) ifar_sim_min=value;          // this condition is used to select the common sim IFAR inferior limit
    value = TMath::MaxElement(sel_entries,V1);
    //cout << "max sim IFAR per year : " << value/YEAR << endl;
    if(value<ifar_sim_max) ifar_sim_max=value;          // this condition is used to select the common sim IFAR superior limit
  } 
 
  if(vIFAR.size()>0) {
    ifar_sim_min/=YEAR;
    ifar_sim_max/=YEAR;
 
    cout << endl;
    cout << "ifar_sim_min : " << ifar_sim_min << " years" << endl;
    cout << "ifar_sim_max : " << ifar_sim_max << " years" << endl;
  }
 
  // ----------------------------------------------------- 
  // sort ifar simulated events
  // ----------------------------------------------------- 
 
  if(vIFAR.size()>0) {
    int size = vIFAR.size();
    int *index = new int[size];
    double *IFAR = new double[size];
    for(int k=0;k<size;k++) IFAR[k]=vIFAR[k];
    TMath::Sort(size,IFAR,index,true);
    SIM_IFAR.resize(2*size-1);
    SIM_NEVT.resize(2*size-1);
    SIM_BKG_NEVT.resize(2*size-1);
    int n=0;
//    vector<double> vWEIGHT(size);       // store weigths used for normalization
//    vWEIGHT[0]=1./vNINJ[0];
//    for(int k=1;k<size;k++) vWEIGHT[k]=vWEIGHT[k-1]+1./vNINJ[k-1]; // init weigths   
    for(int k=size-1;k>=0;k--) {    // invert the ifar order
      SIM_IFAR[n] = IFAR[index[k]]/YEAR;  // sorted ifar simulated events
#ifdef APPLY_MAX_SIM_IFAR 
      if(SIM_IFAR[n]<ifar_sim_min) SIM_IFAR[n]=ifar_sim_min;
      if(SIM_IFAR[n]>ifar_sim_max) SIM_IFAR[n]=ifar_sim_max;
#endif
//      SIM_NEVT[n] = vWEIGHT[k];         // number of simulated events with ifar > SIM_IFAR[k]
      SIM_NEVT[n] = k+1;         // number of simulated events with ifar > SIM_IFAR[k]
      n++;
      if(n>=SIM_NEVT.size()-1) break;
      // add point to make ZL STYLE STEPWISE
      SIM_IFAR[n] = SIM_IFAR[n-1];     
      SIM_NEVT[n] = SIM_NEVT[n-1]-1;            
//      SIM_NEVT[n] = vWEIGHT[k-1];          
      n++;
    }
    delete [] index; 
    delete [] IFAR; 
  }
  SIM_NEVT*=1./gOPT.sfactor;        // expected events from astrophysical reates
//  SIM_NEVT*=gOPT.sfactor;         // expected events from astrophysical reates
  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
 
#ifdef APPLY_MAX_SIM_IFAR 
  // set cut off of zero lag events according to the common sim IFAR inferior limit 
  if(SIM_NEVT.size()>0) {
    for(int i=0;i<ZL_NEVT.size();i++) if(ZL_IFAR[i]>ifar_sim_max) ZL_IFAR[i]=ifar_sim_max;
  }
#endif
 
  // ----------------------------------------------------- 
  // make sim plot 
  // ----------------------------------------------------- 
 
  grsim=NULL;
  grsimbkg=NULL;
  if(SIM_NEVT.size()>0) {
 
// (RATE_TAG)
    // compute the index which contains the maximum common IFAR
    int ifar_cmax_index=0;
    for(int i=0;i<SIM_IFAR.size();i++) if(SIM_IFAR.data[i]<pow(10,ifar_cmax)/YEAR) ifar_cmax_index=i;
 
// (RATE_TAG)
    grsim    = new TGraph(ifar_cmax_index,SIM_IFAR.data,SIM_NEVT.data);                  // expected events from astrophysical reates
    grsimbkg = new TGraph(ifar_cmax_index,SIM_IFAR.data,SIM_BKG_NEVT.data);              // expected events from astrophysical reates + expected from background
  }
 
  // ----------------------------------------------------- 
  // make sim belts at FAP0
  // ----------------------------------------------------- 
 
  if(SIM_NEVT.size()>0) { 
    int M=0;
    for(int n=0;n<IFAR_NBIN;n++) {
      double x = inf+n*difar_bin;
      if(x<(ifar_cmin-difar_bin) && x<inf) continue;  // RATE_TAG
      if(x>ifar_cmax) continue;           // RATE_TAG
      xSIM_BKG_IFAR[M] = pow(10,x)/YEAR;
      if((xSIM_BKG_IFAR[M]<SIM_IFAR[0]) || (xSIM_BKG_IFAR[M]>SIM_IFAR[SIM_IFAR.size()-1])) continue; 
      double mu = grsimbkg->Eval(xSIM_BKG_IFAR[M]);        // expected events in ZERO_LAG_TIME
      ySIM_BKG[M]=mu;
 
      double q[6];
      if(gOPT.belt_style=="continuos") GetPercentiles(mu, q);
 
      double nevt=0;
      nevt = gOPT.belt_style=="continuos" ? q[3] : GetPoissonNsup(mu,FAP0/2.); 
      eySIM_BKG_sup0[M] = (nevt-mu); 
      nevt = gOPT.belt_style=="continuos" ? q[2] : GetPoissonNinf(mu,FAP0/2.); 
      eySIM_BKG_inf0[M] = (mu-nevt); 
      M++;
    }
    xSIM_BKG_IFAR.resize(M);
    ySIM_BKG.resize(M);
    eySIM_BKG_inf0.resize(M);
    eySIM_BKG_sup0.resize(M);
  }
 
 
  // ======================================================================= 
  // MAKE BELTS FOR EXPECTED BACKGROUND
  // ======================================================================= 
 
  // ----------------------------------------------------- 
  // make belts at FAP0,FAP1,FAP2
  // ----------------------------------------------------- 
 
  if(gOPT.belt_style=="step") {
    // add points to make BELT STYLE STEPWISE
 
    xIFAR.resize(2*IFAR_NBIN);
    xFAR.resize(2*IFAR_NBIN);
    bRATE.resize(2*IFAR_NBIN);
    yBKG.resize(2*IFAR_NBIN);
 
    exFAR.resize(2*IFAR_NBIN); exFAR=0;
    eyRATEinf0.resize(2*IFAR_NBIN);
    eyRATEsup0.resize(2*IFAR_NBIN);
    eyRATEinf1.resize(2*IFAR_NBIN);
    eyRATEsup1.resize(2*IFAR_NBIN);
    eyRATEinf2.resize(2*IFAR_NBIN);
    eyRATEsup2.resize(2*IFAR_NBIN);
  }
 
  cout << "Make_PP_IFAR.C - generate belts : be patient, it takes a while ..." << endl;
  int N=0;
  for(int n=0;n<IFAR_NBIN;n++) {
    double x = inf+n*difar_bin;
    //if(x<(ifar_cmin-difar_bin) || x>(ifar_max+2*difar_bin)) continue;  // RATE_TAG
    if(x<(ifar_cmin-difar_bin) && x<inf) continue;  // RATE_TAG
    xIFAR[N] = pow(10,x);
    xFAR[N]  = 1./xIFAR[N];
    bRATE[N] = xFAR[N];
    double mu = xFAR[N]*obsTime;        // expected bakground events in ZERO_LAG_TIME
 
    double q[6];
    if(gOPT.belt_style=="continuos") GetPercentiles(mu, q);
 
    yBKG[N] = gOPT.belt_style=="continuos" ? q[3] : GetPoissonNsup(mu,FAP0/2.); 
    eyRATEsup0[N] = (yBKG[N]-mu)/obsTime; 
    yBKG[N] = gOPT.belt_style=="continuos" ? q[2] : GetPoissonNinf(mu,FAP0/2.); 
    eyRATEinf0[N] = (mu-yBKG[N])/obsTime; 
    yBKG[N] = gOPT.belt_style=="continuos" ? q[4] : GetPoissonNsup(mu,FAP1/2.); 
    eyRATEsup1[N] = (yBKG[N]-mu)/obsTime; 
    yBKG[N] = gOPT.belt_style=="continuos" ? q[1] : GetPoissonNinf(mu,FAP1/2.); 
    eyRATEinf1[N] = (mu-yBKG[N])/obsTime; 
    yBKG[N] = gOPT.belt_style=="continuos" ? q[5] : GetPoissonNsup(mu,FAP2/2.); 
    eyRATEsup2[N] = (yBKG[N]-mu)/obsTime; 
    yBKG[N] = gOPT.belt_style=="continuos" ? q[0] : GetPoissonNinf(mu,FAP2/2.); 
    eyRATEinf2[N] = (mu-yBKG[N])/obsTime; 
 
    N++;
 
    if(gOPT.belt_style=="step") {   // add points to make BELT STYLE STEPWISE
 
      if(N>1) {
 
        xIFAR[N]   = xIFAR[N-1];
        xFAR[N]    = xFAR[N-1];
        bRATE[N]   = bRATE[N-1];
 
        xIFAR[N-1] = xIFAR[N-2];
        xFAR[N-1]  = xFAR[N-2];
        bRATE[N-1] = bRATE[N-2];
 
        yBKG[N-1]  = mu;
        yBKG[N]    = mu;
 
        eyRATEsup0[N] = eyRATEsup0[N-1];
        eyRATEinf0[N] = eyRATEinf0[N-1];
        eyRATEsup1[N] = eyRATEsup1[N-1];
        eyRATEinf1[N] = eyRATEinf1[N-1];
        eyRATEsup2[N] = eyRATEsup2[N-1];
        eyRATEinf2[N] = eyRATEinf2[N-1];
 
        N++;
 
      } else {
        yBKG[N-1]=mu;
      }
    } else {
      yBKG[N-1]=bRATE[N-1]*obsTime;
    }
  }
  xIFAR.resize(N);
  xFAR.resize(N);
 
  // normalize -> per SEC -> per YEAR
  xFAR *=YEAR;  
  xIFAR*=1./YEAR;
 
  // compute number of events in the obsTime period
  bRATE*=obsTime;
  eyRATEsup0*=obsTime;
  eyRATEinf0*=obsTime;
  eyRATEsup1*=obsTime;
  eyRATEinf1*=obsTime;
  eyRATEsup2*=obsTime;
  eyRATEinf2*=obsTime;
 
 
  // ----------------------------------------------------- 
  // print sigma belts
  // ----------------------------------------------------- 
 
  double sigma[3];
  sigma[0] = sqrt(2)*TMath::ErfcInverse(FAP0);
  sigma[1] = sqrt(2)*TMath::ErfcInverse(FAP1);
  sigma[2] = sqrt(2)*TMath::ErfcInverse(FAP2);
 
  cout << endl << "--------------------------------------------------------" << endl << endl;;
  cout << "FAP = " << FAP0 << "\t -> SIGMA = " << sigma[0] << endl;
  cout << "FAP = " << FAP1 << "\t -> SIGMA = " << sigma[1] << endl;
  cout << "FAP = " << FAP2 << "\t -> SIGMA = " << sigma[2] << endl;
  cout << endl << "--------------------------------------------------------" << endl << endl;;
 
 
  // ======================================================================= 
  // DUMP PLOTS AND ASCII FILES
  // ======================================================================= 
 
  // ----------------------------------------------------- 
  // Make Plot 
  // ----------------------------------------------------- 
 
  char title[1024];
  if(gOPT.bbh) {
    sprintf(title,"cWB %s (run=%s : lag=%d : slag=%d : included BBH)",gOPT.search.Data(),gOPT.run.Data(),gOPT.lag,gOPT.slag);
  } else {
    sprintf(title,"cWB %s (run=%s : lag=%d : slag=%d : excluded BBH)",gOPT.search.Data(),gOPT.run.Data(),gOPT.lag,gOPT.slag);
  }
 
  sprintf(title,"");
 
  MakePlot(title, obsTime, ifar_cmax);
 
  // ----------------------------------------------------- 
  // Dump Zero Lag Events to ascii file 
  // ----------------------------------------------------- 
 
  DumpPeriod(liveTot);
  DumpLoudest(obsTime);
  DumpBackground();
  DumpForeground();
  DumpSimulation(ifar_cmax);
 
  if(gOPT.exit && gOPT.pfname!="") exit(0);
}
 
void MakePlot(TString title, double obsTime, double ifar_cmax) {
 
  canvas = new TCanvas("FAR", "FAR", 300,40, 800, 600);
  canvas->SetGridx();
  canvas->SetGridy();
  canvas->SetLogx();
  canvas->SetLogy();
//  canvas->SetFrameBorderMode(0);
//  canvas->SetFrameLineColor(17);
//  canvas->SetFrameBorderMode(0);
 
  gr0 = new TGraph(xIFAR.size(),xIFAR.data,bRATE.data);  
  gr0->SetTitle(title);
  gr0->Draw("LA3");
  gr0->SetLineWidth(2);
  gr0->SetLineStyle(1);
  gr0->SetLineColor(kGray+2);
  gr0->GetHistogram()->GetXaxis()->CenterTitle();
  gr0->GetHistogram()->GetYaxis()->CenterTitle();
  gr0->GetHistogram()->GetXaxis()->SetLabelOffset(0.0001);
  gr0->GetHistogram()->GetXaxis()->SetTickLength(0.02);
  gr0->GetHistogram()->GetXaxis()->SetLabelSize(0.035);
  gr0->GetHistogram()->GetYaxis()->SetLabelSize(0.035);
  gr0->GetHistogram()->GetXaxis()->SetTitleSize(0.035);
  gr0->GetHistogram()->GetYaxis()->SetTitleSize(0.035);
  gr0->GetHistogram()->GetXaxis()->SetTitleOffset(1.7);
  gr0->GetHistogram()->GetYaxis()->SetTitleOffset(1.5);
  gr0->GetHistogram()->GetXaxis()->SetTitleOffset(1.4);
  gr0->GetHistogram()->GetYaxis()->SetTitleOffset(1.2);
  gr0->GetHistogram()->GetXaxis()->SetAxisColor(17);
  gr0->GetHistogram()->GetYaxis()->SetAxisColor(17);
  canvas->RedrawAxis();
  canvas->Update();
 
  // set axis ranges
  if(xIFAR.size()>0) {
    double ymin = 0.7;
    double ymax = gr0->GetHistogram()->GetMaximum();
    if(ZL_NEVT.size()>0) ymax=2*ZL_NEVT[0];
    if(ymax/ymin<10) ymax=10*ymin;
    //if(ymax>gr0->GetHistogram()->GetMaximum()) ymax=gr0->GetHistogram()->GetMaximum();
    if(ymax<1) ymax=2;
    gr0->GetHistogram()->GetYaxis()->SetRangeUser(ymin,ymax);
    gr0->GetHistogram()->GetXaxis()->SetRangeUser(0.01,1e5);
    canvas->RedrawAxis();
    canvas->Update();
    if((xIFAR[0]>0) && (xIFAR[xIFAR.size()-1]/xIFAR[0])<10) {
      canvas->SetLogx(false);
      canvas->SetLogy(false);
    }
  } else {
    canvas->SetLogx(false);
    canvas->SetLogy(false);
  }
  gr0->GetHistogram()->GetXaxis()->SetTitle("Inverse False Alarm Rate (yr)");
  gr0->GetHistogram()->GetYaxis()->SetTitle("Cumulative Number of Events");
  gr0->GetHistogram()->GetXaxis()->SetLabelFont(132);
  gr0->GetHistogram()->GetYaxis()->SetLabelFont(132);
  gr0->GetHistogram()->GetXaxis()->SetTitleFont(132);
  gr0->GetHistogram()->GetYaxis()->SetTitleFont(132);
  gr0->GetHistogram()->GetXaxis()->SetTickLength(0);
 
  // Plot FAP2 
  egr2 = new TGraphAsymmErrors(xIFAR.size(),xIFAR.data,yBKG.data,exFAR.data,exFAR.data,eyRATEinf2.data,eyRATEsup2.data);
  egr2->SetMarkerStyle(20);
  egr2->SetMarkerSize(0);
  egr2->SetLineColor(15);
  egr2->SetFillStyle(1001);
  egr2->SetFillColor(18);  // light gray
  egr2->Draw("3SAME");
 
  // Plot FAP1
  egr1 = new TGraphAsymmErrors(xIFAR.size(),xIFAR.data,yBKG.data,exFAR.data,exFAR.data,eyRATEinf1.data,eyRATEsup1.data);
  egr1->SetMarkerStyle(20);
  egr1->SetMarkerSize(0);
  egr1->SetLineColor(15);
  egr1->SetFillStyle(1001);
  egr1->SetFillColor(17);  // gray
  egr1->Draw("3SAME");
 
  // Plot FAP0
  egr0 = new TGraphAsymmErrors(xIFAR.size(),xIFAR.data,yBKG.data,exFAR.data,exFAR.data,eyRATEinf0.data,eyRATEsup0.data);
  egr0->SetMarkerStyle(20);
  egr0->SetMarkerSize(0);
  egr0->SetLineColor(15);
  egr0->SetFillStyle(1001);
  egr0->SetFillColor(16);  // dark gray
  egr0->Draw("3SAME");
 
// BEG TEST P-ASTRO
  // Plot Simulated
  if(grsimbkg!=NULL) {
 
    // compute the index which contains the maximum common IFAR (RATE_TAG)
    int ifar_cmax_index=0;
    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;
 
    // Plot sim FAP0
    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)
    egrsimbkg0->SetMarkerSize(0);
    egrsimbkg0->SetLineColor(15);
    egrsimbkg0->SetFillColorAlpha(kGreen, 0.1);
    egrsimbkg0->Draw("3SAME");   // belt
 
    grsim->SetLineColor(kGreen);
    grsim->SetLineWidth(1);
    grsim->SetMarkerSize(0);
    grsim->Draw("LSAME");  // expected signal from astrophysical rates
 
    grsimbkg->SetLineColor(kGreen+2);
    grsimbkg->SetLineWidth(1);
    grsimbkg->SetMarkerSize(0);
    grsimbkg->Draw("LSAME");  // expected signal from astrophysical rates + expected from backgrouns
  }
// END TEST P-ASTRO
 
  // Plot Foreground
  if(ZL_NEVT.size()>0) {
    gr = new TGraph(ZL_NEVT.size(),ZL_IFAR.data,ZL_NEVT.data);
//gr->SetLineStyle(2);
//gr->SetLineWidth(2);
    gr->SetLineColor(kRed);
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(23);
    if(gOPT.zl_style=="marker") {
      gr->SetMarkerSize(0.4);
      gr->Draw("PSAME");   // marker
    } else {
      gr->SetMarkerSize(0);
      gr->Draw("SAME");    // stepwise
    }
  } else {
    gr = NULL;
    //canvas->SetLogy(false);
  }
 
#ifdef PLOT_FOREGROUND_NOBBH
 
  wavearray<double> ZL_IFAR_NOBBH(ZL_NEVT.size());    // sorted ifar events
 
  // select event without known BBH
  int nobbh_size=0;
  for(int i=0;i<ZL_NEVT.size();i++) {
    if((GetNameBBH(ZL_GPS[i])=="")) ZL_IFAR_NOBBH[nobbh_size++]=ZL_IFAR[i];
  }
  if(ZL_IFAR_NOBBH[nobbh_size-1]==ZL_IFAR_NOBBH[nobbh_size-2]) nobbh_size-=1;
  ZL_IFAR_NOBBH.resize(nobbh_size);
  wavearray<double> ZL_NEVT_NOBBH(nobbh_size);  // number of nobbh events with ifar > ZL_IFAR[k]
  for(int i=0;i<nobbh_size;i+=2) {
    ZL_NEVT_NOBBH[i]=(nobbh_size-i)/2+1;
    ZL_NEVT_NOBBH[i+1]=(nobbh_size-i)/2;
  }
  ZL_NEVT_NOBBH[0]=(nobbh_size)/2+1;
 
  // Plot Foreground no bbh
  TGraph* gr_nobbh;
  if(ZL_NEVT_NOBBH.size()>0) {
    gr_nobbh = new TGraph(ZL_NEVT_NOBBH.size(),ZL_IFAR_NOBBH.data,ZL_NEVT_NOBBH.data);
    //gr_nobbh->SetLineColor(kRed);
    //gr_nobbh->SetLineColor(kGreen+2);
    //gr_nobbh->SetLineColor(kOrange+7);
    gr_nobbh->SetLineColor(kBlack);
    gr_nobbh->SetLineWidth(1);
    gr_nobbh->SetLineStyle(3);
    gr_nobbh->SetMarkerColor(kRed);
    gr_nobbh->SetMarkerStyle(23);
    if(gOPT.zl_style=="marker") {
      gr_nobbh->SetMarkerSize(0.4);
      gr_nobbh->Draw("PSAME");   // marker
    } else {
      gr_nobbh->SetMarkerSize(0);
      gr_nobbh->Draw("SAME"); // stepwise
    }
  } 
 
#endif
 
  gr0->Draw("LSAME");
 
  // Plot blue marker & name for known BBH
  bool anyBBH=false;
  for(int i=0;i<ZL_NEVT.size();i++) {
    if((GetNameBBH(ZL_GPS[i])!="")&&(i%2==0)) {
      marker = new TMarker(ZL_IFAR[i],ZL_NEVT[i], 20);
      bool selected=false;
      if((GetNameBBH(ZL_GPS[i])=="GW150914")) selected=true;
      if((GetNameBBH(ZL_GPS[i])=="GW170814")) selected=true;
      if((GetNameBBH(ZL_GPS[i])=="GW170104")) selected=true;
      if((GetNameBBH(ZL_GPS[i])=="GW170608")) selected=true;
      if(selected) {
        arrow = new TArrow(ZL_IFAR[i],ZL_NEVT[i],ZL_IFAR[i]*1.5,ZL_NEVT[i],0.05,"|>");
        arrow->SetAngle(60);
        arrow->SetLineWidth(2);
        arrow->SetLineColor(kBlue);
        arrow->SetFillColor(kBlue);
        arrow->SetArrowSize(0.012);
        arrow->Draw();
      }
      marker->SetMarkerSize(1.2);
      marker->SetMarkerColor(kBlue);
      marker->Draw();
 
      latex = new TLatex(ZL_IFAR[i]*(4./3.),ZL_NEVT[i]*1.1, GetNameBBH(ZL_GPS[i]).Data());
      latex->SetTextFont(132);
      latex->SetTextSize(0.028);
      latex->SetTextColor(kBlack);
      latex->Draw();
 
      anyBBH=true;
    }
  }
  // if there no BBH events then a red marker is draw for loudest event
  if((!anyBBH)&&(ZL_NEVT.size()>0)&&(GetNameBBH(ZL_GPS[ZL_NEVT.size()-1])=="")) { 
    marker = new TMarker(ZL_IFAR[ZL_NEVT.size()-1],ZL_NEVT[ZL_NEVT.size()-1], 23);
    marker->SetMarkerSize(1.2);
    marker->SetMarkerColor(kRed);
    marker->Draw();
  }
 
  // draw legend
#ifdef PLOT_FOREGROUND_NOBBH
  leg = new TLegend(0.6,0.6,0.90,0.90,NULL,"brNDC");
#else
  if(grsimbkg!=NULL) {
    leg = new TLegend(0.6,0.630,0.90,0.90,NULL,"brNDC");
  } else {
    leg = new TLegend(0.6,0.6,0.90,0.90,NULL,"brNDC");
  }
#endif
 
  leg->SetBorderSize(1);
  leg->SetTextAlign(22);
  leg->SetTextFont(132);
  leg->SetLineColor(1);
  leg->SetLineStyle(1);
  leg->SetLineWidth(1);
  leg->SetFillColor(0);
  leg->SetFillStyle(1001);
  if(grsimbkg!=NULL) leg->SetTextSize(0.035); else leg->SetTextSize(0.035);
  leg->SetLineColor(kBlack);
  leg->SetFillColor(kWhite);
 
  char legLabel[64];
  sprintf(legLabel,"Foreground"); if(gr!=NULL) leg->AddEntry(gr,legLabel,"lp");
#ifdef PLOT_FOREGROUND_NOBBH
  sprintf(legLabel,"Residual Foreground"); if(gr_nobbh!=NULL) leg->AddEntry(gr_nobbh,legLabel,"lp");
#endif
  sprintf(legLabel,"Expected Background"); leg->AddEntry(gr0,legLabel,"lp");
  if(gOPT.bbh && anyBBH)    leg->AddEntry(marker,"BBH Events","p");
  else if(ZL_NEVT.size()>0) leg->AddEntry(marker,"Loudest","p");
  if(grsimbkg!=NULL) {
    leg->AddEntry(grsim,"Signal Model","lp");;
    leg->AddEntry(grsimbkg,"Noise+Signal Model","lp");;
  } else {
    // add background belts labels
    leg->AddEntry(egr2,"< 3 #sigma","f");
    leg->AddEntry(egr1,"< 2 #sigma","f");
    leg->AddEntry(egr0,"< 1 #sigma","f");
  }
  leg->Draw();
 
  // fix frame lines color
//  double xmax = gr0->GetHistogram()->GetXaxis()->GetXmax();
//  double xmin = gr0->GetHistogram()->GetXaxis()->GetXmin();
  double xmax = pow(10,gPad->GetUxmax());
  double xmin = pow(10,gPad->GetUxmin());
  //double ymin = pow(10,gPad->GetUymin());
  double ymin = 0.0;
  double ymax = pow(10,gPad->GetUymax());
  cout << "xmin=" << xmin << " xmax=" << xmax << endl;
  cout << "ymin=" << ymin << " ymax=" << ymax << endl;
 
  TLine* frameLine[4];
  frameLine[0] = new TLine(xmin,ymin,xmax,ymin);
  frameLine[1] = new TLine(xmin,ymin,xmin,ymax);
  frameLine[2] = new TLine(xmax,ymin,xmax,ymax);
  frameLine[3] = new TLine(0.,ymax,xmax,ymax);
  for(int i=0;i<4;i++) frameLine[i]->Draw();
 
 
  // dump plot to disk
  if(gOPT.pfname) {
    char ofname[1024];
    TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
    if(gOPT.bbh) sprintf(ofname,"%s_bbh_plot.png",PFNAME.Data());
    else         sprintf(ofname,"%s_nobbh_plot.png",PFNAME.Data());
    canvas->Print(ofname);
    if(gOPT.bbh) sprintf(ofname,"%s_bbh_plot.root",PFNAME.Data());
    else         sprintf(ofname,"%s_nobbh_plot.root",PFNAME.Data());
    canvas->Print(ofname);
    if(gOPT.bbh) sprintf(ofname,"%s_bbh_plot.pdf",PFNAME.Data());
    else         sprintf(ofname,"%s_nobbh_plot.pdf",PFNAME.Data());
    canvas->Print(ofname);
  }
}
 
int GetPercentiles(double mu, double* q) {
 
  double max = mu>10 ? 10*mu : 100;
 
  TF1 f("poisson","TMath::Poisson(x,[1])",0,max);
  f.SetParameter(1,mu);
  f.SetNpx(1000);
 
  double probSum[6] = {FAP2/2,FAP1/2,FAP0/2,1-FAP0/2,1-FAP1/2,1-FAP2/2};
  return f.GetQuantiles(6,q,probSum);
}
 
int GetPoissonNsup(double mu, double fap) {
 
  double FAP = 1.;
  for(int j=0;j<1000000;j++) {
    //FAP -= pow(mu,j)*exp(-mu)/TMath::Factorial(j);
    FAP -= TMath::PoissonI(j,mu);
    //cout << " N : " << j << " MU : " << mu << " FAP " << FAP << " fap : " << fap << endl;
    if(FAP<fap) return j==0 ? 0 : j+1;
  }
  return 0;
}
 
int GetPoissonNinf(double mu, double fap) {
 
  double FAP = 0.;
  for(int j=0;j<1000000;j++) {
    //FAP -= pow(mu,j)*exp(-mu)/TMath::Factorial(j);
    FAP += TMath::PoissonI(j,mu);
    //cout << " N : " << j << " MU : " << mu << " FAP " << FAP << " fap : " << fap << endl;
    if(FAP>fap) return j==0 ? 0 : j-1;
  }
  return 0;
}
 
int ReadFileList(TString ifName, TChain** live, TChain** wave, TChain** mdc, int* lag, int* slag, int* chunk, int* bin, TString* run) {
 
  // Read WAVE/LIVE file list
 
  ifstream in;
  in.open(ifName.Data(),ios::in);
  if (!in.good()) {cout << "ReadFileList Error Opening File : " << ifName.Data() << endl;exit(1);}
 
  int size=0;
  char str[1024];
  int fpos=0;
  while(true) {
    in.getline(str,1024);
    if (!in.good()) break;
    if(str[0] != '#') size++;
  }
  in.clear(ios::goodbit);
  in.seekg(0, ios::beg);
  if (size==0) {cout << "ReadFileList Error : File " << ifName.Data() << " is empty" << endl;exit(1);}
  if (size>MAX_LIST_SIZE) {cout << "ReadFileList Error : Files in " << ifName.Data() << " > " << MAX_LIST_SIZE << endl;exit(1);}
 
  char sfile[1024];
  int xlag,xslag,ichunk,ibin;
  char srun[256];
  int nwave=0;
  while(true) {
    fpos=in.tellg();
    in.getline(str,1024);
    if (!in.good()) break;
    if(str[0] == '#' || str[0]=='\0') continue;
    in.seekg(fpos, ios::beg);
    in >> sfile >> xlag >> xslag >> ichunk >> ibin >> srun;
    if(!in.good()) break;
    lag[nwave]=xlag;
    slag[nwave]=xslag;
    TString fwave = sfile;
    if(!CWB::Toolbox::isFileExisting(fwave)) {cout << "ReadFileList Error : File not exist : " << fwave << endl;exit(1);} 
    wave[nwave]->Add(fwave);
    wave[nwave]->SetTitle(fwave);
    if(ibin>=0) { // ibin<=0/ibin>=0 -> simulation/background
      TString flive = sfile;
      flive.ReplaceAll("wave_","live_");
      if(!CWB::Toolbox::isFileExisting(flive)) {cout << "ReadFileList Error : File not exist : " << flive << endl;exit(1);} 
      live[nwave]->Add(flive);
      live[nwave]->SetTitle(flive);
    } else {      // simulation -> read mdc
      TString fmdc = sfile;
      fmdc.ReplaceAll("wave_","mdc_");
      if(!CWB::Toolbox::isFileExisting(fmdc)) {cout << "ReadFileList Error : File not exist : " << fmdc << endl;exit(1);} 
      mdc[nwave]->Add(fmdc);
      mdc[nwave]->SetTitle(fmdc);
    }
    chunk[nwave]=ichunk;
    bin[nwave]=ibin;
    run[nwave]=srun;
    cout << nwave+1 << "\t" << fwave << "\t" << xlag << "\t" << xslag << "\t" << ichunk << "\t" << ibin << "\t" << srun << endl;
    nwave++;
  }
  in.close();
 
  return nwave;
}
 
void GetLiveTime(int nwave, TChain** live, int* lag, int* slag, int* bin, TString* run, double* liveZero, double* liveTot) {
 
  CWB::Toolbox TB;
 
  wavearray<double> Trun(500000); Trun = 0.;
  wavearray<double> Wlag[NIFO_MAX+1];
  wavearray<double> Wslag[NIFO_MAX+1];
  wavearray<double> Tdlag;
  wavearray<double> Tlag;
  wavearray<double> Rlag;
  wavearray<double> Elag;
  wavearray<double> Olag;
 
  cout << "Start CWB::Toolbox::getLiveTime : be patient, it takes a while ..." << endl;
  for(int m=0;m<nwave;m++) {
    if(bin[m]<0) continue;    // simulation
    cout << "Compute liveTot : Read : " << live[m]->GetTitle() << endl;
    if(run[m]=="O1") liveTot[0]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,lag[m],slag[m]);
    if(run[m]=="O2") liveTot[1]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,lag[m],slag[m]);
    if(run[m]=="O3") liveTot[2]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,lag[m],slag[m]);
  }
 
  // get live time zero (when is not included into Tlag array)
  for(int m=0;m<nwave;m++) {
    if(bin[m]<0) continue;    // simulation
    if((lag[m]!=0)&&(slag[m]!=0)) {
      cout << "Compute liveZero : Read : " << live[m]->GetTitle() << endl;
      if(run[m]=="O1") liveZero[0]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,0,0);
      if(run[m]=="O2") liveZero[1]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,0,0);
      if(run[m]=="O3") liveZero[2]+=TB.getLiveTime(gOPT.nifo,*live[m],Trun,Wlag,Wslag,Tlag,Tdlag,0,0);
      //liveZero+=TB.getZeroLiveTime(gOPT.nifo,*live[m]);   // crash in ATLAS ???
      if(run[m]=="O1") printf("live time: zero lags = %10.1f, %10.1f days \n",liveZero[0],liveZero[0]/DAY);
      if(run[m]=="O2") printf("live time: zero lags = %10.1f, %10.1f days \n",liveZero[1],liveZero[1]/DAY);
      if(run[m]=="O3") printf("live time: zero lags = %10.1f, %10.1f days \n",liveZero[2],liveZero[1]/DAY);
    }
  }
}
 
void InitBBH() {
 
// ---------------------------------------------------------------------------------
// INIT KNOWN O1,O2,O3 BBH
// ---------------------------------------------------------------------------------
 
  // O1 known BBH events
  O1_BBH_NAME.resize(3);
  O1_BBH_GPS.resize(3);
 
  O1_BBH_NAME[0]="GW150914";  O1_BBH_GPS[0]=1126259462.40;
  O1_BBH_NAME[1]="LVT151012"; O1_BBH_GPS[1]=1128678900.55;
  O1_BBH_NAME[2]="GW151226";  O1_BBH_GPS[2]=1135136350.56;
 
  // O2 known BBH events
  O2_BBH_NAME.resize(7);
  O2_BBH_GPS.resize(7);
 
  O2_BBH_NAME[0]="GW170104";  O2_BBH_GPS[0]=1167559936.59;
  O2_BBH_NAME[1]="GW170608";  O2_BBH_GPS[1]=1180922494.38;
  O2_BBH_NAME[2]="GW170729";  O2_BBH_GPS[2]=1185389807.33;
  O2_BBH_NAME[3]="GW170809";  O2_BBH_GPS[3]=1186302519.73;
  O2_BBH_NAME[4]="GW170814";  O2_BBH_GPS[4]=1186741861.50;
  O2_BBH_NAME[5]="GW170818";  O2_BBH_GPS[5]=1187058327.07;
  O2_BBH_NAME[6]="GW170823";  O2_BBH_GPS[6]=1187529256.51;
 
  // O3 known BBH events
  ReadPublicAlerts(); 
 
}
 
TString GetCutBBH() {
 
  TString CutBBH="";
 
  if(gOPT.run.Contains("O1")) {
    CutBBH += TString::Format("   !(abs(time[0]-%f)<0.2) ",O1_BBH_GPS[0]);
  }
  if(gOPT.run.Contains("O1")) {
    for(int i=1;i<O1_BBH_GPS.size();i++) CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.2) ",O1_BBH_GPS[i]);
  }
 
  if(gOPT.run=="O1O2") {
    CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.2) ",O2_BBH_GPS[0]);
  } else if(gOPT.run=="O2") {
    CutBBH += TString::Format("   !(abs(time[0]-%f)<0.2) ",O2_BBH_GPS[0]);
  }
  if(gOPT.run.Contains("O2")) {
    for(int i=1;i<O2_BBH_GPS.size();i++) CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.2) ",O2_BBH_GPS[i]);
  }
 
  if(gOPT.run.Contains("O3")) {
    CutBBH += TString::Format("   !(abs(time[0]-%f)<0.2) ",O3_BBH_GPS[0]);
  }
  if(gOPT.run.Contains("O3")) {
    for(int i=1;i<O3_BBH_GPS.size();i++) CutBBH += TString::Format("&& !(abs(time[0]-%f)<0.5) ",O3_BBH_GPS[i]);
  }
 
  return CutBBH;
}
 
TString GetNameBBH(double gps) {
 
  // known O1 BBH
  for(int i=0;i<O1_BBH_GPS.size();i++) if(fabs(gps-O1_BBH_GPS[i])<0.2) return O1_BBH_NAME[i]; 
 
  // known O2 BBH
  for(int i=0;i<O2_BBH_GPS.size();i++) if(fabs(gps-O2_BBH_GPS[i])<0.2) return O2_BBH_NAME[i]; 
 
  // known O3 BBH
  for(int i=0;i<O3_BBH_GPS.size();i++) if(fabs(gps-O3_BBH_GPS[i])<0.5) return O3_BBH_NAME[i]; 
 
  return "";
}
 
double GetGPSBBH(TString name) {
 
  // known O1 BBH
  for(int i=0;i<O1_BBH_NAME.size();i++) if(O1_BBH_NAME[i]==name) return O1_BBH_GPS[i]; 
 
  // known O2 BBH
  for(int i=0;i<O2_BBH_NAME.size();i++) if(O2_BBH_NAME[i]==name) return O2_BBH_GPS[i]; 
 
  // known O3 BBH
  for(int i=0;i<O3_BBH_NAME.size();i++) if(O3_BBH_NAME[i]==name) return O3_BBH_GPS[i]; 
 
  return 0;
}
 
void ResetUserOptions() {
 
    gOPT.run            = PP_RUN;
    gOPT.search         = PP_SEARCH;
 
    gOPT.nifo           = PP_NIFO;
    gOPT.lag         = PP_LAG;
    gOPT.slag        = PP_SLAG;
    gOPT.chunk       = PP_CHUNK;
 
    gOPT.trials         = PP_TRIALS;
 
    gOPT.pfname         = PP_PFNAME;
    gOPT.zl_style       = PP_ZL_STYLE;
    gOPT.belt_style     = PP_BELT_STYLE;
 
    gOPT.ifar_max       = PP_IFAR_MAX;
    gOPT.rho_min        = PP_RHO_MIN;
 
    gOPT.sfactor        = PP_SFACTOR;
 
    gOPT.bbh         = PP_BBH;
    gOPT.exit        = PP_EXIT;
 
    gOPT.gw151226_ifar  = PP_GW151226_IFAR;
    gOPT.zl_max_ifar    = PP_ZL_MAX_IFAR;
}
 
void PrintUserOptions() {
 
    cout.precision(14);
 
    TString sexit = gOPT.exit?"true":"false";
    TString bbh   = gOPT.bbh?"true":"false";
 
    cout << endl;
 
    cout << "-----------------------------------------"     << endl;
    cout << "PP IFAR config options                   "     << endl;
    cout << "-----------------------------------------"     << endl << endl;
 
    cout << "PP_RUN              " << gOPT.run              << endl;
    cout << "PP_SEARCH           " << gOPT.search           << endl;
 
    cout << "PP_NIFO             " << gOPT.nifo             << endl;
    cout << "PP_LAG              " << gOPT.lag              << endl;
    cout << "PP_SLAG             " << gOPT.slag             << endl;
    cout << "PP_CHUNK            " << gOPT.chunk            << endl;
 
    cout << "PP_TRIALS           " << gOPT.trials           << endl;
 
    cout << "PP_PFNAME           " << gOPT.pfname           << endl;
    cout << "PP_ZL_STYLE         " << gOPT.zl_style         << endl;
    cout << "PP_BELT_STYLE       " << gOPT.belt_style       << endl;
 
    cout << "PP_IFAR_MAX         " << gOPT.ifar_max << " yr"<< endl;
    cout << "PP_RHO_MIN          " << gOPT.rho_min          << endl;
 
    cout << "PP_SFACTOR          " << gOPT.sfactor          << endl;
 
    cout << "PP_BBH              " << bbh                   << endl;
    cout << "PP_EXIT             " << sexit                 << endl;
 
    cout << "PP_GW151226_IFAR    " << gOPT.gw151226_ifar    << endl;
    cout << "PP_ZL_MAX_IFAR      " << gOPT.zl_max_ifar      << endl;
 
    cout << endl;
}
 
void ReadUserOptions(TString options) {
 
  // get plugin options 
 
  if(TString(options)!="") {
 
    //cout << "PP_IFAR options : " << options << endl;
    TObjArray* token = TString(options).Tokenize(TString(' '));
    for(int j=0;j<token->GetEntries();j++) {
 
      TObjString* tok = (TObjString*)token->At(j);
      TString stok = tok->GetString();
 
      if(stok.Contains("run=")) {
        TString pp_run=stok;
        pp_run.Remove(0,pp_run.Last('=')+1);
        gOPT.run=pp_run;
      }
 
      if(stok.Contains("search=")) {
        TString pp_search=stok;
        pp_search.Remove(0,pp_search.Last('=')+1);
        gOPT.search=pp_search;
      }
 
      if(stok.Contains("nifo=")) {
        TString pp_nifo=stok;
        pp_nifo.Remove(0,pp_nifo.Last('=')+1);
        if(pp_nifo.IsDigit()) gOPT.nifo=pp_nifo.Atoi();
      }
 
      if(stok.Contains("lag=") && !stok.Contains("slag=")) {
        TString pp_lag=stok;
        pp_lag.Remove(0,pp_lag.Last('=')+1);
        if(pp_lag.IsDigit()) gOPT.lag=pp_lag.Atoi();
      }
 
      if(stok.Contains("slag=")) {
        TString pp_slag=stok;
        pp_slag.Remove(0,pp_slag.Last('=')+1);
        if(pp_slag.IsDigit()) gOPT.slag=pp_slag.Atoi();
      }
 
      if(stok.Contains("chunk=")) {
        TString pp_chunk=stok;
        pp_chunk.Remove(0,pp_chunk.Last('=')+1);
        if(pp_chunk.IsDigit()) gOPT.chunk=pp_chunk.Atoi();
      }
 
      if(stok.Contains("trials=")) {
        TString pp_trials=stok;
        pp_trials.Remove(0,pp_trials.Last('=')+1);
        if(pp_trials.IsDigit()) gOPT.trials=pp_trials.Atoi();
      }
 
      if(stok.Contains("pfname=")) {
        TString pp_pfname=stok;
        pp_pfname.Remove(0,pp_pfname.Last('=')+1);
        gOPT.pfname=pp_pfname;
      }
 
      if(stok.Contains("zl_style=")) {
        TString pp_zl_style=stok;
        pp_zl_style.Remove(0,pp_zl_style.Last('=')+1);
        gOPT.zl_style=pp_zl_style;
      }
 
      if(stok.Contains("belt_style=")) {
        TString pp_belt_style=stok;
        pp_belt_style.Remove(0,pp_belt_style.Last('=')+1);
        gOPT.belt_style=pp_belt_style;
      }
 
      if(stok.Contains("ifar_max=")) {
        TString pp_ifar_max=stok;
        pp_ifar_max.Remove(0,pp_ifar_max.Last('=')+1);
        if(pp_ifar_max.IsFloat()) gOPT.ifar_max=pp_ifar_max.Atof();
      }
 
      if(stok.Contains("rho_min=")) {
        TString pp_rho_min=stok;
        pp_rho_min.Remove(0,pp_rho_min.Last('=')+1);
        if(pp_rho_min.IsFloat()) gOPT.rho_min=pp_rho_min.Atof();
      }
 
      if(stok.Contains("sfactor=")) {
        TString pp_sfactor=stok;
        pp_sfactor.Remove(0,pp_sfactor.Last('=')+1);
        if(pp_sfactor.IsFloat()) gOPT.sfactor=pp_sfactor.Atof();
      }
 
      if(stok.Contains("bbh=")) {
        TString pp_bbh=stok;
        pp_bbh.Remove(0,pp_bbh.Last('=')+1);
        if(pp_bbh=="true")  gOPT.bbh=true;
        if(pp_bbh=="false") gOPT.bbh=false;
      }
 
      if(stok.Contains("exit=")) {
        TString pp_exit=stok;
        pp_exit.Remove(0,pp_exit.Last('=')+1);
        if(pp_exit=="true")  gOPT.exit=true;
        if(pp_exit=="false") gOPT.exit=false;
      }
 
      if(stok.Contains("gw151226_ifar=")) {
        TString pp_gw151226_ifar=stok;
        pp_gw151226_ifar.Remove(0,pp_gw151226_ifar.Last('=')+1);
        if(pp_gw151226_ifar.IsFloat()) gOPT.gw151226_ifar=pp_gw151226_ifar.Atof();
      }
 
      if(stok.Contains("zl_max_ifar=")) {
        TString pp_zl_max_ifar=stok;
        pp_zl_max_ifar.Remove(0,pp_zl_max_ifar.Last('=')+1);
        if(pp_zl_max_ifar.IsFloat()) gOPT.zl_max_ifar=pp_zl_max_ifar.Atof();
      }
    }
  }
}
 
void DumpLoudest(double obsTime) {
 
  if(gOPT.pfname=="") return;
 
  char ofname[1024];
  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
  if(gOPT.bbh) sprintf(ofname,"%s_bbh_loudest.txt",PFNAME.Data());
  else         sprintf(ofname,"%s_nobbh_loudest.txt",PFNAME.Data());
 
  char sout[1024];
 
  ofstream out;
  out.open(ofname,ios::out);
  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
  cout << "Create Output File : " << ofname << endl;
  out.precision(19);
 
  // write header
  sprintf(sout,"%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s",
               4, "#run", 6, "chunk", 18, "gps_time", 10, "known.BBH", 16, "ifar(sec)", 16, "ifar(year)", 16, "obs_time(sec)", 
               16, "obs_time(days)", 16, "expected", 10, "observed", 11, "cumul.FAP", 8, "sigma");
  out << sout << endl;
 
  // write data
  int observed=1;
  for (int i=ZL_IFAR.size()-1; i>=0; i--) {
 
    double expected = obsTime/(ZL_IFAR[i]*YEAR); 
 
    // compute FAP
    double FAP = 1.;
    for(int j=0;j<observed;j++) {
      FAP -= TMath::PoissonI(j,expected);
      //cout << " N : " << j << " MU : " << expected << " fap : " << fap << endl;
    }
 
    // compute sigma
    double sigma = sqrt(2)*TMath::ErfcInverse(FAP);
 
    TString bbh_name = GetNameBBH(ZL_GPS[i]);
    if(bbh_name=="") bbh_name=".";
 
    TString run;
    int chunkID = GetChunkID(ZL_GPS[i],run);
 
    sprintf(sout,"%*s %*d %*.6f %*s %*.1f %*.4f %*f %*.1f %*.9f %*d %*.3e %*.2f",
                 4, run.Data(), 6, chunkID, 18, ZL_GPS[i], 10, bbh_name.Data(), 16, 
                 ZL_IFAR[i]*YEAR, 16, ZL_IFAR[i], 16, obsTime, 
                 16, obsTime/DAY, 16, expected, 10,  observed, 11,  FAP, 8,  sigma);
    if(gOPT.zl_style=="marker") {
      out << sout << endl;
      observed++;
    } else {
      if(i%2==0) {out << sout << endl; observed++;}
    }
  }
 
  out.close();
}
 
int  GetChunkID(double gps, TString& run) {
 
  for(int n=0;n<MAX_RUNS;n++) {  
    for(int i=0;i<chunk_size[n];i++) {
      if(gps>chunk_start[n][i] && gps<chunk_stop[n][i]) {
        if(n==0) run="O1";
        if(n==1) run="O2";
        if(n==2) run="O3";
        return chunk_id[n][i];
      }
    }
  }
 
  return -1;
}
 
void DumpPeriod(double* obsTime) {
 
  if(gOPT.pfname=="") return;
 
  char ofname[1024];
  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
  if(gOPT.bbh) sprintf(ofname,"%s_bbh_period.txt",PFNAME.Data());
  else         sprintf(ofname,"%s_nobbh_period.txt",PFNAME.Data());
 
  char sout[1024];
 
  ofstream out;
  out.open(ofname,ios::out);
  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
  cout << "Create Output File : " << ofname << endl;
  out.precision(19);
 
  // write header
  sprintf(sout,"%*s %*s %*s %*s %*s %*s %*s %*s",
               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)"); 
  out << sout << endl;
 
  for(int n=0;n<MAX_RUNS;n++) {  
 
    TString RUN;
 
    if(n==0 && gOPT.run.Contains("O1")) RUN="O1";
    else
    if(n==1 && gOPT.run.Contains("O2")) RUN="O2";
    else
    if(n==2 && gOPT.run.Contains("O3")) RUN="O3";
    else continue;
 
    wat::Time beg_date(chunk_start[n][0]);
    wat::Time end_date(chunk_stop[n][chunk_size[n]-1]);
 
    double interval = (end_date.GetGPS()-beg_date.GetGPS())/DAY;
 
    TString sbeg_date = beg_date.GetDateString();sbeg_date.Resize(19);
    TString send_date = end_date.GetDateString();send_date.Resize(19);
 
    sbeg_date.ReplaceAll(" ","-");
    send_date.ReplaceAll(" ","-");
 
    sprintf(sout,"%*s %*.0d %*s %*.0d %*s %*.2f %*f %*.1f", 4, RUN.Data(), 20, beg_date.GetGPS(), 20, sbeg_date.Data(), 
                                20, end_date.GetGPS(), 20, send_date.Data(), 20, interval, 20, obsTime[n]/gOPT.trials, 20, obsTime[n]/DAY/gOPT.trials);
    out << sout << endl;
  }
 
  out.close();
}
 
void DumpUserOptions() {
 
  if(gOPT.pfname=="") return;
 
  char ofname[1024];
  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
  if(gOPT.bbh) sprintf(ofname,"%s_bbh_config.txt",PFNAME.Data());
  else         sprintf(ofname,"%s_nobbh_config.txt",PFNAME.Data());
 
  ofstream out;
  out.open(ofname,ios::out);
  if(!out.good()) {cout << "DumpUserOptions : Error Opening File : " << ofname << endl;exit(1);}
  cout << "dump: " << ofname << endl;
 
  TString info="";
  TString tabs="\t\t\t\t";
 
  // cWB library
  char version[128];
  sprintf(version,"WAT Version : %s - SVN Revision : %s - Tag/Branch : %s",watversion('f'),watversion('r'),watversion('b'));
 
  // cWB config
 
  // get CWB_CONFIG
  char cwb_config_env[1024] = "";
  if(gSystem->Getenv("CWB_CONFIG")!=NULL) {
    strcpy(cwb_config_env,TString(gSystem->Getenv("CWB_CONFIG")).Data());
  } 
  char cfg_version[256];
  TString cfg_branch = GetGitInfos("branch",cwb_config_env);
  TString cfg_tag    = GetGitInfos("tag",cwb_config_env);
  TString cfg_diff   = GetGitInfos("diff",cwb_config_env);
  if(cfg_branch!="") {
    if(cfg_diff!="") cfg_branch=cfg_branch+"/M";
    sprintf(cfg_version,"cWB Config  Branch\t%s",cfg_branch.Data());
  } else if(cfg_tag!="") {
    if(cfg_diff!="") cfg_branch=cfg_branch+"/M";
    sprintf(cfg_version,"cWB Config  Tag\t%s",cfg_tag.Data());
  } else {
    sprintf(cfg_version,"cWB Config\t%s","Undefined");
  }
 
  out << endl;
  out << "--------------------------------"     << endl;
  out << "PP IFAR version                 "     << endl;
  out << "--------------------------------"     << endl;
  out << endl;
  out << version << endl;
  out << endl;
  out<< cfg_version <<endl;
  out<< "cWB Config  Hash\t" << GetGitInfos("hash",cwb_config_env) <<endl;
  out << endl;
 
  out << "--------------------------------"     << endl;
  out << "PP IFAR config options          "     << endl;
  out << "--------------------------------"     << endl;
  out << endl;
 
  out << "pp_run                          " << gOPT.run << endl;
  info = "// run type: allowed types are O1/O2/O1O2/O3 (Default: O3)";
  out << tabs << info << endl << endl;
 
  out << "pp_search                       " << gOPT.search << endl;
  info = "// search type: allowed types are BBH,BurstLF,BurstHF,BurstLD,IMBHB + bin1/bin2/... (Default: none)";
  out << tabs << info << endl << endl;
 
  out << "pp_nifo                         " << gOPT.nifo << endl;
  info = "// number of detectors  (Default: 2)";
  out << tabs << info << endl << endl;
 
  out << "pp_lag                          " << gOPT.lag << endl;
  info = "// lag   -> used only when the input file name is a root file (Default: -1)";
  out << tabs << info << endl << endl;
 
  out << "pp_slag                         " << gOPT.slag << endl;
  info = "// slag  -> used only when the input file name is a root file (Default: -1)";
  out << tabs << info << endl << endl;
 
  out << "pp_chunk                        " << gOPT.chunk << endl;
  info = "// chunk -> used only when the input file name is a root file (Default: -1)";
  out << tabs << info << endl << endl;
 
  out << "pp_trials                       " << gOPT.trials << endl;
  info = "// trials factor (Default: 1)";
  out << tabs << info << endl << endl;
 
  out << "pp_pfname                       " << gOPT.pfname << endl;
  info = "// output plot file name, included the directory (Default: none)";
  out << tabs << info << endl << endl;
 
  out << "pp_zl_style                     " << gOPT.zl_style << endl;
  info = "// zero lag plot style: step/marker -> zero lag is displayed with step/marker (Default: step)";
  out << tabs << info << endl << endl;
 
  out << "pp_belt_style                   " << gOPT.belt_style << endl;
  info = "// belt plot style: step/continous -> belts are displayed with step/continuos poisson distribution (Default: continuos)";
  out << tabs << info << endl << endl;
 
  out << "pp_ifar_max                     " << gOPT.ifar_max << endl;
  info = "// max ifar (year) used for plot, if =-1 then it is the maximum ifar from root files (Default: -1)";
  out << tabs << info << endl << endl;
 
  out << "pp_rho_min                      " << gOPT.rho_min << endl;
  info = "// min rho[0/1] read from output root files (Default: 4.5)";
  out << tabs << info << endl << endl;
 
  out << "pp_sfactor                      " << gOPT.sfactor << endl;
  info = "// factor used to normalize the simulation events for the computation of astrophysical rates (Default: 1)";
  out << tabs << info << endl << endl;
 
  out << "bbh                             " << gOPT.bbh << endl;
  info = "// false/true -> if false the known bbh events are excluded from analysis (Default: false)";
  out << tabs << info << endl << endl;
 
  out << "exit                            " << gOPT.exit << endl;
  info = "// false/true -> if true program exit ayt the end of the execution (Default: true)";
  out << tabs << info << endl << endl;
 
  out << "pp_gw151226_ifar                " << gOPT.gw151226_ifar << endl;
  info = "// GW151226 IFAR(YEARS): if > 0 the GW151226 event is added to the foreground (Default: 0, not added)";
  out << tabs << info << endl << endl;
 
  out << "pp_zl_max_ifar                  " << gOPT.zl_max_ifar << endl;
  info = "// ZL_MAX_IFAR(YEARS): if > 0 all foreground with IFAR>zl_max_ifar is set to zl_max_ifar (Default: 0, not used)";
  out << tabs << info << endl << endl;
 
}
 
// Dumps Background & Belts in ASCII 
void DumpBackground() {
 
  if(gOPT.pfname=="") return;
 
  char ofname[1024];
  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
  if(gOPT.bbh) sprintf(ofname,"%s_bbh_background.txt",PFNAME.Data());
  else         sprintf(ofname,"%s_nobbh_background.txt",PFNAME.Data());
 
  ofstream out;
  out.open(ofname,ios::out);
  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
  cout << "Create Output File : " << ofname << endl;
  out.precision(19);
 
  // write header
  out << "#background data : ifar(years), expected, lower_sigma1, lower_sigma2, lower_sigma3, upper_sigma1, upper_sigma2, upper_sigma3" << endl;
 
  for(int i=0;i<xIFAR.size();i++) {
 
    double med = bRATE[i];
 
    double l900 = yBKG[i]-eyRATEinf0[i];
    double u900 = yBKG[i]-eyRATEsup0[i];
 
    double l990 = yBKG[i]-eyRATEinf1[i];
    double u990 = yBKG[i]+eyRATEsup1[i];
 
    double l999 = yBKG[i]-eyRATEinf2[i];
    double u999 = yBKG[i]+eyRATEsup2[i];
 
    out << xIFAR[i] << " " << med 
        << " " << l900 << " " << l990 << " " << l999
        << " " << u900 << " " << u990 << " " << u999
        << endl;
  }
 
  out.close();
}
 
// Dumps Foreground in ASCII format 
void DumpForeground() {
 
  if(gOPT.pfname=="") return;
 
  char ofname[1024];
  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
  if(gOPT.bbh) sprintf(ofname,"%s_bbh_foreground.txt",PFNAME.Data());
  else         sprintf(ofname,"%s_nobbh_foreground.txt",PFNAME.Data());
 
  ofstream out;
  out.open(ofname,ios::out);
  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
  cout << "Create Output File : " << ofname << endl;
  out.precision(19);
 
  // write header
  out << "#foreground data : ifar(years), observed" << endl;
 
  for(int i=0;i<ZL_NEVT.size();i++) {
 
    out << ZL_IFAR[i] << " " << ZL_NEVT[i] << endl;
  }
 
  out.close();
}
 
// Dumps Simulation & Belts in ASCII format 
void DumpSimulation(double ifar_cmax) {
 
  if(SIM_NEVT.size()<=0) return;
  if(gOPT.pfname=="") return;
  if(xSIM_BKG_IFAR.size()<=0) return;
 
  char ofname[1024];
  TString PFNAME = gOPT.pfname(0,gOPT.pfname.Last('.'));
  if(gOPT.bbh) sprintf(ofname,"%s_bbh_simulation.txt",PFNAME.Data());
  else         sprintf(ofname,"%s_nobbh_simulation.txt",PFNAME.Data());
 
  ofstream out;
  out.open(ofname,ios::out);
  if (!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
  cout << "Create Output File : " << ofname << endl;
  out.precision(19);
 
  // write header
  out << "#simulation data : ifar(years), estimated+background, lower_90.0, upper_90.0" << endl;
 
  // compute the index which contains the maximum common IFAR (RATE_TAG)
  int ifar_cmax_index=0;
  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;
 
  for(int i=0;i<ifar_cmax_index;i++) {
 
    double med = xSIM_BKG_IFAR[i];
 
    double l90 = med-eySIM_BKG_inf0[i];
    double u90 = med+eySIM_BKG_sup0[i];
 
    out << SIM_IFAR[i] << " " << med << " " << l90 << " " << u90 << endl;
  }
 
  out.close();
}
 
int ReadPublicAlerts() {
 
  if(gOPT.run!="O3") return 0;   // this function is implemente only for O3
 
  // ----------------------------------------------------- 
  // Read Public Alerts Name & GPS times
  // ----------------------------------------------------- 
 
  // get CWB_CONFIG
  char cwb_config_env[1024] = "";
  if(gSystem->Getenv("CWB_CONFIG")!=NULL) {
    strcpy(cwb_config_env,TString(gSystem->Getenv("CWB_CONFIG")).Data());
  } 
 
  char pa_file_list[1024];
  sprintf(pa_file_list,"%s/%s/PUBLIC_ALERTS/%s",cwb_config_env,gOPT.run.Data(),PA_FILE_LIST);
  cout << pa_file_list << endl;
 
  CWB::Toolbox::checkFile(pa_file_list);
 
  // Open chunk list
  ifstream in;
  in.open(pa_file_list,ios::in);
  if (!in.good()) {cout << "Error Opening File : " << pa_file_list << endl;exit(1);}
 
  int isize=0;
  char str[1024];
  int fpos=0;
  while(true) {
    in.getline(str,1024);
    if (!in.good()) break;
    if(str[0] != '#') isize++;
  }
  in.clear(ios::goodbit);
  in.seekg(0, ios::beg);
  if(isize==0)  {cout << "Error : File " << pa_file_list << " is empty" << endl;exit(1);}
 
  // O3 known PA events
  if(gOPT.run=="O3") {
    O3_BBH_NAME.resize(isize);
    O3_BBH_GPS.resize(isize);
  }
 
  int k=0;
  char name[256];
  double time;
  while(true) {
    // IN  -> S190326b  1237605877.39
    in >> name >> time;
    if(!in.good()) break;
    if(name[0]=='#') continue;
    wat::Time sdate(time);
 
    TString month[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
 
    // OUT -> S190326b          - Mar 26 03:24
    char sout_txt[1024];
    sprintf(sout_txt,"%s                - %s %02d %02d:%02d",name,
            month[sdate.GetMonth()-1].Data(),sdate.GetDay(),sdate.GetHour(),sdate.GetMinute());
 
    // O3 known PA events
    if(gOPT.run=="O3") {
      O3_BBH_NAME[k]=name;    
      O3_BBH_GPS[k]=time; 
    }
 
    //cout << sout_txt << endl;
 
    k++;
  }
  in.close();
 
  cout.precision(12);
  cout << endl;
  for(int i=0;i<k;i++) cout << O3_BBH_NAME[i] << "\t" << O3_BBH_GPS[i] << endl;
  cout << endl;
 
  return k;
}