source: trunk/MagicSoft/Mars/mtemp/mifae/macros/cleanComp.C@ 6196

Last change on this file since 6196 was 4822, checked in by rico, 21 years ago
*** empty log message ***
File size: 14.2 KB
Line 
1/* ======================================================================== *\!
2! Author(s): Javier Rico
3\* ======================================================================== */
4
5void cleanComp()
6{
7 // general settings
8 gROOT->Reset();
9 gStyle->SetCanvasColor(0);
10 gStyle->SetCanvasBorderMode(0);
11 gStyle->SetPadBorderMode(0);
12 gStyle->SetFrameBorderMode(0);
13 gStyle->SetOptTitle(1);
14 // gStyle->SetTitleOffset(1.7,"y");
15 // gStyle->SetPadLeftMargin(0.15);
16 gStyle->SetOptStat(111110);
17 gStyle->SetStatColor(0);
18 gStyle->SetStatBorderSize(1);
19 gStyle->SetStatW(0.2);
20 gStyle->SetStatH(0.1);
21 gStyle->SetStatX(0.9);
22 gStyle->SetStatY(0.9);
23
24 // define files, chain, trees, canvas...
25 TChain* ton;
26 ton= new TChain("Events");
27
28 // data file(s)
29 ton->Add("/local_disk/jrico/mc/cleanComp0-3.root");
30
31 // define aliases
32 ton->SetAlias("run","MRawRunHeader.fRunNumber");
33 ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
34 ton->SetAlias("energy","MMcEvt.fEnergy");
35
36 // alpha0 is for the file without noise
37 // alphaNOCL is for the file with noise and no cleaning
38 ton->SetAlias("alpha0","MHillasSrc0.fAlpha");
39 ton->SetAlias("alphaNOCL","MHillasSrcnoclean.fAlpha");
40 ton->SetAlias("alpha3015","MHillasSrc3015.fAlpha");
41 ton->SetAlias("alpha4025","MHillasSrc4025.fAlpha");
42 ton->SetAlias("alpha5035","MHillasSrc5035.fAlpha");
43 ton->SetAlias("alpha6045","MHillasSrc6045.fAlpha");
44
45 // ton->Scan("run:event","alpha0<-30 && energy>2500 && energy<5000");
46
47 // plotAlphaWidth(ton);
48 plotAlphaWidthVsEnergy(ton);
49 // plotEfficiencies(ton);
50}
51
52void plotAlphaWidth(TChain* ton)
53{
54 // define and fill histos
55 const Int_t nbins=90;
56 TH1F* halpha3015 = new TH1F("Alpha3015" ,"Alpha3015",nbins,-180,180);
57 TH1F* halpha4025 = new TH1F("Alpha4025" ,"Alpha4025",nbins,-180,180);
58 TH1F* halpha5035 = new TH1F("Alpha5035" ,"Alpha5035",nbins,-180,180);
59 TH1F* halpha6045 = new TH1F("Alpha6045" ,"Alpha6045",nbins,-180,180);
60
61 halpha3015->SetXTitle("alpha-alpha0 (deg)");
62 halpha3015->SetYTitle("Entries");
63
64 halpha3015->SetLineColor(1);
65 halpha4025->SetLineColor(2);
66 halpha5035->SetLineColor(3);
67 halpha6045->SetLineColor(4);
68
69 ton->Draw("alpha3015-alpha0>>Alpha3015","alpha3015!=0");
70 ton->Draw("alpha4025-alpha0>>Alpha4025","alpha4025!=0");
71 ton->Draw("alpha5035-alpha0>>Alpha5035","alpha5035!=0");
72 ton->Draw("alpha6045-alpha0>>Alpha6045","alpha6045!=0");
73
74 // plot
75 TCanvas *c1 = new TCanvas("c1","c1",800,500);
76 c1->cd(1);
77
78 halpha3015->Draw("HIST");
79 halpha4025->Draw("HISTSAMES");
80 halpha5035->Draw("HISTSAMES");
81 halpha6045->Draw("HISTSAMES");
82
83 gPad->Update();
84 TPaveStats* pvstat4025 = (TPaveStats*) halpha4025->GetListOfFunctions()->FindObject("stats");
85 TPaveStats* pvstat5035 = (TPaveStats*) halpha5035->GetListOfFunctions()->FindObject("stats");
86 TPaveStats* pvstat6045 = (TPaveStats*) halpha6045->GetListOfFunctions()->FindObject("stats");
87
88 Float_t shifty = pvstat4025->GetY2NDC()-pvstat4025->GetY1NDC();
89
90 pvstat4025->SetY1NDC(pvstat4025->GetY1NDC()-shifty);
91 pvstat4025->SetY2NDC(pvstat4025->GetY2NDC()-shifty);
92 pvstat5035->SetY1NDC(pvstat5035->GetY1NDC()-2*shifty);
93 pvstat5035->SetY2NDC(pvstat5035->GetY2NDC()-2*shifty);
94 pvstat6045->SetY1NDC(pvstat6045->GetY1NDC()-3*shifty);
95 pvstat6045->SetY2NDC(pvstat6045->GetY2NDC()-3*shifty);
96
97 gPad->Modified();
98 gPad->Update();
99
100 c1->Print("cleanComp.ps");
101}
102
103/*******************************************************************************/
104void plotAlphaWidthVsEnergy(TChain* ton)
105{
106 const Int_t nbins = 8;
107 Float_t bound[nbins+1] ={0,30,50,80,150,300,600,1200,3000}; // boundary of considered size bins
108
109 // define the arrays that will feed the graphs
110 // value of rms of the alpha-alpha0 distribution, and errors
111 Float_t alpha0[nbins];
112 Float_t alpha3015[nbins];
113 Float_t alpha4025[nbins];
114 Float_t alpha5035[nbins];
115 Float_t alpha6045[nbins];
116
117 Float_t ealpha0[nbins];
118 Float_t ealpha3015[nbins];
119 Float_t ealpha4025[nbins];
120 Float_t ealpha5035[nbins];
121 Float_t ealpha6045[nbins];
122
123 // x axis and error
124 Float_t x[nbins];
125 Float_t ex[nbins];
126
127 TH1F* hp[nbins];
128
129 // define the function that will be used to fit the histograms
130 TF1 f2("f2","[0]+[1]/[3]*exp(-(x-[2])*(x-[2])/(2*[3]*[3]))",-180.,180.);
131 f2.SetParameters(0,100,0,02);
132 f2.FixParameter(0,0);
133 f2.SetLineWidth(1);
134
135 // fill the arrays
136 for(Int_t i=0;i<nbins;i++)
137 {
138 // cout << "Bin " << i << endl;
139 TH1F dummy0("dummy0" ,"dummy0",100,-90,90);
140 TH1F dummy3015("dummy3015" ,"dummy3015",100,-180,180);
141 TH1F dummy4025("dummy4025" ,"dummy4025",100,-180,180);
142 TH1F dummy5035("dummy5035" ,"dummy5035",100,-180,180);
143 TH1F dummy6045("dummy6045" ,"dummy6045",100,-180,180);
144
145 TH1F xdummy("xdummy" ,"xdummy",100,bound[i],bound[i+1]);
146
147 // define the cuts for each of the cleaning levels
148 char cut[100];
149 sprintf(cut,"energy>%f && energy<%f && alpha0!=0",bound[i],bound[i+1]);
150 char cut3015[200];
151 sprintf(cut3015,"%s && alpha3015!=0",cut);
152 char cut4025[200];
153 sprintf(cut4025,"%s && alpha4025!=0",cut);
154 char cut5035[200];
155 sprintf(cut5035,"%s && alpha5035!=0",cut);
156 char cut6045[200];
157 sprintf(cut6045,"%s && alpha6045!=0",cut);
158
159 TF1 f1("f1","[0]+[1]/[3]*exp(-(x-[2])*(x-[2])/(2*[3]*[3]))",-180.,180.);
160 f1.SetParameters(0,100,0,02);
161 f1.FixParameter(0,0);
162
163
164 // fill the histograms to extract later the information from rms, entries, means, etc
165 ton->Draw("alpha0>>dummy0",cut);
166 ton->Draw("alpha3015-alpha0>>dummy3015",cut3015);
167 ton->Draw("alpha4025-alpha0>>dummy4025",cut4025);
168 ton->Draw("alpha5035-alpha0>>dummy5035",cut5035);
169 ton->Draw("alpha6045-alpha0>>dummy6045",cut6045);
170
171 ton->Draw("energy>>xdummy",cut);
172
173 // fit to a gaussian and get the sigma to be plotted later
174 dummy0.Fit("f1");
175 alpha0[i] =f1.GetParameter(3);
176 dummy3015.Fit("f1");
177 alpha3015[i]=f1.GetParameter(3);
178 dummy4025.Fit("f1");
179 alpha4025[i]=f1.GetParameter(3);
180 dummy5035.Fit("f1");
181 alpha5035[i]=f1.GetParameter(3);
182 dummy6045.Fit("f1");
183 alpha6045[i]=f1.GetParameter(3);
184
185 Float_t nevt0 = dummy0.GetEntries();
186 Float_t nevt3015 = dummy3015.GetEntries();
187 Float_t nevt4025 = dummy4025.GetEntries();
188 Float_t nevt5035 = dummy5035.GetEntries();
189 Float_t nevt6045 = dummy6045.GetEntries();
190
191 // since gaussian is not so good aproximation, compute error like follows
192 ealpha0[i] = alpha0[i]/TMath::Sqrt(2*nevt0);
193 ealpha3015[i]= alpha3015[i]/TMath::Sqrt(2*nevt3015);
194 ealpha4025[i]= alpha4025[i]/TMath::Sqrt(2*nevt4025);
195 ealpha5035[i]= alpha5035[i]/TMath::Sqrt(2*nevt5035);
196 ealpha6045[i]= alpha6045[i]/TMath::Sqrt(2*nevt6045);
197
198 // copy some histos to be plotted later
199 Char_t nom[10];
200 sprintf(nom,"hp[%02d]",i);
201 hp[i] = new TH1F(nom,nom,100,-180,180);
202 hp[i]->GetXaxis()->SetTitle("alpha (deg)");
203 hp[i]->Add(&dummy6045);
204 hp[i]->Fit("f2");
205
206 // cout << xdummy.GetMean() << " " << xdummy.GetRMS()<< " " << xdummy.GetEntries() << endl;
207 x[i] = xdummy.GetMean();
208 // Float_t rms = xdummy.GetRMS();
209 ex[i] = 0;//rms*rms/xdummy.GetEntries();
210 }
211
212 //#if 0
213 TCanvas* cdum = new TCanvas("cdum","cdum",600,900);
214 cdum->Divide(2,nbins/2+1);
215 for(Int_t i=0;i<nbins;i++)
216 {
217 cdum->cd(i+1);
218 hp[i]->Draw();
219 gPad->Modified();
220 gPad->Update();
221 }
222 //#endif
223
224 // create the graphs
225 TGraphErrors* gr0 = new TGraphErrors(nbins,x,alpha0,ex,ealpha0);
226 TGraphErrors* gr3015 = new TGraphErrors(nbins,x,alpha3015,ex,ealpha3015);
227 TGraphErrors* gr4025 = new TGraphErrors(nbins,x,alpha4025,ex,ealpha4025);
228 TGraphErrors* gr5035 = new TGraphErrors(nbins,x,alpha5035,ex,ealpha5035);
229 TGraphErrors* gr6045 = new TGraphErrors(nbins,x,alpha6045,ex,ealpha6045);
230
231 // marker/line colors
232 gr0->SetMarkerStyle(20);
233 gr3015->SetMarkerStyle(24);
234 gr4025->SetMarkerStyle(21);
235 gr5035->SetMarkerStyle(25);
236 gr6045->SetMarkerStyle(22);
237
238 gr0->SetMarkerColor(1);
239 gr3015->SetMarkerColor(2);
240 gr4025->SetMarkerColor(3);
241 gr5035->SetMarkerColor(4);
242 gr6045->SetMarkerColor(6);
243
244 gr0->SetLineColor(1);
245 gr3015->SetLineColor(2);
246 gr4025->SetLineColor(3);
247 gr5035->SetLineColor(4);
248 gr6045->SetLineColor(6);
249
250 //////////////////////
251 // plot alpha width //
252 //////////////////////
253 TCanvas *c1 = new TCanvas("c1","c1",800,500);
254 c1->cd(1);
255 gPad->SetGridx();
256 gPad->SetGridy();
257 gPad->SetLogx();
258
259 TGraphErrors* grfirst = gr4025;
260 grfirst->Draw("APL");
261 gr0->Draw("PL");
262 gr3015->Draw("PL");
263 gr4025->Draw("PL");
264 gr5035->Draw("PL");
265 gr6045->Draw("PL");
266
267 grfirst->GetHistogram()->SetTitle("Cleaning comparison (alpha)");
268 grfirst->GetHistogram()->GetXaxis()->SetTitle("E_{MC} [GeV]");
269 grfirst->GetHistogram()->GetYaxis()->SetTitle("#sigma (alpha-alpha_{0}) [deg]");
270 grfirst->GetHistogram()->SetMaximum(80);
271
272 TPaveLabel* label = new TPaveLabel(0.25,0.75,0.6,0.85,"MC Gammas - zbins 0-3","NDC");
273 label->SetBorderSize(0);
274 label->SetFillColor(0);
275 label->Draw();
276
277 // legend
278 TLegend* leg = new TLegend(.6,.5,.8,.8);
279 leg->SetHeader("");
280 leg->SetFillColor(0);
281 leg->SetLineColor(0);
282 leg->SetBorderSize(0);
283 leg->AddEntry(gr0,"No noise","p");
284 leg->AddEntry(gr3015,"(3.0,1.5)","p");
285 leg->AddEntry(gr4025,"(4.0,2.5)","p");
286 leg->AddEntry(gr5035,"(5.0,3.5)","p");
287 leg->AddEntry(gr6045,"(6.0,4.5)","p");
288 leg->Draw();
289
290 gPad->Modified();
291 gPad->Update();
292 // print to file
293 c1->Print("cleanAlphaWidthVSEnergy.ps");
294}
295
296/*******************************************************************************/
297void plotEfficiencies(TChain* ton)
298{
299 const Int_t nbins = 8;
300 Float_t bound[nbins+1] ={0,30,50,80,150,300,600,1200,3000}; // boundary of considered size bins
301
302 // define the arrays that will feed the graphs
303 // gamma detection efficiency, and errors
304 Float_t effNOCL[nbins];
305 Float_t eff3015[nbins];
306 Float_t eff4025[nbins];
307 Float_t eff5035[nbins];
308 Float_t eff6045[nbins];
309
310 Float_t eeffNOCL[nbins];
311 Float_t eeff3015[nbins];
312 Float_t eeff4025[nbins];
313 Float_t eeff5035[nbins];
314 Float_t eeff6045[nbins];
315
316 // x axis and error
317 Float_t x[nbins];
318 Float_t ex[nbins];
319
320 // fill the arrays
321 for(Int_t i=0;i<nbins;i++)
322 {
323 // cout << "Bin " << i << endl;
324 TH1F dummyNOCL("dummyNOCL" ,"dummyNOCL",100,-90,90);
325 TH1F dummy3015("dummy3015" ,"dummy3015",100,-90,90);
326 TH1F dummy4025("dummy4025" ,"dummy4025",100,-90,90);
327 TH1F dummy5035("dummy5035" ,"dummy5035",100,-90,90);
328 TH1F dummy6045("dummy6045" ,"dummy6045",100,-90,90);
329
330 TH1F xdummy("xdummy" ,"xdummy",100,bound[i],bound[i+1]);
331
332 // define the cuts for each of the cleaning levels
333 char cut[100];
334 sprintf(cut,"energy>%f && energy<%f && alpha0!=0",bound[i],bound[i+1]);
335 char cutNOCL[200];
336 sprintf(cutNOCL,"%s && alphaNOCL!=0",cut);
337 char cut3015[200];
338 sprintf(cut3015,"%s && alpha3015!=0",cut);
339 char cut4025[200];
340 sprintf(cut4025,"%s && alpha4025!=0",cut);
341 char cut5035[200];
342 sprintf(cut5035,"%s && alpha5035!=0",cut);
343 char cut6045[200];
344 sprintf(cut6045,"%s && alpha6045!=0",cut);
345
346
347 // fill the histograms to extract later the information from rms, entries, means, etc
348 ton->Draw("alpha0>>dummyNOCL",cutNOCL);
349 ton->Draw("alpha0>>dummy3015",cut3015);
350 ton->Draw("alpha0>>dummy4025",cut4025);
351 ton->Draw("alpha0>>dummy5035",cut5035);
352 ton->Draw("alpha0>>dummy6045",cut6045);
353
354 ton->Draw("energy>>xdummy",cut);
355
356 Float_t nevt0 = dummyNOCL.GetEntries();
357 Float_t nevt3015 = dummy3015.GetEntries();
358 Float_t nevt4025 = dummy4025.GetEntries();
359 Float_t nevt5035 = dummy5035.GetEntries();
360 Float_t nevt6045 = dummy6045.GetEntries();
361
362 // fill the array of efficiencies
363 eff3015[i] = nevt3015/nevt0*100;
364 eff4025[i] = nevt4025/nevt0*100;
365 eff5035[i] = nevt5035/nevt0*100;
366 eff6045[i] = nevt6045/nevt0*100;
367
368 eeff3015[i] = TMath::Sqrt(nevt3015)/nevt0*100;
369 eeff4025[i] = TMath::Sqrt(nevt4025)/nevt0*100;
370 eeff5035[i] = TMath::Sqrt(nevt5035)/nevt0*100;
371 eeff6045[i] = TMath::Sqrt(nevt6045)/nevt0*100;
372
373 // cout << xdummy.GetMean() << " " << xdummy.GetRMS()<< " " << xdummy.GetEntries() << endl;
374 x[i] = xdummy.GetMean();
375 // Float_t rms = xdummy.GetRMS();
376 ex[i] = 0;//rms*rms/xdummy.GetEntries();
377 }
378
379 // create the graphs
380 TGraphErrors* efgr3015 = new TGraphErrors(nbins,x,eff3015,ex,eeff3015);
381 TGraphErrors* efgr4025 = new TGraphErrors(nbins,x,eff4025,ex,eeff4025);
382 TGraphErrors* efgr5035 = new TGraphErrors(nbins,x,eff5035,ex,eeff5035);
383 TGraphErrors* efgr6045 = new TGraphErrors(nbins,x,eff6045,ex,eeff6045);
384
385 // marker/line colors
386 efgr3015->SetMarkerStyle(24);
387 efgr4025->SetMarkerStyle(21);
388 efgr5035->SetMarkerStyle(25);
389 efgr6045->SetMarkerStyle(22);
390
391 efgr3015->SetMarkerColor(2);
392 efgr4025->SetMarkerColor(3);
393 efgr5035->SetMarkerColor(4);
394 efgr6045->SetMarkerColor(6);
395
396 efgr3015->SetLineColor(2);
397 efgr4025->SetLineColor(3);
398 efgr5035->SetLineColor(4);
399 efgr6045->SetLineColor(6);
400
401 ///////////////////////
402 // plot efficiencies //
403 ///////////////////////
404 TCanvas *c1 = new TCanvas("c1","c1",800,500);
405 c1->cd(1);
406 gPad->SetGridx();
407 gPad->SetGridy();
408 gPad->SetLogx();
409
410 TGraphErrors* efgrfirst = efgr6045;
411 efgrfirst->Draw("APL");
412 efgr3015->Draw("PL");
413 efgr4025->Draw("PL");
414 efgr5035->Draw("PL");
415 efgr6045->Draw("PL");
416
417 efgrfirst->GetHistogram()->SetTitle("Cleaning comparison (efficiency)");
418 efgrfirst->GetHistogram()->GetXaxis()->SetTitle("E_{MC} [GeV]");
419 efgrfirst->GetHistogram()->GetYaxis()->SetTitle("Detection efficiency (%)");
420
421 // label
422 TPaveLabel* eflabel = new TPaveLabel(0.5,0.49,0.85,0.59,"MC Gammas - zbins 0-3","NDC");
423 eflabel->SetBorderSize(0);
424 eflabel->SetFillColor(0);
425 eflabel->Draw();
426
427 // legend
428 TLegend* efleg = new TLegend(.6,.2,.8,.5);
429 efleg->SetHeader("");
430 efleg->SetFillColor(0);
431 efleg->SetLineColor(0);
432 efleg->SetBorderSize(0);
433 efleg->AddEntry(efgr3015,"(3.0,1.5)","p");
434 efleg->AddEntry(efgr4025,"(4.0,2.5)","p");
435 efleg->AddEntry(efgr5035,"(5.0,3.5)","p");
436 efleg->AddEntry(efgr6045,"(6.0,4.5)","p");
437 efleg->Draw();
438
439 gPad->Modified();
440 gPad->Update();
441
442 // print to file
443 c1->Print("cleanEfficiency.ps");
444}
445
Note: See TracBrowser for help on using the repository browser.