Logo coherent WaveBurst  
Library Reference Guide
Logo
DrawPvalueHM.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 // this macro produces the PE HM plot
19 //
20 // How to use: root 'DrawPvalueHM("roc_config.txt")'
21 //
22 // roc_config.txt format
23 //
24 // signal_statistic_1.txt noise_statistic_1.txt roc#1
25 // signal_statistic_2.txt noise_statistic_2.txt roc#2
26 //
27 // the files *_statistic_*.txt are produced in the final output directory by the cwb_pereport command
28 // - */distributions/pestat_ResidualEnergy_statistic.txt
29 //
30 // roc#* are the labels displayed in the output plot legend
31 //
32 
33 
34 #define nHM_MAX 10
35 
36 //#define SHOW_RESIDUALS
37 #define DUMP_PVALUES
38 
39 int readPvalue(TString fname, vector<float>& valpha, vector<float>& vpvalue, vector<float>& vre);
40 int readConfig(TString fname, vector<TString>& vfpvalue, vector<TString>& vname);
41 void PlotPvalueHM(int nroc, TString ofname="");
42 
44 
49 
51 
52 void DrawPvalueHM(TString ifname, TString ofname="") {
53 
54  // init blind colors
55  gCOLOR[0] = CWB::Toolbox::getTableau10BlindColor("DeepSkyBlue4");
56  gCOLOR[1] = CWB::Toolbox::getTableau10BlindColor("DarkOrange1");
65 
66  vector<TString> vfpvalue;
67  vector<TString> vname;
68  int nHM = readConfig(ifname, vfpvalue, vname);
69 
70  for(int n=0;n<nHM;n++) {
71 
72  gNAME[n] = vname[n];
73  gNAME[n].ReplaceAll("*"," ");
74 
75  vector<float> valpha,vpvalue,vre;
76 
77  int size = readPvalue(vfpvalue[n], valpha, vpvalue, vre);
78  wSIZE[n]=size;
79 
80  wALPHA[n].resize(size);
81  wPVALUE[n].resize(size);
82  wRE[n].resize(size);
83  int j=0;
84  for(int i=0;i<size;i++) {
85 #ifdef DUMP_PVALUES
86  if(valpha[i]>=0.2 && valpha[i]<3.5) {
87  cout << valpha[i] << "\t" << vpvalue[i] << endl;
88  }
89 #endif
90  wALPHA[n][j] = valpha[i];
91  wPVALUE[n][j] = vpvalue[i];
92  wRE[n][j] = vre[i];
93  j++;
94  }
95  wSIZE[n]=j;
96  wALPHA[n].resize(j);
97  wPVALUE[n].resize(j);
98  wRE[n].resize(j);
99  }
100 
101  PlotPvalueHM(nHM, ofname);
102 }
103 
104 void PlotPvalueHM(int nHM, TString ofname) {
105 
106  // create plots
107  gStyle->SetFrameBorderMode(0); // remove the red box around canvas
108  gROOT->ForceStyle();
109 
110  gStyle->SetMarkerColor(50);
111  gStyle->SetLineColor(kWhite);
112  gStyle->SetTitleW(0.98);
113  gStyle->SetTitleH(0.08);
114  gStyle->SetTitleY(0.98);
115  gStyle->SetFillColor(kWhite);
116  gStyle->SetLineColor(kWhite);
117  gStyle->SetTitleFont(12,"D");
118 // gStyle->SetTitleSize(0.18);
119 
120  TCanvas *canvas = new TCanvas("Higher Modes", "Higher Modes", 300, 40, 800, 500);
121  canvas->Clear();
122  canvas->ToggleEventStatus();
123 // canvas->SetLogx();
124  canvas->SetLogy();
125  canvas->SetGridx();
126  canvas->SetGridy();
127  canvas->SetFillColor(kWhite);
128 
129  double xmin=2e20;
130  double xmax=0;
131  double ymin=2e20;
132  double ymax=0;
133  TGraph* gr[nHM_MAX];
134  for(int n=0;n<nHM;n++) {
135 #ifdef SHOW_RESIDUALS
136  gr[n] = new TGraph(wSIZE[n],wALPHA[n].data,wRE[n].data);
137 #else
138  gr[n] = new TGraph(wSIZE[n],wALPHA[n].data,wPVALUE[n].data);
139 #endif
140  for(int i=0;i<wSIZE[n];i++) {
141  if(wALPHA[n][i]<xmin) xmin=wALPHA[n][i];
142  if(wALPHA[n][i]<xmax) xmax=wALPHA[n][i];
143  if(wPVALUE[n][i]<ymin) ymin=wPVALUE[n][i];
144  if(wPVALUE[n][i]<ymax) ymax=wPVALUE[n][i];
145  }
146  }
147 
148  for(int n=0;n<nHM;n++) {
149  gr[n]->SetLineWidth(2);
150  gr[n]->SetMarkerColor(gCOLOR[n]);
151  gr[n]->SetMarkerStyle(20);
152  gr[n]->SetMarkerSize(0.8);
153  gr[n]->SetLineColor(gCOLOR[n]);
154  gr[n]->SetLineStyle(1);
155  }
156 
157  TMultiGraph* mg = new TMultiGraph();
158  char gTitle[256];
159  sprintf(gTitle,"Higher Modes");
160  mg->SetTitle(gTitle);
161  for(int n=0;n<nHM;n++) mg->Add(gr[n]);
162  //mg->Paint("APL");
163  mg->Paint("AP");
164 
165  mg->GetHistogram()->GetXaxis()->SetLabelSize(0.04);
166  mg->GetHistogram()->GetYaxis()->SetLabelSize(0.04);
167  mg->GetHistogram()->GetXaxis()->SetTitleSize(0.04);
168  mg->GetHistogram()->GetYaxis()->SetTitleSize(0.04);
169  mg->GetHistogram()->GetXaxis()->SetLabelFont(42);
170  mg->GetHistogram()->GetYaxis()->SetLabelFont(42);
171  mg->GetHistogram()->GetYaxis()->SetLabelOffset(0.01);
172  mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.5);
173 
174  mg->GetHistogram()->GetXaxis()->SetRangeUser(0.001,7.5);
175 #ifndef SHOW_RESIDUALS
176  mg->GetHistogram()->GetYaxis()->SetRangeUser(0.001,1);
177 #endif
178 
179  mg->GetHistogram()->GetXaxis()->SetNdivisions(509);
180  mg->GetHistogram()->GetYaxis()->SetNdivisions(509);
181 
182  mg->GetXaxis()->SetTitle(gr[0]->GetXaxis()->GetTitle());
183  mg->GetXaxis()->SetLabelFont(42);
184  mg->GetYaxis()->SetLabelFont(42);
185  mg->GetXaxis()->SetTitleFont(42);
186  mg->GetYaxis()->SetTitleFont(42);
187  mg->GetXaxis()->SetTitleOffset(0.80);
188  mg->GetYaxis()->SetTitleOffset(0.90);
189  mg->GetXaxis()->SetTitleSize(0.05);
190  mg->GetYaxis()->SetTitleSize(0.05);
191  mg->GetXaxis()->SetTitle("alpha");
192 #ifdef SHOW_RESIDUALS
193  mg->GetYaxis()->SetTitle("residuals");
194 #else
195  mg->GetYaxis()->SetTitle("p-value");
196 #endif
197  mg->GetXaxis()->CenterTitle(true);
198  mg->GetYaxis()->CenterTitle(true);
199 
200  mg->Draw("ALP");
201 
202  // draw the legend
203 
204  TLegend* leg;
205  double hleg = 0.18+nHM*0.06;
206  leg = new TLegend(0.593985,0.128692,0.8834586,hleg,NULL,"brNDC");
207 
208  leg->SetBorderSize(1);
209  leg->SetTextAlign(22);
210  leg->SetTextFont(12);
211  leg->SetLineColor(1);
212  leg->SetLineStyle(1);
213  leg->SetLineWidth(1);
214  leg->SetFillColor(0);
215  leg->SetFillStyle(1001);
216  leg->SetTextSize(0.05);
217  leg->SetLineColor(kBlack);
218  leg->SetFillColor(kWhite);
219 
220  for(int n=0;n<nHM;n++) {
221  char legLabel[256];
222  sprintf(legLabel,"%s",gNAME[n].Data());
223  leg->AddEntry(gr[n],legLabel,"lp");
224  }
225  leg->Draw();
226 
227  // save plot
228  if(ofname!="") canvas->Print(ofname);
229 
230  return;
231 }
232 
233 int readPvalue(TString fname, vector<float>& valpha, vector<float>& vpvalue, vector<float>& vre) {
234 
235  float a1,a2,pvalue,re;
236 
237  ifstream in;
238  in.open(fname.Data(),ios::in);
239  if (!in.good()) {cout << "Error Opening File : " << fname.Data() << endl;exit(1);}
240 
241 // S190814bv : ResidualEnergy 2.0572 pvalue 0.0537 ccut_wdm_fres 16 ccut_bchirp 2.2500 ccut_uchirp 2.5500 ccut_ltime 0.5000 ccut_rtime 0.0000
242 
243  char dy[32];
244 
245  while (1) {
246  in >> dy >> dy >> dy >> re >> dy >> pvalue >> dy >> dy >> dy >> a1 >> dy >> a2 >> dy >> dy >> dy >> dy;
247  if (!in.good()) break;
248  valpha.push_back((a1+a2)/2);
249  vpvalue.push_back(pvalue);
250  vre.push_back(re);
251  }
252  in.close();
253 
254  return vpvalue.size();
255 }
256 
257 int readConfig(TString fname, vector<TString>& vfpvalue, vector<TString>& vname) {
258 
259  char fpvalue[1024];
260  char name[1024];
261 
262  ifstream in;
263  in.open(fname.Data(),ios::in);
264  if (!in.good()) {cout << "Error Opening File : " << fname.Data() << endl;exit(1);}
265 
266  while (1) {
267  in >> fpvalue >> name;
268  if (!in.good()) break;
269  if(TString(fpvalue[0])=="#") continue;
270  vfpvalue.push_back(fpvalue);
271  vname.push_back(name);
272  }
273  in.close();
274 
275  return vname.size();
276 }
277 
TString gNAME[nHM_MAX]
Definition: DrawPvalueHM.C:50
wavearray< double > wALPHA[nHM_MAX]
Definition: DrawPvalueHM.C:45
int gCOLOR[nHM_MAX]
Definition: DrawPvalueHM.C:43
par [0] name
int n
Definition: cwb_net.C:28
TString("c")
wavearray< double > wPVALUE[nHM_MAX]
Definition: DrawPvalueHM.C:46
return wmap canvas
Long_t size
int j
Definition: cwb_net.C:28
i drho i
wavearray< double > wRE[nHM_MAX]
Definition: DrawPvalueHM.C:47
#define nHM_MAX
Definition: DrawPvalueHM.C:34
TGraph * gr
cout<<"Number of Entries: "<< num<< endl;double *slag1=new double[slag_entries];slag1=wave.GetV1();double *slag2=new double[slag_entries];slag2=wave.GetV2();char mytitle[256];double SlagMax=wave.GetMaximum("slag")+segLen/2.;double SlagMin=wave.GetMinimum("slag") -segLen/2.;int NSlag=TMath::FloorNint((SlagMax-SlagMin)/segLen);cout<< "SLAG MAX : "<< wave.GetMaximum("slag")<< " s SLAG MIN : "<< wave.GetMinimum("slag")<< " s #SLAGS : "<< NSlag-1<< endl;if(NSlag==1){cout<<"Just one slag....Skipping further execution!"<< endl;exit(0);} sprintf(mytitle,"FAR distribution over slags (post cat3 & rho>%f)", T_cut);TH2F *Slag=new TH2F("SLAG", mytitle, NSlag, SlagMin/86400., SlagMax/86400., NSlag, SlagMin/86400., SlagMax/86400.);Slag-> GetXaxis() -> SetTitle("slag[1] shift [day]")
char fname[1024]
void DrawPvalueHM(TString ifname, TString ofname="")
Definition: DrawPvalueHM.C:52
double vpvalue[NTRIALS][NPVALUES]
ifstream in
static double a2
Definition: geodesics.cc:26
void PlotPvalueHM(int nroc, TString ofname="")
Definition: DrawPvalueHM.C:104
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
int wSIZE[nHM_MAX]
Definition: DrawPvalueHM.C:48
int readPvalue(TString fname, vector< float > &valpha, vector< float > &vpvalue, vector< float > &vre)
Definition: DrawPvalueHM.C:233
static Int_t getTableau10BlindColor(Int_t index)
Definition: Toolbox.cc:7275
int readConfig(TString fname, vector< TString > &vfpvalue, vector< TString > &vname)
Definition: DrawPvalueHM.C:257
virtual void resize(unsigned int)
Definition: wavearray.cc:463
vector< TString > vname
Definition: cwb_report_pe.C:94
TMultiGraph * mg
exit(0)