Logo coherent WaveBurst  
Library Reference Guide
Logo
cwb_report_slags.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 {
20 
21  cout<<"slag starts..."<<endl;
22 
23  if(nIFO>3) {
24  cout << "cwb_report_slag.C : Skipping plot production - wave tree should come from a 2- or 3-fold network : "
25  << net_file_name << endl;
26  exit(1);
27  }
28 
30 
31  TB.checkFile(gSystem->Getenv("CWB_ROOTLOGON_FILE"));
32  TB.checkFile(gSystem->Getenv("CWB_PARAMETERS_FILE"));
33  TB.checkFile(gSystem->Getenv("CWB_UPARAMETERS_FILE"));
34  TB.checkFile(gSystem->Getenv("CWB_PPARAMETERS_FILE"));
35  TB.checkFile(gSystem->Getenv("CWB_UPPARAMETERS_FILE"));
36  TB.checkFile(gSystem->Getenv("CWB_EPPARAMETERS_FILE"));
37 
38  gStyle->SetTitleFillColor(kWhite);
39  gStyle->SetLineColor(kWhite);
40  gStyle->SetNumberContours(256);
41  gStyle->SetCanvasColor(kWhite);
42  gStyle->SetStatBorderSize(1);
43  gStyle->SetOptStat(kFALSE);
44 
45  // remove the red box around canvas
46  gStyle->SetFrameBorderMode(0);
47  gROOT->ForceStyle();
48 
49 
50  // ----------------------------------------------------------
51  // canvas
52  // ----------------------------------------------------------
53  TCanvas *c1 = new TCanvas("Slag", "Slag",32,55,750,502);
54  c1->Range(-0.75,-1.23433,6.75,8.17182);
55  c1->SetBorderSize(1);
56  c1->SetFillColor(0);
57  c1->SetGridx();
58  c1->SetGridy();
59  c1->SetFrameFillColor(0);
60  gStyle->SetPalette(1,0);
61 
62  TChain wave("waveburst");
63  TChain live("liveTime");
64 
65  wave.Add(net_file_name);
66  live.Add(liv_file_name);
67 
68  // check if slag is presents in liv tree
69  TBranch* branch;
70  bool slagFound=false;
71  TIter mynext(wave.GetListOfBranches());
72  while ((branch=(TBranch*)mynext())) {
73  if (TString("slag").CompareTo(branch->GetName())==0) slagFound=true;
74  }
75  mynext.Reset();
76  if(!slagFound) {
77  cout << "cwb_report_slags.C : Error - wave tree does not contain slag leaf : "
78  << net_file_name << endl;
79  exit(1);
80  }
81 
82  char fname[1024];
83  char cut[512];
84  char cut2[256];
85 
86  sprintf(cut,"rho[1]>%f",T_cut);
87  //test if vetoes are applied
89  if(!test_veto.IsNull())sprintf(cut,"%s && %s",cut,veto_not_vetoed);
90  sprintf(cut2,"!(lag[%d]==0&&slag[%d]==0)",nIFO,nIFO);
91  long int num = wave.GetEntries();
92  wave.SetEstimate(num);
93  TString sel("slag[1]:slag[2]");
94  if(nIFO==2) {sel.ReplaceAll("1","0");sel.ReplaceAll("2","1");}
95  num=wave.Draw(sel,cut,"goff");
96  const int slag_entries = num;
97  cout<<"Number of Entries: "<<num<<endl;
98  double* slag1 = new double[slag_entries];
99  slag1 = wave.GetV1();
100  double* slag2 = new double[slag_entries];
101  slag2 = wave.GetV2();
102  char mytitle[256];
103  double SlagMax = wave.GetMaximum("slag")+segLen/2.;
104  double SlagMin = wave.GetMinimum("slag")-segLen/2.;
105  int NSlag = TMath::FloorNint((SlagMax-SlagMin)/segLen);
106  cout << "SLAG MAX : "<< wave.GetMaximum("slag")<< " s SLAG MIN : "<< wave.GetMinimum("slag")<< " s #SLAGS : "<<NSlag-1<<endl;
107 
108  if(NSlag==1){cout <<"Just one slag....Skipping further execution!"<<endl; exit(0);}
109  sprintf(mytitle,"FAR distribution over slags (post cat3 & rho>%f)",T_cut);
110  TH2F* Slag = new TH2F("SLAG",mytitle,NSlag,SlagMin/86400.,SlagMax/86400.,NSlag,SlagMin/86400.,SlagMax/86400.);
111  Slag->GetXaxis()->SetTitle("slag[1] shift [day]");
112  Slag->GetYaxis()->SetTitle("slag[2] shift [day]");
113  Slag->SetStats(kFALSE);
114 
115  TH2F* lSlag = new TH2F("LSLAG","Live time distribution over slags",NSlag,SlagMin/86400.,SlagMax/86400.,NSlag,SlagMin/86400.,SlagMax/86400.);
116  for(long int l=0; l<num; l++){
117  Slag->Fill(slag1[l]/86400.,slag2[l]/86400.);
118  }
119 
120  float xlag[NIFO_MAX+1];
121  float xslag[NIFO_MAX+1];
122  double xlive;
123  live.SetBranchAddress("slag",xslag);
124  live.SetBranchAddress("lag",xlag);
125  live.SetBranchAddress("live",&xlive);
126  live.SetBranchStatus("gps",false);
127 
128  int ntrg = live.GetEntries();
129  for(int i=0; i<ntrg; i++) {
130  live.GetEntry(i);
131  lSlag->Fill(xslag[0]/86400.,xslag[1]/86400.,xlive);
132  }
133 
134  Slag->Divide(lSlag);
135 
136  TH1D* lSlag2 = lSlag->ProjectionY();
137  lSlag2->GetXaxis()->SetTitle("slag shift [day]");
138  lSlag2->GetYaxis()->SetTitle("Live time [day]");
139  //lSlag2->GetYaxis()->SetNdivisions(10,kFALSE);
140  lSlag2->SetMarkerStyle(20);
141  lSlag2->Scale(1./86400); //Scaling to get the live time in days
142  sprintf(mytitle,"Live time (%zu lags) over slag shift",lagSize);
143  lSlag2->SetTitle(mytitle);
144  lSlag2->Draw("P");
145  sprintf(fname,"%s/Live_vs_slagtime.png",netdir);
146  c1->Update(); c1->SaveAs(fname);
147 
148  TAxis *xaxis = Slag->GetXaxis();
149  TAxis *yaxis = Slag->GetYaxis();
150  int mycount = 0;
151  double FAR[50000];
152  double LSlag[50000];
153  double eFAR[50000];
154  double D[50000];
155  double Dx[50000];
156  double Dy[50000];
157  for(int i =1 ;i<NSlag+1;i++){
158  for(int j =1 ;j<NSlag+1;j++){
159  if(Slag->GetBinContent(i,j)>=1.E-15){
160  FAR[mycount] = Slag->GetBinContent(i,j)*86400.*365;
161  LSlag[mycount] = lSlag->GetBinContent(i,j);
162  LSlag[mycount] = lSlag->GetBinContent(i,j);
163  Dx[mycount] = xaxis->GetBinCenter(i);
164  Dy[mycount] = yaxis->GetBinCenter(j);
165  D[mycount] = TMath::Sqrt(pow(Dx[mycount],2)+pow(Dy[mycount],2));
166  eFAR[mycount] = Slag->GetBinError(i,j)/lSlag->GetBinError(i,j)*86400.*365;
167  mycount++;
168  if(i*j%100==0) {
169  cout << "slag[1] :"<<xaxis->GetBinCenter(i)<<" slag[2] :"<<yaxis->GetBinCenter(j)
170  <<" FAR :"<<FAR[mycount-1] << "+/-" <<eFAR[mycount-1] << " " <<endl;}
171  }
172  }
173  }
174 
175  cout << "Number of slags : " << mycount << endl;
176  double liveTot = lSlag->GetSumOfWeights();
177  double liveMax = lSlag->GetMaximum();
178  cout << "Total BKG live time : " << liveTot << endl;
179  double* eN = new double[mycount];
180  for(int i =0 ;i<mycount;i++) eN[i]=0.0;
181 
182  TCanvas *c2 = new TCanvas("Slag2", "Slag2",32,55,1450,502);
183  c2->SetBorderSize(1);
184  c2->SetFillColor(0);
185  c2->SetGridx();
186  c2->SetGridy();
187  c2->SetFrameFillColor(0);
188 
189  sprintf(fname,"%s/FAR_vs_slag.txt",netdir);
190  FILE* frate = fopen(fname,"w");
191  fprintf(frate,"#FAR[years^-1] eFAR Tshift[days]\n");
192  for(int i=0; i<mycount;i++){fprintf(frate,"%.4f %.4f %.4f\n",FAR[i],eFAR[i],Dy[i]);}
193  fclose(frate);
194 
195  TGraphErrors* FarPlot= new TGraphErrors(mycount, Dy, FAR, eN, eFAR);
196  FarPlot->SetLineColor(kBlue);
197  FarPlot->SetMarkerColor(kBlue);
198  FarPlot->SetTitle("");
199  FarPlot->SetMarkerStyle(20);
200  FarPlot->SetMarkerSize(1);
201  FarPlot->SetLineStyle(3);
202  FarPlot->GetHistogram()->GetXaxis()->SetRangeUser(TMath::FloorNint(SlagMin/86400.),TMath::CeilNint(SlagMax/86400.));
203  FarPlot->GetHistogram()->SetXTitle("time shift [day]");
204  FarPlot->GetHistogram()->SetYTitle("FAR [1/year]");
205  FarPlot->GetHistogram()->GetXaxis()->SetTitleSize(0.05);
206  FarPlot->GetHistogram()->GetYaxis()->SetTitleSize(0.05);
207  FarPlot->GetHistogram()->GetXaxis()->SetLabelSize(0.05);
208  FarPlot->GetHistogram()->GetYaxis()->SetLabelSize(0.05);
209  c2->Clear();
210  c2->SetLogy(kFALSE);
211 
212  FarPlot->Draw("AP");
213  TF1 *g1 = new TF1("g1","pol1",0,SlagMax/86400.);
214  g1->SetLineColor(2);
215  g1->SetLineWidth(5);
216  TF1 *g2 = new TF1("g2","pol1",-SlagMax/86400.,0);
217  g2->SetLineColor(3);
218  g2->SetLineWidth(5);
219  FarPlot->Fit(g1,"R");
220  FarPlot->Fit(g2,"R+","SAME");
221  legr = new TLegend(0.6,0.8,0.885,0.98,"","brNDC");
222  legl = new TLegend(0.15,0.8,0.4,0.98,"","brNDC");
223  char lab[256];
224  sprintf(lab,"#chi^{2}/ndf(right) : %3.3g / %d",g1->GetChisquare(),g1->GetNDF());
225  legr->AddEntry(g1, lab, "l");
226  sprintf(lab,"P0(right) : %.2f +/- %.2f",g1->GetParameters()[0],g1->GetParErrors()[0]);
227  legr->AddEntry("", lab, "a");
228  sprintf(lab,"P1(right) : %.2g +/- %.2g",g1->GetParameters()[1],g1->GetParErrors()[1]);
229  legr->AddEntry("", lab, "a");
230  legr->SetLineColor(1);
231 
232  sprintf(lab,"#chi^{2}/ndf(left) : %3.3g / %d",g2->GetChisquare(),g2->GetNDF());
233  legl->AddEntry(g2, lab, "l");
234  sprintf(lab,"P0(left) : %.2f +/- %.2f",g2->GetParameters()[0],g2->GetParErrors()[0]);
235  legl->AddEntry("", lab, "a");
236  sprintf(lab,"P1(left) : %.2g +/- %.2g",g2->GetParameters()[1],g2->GetParErrors()[1]);
237  legl->AddEntry("", lab, "a");
238  legl->SetLineColor(1);
239 
240  legr->Draw();
241  legl->Draw();
242  g1->Draw("same");
243  g2->Draw("same");
244  sprintf(fname,"%s/FAR_vs_slag.png",netdir);
245  c2->Update(); c2->SaveAs(fname);
246 
247  exit(0);
248 }
float xlag[NIFO_MAX+1]
double liveMax
bool slagFound
TF1 * g2
TCanvas * c2
TString("c")
TBranch * branch
CWB::Toolbox TB
double Dy[50000]
TH1D * lSlag2
int j
Definition: cwb_net.C:28
i drho i
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:4670
s veto_not_vetoed
char cut2[256]
#define nIFO
TF1 * g1
TChain live("liveTime")
TH2F * lSlag
double FAR[50000]
double LSlag[50000]
double eFAR[50000]
char netdir[1024]
const int NIFO_MAX
Definition: wat.hh:22
const int slag_entries
char lab[256]
double D[50000]
char cut[512]
int l
double xlive
segLen
Definition: cwb_eced.C:24
char fname[1024]
double liveTot
TString sel("slag[1]:slag[2]")
int mycount
exit(0)
fprintf(frate,"#FAR[years^-1] eFAR Tshift[days]\)
char net_file_name[256]
char liv_file_name[256]
FILE * frate
TAxis * yaxis
TGraphErrors * FarPlot
TIter mynext(wave.GetListOfBranches())
float xslag[NIFO_MAX+1]
long int num
double T_cut
fclose(frate)
TString test_veto
TChain wave("waveburst")
TAxis * xaxis
sprintf(cut,"rho[1]>%f", T_cut)
TCanvas * c1
int ntrg
double Dx[50000]
size_t lagSize
Definition: test_config1.C:52