source: trunk/Mars/mmontecarlo/MMcThresholdCalc.cc@ 17689

Last change on this file since 17689 was 2207, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 5.1 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): Javier Lopez 05/2001 <mailto:jlopez@ifae.es>
19! Author(s): Thomas Bretz 06/2001 <mailto:tbretz@uni-sw.gwdg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26///////////////////////////////////////////////////////////////////////////
27//
28// MMcThresholdCalc
29//
30// Input Containers:
31// MMcEvt, MMcTrig;*
32//
33// Output Containers:
34// MHMcEnergies
35//
36/////////////////////////////////////////////////////////////////////////////
37
38#include "MMcThresholdCalc.h"
39
40#include <math.h>
41
42#include "MParList.h"
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MMcEvt.hxx"
48#include "MMcTrig.hxx"
49
50#include "MHMcEnergy.h"
51
52ClassImp(MMcThresholdCalc);
53
54using namespace std;
55
56const Float_t MMcThresholdCalc::fSqrt2 = sqrt(2.);
57
58// --------------------------------------------------------------------------
59//
60// Default Constructor.
61//
62// Specify the number of trigger conditions you want to use.
63// The default is 0.
64//
65// dim < 0: use only condition number dim (eg "MMcTrig;3")
66// dim = 0: use only condition without a number ("MMcTrig")
67// dim > 0: use conditions up to dim (from "MMcTrig;1" to "MMcTrig;dim")
68//
69MMcThresholdCalc::MMcThresholdCalc(const Int_t dim, const char* name,
70 const char* title)
71 : fMcTrig(NULL), fEnergy(NULL)
72{
73 fName = name ? name : "MMcThresholdCalc";
74 fTitle = title ? title : "Task to calculate the energy threshold from Monte Carlo";
75
76 fFirst = dim>0 ? 1 : -dim;
77 fLast = dim>0 ? dim : -dim;
78
79 fNum = fLast-fFirst+1;
80
81 AddToBranchList("MMcEvt.fEnergy");
82 AddToBranchList("MMcTrig", "fNumFirstLevel", fFirst, fLast);
83}
84
85// -------------------------------------------------------------------------
86//
87// Destructor.
88//
89MMcThresholdCalc::~MMcThresholdCalc()
90{
91 if (fMcTrig)
92 delete fMcTrig;
93
94 if (fEnergy)
95 delete fEnergy;
96}
97
98// --------------------------------------------------------------------------
99//
100// connect Monte Carlo data with this task
101//
102Int_t MMcThresholdCalc::PreProcess(MParList* pList)
103{
104 //
105 // This task takes into accout if the root file has one trigger
106 // condition (MMcTrig) or severl of them (MMcTrig;#.)
107
108 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
109 if (!fMcEvt)
110 {
111 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
112 return kFALSE;
113 }
114
115 UInt_t num;
116
117 fMcTrig = new TObjArray(pList->FindObjectList("MMcTrig", fFirst, fLast));
118 num = fMcTrig->GetEntriesFast();
119 if (num != fNum)
120 {
121 *fLog << err << dbginf << fNum << " MMcTrig objects requested, ";
122 *fLog << num << " are available... aborting." << endl;
123 return kFALSE;
124 }
125
126 fEnergy = new TObjArray(pList->FindCreateObjList("MHMcEnergy", fFirst, fLast));
127 num = fMcTrig->GetEntriesFast();
128 if (num != fNum)
129 {
130 *fLog << err << dbginf << fNum << " MHMcEnergy objects requested, ";
131 *fLog << num << " are available... aborting." << endl;
132 return kFALSE;
133 }
134
135 return kTRUE;
136}
137
138// --------------------------------------------------------------------------
139//
140// The histograms are filled with log10 of the energy for triggered
141// events and weighted with 1/E because it is needed the dN/dE vs. logE
142// distribution to get the energy threshold.
143//
144Int_t MMcThresholdCalc::Process()
145{
146 const Float_t energy = fMcEvt->GetEnergy();
147 const Float_t lg10 = log10(energy);
148 const Float_t reciproc = 1./energy;
149
150 for (UInt_t i=0; i<fNum; i++)
151 if (GetTrig(i)->GetFirstLevel()>0)
152 GetHEnergy(i)->Fill(lg10, reciproc);
153
154 return kTRUE;
155}
156
157// --------------------------------------------------------------------------
158//
159// fit the energy distribution to get the threshold
160// Some iterations are done to be sure the fit parameters converge.
161//
162Int_t MMcThresholdCalc::PostProcess()
163{
164 for (UInt_t i=0; i<fNum; i++)
165 {
166 MHMcEnergy &hist = *GetHEnergy(i);
167
168 Float_t peak;
169 Float_t sigma;
170
171 hist.Fit(1, 3);
172
173 peak = hist.GetGaussPeak();
174 sigma = hist.GetGaussSigma();
175 hist.Fit(peak - 2. *sigma, peak + 2. *sigma);
176
177 peak = hist.GetGaussPeak();
178 sigma = hist.GetGaussSigma();
179 hist.Fit(peak - fSqrt2*sigma, peak + fSqrt2*sigma);
180
181 hist.SetReadyToSave();
182 }
183 return kTRUE;
184}
185
Note: See TracBrowser for help on using the repository browser.