Draw_FARvsRHO.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/>.
*/
 
 
 
#define nFAR_MAX        32
 
#define RHO_MIN      4.5
//#define RHO_MIN    0.0
 
int nFAR;
TString FAR_PLOT_NAME;
TString FAR_TITLE;
TString FAR_NAME[nFAR_MAX];
int LINE_COLOR[nFAR_MAX];
int LINE_STYLE[nFAR_MAX];
int LINE_MARKER[nFAR_MAX];
double OBS_TIME[nFAR_MAX];
TString FAR_PATH[nFAR_MAX];
double ZERO_LAG_TIME[nFAR_MAX];
 
int readParameters(TString fname, wavearray<double>& RHO, wavearray<double>& eRHO, wavearray<double>& FAR, wavearray<double>& eFAR);
void writeParameters(TString fname, wavearray<double>& RHO, wavearray<double>& FAR);
void PlotFAR(int nfar, TString far_plot_name, double rho_min=-1, double rho_max=-1);
 
TString PLOT_PATH = "./";
 
void 
Draw_FARvsRHO(TString ifile_bkg, TString ifile_fgd, TString ofile="") {
 
 
  Style_t markers[32]= {20,21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,
                        21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,20 };
 
  Color_t colors[32] = { 8, 0, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7,
                         6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7 };
 
  for(int i=0;i<nFAR_MAX;i++) {
    LINE_COLOR[i]  = colors[i];
    LINE_STYLE[i]  = 1;
    LINE_MARKER[i] = markers[i];
    ZERO_LAG_TIME[i] = -1; 
    OBS_TIME[i]  = 1;
  }
 
  FAR_TITLE="";
  FAR_PLOT_NAME="";
 
  nFAR = 2;
  FAR_PATH[0]=ifile_bkg;
  FAR_PATH[1]=ifile_fgd;
 
  PlotFAR(nFAR, ofile, RHO_MIN, -1);
 
}
 
void PlotFAR(int nfar, TString far_plot_name, double rho_min, double rho_max) {
 
  // create plots
  gStyle->SetFrameBorderMode(0);     // remove the red box around canvas
  gROOT->ForceStyle();               
 
  gStyle->SetTitleFont(72);
  gStyle->SetMarkerColor(50);
  gStyle->SetLineColor(kWhite);
  gStyle->SetTitleW(0.98);     
  gStyle->SetTitleH(0.05);     
  gStyle->SetTitleY(0.98);     
  gStyle->SetFillColor(kWhite);
  gStyle->SetLineColor(kWhite);
  gStyle->SetTitleFont(12,"D");
 
  TCanvas *canvas = new TCanvas("roc", "roc", 300,40, 800, 600);
  canvas->Clear();                                                       
  canvas->ToggleEventStatus();                                           
  canvas->SetLogy();                                                     
  canvas->SetGridx();                                                    
  canvas->SetGridy();                                                    
  canvas->SetFillColor(kWhite);                                          
 
  double far_min=1e10;
  double rmin=1e20;
  double rmax=0;
  char fname[1024];
  int nGR=0;
  int line_style[3*nFAR_MAX];
  int line_marker[3*nFAR_MAX];                                    
  int line_color[3*nFAR_MAX];                                    
  TString FAR_name[3*nFAR_MAX];
  bool sigma_lines[3*nFAR_MAX];
  TGraphErrors* gr[3*nFAR_MAX];
  for(int n=0;n<nFAR;n++) {
    cout << n << " OBS_TIME : " << OBS_TIME[n] << endl;
    if(OBS_TIME[n]>0) if(far_min>1./OBS_TIME[n]) far_min=1./OBS_TIME[n];
    vector<double> rho_far;
    vector<double> far;
    // read FAR
    sprintf(fname,"%s",FAR_PATH[n].Data());
    cout << "rate_threshold : " << fname << endl;
    wavearray<double> RHO;
    wavearray<double> FAR;
    wavearray<double> eRHO;
    wavearray<double> eFAR;
    int far_size = readParameters(fname, RHO, eRHO, FAR, eFAR);
    gr[nGR] = NULL;  
    if(far_size==0) continue; 
    if(rmin<RHO[0]) rmin=RHO[0];
    if(rmax<RHO[far_size-1]) rmax=RHO[far_size-1];
    sigma_lines[nGR] = false;
    line_style[nGR]=LINE_STYLE[n];
    line_marker[nGR]=LINE_MARKER[n];                                    
    line_color[nGR]=LINE_COLOR[n];                                    
    FAR_name[nGR]=FAR_NAME[n];
    gr[nGR++] = new TGraphErrors(far_size,RHO.data,FAR.data,eRHO.data,eFAR.data);
    gr[nGR-1]->SetMarkerStyle(20);
    //gr[nGR-1]->SetMarkerSize(0.35);
    gr[nGR-1]->SetMarkerSize(0.55);
    if(n==0) gr[nGR-1]->SetMarkerColor(kBlack);
    if(n==1) {
      gr[nGR-1]->SetMarkerColor(kRed);
      gr[nGR-1]->SetMarkerStyle(22);
      gr[nGR-1]->SetMarkerSize(1.3);
    }
    if(n==2) gr[nGR-1]->SetMarkerColor(kRed);
    gr[nGR-1]->SetName(TString::Format("gr%d",nGR-1));
  }
 
  TMultiGraph* mg = new TMultiGraph();
  char gTitle[256]; 
  if(FAR_TITLE!="") strcpy(gTitle,FAR_TITLE.Data());
  else              sprintf(gTitle,"FAR Comparison");
  if(FAR_PATH[1].Contains("lag0_slag0")) sprintf(gTitle,"FAR vs Rank");
  else                                   sprintf(gTitle,"FAR vs Rank (fake lag)");
  mg->SetName("mg");
  mg->SetTitle(gTitle); 
  for(int n=0;n<nGR;n++) if(gr[n]!=NULL) mg->Add(gr[n]);   
  //mg->Paint("APL");                        
  //mg->Paint("AL");                        
  mg->Paint("AP");                        
 
  mg->GetHistogram()->GetXaxis()->SetLabelSize(0.03);
  mg->GetHistogram()->GetYaxis()->SetLabelSize(0.03);
  mg->GetHistogram()->GetXaxis()->SetTitleSize(0.03);
  mg->GetHistogram()->GetYaxis()->SetTitleSize(0.03);
  mg->GetHistogram()->GetXaxis()->SetLabelFont(42);  
  mg->GetHistogram()->GetYaxis()->SetLabelFont(42);  
  mg->GetHistogram()->GetYaxis()->SetLabelOffset(0.01);
  mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.5); 
//  mg->GetHistogram()->SetMinimum(far_min/2.); 
/*
  if(rho_min>=0) rmin=rho_min*sqrt(2);
  if(rho_max>=0) rmax=rho_max*sqrt(2);
  rmax=16*sqrt(2);
*/
  if(rho_min>=0) rmin=rho_min;
  if(rho_max>=0) rmax=rho_max;
//  if(rmin<RHO_MIN) rmin=RHO_MIN;
//  mg->GetHistogram()->GetXaxis()->SetRangeUser(rmin,rmax);
 
  mg->GetXaxis()->SetTitle(gr[0]->GetXaxis()->GetTitle());
  mg->GetXaxis()->SetLabelFont(42);                       
  mg->GetYaxis()->SetLabelFont(42);                       
  mg->GetXaxis()->SetTitleFont(42);                       
  mg->GetYaxis()->SetTitleFont(42);                       
  mg->GetXaxis()->SetTitleOffset(1.20);                   
  mg->GetYaxis()->SetTitleOffset(1.20);                   
  mg->GetXaxis()->SetTitleSize(0.03);                     
  mg->GetYaxis()->SetTitleSize(0.03);                     
  mg->GetXaxis()->CenterTitle(true);            //PAPER
  mg->GetYaxis()->CenterTitle(true);            //PAPER
  mg->GetXaxis()->SetTitle("Coherent Network SNR ( #rho )");
  //mg->GetXaxis()->SetTitle("#eta_{c}");
  mg->GetYaxis()->SetTitle("FAR ( yr^{-1} )");
 
  mg->Draw("ALP");
 
  // draw the legend
 
  TLegend* leg;                
//  leg = new TLegend(0.62,0.72,0.90,0.90,NULL,"brNDC");
//  leg = new TLegend(0.72,0.85,1.0,1.0,NULL,"brNDC");  
  leg = new TLegend(0.72,0.87,1.0,1.0,NULL,"brNDC");  
 
  leg->SetBorderSize(1);
  leg->SetTextAlign(22);
  leg->SetTextFont(12); 
  leg->SetLineColor(1); 
  leg->SetLineStyle(1); 
  leg->SetLineWidth(1); 
  leg->SetFillColor(0); 
  leg->SetFillStyle(1001);
  leg->SetTextSize(0.03); 
  leg->SetLineColor(kBlack);
  leg->SetFillColor(kWhite);
 
  char legLabel[256];
  for(int n=0;n<nGR;n++) {
    if(n==0) strcpy(legLabel,"Background");
    if(n==1) strcpy(legLabel,"Foreground");
    if(gr[n]!=NULL) TLegendEntry* eleg = leg->AddEntry(gr[n],legLabel,"lp");  
  }
  leg->SetBorderSize(1);
  leg->SetTextAlign(22);
  leg->SetTextFont(12);
  leg->SetLineColor(1);
  leg->SetLineStyle(1);
  leg->SetLineWidth(1);
  leg->SetFillColor(0);
  leg->SetFillStyle(1001);
  leg->SetTextSize(0.05);
  leg->SetLineColor(kBlack);
  leg->Draw();
 
  // save plot
  if(far_plot_name!="") {
    char gfileName[1024];
    sprintf(gfileName,"%s.gif",far_plot_name.Data());
    canvas->Print(gfileName);
    TString pfileName=gfileName;
    pfileName.ReplaceAll(".gif",".png");
    char cmd[1024];
    sprintf(cmd,"convert %s %s",gfileName,pfileName.Data());
    cout << cmd << endl;
    gSystem->Exec(cmd);
    sprintf(cmd,"rm %s",gfileName);
    cout << cmd << endl;
    gSystem->Exec(cmd);
    gSystem->Exit(0);
  }
 
  return;
}
 
void writeParameters(TString fname, wavearray<double>& RHO, wavearray<double>& FAR) {
 
  ofstream out;
  out.open(fname.Data(),ios::out);
  if (!out.good()) {cout << "Error Opening File : " << fname.Data() << endl;exit(1);}
 
  for(int i=0;i<RHO.size();i++) {
    out << RHO[i] << "\t" << FAR[i] << endl; 
  }
 
  out.close();
}
 
 
 
int readParameters(TString fname, wavearray<double>& RHO, wavearray<double>& eRHO, wavearray<double>& FAR, wavearray<double>& eFAR) {
 
  double year = (24.*3600.*365.);
 
  RHO.resize(1000000);
  eRHO.resize(1000000);
  FAR.resize(1000000);
  eFAR.resize(1000000);
 
  RHO=0.;
  FAR=0.;
 
  ifstream in;
  in.open(fname.Data(),ios::in);
  if (!in.good()) {cout << "Error Opening File : " << fname.Data() << endl;exit(1);}
  cout << fname.Data() << endl;
 
  wavearray<double> rho(1000000);
  wavearray<double> far(1000000);
  int m=0;
  while (1) {
    in >> rho[m] >> far[m];
    if (!in.good()) break;
    if(rho[m]>=RHO_MIN) m++;
  }
  in.close();
  rho.resize(m);
  far.resize(m);
 
  if(rho.size()==0) {
    RHO.resize(0); 
    eRHO.resize(0); 
    FAR.resize(0); 
    eFAR.resize(0); 
    return 0;
  }
 
  double pifar=far[rho.size()-1];
  RHO[0] = rho[rho.size()-1];
  eRHO[0] = 0;
  FAR[0] = far[rho.size()-1];
  eFAR[0] = 0;
  int n=1;
  for(int i=rho.size()-2;i>=0;i--) { 
    if(far[i]>pifar) {
      RHO[n] = rho[i];
      eRHO[n] = 0;
      FAR[n] = far[i];
      eFAR[n] = 0;
      pifar=far[i];
      n++;
    }
  }   
 
  RHO.resize(n); 
  eRHO.resize(n); 
  FAR.resize(n); 
  eFAR.resize(n); 
 
  wavearray<double> xRHO=RHO;
  wavearray<double> xFAR=FAR;
  for(int i=0;i<RHO.size();i++) { 
    RHO[i]=xRHO[RHO.size()-i-1];
    FAR[i]=xFAR[FAR.size()-i-1];
    FAR[i]*=year;
  }
  for(int i=0;i<RHO.size();i++) { 
//    cout << i << " " << RHO[i] << " " << FAR[i] << endl;
  }
 
  return RHO.size(); 
}