source: trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.cc@ 6977

Last change on this file since 6977 was 6967, checked in by marcos, 21 years ago
*** empty log message ***
File size: 10.8 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Marcos Lopez 12/2004 <mailto:marcos@gae.ucm.es>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// MHExcessEnergyThetaONOFF //
28// //
29// 3D-histogram in alpha vs. E-est and Theta //
30// //
31//////////////////////////////////////////////////////////////////////////////
32
33#include "MHExcessEnergyThetaONOFF.h"
34
35#include <TCanvas.h>
36#include <THStack.h>
37#include <TLegend.h>
38#include <TStyle.h>
39
40#include "MMcEvt.hxx"
41#include "MHillasSrc.h"
42#include "MEnergyEst.h"
43#include "MPointingPos.h"
44#include "MRawRunHeader.h"
45#include "MHAlphaEnergyTheta.h"
46#include "MHFindSignificance.h"
47
48#include "MBinning.h"
49#include "MParList.h"
50
51#include "MLog.h"
52#include "MLogManip.h"
53
54#include <TGraphErrors.h>
55
56#include "MHMyFindSignificanceONOFF.h"
57
58
59ClassImp(MHExcessEnergyThetaONOFF);
60
61using namespace std;
62
63// --------------------------------------------------------------------------
64//
65// Default Constructor. It sets name and title of the histogram.
66//
67MHExcessEnergyThetaONOFF::MHExcessEnergyThetaONOFF(const char *name, const char *title)
68 : fHist("","",10,0,1, 10,0,1)
69{
70 //
71 // set the name and title of this object
72 //
73 fName = name ? name : "MHExcessEnergyThetaONOFF";
74 fTitle = title ? title : "# Excess events vs. E and Theta";
75
76
77 fHist.SetDirectory(NULL);
78
79 fHist.SetTitle("# Excess events vs. E and Theta");
80 fHist.SetXTitle("E [GeV]");
81 fHist.SetYTitle("\\Theta [\\circ]");
82 fHist.SetZTitle("# Excess events");
83
84
85
86}
87
88// --------------------------------------------------------------------------
89//
90// Set binnings and prepare filling of the histogram
91//
92Bool_t MHExcessEnergyThetaONOFF::SetupFill(const MParList *plist)
93{
94
95 // fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst");
96// if (!fEnergy)
97// {
98// *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
99// return kFALSE;
100// }
101
102
103// //
104// // Binning
105// //
106// MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
107// MBinning* binsalphaflux = (MBinning*)plist->FindObject("BinningAlphaFlux");
108// MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
109// if (!binsenergy || !binsalphaflux || !binstheta)
110// {
111// *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
112// return kFALSE;
113// }
114
115// // SetBinning(&fHist, binsalphaflux, binsenergy, binstheta);
116// SetBinning(&fHist, binsenergy, binstheta, binsalphaflux);
117
118// fHist.Sumw2();
119
120
121 return kTRUE;
122}
123
124// --------------------------------------------------------------------------
125//
126// Look in the parlist for MMcEvt or MPointingPos depending on the run type.
127//
128Bool_t MHExcessEnergyThetaONOFF::ReInit(MParList *pList)
129{
130
131// // This must be done in ReInit because in PreProcess the
132// // runheaders are not available
133// const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
134// if (!runheader)
135// {
136// *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
137// *fLog << warn << dbginf << "Warning - Assuming data type file...searching for MPointingPos..." << endl;
138// }
139
140
141// if (runheader && runheader->IsMonteCarloRun())
142// {
143// fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
144// if (!fMcEvt)
145// {
146// *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
147// return kFALSE;
148// }
149// }
150// else
151// {
152// fPointingPos = (MPointingPos*)pList->FindObject("MPointingPos");
153// if (!fPointingPos)
154// {
155// *fLog << err << dbginf << "PointingPos not found... aborting." << endl;
156// return kFALSE;
157// }
158// }
159
160 return kTRUE;
161}
162
163
164
165// --------------------------------------------------------------------------
166//
167// Fill the histogram
168//
169Bool_t MHExcessEnergyThetaONOFF::Fill(const MParContainer *par, const Stat_t w)
170{
171 // MHillasSrc &hil = *(MHillasSrc*)par;
172
173// // fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(),
174// // fMcEvt->GetTelescopeTheta()*kRad2Deg, w);
175
176
177// const Double_t Zd = (fMcEvt) ? fMcEvt->GetTelescopeTheta()*kRad2Deg : fPointingPos->GetZd() ;
178
179// fHist.Fill(fEnergy->GetEnergy(), Zd, hil.GetAlpha(), w);
180
181// // *fLog<< Zd << " " << fEnergy->GetEnergy() << " " << hil.GetAlpha() << endl;
182
183 return kTRUE;
184}
185
186
187
188// --------------------------------------------------------------------------
189//
190// Get the Alpha distribution for each bin in energy and theta,
191// fit it with MHFindSignificance to get the number of Excess events
192// and it error.
193//
194void MHExcessEnergyThetaONOFF::Calc(MHAlphaEnergyTheta* hAlphaON, MHAlphaEnergyTheta* hAlphaOFF)
195{
196 *fLog << endl;
197 fLog->Separator("Excess Events Calculation");
198
199
200 TH3D* aetON = (TH3D*)hAlphaON->GetHist();
201 TH3D* aetOFF = (TH3D*)hAlphaOFF->GetHist();
202
203 const TAxis* axisEnergy = aetON->GetXaxis();
204 const TAxis* axisTheta = aetON->GetYaxis();
205 const Int_t energyBins = aetON->GetXaxis()->GetNbins();
206 const Int_t thetaBins = aetON->GetYaxis()->GetNbins();;
207
208 MH::SetBinning(&fHist, axisEnergy, axisTheta);
209
210
211
212 //
213 // Loop over Theta bins
214 //
215 for(Int_t iy=1; iy<=thetaBins; iy++)
216 {
217 TCanvas* c= new TCanvas(Form("zbin %d",iy),Form("zbin %d",iy));
218 c->Divide(4,3);
219
220 //
221 // Loop over Energy bins
222 //
223 for(Int_t ix=1; ix<=energyBins; ix++)
224 {
225 const Double_t e = aetON->GetXaxis()->GetBinCenter(ix);
226 const Double_t e1 = aetON->GetXaxis()->GetBinLowEdge(ix);
227 const Double_t e2 = aetON->GetXaxis()->GetBinLowEdge(ix+1);
228
229 //
230 // Get Alpha plot for that Theta & Energy bin
231 //
232 TH1* halphaON = aetON->ProjectionZ("AlphaTempON",ix,ix,iy,iy);
233 halphaON->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2));
234 halphaON->SetDirectory(NULL);
235
236 TH1* halphaOFF = aetOFF->ProjectionZ("AlphaTempOFF",ix,ix,iy,iy);
237 halphaOFF->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2));
238 halphaOFF->SetDirectory(NULL);
239
240 TCanvas *c = new TCanvas;
241 halphaOFF->Draw();
242
243
244 //
245 // Fit Alpha plot to get Nex
246 //
247 MHMyFindSignificanceONOFF findsig;
248 findsig.SetRebin(kFALSE);
249 //findsig.SetRebin(kTRUE);
250 double alphasig=25;
251
252 // if(ix<5)
253// alphasig = 20;
254// else
255// alphasig = 20;
256
257// switch (ix)
258// {
259// case 1:
260// alphasig = 20;
261// break;
262// case 2:
263// alphasig = 20;
264// break;
265// case 3:
266// alphasig = 20;
267// break;
268// case 4:
269// alphasig = 20;
270// break;
271// case 5:
272// alphasig = 20;
273// break;
274// case 6:
275// alphasig = 10;
276// break;
277// case 7:
278// alphasig = 10;
279// break;
280// case 8:
281// alphasig = 10;
282// break;
283// case 9:
284// alphasig = 10;
285// break;
286// case 10:
287
288// alphasig = 10;
289// break;
290// }
291
292 const Double_t AlphaMinON = 30;
293 const Double_t AlphaMaxON = 90;
294 const Double_t AlphaMinOFF = 0;
295 const Double_t AlphaMaxOFF = AlphaMaxON;
296 const Double_t Degree = 4;
297
298 fLog->SetNullOutput(kTRUE);
299
300
301 Bool_t rc = findsig.FindSigmaONOFF(halphaON, halphaOFF, alphasig, AlphaMinON,AlphaMaxON,Degree);
302
303 // Double_t SigmaGauss = findsig.GetSigmaGauss();
304
305 // alphasig = SigmaGauss*3;
306// rc = findsig.FindSigma(halpha, 25,90, degree, alphasig, 0,1,0);
307
308 fLog->SetNullOutput(kFALSE);
309
310
311
312 if(ix<=20)
313 {
314 c->cd(ix);
315 gPad->SetGridx();
316 gPad->SetGridy();
317
318 findsig.Draw("all");
319 }
320
321 Double_t Nex = 0;
322 Double_t dNex = 0;
323 Double_t Sig = 0;
324
325 if(!rc)
326 {
327 cout << "ERROR: Fit not possible. " << endl;
328 }
329 // else
330 {
331 //SigmaGauss = findsig.GetSigmaGauss();
332 Nex = findsig.GetNex();
333 dNex = findsig.GetdNex();
334 Sig = findsig.GetSigLiMa();
335 }
336 // Avoid to have a negative number of excess events, which
337 // will give a negative flux which doesn't make sense
338 if(Nex<0)
339 {
340 Nex = 0;
341 dNex = 0;
342 Sig = 0;
343 }
344
345 fHist.SetBinContent(ix,iy, Nex);
346 fHist.SetBinError(ix,iy, dNex);
347
348
349 // Print
350 *fLog<< endl
351 << " Theta Bin = " << iy << endl
352 << " Energy bin width E = " << ix << endl
353 << " Alphasig = " <<alphasig << endl
354 << " N = " << Nex <<" +- " << dNex << endl
355 << " Significane = " << Sig << endl;
356 }
357 }
358
359
360 *fLog << "Total Number of Excess evts = " << fHist.Integral() << endl;
361
362
363}
364
365
366
367// --------------------------------------------------------------------------
368//
369// Draw the histogram
370//
371void MHExcessEnergyThetaONOFF::Draw(Option_t *opt)
372{
373 TCanvas *c1 = new TCanvas("# Excess events vs. E and Theta","# Excess events vs. E and Theta");
374 c1->SetLogx();
375 c1->SetLogz();
376
377 fHist.SetStats(0);
378 fHist.Draw("lego");
379
380 TCanvas *c2 = new TCanvas;
381 c2->SetLogx();
382
383 TH1* h= fHist.ProjectionX();
384 h->SetStats(0);
385 h->SetTitle("No. Excess events Vs. Energy");
386 h->SetXTitle("E [GeV]");
387 h->SetYTitle("# Excess events");
388
389 h->Draw("e1");
390
391
392 /////////////
393 TCanvas *c3 = new TCanvas;
394
395
396 Double_t res[10] ={0.224281,0.30083, 0.397355, 0.455888, 0.482985, 0.479416, 0.441429, 0.424654, 0.366658, 0.287384};
397
398 TGraphErrors* graphIntegral = new TGraphErrors(h);
399
400 for(int i=0; i<graphIntegral->GetN(); i++)
401 {
402 Double_t errorX = graphIntegral->GetErrorX(i);
403 Double_t errorY = graphIntegral->GetErrorY(i);
404 errorX *= res[i];
405
406 graphIntegral->SetPointError(i,errorX,errorY);
407 }
408
409
410 c3->SetLogx();
411 c3->SetLogy();
412 c3->SetGridx();
413 c3->SetGridy();
414
415 graphIntegral->SetMarkerStyle(8);
416 graphIntegral->GetHistogram() ->SetXTitle("E [GeV]");
417 graphIntegral->GetHistogram() ->SetYTitle("# Excess events");
418 graphIntegral->Draw("Ap");
419
420}
Note: See TracBrowser for help on using the repository browser.