source: trunk/MagicSoft/Mars/mtemp/meth/macros/analysis_read.C@ 5138

Last change on this file since 5138 was 4354, checked in by commichau, 21 years ago
*** empty log message ***
File size: 16.3 KB
Line 
1
2/////////////////////////////////////////////////////////////////////////////
3// //
4// analysis_read.C //
5// //
6// Author(s): S.C. Commichau, 3/2004 //
7// //
8/////////////////////////////////////////////////////////////////////////////
9
10
11// This macro reads the output of analysis.C and plots some of the
12// image parameters. It also produces an ALPHA plot and calculates
13// the significance
14
15void analysis_read(){
16
17 gROOT->Reset();
18
19 // Some parameters for the ALPHA-plot
20 const Int_t inibin = 9;
21 const Int_t nbins = 18;
22
23 // ON data chain
24 TChain* CON = new TChain("Parameters");
25
26 // OFF data chain
27 TChain* COFF = new TChain("Parameters");
28
29 // Add the root file(s) to the chain and treat it like a tree
30 CON->Add("~/Mars/on1.root");
31 COFF->Add("~/Mars/off1.root");
32
33 //CON->Add("~/no_back/hillas_on_all23.root");
34 //COFF->Add("~/no_back/hillas_off_all23.root");
35
36 TString title = "Mrk421 - 04/23/2004";
37
38 // Set some global style options
39 //gStyle->SetOptDate(10);
40 //gStyle->SetDateX(.1);
41 //gStyle->SetDateY(.905);
42 gStyle->SetOptFit(0);
43 gStyle->SetOptStat(110011);
44 gStyle->SetFrameBorderMode(0);
45 //gStyle->SetPalette(1);
46 gStyle->SetPalette(1,0);
47 gStyle->SetFrameBorderSize(0);
48 gStyle->SetCanvasColor(0);
49 gStyle->SetFrameFillColor(0);
50 gStyle->SetTitleFillColor(0);
51 gStyle->SetTitleBorderSize(0);
52 gStyle->SetStatColor(0);
53 gStyle->SetStatBorderSize(1);
54 // Put the statistics box into the right corner!
55 gStyle->SetStatX(.9);
56 gStyle->SetStatY(.9);
57 //gStyle->SetStatTextColor(2);
58
59 // Set alias for each quantity to simplify the handling
60 CON->SetAlias("length","MHillas.fLength*0.6/189.0");
61 CON->SetAlias("width","MHillas.fWidth*0.6/189.0");
62 CON->SetAlias("meanx","MHillas.fMeanX*0.6/189.0");
63 CON->SetAlias("meany","MHillas.fMeanY*0.6/189.0");
64 CON->SetAlias("size","MHillas.fSize");
65 CON->SetAlias("dist","MHillasSrc.fDist*0.6/189.0");
66 CON->SetAlias("alpha","MHillasSrc.fAlpha");
67 CON->SetAlias("delta","abs(MHillas.fDelta)*180/3.1415");
68 //CON->SetAlias("delta","MDCA.fDelta1*180/3.1415");
69 //CON->SetAlias("dca","MDCA.fDCA*0.6/189.0");
70
71 COFF->SetAlias("length","MHillas.fLength*0.6/189.0");
72 COFF->SetAlias("width","MHillas.fWidth*0.6/189.0");
73 COFF->SetAlias("meanx","MHillas.fMeanX*0.6/189.0");
74 COFF->SetAlias("meany","MHillas.fMeanY*0.6/189.0");
75 COFF->SetAlias("size","MHillas.fSize");
76 COFF->SetAlias("dist","MHillasSrc.fDist*0.6/189.0");
77 COFF->SetAlias("alpha","MHillasSrc.fAlpha");
78 COFF->SetAlias("delta","abs(MHillas.fDelta)*180/3.1415");
79
80
81 // Cuts...
82 // Cut for ALPHA
83 TString cut1 = "(dist>0.2) && (dist<.7) && (width>0.04) && (width<0.14) && (length>0.14) && (length<0.26) && (size>500)";
84
85 TString cut2 = "";//size>1000 && size<2000";
86
87 TString cut3 = "";
88
89 TString cut4 = "";
90
91 TString cut5 = "";
92
93
94 // And now the draw area
95
96 //************ Width vs Length ************
97
98 TCanvas* CWL = new TCanvas("CWL","Canvas",0,0,300,300);
99 CWL->cd();
100 CWL->SetBorderMode(0);
101 gPad->SetBorderMode(0);
102 gStyle->SetOptStat(0);
103
104 TH2F* hWL = new TH2F("hWL","",200,0,0.7,200,0,1);
105
106 // Leaf, Cuts, Draw option, use >>+ to append data to an ex. histogram
107 CON->Draw("length:width>>hWL",cut2);
108
109 hWL->SetTitle(title);
110 hWL->GetXaxis()->SetTitle("WIDTH [#circ]");
111 hWL->GetYaxis()->SetTitle("LENGTH [#circ]");
112 hWL->GetYaxis()->SetTitleOffset(1.2);
113 hWL->SetMarkerColor(2);
114
115 hWL->SetMaximum(150);
116 hWL->Draw("COLZ");//CONT
117
118
119 //************ Width and Length ************
120
121 // Create Canvas and set Canvas some options
122 TCanvas* Cwl = new TCanvas("Cwl","Canvas",325,0,300,300);
123 Cwl->cd();
124 Cwl->SetBorderMode(0);
125 gPad->SetBorderMode(0);
126 gStyle->SetOptStat(110011);
127
128 TH1F* hLength = new TH1F("hLength","",100,0,1.1);
129
130 CON->Draw("length>>hLength",cut5,"");
131
132 hLength->SetTitle(title);
133 hLength->GetXaxis()->SetTitle("LENGTH and WIDTH [#circ]");
134 hLength->GetYaxis()->SetTitle("Entries");
135 hLength->GetYaxis()->SetTitleOffset(1.6);
136 hLength->SetFillColor(50);
137 hLength->SetLineColor(50);
138 hLength->SetLineWidth(2);
139 hLength->SetFillStyle(3005);
140 //hLength->SetMaximum(55000);
141 hLength->SetName("LENGTH");
142
143 TH1F* hWidth = new TH1F("hWidth","",100,0,1.1);
144
145 CON->Draw("width>>hWidth",cut5,"");
146
147 gStyle->SetStatTextColor(4);
148 hWidth->SetFillColor(4);
149 hWidth->SetLineColor(4);
150 hWidth->SetFillStyle(3004);
151 hWidth->SetLineWidth(2);
152 hWidth->SetName("WIDTH");
153
154 TLegend* leg = new TLegend(.59,.75,.8,.84);
155 leg->AddEntry(hLength,"LENGTH","F");
156 leg->AddEntry(hWidth,"WIDTH","F");
157 leg->SetFillColor(0);
158 leg->SetLineColor(1);
159 leg->SetBorderSize(0);
160
161 hLength->Draw();//DrawNormalized or DrawCopy are possible options
162 hWidth->Draw("sames");
163 //leg->Draw("same");
164
165 // Beautify the statistic boxes!
166 gPad->Update();
167 TPaveStats* pavstat = (TPaveStats*) hLength->GetListOfFunctions()->FindObject("stats");
168 if(pavstat)
169 {
170 Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
171 pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
172 pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
173 pavstat->SetTextColor(50);
174 }
175
176 gPad->Modified();
177 gPad->Update();
178
179
180 //************ Dist ************
181
182 TCanvas* CDist = new TCanvas("CDist","Canvas",650,0,300,300);
183 CDist->cd();
184 CDist->SetBorderMode(0);
185 gStyle->SetOptStat(110011);
186 //CDist->SetLogy();
187
188 TH1F* hDist = new TH1F("hDist","",90,0,1.8);
189
190 CON->Draw("dist>>hDist");
191
192 hDist->GetXaxis()->SetTitle("DIST [#circ]");
193 hDist->GetYaxis()->SetTitle("Entries");
194 hDist->GetYaxis()->SetTitleOffset(1.6);
195 //hDist->GetXaxis()->SetTitleOffset(1.3);
196 hDist->SetLineColor(4);
197 hDist->SetMarkerColor(4);
198 hDist->SetFillColor(4);
199 hDist->SetFillStyle(3004);
200 hDist->SetLineWidth(2);
201 hDist->SetTitle(title);
202 hDist->SetName("DIST");
203
204 hDist->Draw();
205
206
207 //************ MeanY vs MeanX ************
208
209 TCanvas* CXY = new TCanvas("CXY","Canvas",975,0,300,300);
210 CXY->cd();
211 CXY->SetBorderMode(0);
212 gStyle->SetOptStat(0);
213
214 TH2F* hXY = new TH2F("hXY","",500,-1.5,1.5,500,-1.5,1.5);
215
216 CON->Draw("meany:meanx>>hXY",cut3,"");
217
218 TString titlexy = TString(title);
219
220 hXY->SetTitle(titlexy);
221 hXY->GetXaxis()->SetTitle("MeanX [#circ]");
222 hXY->GetYaxis()->SetTitle("MeanY [#circ]");
223 hXY->GetYaxis()->SetTitleOffset(1.2);
224 hXY->SetMaximum(15);
225
226 hXY->Draw("COLZ");
227
228
229 //************ DCA ************
230 /*
231 TCanvas* CDCA = new TCanvas("CDCA","Canvas",10,400,300,300);
232 CDCA->cd();
233 CDCA->SetBorderMode(0);
234 gStyle->SetOptStat(110011);
235
236 TH1F* hDCA = new TH1F("hDCA","",100,-1.5,1.5);
237
238 CON->Draw("dca>>hDCA");
239
240 hDCA->GetXaxis()->SetTitle("DCA [#circ]");
241 hDCA->GetYaxis()->SetTitle("Entries");
242 hDCA->GetYaxis()->SetTitleOffset(2);
243 hDCA->SetLineColor(4);
244 hDCA->SetFillColor(4);
245 hDCA->SetFillStyle(3004);
246 hDCA->SetTitle(title);
247 hDCA->SetName("DCA");
248
249 hDCA->Draw();
250 */
251
252
253 //************ Delta ************
254
255 TCanvas* CDelta = new TCanvas("CDelta","Canvas",0,335,300,300);
256 CDelta->cd();
257 CDelta->SetBorderMode(0);
258 gStyle->SetOptStat(110011);
259
260 TH1F* hDelta = new TH1F("hDelta","",100,0,90);
261
262 CON->Draw("delta>>hDelta",cut3);
263
264 hDelta->GetXaxis()->SetTitle("#delta [#circ]");
265 hDelta->GetYaxis()->SetTitle("Entries");
266 hDelta->GetYaxis()->SetTitleOffset(1.4);
267 hDelta->SetLineColor(4);
268 hDelta->SetFillColor(4);
269 hDelta->SetFillStyle(3004);
270 hDelta->SetLineWidth(2);
271 hDelta->SetTitle(title);
272 hDelta->SetName("Delta");
273
274 hDelta->Draw();
275
276
277 //************ Size ************
278
279 TCanvas* CSize = new TCanvas("CSize","Canvas",325,335,300,300);
280 CSize->cd();
281 CSize->SetBorderMode(0);
282 gStyle->SetOptStat(110011);
283
284 TH1F* hSize = new TH1F("hSize","",100,0.1,5000);
285
286 CON->Draw("size>>hSize");
287
288 hSize->GetXaxis()->SetTitle("Size [Photons]");
289 hSize->GetYaxis()->SetTitle("Entries");
290 hSize->GetYaxis()->SetTitleOffset(1.6);
291 hSize->GetXaxis()->SetTitleOffset(1.2);
292 hSize->SetLineColor(4);
293 hSize->SetFillColor(4);
294 hSize->SetLineWidth(2);
295 hSize->SetFillStyle(3004);
296 hSize->SetTitle(title);
297 hSize->SetName("SIZE");
298
299 hSize->Draw();
300
301
302 //************ Alpha ************
303
304 TCanvas* CAlpha = new TCanvas("CAlpha","Canvas",650,335,300,300);
305 CAlpha->cd();
306 CAlpha->SetBorderMode(0);
307 CAlpha->SetGridx();
308 CAlpha->SetGridy();
309 gStyle->SetOptStat(11);
310
311 TH1F *hAlpha = (TH1F*)new TH1F("hAlpha","Alpha ON",nbins,0,90);
312 TH1F *hAlphaoff = (TH1F*)new TH1F("hAlphaoff","Alpha OFF",nbins,0,90);
313
314 CON->Draw("abs(alpha)>>hAlpha",cut1);
315
316 COFF->Draw("abs(alpha)>>hAlphaoff",cut1);
317
318
319 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
320 // Normalize ON and OFF data
321 cout << "******************************************************" << endl;
322 cout << " Bin error before scaling: " << hAlphaoff->GetBinError(2) << endl;
323
324 // If TH1::Sumw2 has been called before filling, the sum of squares of
325 // weights is also stored - I do not understand this!!
326 hAlphaoff->Sumw2();
327
328 Float_t level = 0, leveloff = 0;
329
330 // Calculate the mean of entries beside the signal region
331 for(Int_t ibin = inibin; ibin<=nbins; ibin++)
332 level += hAlpha->GetBinContent(ibin);
333 level /= (nbins-inibin+1);
334
335 cout << " Mean beside signal: " << level << endl;
336
337 for(Int_t ibin = inibin; ibin<=nbins; ibin++)
338 leveloff += hAlphaoff->GetBinContent(ibin);
339 leveloff /= (nbins-inibin+1);
340
341 // Compare the mean values of both histograms
342 if(leveloff>0)
343 const Float_t norm = level/leveloff;//*nbins/hAlphaoff->GetEntries();
344 cout << " Normalizing by factor " << norm <<endl;
345 hAlphaoff->Scale(norm);
346
347 cout << " Bin error after scaling: " << hAlphaoff->GetBinError(2) << endl;
348
349 cout << " Cut on (ALPHA): " << cut1 << endl;
350 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
351
352 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
353 // Significance calculation
354 Float_t significance = 0;
355 const Int_t signbins = 3;
356
357 Int_t Non = 0;
358 Int_t Noff = 0;
359
360 cout << " Significance Bins: " << signbins << endl;
361
362 // Determine number of excess events and off events in the same region
363 for(Int_t ibin = 1; ibin<=signbins;ibin++)
364 {
365 Non +=hAlpha->GetBinContent(ibin);
366 Noff+=hAlphaoff->GetBinContent(ibin);
367 }
368
369 cout << " Non = " << Non << " Noff = " << Noff << endl;
370
371 significance = Significance(Non, Noff, 1.0, 0); //norm
372
373 cout << " Significance: " << significance << endl;
374 cout << "******************************************************" << endl;
375 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
376
377 // Set some draw options and draw the histograms
378 hAlphaoff->GetXaxis()->SetTitle("ALPHA [#circ]");
379 hAlphaoff->GetYaxis()->SetTitle("Entries");
380 hAlphaoff->GetYaxis()->SetTitleOffset(1.4);
381 hAlphaoff->GetXaxis()->SetTitleOffset(1.1);
382
383 char legentry[128];
384 sprintf(legentry," (Sign.: %2.3f #sigma, Norm.: %2.3f, |ALPHA| #leq %d#circ)", significance, norm, 90/nbins*signbins);
385
386 hAlphaoff->SetTitle(title+legentry);
387 hAlphaoff->SetName("OFF source");
388 hAlphaoff->SetLineColor(4);
389 hAlphaoff->SetLineWidth(2);
390 hAlphaoff->SetMarkerColor(8);
391 hAlphaoff->SetMarkerStyle(21);
392 hAlphaoff->SetMarkerSize(0);
393 hAlphaoff->SetFillColor(18);
394 hAlphaoff->SetMaximum(hAlpha->GetMaximum()+200);
395
396 hAlpha->SetLineColor(1);
397 hAlpha->SetLineWidth(2);
398 hAlpha->SetMarkerColor(50);
399 hAlpha->SetMarkerStyle(20);
400 hAlpha->SetMarkerSize(0.9);
401 hAlpha->SetName("ON source");
402 //hAlpha->SetMaximum(hAlpha->GetMaximum());
403
404 hAlphaoff->Draw("eH");
405 hAlpha->Draw("e1sames");
406
407 // Beautify the statistic boxes!
408 gPad->Update();
409 TPaveStats* pavstat = (TPaveStats*) hAlphaoff->GetListOfFunctions()->FindObject("stats");
410 if(pavstat)
411 {
412 Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
413 pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
414 pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
415 pavstat->SetTextColor(4);
416 }
417
418 TPaveStats* pavstaton = (TPaveStats*) hAlpha->GetListOfFunctions()->FindObject("stats");
419 if(pavstaton)
420 {
421 pavstaton->SetTextColor(50);
422 }
423
424 gPad->Modified();
425 gPad->Update();
426
427
428 //************ Length/Size ************
429
430 TCanvas* CLSiz = new TCanvas("CLSiz","Canvas",975,335,300,300);
431 CLSiz->cd();
432 CLSiz->SetBorderMode(0);
433 gStyle->SetOptStat(110011);
434
435 TH1F* hLenSiz = new TH1F("hLenSiz","",100,0,0.005);
436
437 CON->Draw("length/size>>hLenSiz",cut2,"");
438
439 hLenSiz->GetXaxis()->SetTitle("LENGTH/SIZE [#circ/Photons]");
440 hLenSiz->GetYaxis()->SetTitle("Entries");
441 hLenSiz->GetYaxis()->SetTitleOffset(1.5);
442 hLenSiz->GetXaxis()->SetTitleOffset(1.1);
443 hLenSiz->SetLineColor(4);
444 hLenSiz->SetFillColor(4);
445 hLenSiz->SetLineWidth(2);
446 hLenSiz->SetFillStyle(3004);
447 hLenSiz->SetNdivisions(505);
448
449 hLenSiz->SetTitle(title);
450 hLenSiz->SetName("LENGTH/SIZE");
451
452 hLenSiz->Draw();
453
454
455 // Uncomment the following lines for MC only!
456
457 //************ Energy ************
458 /*
459
460 TH1F *hEnergy = (TH1F*)gDirectory->Get("hEnergy");
461
462 TCanvas* CEnergy = new TCanvas("CEnergy","Canvas",0,670,300,300);
463 CEnergy->cd();
464 CEnergy->SetBorderMode(0);
465 //CEnergy->SetLogy();
466 gStyle->SetOptStat(110011);
467
468 TH1F* hEnergy = new TH1F("hEnergy","",100,0,1000);
469
470 CON->Draw("energy>>hEnergy","","");
471
472 hEnergy->GetXaxis()->SetTitle("Energy [GeV]");
473 hEnergy->GetYaxis()->SetTitle("Entries");
474 hEnergy->GetYaxis()->SetTitleOffset(1.2);
475 hEnergy->GetXaxis()->SetTitleOffset(1.1);
476 hEnergy->SetLineColor(4);
477 hEnergy->SetFillColor(4);
478 hEnergy->SetLineWidth(2);
479 hEnergy->SetFillStyle(3004);
480 hEnergy->SetNdivisions(505);
481
482 hEnergy->SetTitle(title);
483 hEnergy->SetName("Energy");
484
485 hEnergy->Draw();
486 */
487
488
489 //************ MeanX and MeanY ************
490
491 // Create Canvas and set Canvas some options
492 TCanvas* Cxy = new TCanvas("Cxy","Canvas",325,670,300,300);
493 Cxy->cd();
494 Cxy->SetBorderMode(0);
495 gPad->SetBorderMode(0);
496 gStyle->SetOptStat(111);
497
498 TH1F* hMeanX = new TH1F("hMeanX","",100,-1.5,1.5);
499
500 CON->Draw("meanx>>hMeanX",cut3,"");
501
502 TH1F* hMeanY = new TH1F("hMeanY","",100,-1.5,1.5);
503
504 CON->Draw("meany>>hMeanY",cut3,"");
505
506 hMeanX->SetTitle(title);
507 hMeanX->GetXaxis()->SetTitle("MeanX and MeanY [#circ]");
508 hMeanX->GetYaxis()->SetTitle("Entries");
509 hMeanX->GetYaxis()->SetTitleOffset(1.6);
510 hMeanX->SetFillColor(50);
511 hMeanX->SetLineColor(50);
512 hMeanX->SetLineWidth(2);
513 hMeanX->SetFillStyle(3005);
514 //hMeanX->SetMaximum(16000);
515 //hMeanX->SetMaximum(6500);
516 hMeanX->SetName("MeanX");
517
518 gStyle->SetStatTextColor(4);
519 hMeanY->SetFillColor(4);
520 hMeanY->SetLineColor(4);
521 hMeanY->SetFillStyle(3004);
522 hMeanY->SetLineWidth(2);
523 hMeanY->SetName("MeanY");
524
525 hMeanX->Draw();//DrawNormalized or DrawCopy are possible options
526 hMeanY->Draw("sames");//esames
527
528 // Beautify the statistic boxes!
529 gPad->Update();
530 TPaveStats* pavstat = (TPaveStats*) hMeanX->GetListOfFunctions()->FindObject("stats");
531 if(pavstat)
532 {
533 Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
534 pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
535 pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
536 pavstat->SetTextColor(50);
537 }
538
539 gPad->Modified();
540 gPad->Update();
541
542 //delete hWL;
543 //delete hXY;
544}
545
546
547Double_t Significance(Float_t Non, Float_t Noff, Float_t alpha, Int_t type = 0)
548{
549 // Calculate significance after Li and Ma
550
551 // Formula Nr. 17
552 if(type == 0)
553 return TMath::Sqrt(2*(Non*TMath::Log(((1+alpha)/alpha)*(Non/(Non+Noff))) + Noff*TMath::Log((1+alpha)*Noff/(Non+Noff))));
554
555 // Formula Nr. 9
556 if(type == 1)
557 return (Non-alpha*Noff)/TMath::Sqrt(alpha*(Non+Noff));
558
559}
560
561
562
563
564
565
566
567
568
Note: See TracBrowser for help on using the repository browser.