DrawSearchAreaPRC.C

Go to the documentation of this file.

//
// Draw PRC Search Area
// Note : this macro is used to generate the PRC report (cwb_report merge_label prc)
// Author : Gabriele Vedovato

#include "TROOT.h"                               
#include "TSystem.h"                             
#include "TFile.h"                             
#include "TTree.h"                             
#include <fstream>                               
#include <iostream>                              
#include "TGraph.h"                              
#include "TMultiGraph.h"                         
#include "TCanvas.h"                             
#include "TH1F.h"                                
#include "TMath.h"                               
#include "TH1F.h"                                
#include "TH2F.h"                                
#include <TComplex.h>                            
#include <TStyle.h>                              
#include <TRandom.h>                             
#include "TVector3.h"                            
#include "TRotation.h"                           
#include "Math/Vector3Dfwd.h"                    
#include "Math/Rotation3D.h"                     

TH1F* sa_hist;
TCanvas* sa_canvas;

void
DrawSearchAreaPRC(TString data_label, TString odir, TString merge_label, int nIFO, float T_win, 
                  int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar) {

  using namespace ROOT::Math;

  char wave_file_name[1024];
  sprintf(wave_file_name,"merge/wave_%s.%s.root",data_label.Data(),merge_label.Data());

  TFile* ifile = TFile::Open(wave_file_name);
  TTree* itree = (TTree *) gROOT->FindObject("waveburst");   

  char cut[1024];
  char tmp[1024];
  sprintf(cut,"abs(time[0]-time[%d])<%f && netcc[%d]>%f && rho[%d]>%f",
          nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
  if(T_vED>0)  {strcpy(tmp,cut);sprintf(cut,"%s && neted[0]/ecor<%f",tmp,T_vED);}
  if(T_pen>0)  {strcpy(tmp,cut);sprintf(cut,"%s && penalty>%f",tmp,T_pen);}
  if(T_ifar>0) {strcpy(tmp,cut);sprintf(cut,"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}

  itree->Draw("erA[0]",cut,"goff");
  int size = itree->GetSelectedRows();
  cout << "detected events : " << size << endl;

  double* erA0 = itree->GetV1();
  Int_t *index = new Int_t[size];

  TMath::Sort(size,itree->GetV1(),index,false);

  double area_min = pow(erA0[index[0]],2);
  double area_max = pow(erA0[index[size-1]],2);
  area_min = area_min ? area_min : 0.1;
  int nbin = TMath::Ceil(area_max/area_min);

  sa_hist = new TH1F("sa_hist","sa_hist",nbin,area_min,area_max);

  for(int i=0;i<size;i++) {
    double sarea = erA0[index[i]]*erA0[index[i]];
    sa_hist->Fill(sarea);
  }

  // Normalization
  double integral = 0;
  for (int i=0;i<=sa_hist->GetNbinsX();i++) integral+=sa_hist->GetBinContent(i);
  for (int i=0;i<=sa_hist->GetNbinsX();i++) sa_hist->SetBinContent(i,sa_hist->GetBinContent(i)/integral);
  for (int i=2;i<=sa_hist->GetNbinsX();i++) sa_hist->SetBinContent(i,sa_hist->GetBinContent(i)+sa_hist->GetBinContent(i-1));
  // find area at 50%
  double search_area50=0;
  for (int i=0;i<=sa_hist->GetNbinsX();i++) if(sa_hist->GetBinContent(i)>0.5) {search_area50=sa_hist->GetBinCenter(i);break;}
  cout << "search_area50 : " << search_area50 << endl;
  // find area at 90%
  double search_area90=0;
  for (int i=0;i<=sa_hist->GetNbinsX();i++) if(sa_hist->GetBinContent(i)>0.9) {search_area90=sa_hist->GetBinCenter(i);break;}
  cout << "search_area90 : " << search_area90 << endl;

  char title[256];
  sprintf(title,"Ndet = %d",size); 
  sa_hist->SetTitle(title);
  sa_hist->GetXaxis()->SetTitle("search area (deg^{2})");
  sa_hist->GetYaxis()->SetTitle("cumulative probability");
  sa_hist->GetXaxis()->SetRangeUser(area_min,area_max);
  sa_hist->GetYaxis()->SetRangeUser(0,1);

  sa_hist->SetStats(kFALSE);
  sa_canvas = new TCanvas("fom", "PRC", 300,40, 600, 600);
  sa_canvas->SetGridx();
  sa_canvas->SetGridy();
  sa_canvas->SetLogx();
  //sa_canvas->SetLogy();

  gStyle->SetTitleH(0.062);
  gStyle->SetTitleW(0.98);
  gStyle->SetTitleY(0.98);
  gStyle->SetTitleFont(72);
  gStyle->SetLineColor(kWhite);

  sa_hist->Draw();  
  sa_hist->GetXaxis()->SetTitleOffset(1.35);
  sa_hist->GetYaxis()->SetTitleOffset(1.35);
  sa_hist->GetXaxis()->CenterTitle(true);
  sa_hist->GetYaxis()->CenterTitle(true);
  sa_hist->GetXaxis()->SetLabelFont(42);
  sa_hist->GetXaxis()->SetTitleFont(42);
  sa_hist->GetYaxis()->SetLabelFont(42);
  sa_hist->GetYaxis()->SetTitleFont(42);

  // draw the legend
  TLegendEntry* entry;
  TLegend *legend = new TLegend(0.3708054,0.174216,0.8892617,0.2770035,NULL,"brNDC");
  legend->SetTextFont(42);
  legend->SetTextSize(0.03);
  legend->SetLineColor(kBlack);
  legend->SetFillColor(kWhite);
  char label[256];
  sprintf(label,"search area 50%% : %0.2f (deg^{2})",search_area50);
  legend->AddEntry(sa_hist,label,"lp");
  sprintf(label,"search area 90%% : %0.2f (deg^{2})",search_area90);
  legend->AddEntry(sa_hist,label,"lp");
  legend->Draw();

  char ofname[1024];
  sprintf(ofname,"%s/search_area.gif",odir.Data());
  sa_canvas->Print(ofname);
  char cmd[1024];
  TString pfname(ofname);
  pfname.ReplaceAll(".gif",".png");
  sprintf(cmd,"convert %s %s",ofname,pfname.Data());
  cout << cmd << endl;
  gSystem->Exec(cmd);
  sprintf(cmd,"rm %s",ofname);
  cout << cmd << endl;
  gSystem->Exec(cmd);
  //exit(0);

}