Logo coherent WaveBurst  
Library Reference Guide
Logo
DrawLarsHistogramPRC.C
Go to the documentation of this file.
1 //
2 // Draw Cumulative fraction of the sky as a function of the perc of error region
3 // Author : Gabriele Vedovato
4 // Ex : broot 'DrawLarsHistogramPRC.C("skyloc.lst","MEDIAN50","ADV_SIM_BRST_LF_L1H1V1_2G.M1.png")'
5 // ilist_file format : skymap.root label
6 // skymap.root is a root file produced by the macro DrawSkyMapPRC.C
7 //
8 
9 #define APPLY_FIT
10 
11 #define MAX_FILES 20
12 
13 
14 void DrawLarsHistogramPRC(TString ilist_file, TString pType="", TString ofName="") {
15 
16  if((ofName!="")&&(!ofName.Contains(".png"))) {
17  cout << "Error Outpu File : " << ilist_file.Data() << endl;
18  cout << "Must have .png extension" << endl;
19  exit(1);
20  } else {
21  ofName.ReplaceAll(".png",".gif");
22  }
23 
24  TString ptitle="Cumulative fraction of the sky as a function of the ";
25 
26  if((pType!="MEDIAN50")&&(pType!="MEDIAN90")&&(pType!="WRC50")&&(pType!="EFFICIENCY")&&(pType!="")) {
27  ptitle=pType;
28  } else {
29  if(pType=="") pType="MEDIAN50";
30  if(pType=="MEDIAN50") ptitle=ptitle+"50% error region";
31  if(pType=="MEDIAN90") ptitle=ptitle+"90% error region";
32  if(pType=="WRC50") ptitle=ptitle+"50% normalized residual energy";
33  if(pType=="EFFICIENCY") ptitle=ptitle+"efficiency";
34  }
35 
38 
39  // Open list
40  ifstream in;
41  in.open(ilist_file.Data(),ios::in);
42  if (!in.good()) {cout << "Error Opening File : " << ilist_file.Data() << endl;exit(1);}
43 
44  int size=0;
45  char str[1024];
46  int fpos=0;
47  while(true) {
48  in.getline(str,1024);
49  if (!in.good()) break;
50  if(str[0] != '#') size++;
51  }
52  cout << "size " << size << endl;
53  in.clear(ios::goodbit);
54  in.seekg(0, ios::beg);
55  if (size==0) {cout << "Error : File " << ilist_file.Data() << " is empty" << endl;exit(1);}
56 
57  char sfile[256];
58  char stitle[256];
59  int NFILE=0;
60  while(true) {
61  in >> sfile >> stitle;
62  if(!in.good()) break;
63  if(sfile[0]=='#') continue;
64  ifile[NFILE]=sfile;
65  title[NFILE]=stitle;
66  cout << NFILE+1 << " " << ifile[NFILE] << " " << title[NFILE] << endl;
67  NFILE++;
68  }
69  in.close();
70 
71 #ifdef APPLY_FIT
72  // sigmoid fit function
73  TF1 *f1;
74  f1=new TF1("logNfit",logNfit,2,1024,5);
75  f1->SetParameters(16,0.7,1.,1.,0);
76  f1->SetParNames("ptype","sigma","betam","betap");
77  f1->SetParLimits(0,log10(1.),log10(1024.));
78  f1->SetParLimits(1,0.1,10.);
79  f1->SetParLimits(2,0.1,4.);
80  f1->SetParLimits(3,0.1,4.);
81 #endif
82 
83  TCanvas* canvas;
84  canvas = new TCanvas("lars histogram", "canvas", 35, 46, 817, 472);
85  canvas->Clear();
86  canvas->ToggleEventStatus();
87  canvas->SetGridx(false);
88  canvas->SetGridx(false);
89  canvas->SetFillColor(kWhite);
90  canvas->SetLeftMargin(0.08);
91  canvas->SetBottomMargin(0.13);
92  canvas->SetBorderMode(0);
93 
94  // remove the red box around canvas
95  gStyle->SetFrameBorderMode(0);
96  gROOT->ForceStyle();
97 
98  gROOT->SetStyle("Plain");
99  gPad->UseCurrentStyle();
100  gROOT->ForceStyle();
101 
102  gStyle->SetTitleFont(12,"D");
103  gStyle->SetTextFont(12);
104  gStyle->SetTitleFillColor(kWhite);
105  gStyle->SetLineColor(kWhite);
106  //gStyle->SetTitleAlign(11);
107  gStyle->SetTitleW(1.0);
108  gStyle->SetTitleH(0.050);
109  gStyle->SetTitleY(0.98);
110 
113  for(int n=0;n<NFILE;n++) {
114  cout << "Process File : " << ifile[n].Data() << endl;
115  gNET[n].LoadObject((char*)ifile[n].Data());
116  gSM[n] = gNET[n].GetGskymap();
117  }
118 
119  int L = gSM[0]->size();
120 
121  double pi = TMath::Pi();
122  double sr=0;
123  sr = cos(gSM[0]->theta_1*pi/180.)-cos(gSM[0]->theta_2*pi/180.);
124  sr*= (gSM[0]->phi_2-gSM[0]->phi_1)*180/pi/L;
125  sr = double(int(1000*sr))/1000.;
126  cout << "sky resolution (deg^2) : " << sr << endl;
127 
128  double mean=0;
129  for(int l=0;l<L;l++) mean+=gSM[0]->get(l);
130  mean/=L;
131  cout << mean << endl;
132  double integral = 0;
133 
134  Color_t colors[32] = { 6, 3, 2, 8,43, 7, 8, 4, 797, 2,43, 1, 3, 1, 6, 7,
135  6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7 };
136 
137  int lstyle[32] = {2, 0, 0, 0, 9, 9, 9, 2, 2, 2, 2,43, 1, 3, 1, 6, 7,
138  6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6 };
139 
140  double ptype[MAX_FILES];
141  TH1F* h[MAX_FILES];
142  for(int n=0;n<NFILE;n++) {
143  cout << "Process : " << title[n].Data() << endl;
144  if(pType=="EFFICIENCY") {
145  h[n] = new TH1F(title[n].Data(),title[n].Data(),20,0,1.0001); // NOTE! 1.0001 fix hist when filled with value=1
146  } else if(pType=="WRC50") {
147  h[n] = new TH1F(title[n].Data(),title[n].Data(),20,0,0.3);
148  } else {
149  h[n] = new TH1F(title[n].Data(),title[n].Data(),20,0,log(1024)/log(2));
150  }
151  //h[n]->SetTitle(title[n].Data());
152  h[n]->SetMarkerStyle(20);
153  h[n]->SetMarkerSize(1);
154  h[n]->SetMarkerColor(colors[n]);
155  h[n]->SetLineColor(colors[n]);
156  // h[n]->SetLineStyle(lstyle[n]);
157  h[n]->SetLineWidth(2);
158  h[n]->SetStats(kFALSE);
159  h[n]->GetXaxis()->SetLabelFont(42);
160  h[n]->GetYaxis()->SetLabelFont(42);
161  h[n]->GetYaxis()->SetRangeUser(0,1.0);
162  h[n]->GetXaxis()->SetTitleFont(42);
163  h[n]->GetYaxis()->SetTitleFont(42);
164  if(pType=="EFFICIENCY") {
165  h[n]->GetXaxis()->SetTitle("efficiency ");
166  } else if(pType=="WRC50") {
167  h[n]->GetXaxis()->SetTitle("normalized residual energy");
168  } else {
169  h[n]->GetXaxis()->SetTitle("#sigma^{2} (deg^{2}) ");
170  }
171  h[n]->GetYaxis()->SetTitle("#frac{#Omega}{4#pi}");
172  if(pType=="EFFICIENCY") {
173  h[n]->GetXaxis()->SetRangeUser(0.5,1.0);
174  for(int l=0;l<L;l++) {
175  double sigma2 = gSM[n]->get(l);
176  h[n]->Fill(sigma2);
177  }
178  } else if(pType=="WRC50") {
179  h[n]->GetXaxis()->SetRangeUser(0,0.3);
180  for(int l=0;l<L;l++) {
181  double sigma2 = gSM[n]->get(l);
182  h[n]->Fill(sigma2);
183  }
184  } else {
185  h[n]->GetXaxis()->SetRangeUser(0,log(1024)/log(2));
186  for(int l=0;l<L;l++) {
187  double sigma2 = log(pow(gSM[n]->get(l),2))/log(2);
188  if(sigma2<0) sigma2=0.25;
189  if(gSM[n]->get(l)==0) sigma2=-0.25;
190  h[n]->Fill(sigma2);
191  }
192  }
193 
194  // Normalization
195  integral = 0;
196  for (int i=0;i<=h[n]->GetNbinsX();i++) integral+=h[n]->GetBinContent(i);
197  for (int i=0;i<=h[n]->GetNbinsX();i++) h[n]->SetBinContent(i,h[n]->GetBinContent(i)/integral);
198  for (int i=2;i<=h[n]->GetNbinsX();i++) h[n]->SetBinContent(i,h[n]->GetBinContent(i)+h[n]->GetBinContent(i-1));
199 
200  if(n==0) {
201  if(pType=="EFFICIENCY" || pType=="WRC50") {
202  h[0]->GetXaxis()->SetLabelFont(42);
203  h[0]->GetXaxis()->SetLabelOffset(0.01);
204  h[0]->GetXaxis()->LabelsOption("h");
205  } else {
206  char xlabel[16];
207  for (int i=0;i<=h[0]->GetNbinsX();i++) {
208  double xvalue = h[0]->GetXaxis()->GetBinCenter(i)+h[0]->GetXaxis()->GetBinWidth(i)/2.;
209  if(i%4==0) {
210  sprintf(xlabel,"%2.0f",pow(2,xvalue));
211  } else {
212  sprintf(xlabel,"");
213  }
214  h[0]->GetXaxis()->SetBinLabel(i,xlabel);
215  h[0]->GetXaxis()->SetLabelSize(0.07);
216  h[0]->GetXaxis()->SetLabelFont(42);
217  h[0]->GetXaxis()->SetLabelOffset(0.01);
218  h[0]->GetXaxis()->LabelsOption("h");
219  }
220  }
221  h[0]->SetTitle(ptitle);
222  h[0]->Draw("lp");
223  } else h[n]->Draw("lpsame");
224 
225 #ifdef APPLY_FIT
226  int size = h[n]->GetNbinsX();
227  double* x = new double[size];
228  double* y = new double[size];
229  if(pType=="EFFICIENCY" || pType=="WRC50") {
230  for (int i=1;i<=size;i++) {
231  x[i-1] = h[n]->GetXaxis()->GetBinCenter(i);
232  y[i-1] = h[n]->GetBinContent(i);
233  }
234  } else {
235  for (int i=1;i<=size;i++) {
236  x[i-1] = h[n]->GetXaxis()->GetBinCenter(i)+h[n]->GetXaxis()->GetBinWidth(i)/2.;
237  x[i-1] = pow(2,x[i-1]);
238  y[i-1] = h[n]->GetBinContent(i);
239  //cout << i-1 << " : " << x[i-1] << " " << y[i-1] << endl;
240  }
241  }
242  TGraph* gr = new TGraph(size,x,y);
243  // Fit with sigmoid
244  if(pType=="EFFICIENCY" || pType=="WRC50") {
245  TGraphSmooth* gs = new TGraphSmooth("normal");
246  wavearray<double> xs(1000);for(int i=0;i<xs.size();i++) xs[i]=double(i)/double(xs.size());
247  TGraph* grs = gs->Approx(gr,"linear", xs.size(), xs.data);
248  // find ptype[n] for ys=0.5
249  for(int i=0;i<xs.size();i++) {
250  double ys;
251  grs->GetPoint(i,ptype[n],ys);
252  if(ys>0.5) break;
253  }
254  } else {
255  f1->SetLineColor(colors[n]);
256  f1->SetNpx(10000);
257  gr->Fit("logNfit");
258  ptype[n]=pow(10.,f1->GetParameter(0));
259  double chi2=f1->GetChisquare();
260  double sigma=f1->GetParameter(1);
261  double betam=f1->GetParameter(2);
262  double betap=f1->GetParameter(3);
263  cout << "Fit ptype% " << ptype[n] << endl;
264  //f1->Draw("same");
265  }
266  delete gr;
267 #endif
268  }
269 
270  // draw the legend
271  double hleg = 0.20+NFILE*0.05;
272  TLegend* leg;
273  if(pType=="EFFICIENCY" || pType=="MEDIAN90") {
274  double hleg = 0.8-NFILE*0.05;
275  leg = new TLegend(0.1291513,hleg,0.6482165,0.8738739,NULL,"brNDC");
276  } else {
277  double hleg = 0.20+NFILE*0.05;
278  leg = new TLegend(0.5202952,0.1390135,0.9901599,hleg,NULL,"brNDC");
279  }
280 
281  leg->SetBorderSize(1);
282  leg->SetTextAlign(22);
283  leg->SetTextFont(12);
284  leg->SetLineColor(1);
285  leg->SetLineStyle(1);
286  leg->SetLineWidth(1);
287  leg->SetFillColor(0);
288  leg->SetFillStyle(1001);
289  leg->SetTextSize(0.04);
290  leg->SetLineColor(kBlack);
291  leg->SetFillColor(kWhite);
292 
293  for(int n=0;n<NFILE;n++) {
294  char legLabel[256];
295 #ifdef APPLY_FIT
296  sprintf(legLabel,"%s [ %2.2f ]",title[n].Data(),ptype[n]);
297 #else
298  sprintf(legLabel,"%s",title[n].Data());
299 #endif
300  leg->AddEntry(h[n],legLabel,"lp");
301  }
302  leg->Draw();
303 
304  canvas->SetGridy(true);
305 
306  if(ofName!="") {
307  TString gfileName=ofName;
308  canvas->Print(gfileName);
309  TString pfileName=gfileName;
310  pfileName.ReplaceAll(".gif",".png");
311  char cmd[1024];
312  sprintf(cmd,"convert %s %s",gfileName.Data(),pfileName.Data());
313  cout << cmd << endl;
314  gSystem->Exec(cmd);
315  sprintf(cmd,"rm %s",gfileName.Data());
316  cout << cmd << endl;
317  gSystem->Exec(cmd);
318  }
319 
320  gSystem->Exit(0);
321 }
322 
gskymap * gSM
TString ofName
gnetwork * gNET
TString ptype[nPLOT]
Definition: cwb_mkfad.C:112
int n
Definition: cwb_net.C:28
TString("c")
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
return wmap canvas
Long_t size
Color_t colors[16]
i drho i
double pi
Definition: TestChirp.C:18
virtual size_t size() const
Definition: wavearray.hh:145
char str[1024]
TGraph * gr
wavearray< double > h
Definition: Regression_H1.C:25
double Pi
bool log
Definition: WaveMDC.C:41
int l
void DrawLarsHistogramPRC(TString ilist_file, TString pType="", TString ofName="")
void LoadObject(char *file)
Definition: gnetwork.cc:1506
double sigma
TFile * ifile
gskymap * GetGskymap()
Definition: gnetwork.hh:44
char title[256]
Definition: SSeriesExample.C:1
ifstream in
#define MAX_FILES
char cmd[1024]
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double phi_2
Definition: skymap.hh:324
DataType_t * data
Definition: wavearray.hh:319
double phi_1
Definition: skymap.hh:323
double get(size_t i)
param: sky index
Definition: skymap.cc:699
size_t size()
Definition: skymap.hh:136
wavearray< double > y
Definition: Test10.C:31
Double_t logNfit(Double_t *x, Double_t *par)
Definition: Toolfun.hh:42
exit(0)