source: branches/removing_cpp11_features/mmontecarlo/MMcUnfoldCoeffCalc.cc@ 19796

Last change on this file since 19796 was 6557, checked in by rico, 20 years ago
*** empty log message ***
File size: 6.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! Author(s): Javier Rico 2/2005 <mailto:jrico@ifae.es>
18!
19! Copyright: MAGIC Software Development, 2000-2005
20!
21!
22\* ======================================================================== */
23
24//////////////////////////////////////////////////////////////////////////////
25//
26// MMcUnfoldCoeffCalc
27// Task to compute the coefficients for energy unfolding (to be used
28// typically in an iterative process).
29//
30// Input containers:
31// MMcEvt
32// MMcEvtBasic
33// MEnergyEst
34// binsTheta [MBinning]
35// binsEnergy [MBinning]
36//
37// Output containers:
38// MHMcUnfoldCoeff
39//
40// The binning for energy and theta are looked in the parameter list
41// The tentative spectrum MUST be provided using SetSpectrum, otherwise
42// it is used for default pow(energy,-1.6) [~Crab at 1 TeV]
43//
44// This task is supposed to be in a two-looped macro/program:
45// In the first one MMcEvtBasic is read in order to compute the
46// weights to convert from original to tentative spectrum.
47// In the second one MMcEvt and MEnergyEst are read and fill the
48// MC and estimated energy that are divided in the postprocess to
49// compute the coefficients
50//
51//////////////////////////////////////////////////////////////////////////////
52
53#include "TF1.h"
54
55#include "MMcUnfoldCoeffCalc.h"
56#include "MHMcUnfoldCoeff.h"
57
58#include "MParList.h"
59
60#include "MLog.h"
61#include "MLogManip.h"
62
63#include "MMcEvt.hxx"
64#include "MMcEvtBasic.h"
65#include "MEnergyEst.h"
66
67ClassImp(MMcUnfoldCoeffCalc);
68
69using namespace std;
70
71// -------------------------------------------------------------------------
72//
73// Constructor
74//
75MMcUnfoldCoeffCalc::MMcUnfoldCoeffCalc(const char* name, const char* title) :
76 fSpectrum(NULL)
77{
78 fName = name ? name : "MMcUnfoldCoeffCalc";
79 fTitle = title ? title : "Task to calculate the unfolding coefficients";
80}
81
82// -------------------------------------------------------------------------
83//
84// PreProcess.
85// Create MHMcUnfoldCoeff object if necessary. Search for other necessary
86// containers: if MMcEvt is found, it means we are in the loop over the Events
87// tree, and so we must fill the histogram MHMcUnfoldCoeff::fHistMcE and
88// MHMcUnfoldCoeff::fHistEstE (using MHMcUnfoldCoeff::FillSel()).
89// If MMcEvt is not found, it means that we are in the loop over the
90// "OriginalMC" tree, containing all the original Corsika events
91// produced, and hence we must fill the histogram fHistAll through
92// MHMcUnfoldCoeff::FillAll()
93//
94Int_t MMcUnfoldCoeffCalc::PreProcess (MParList *pList)
95{
96 fUnfoldCoeff = (MHMcUnfoldCoeff*)pList->FindCreateObj("MHMcUnfoldCoeff");
97
98 if (!fBinsTheta)
99 {
100 fBinsTheta = (MBinning*) pList->FindObject("binsTheta");
101 if (!fBinsTheta)
102 {
103 *fLog << err << "Coarse Theta binning not found... Aborting."
104 << endl;
105 return kFALSE;
106 }
107 }
108
109 if (!fBinsEnergy)
110 {
111 fBinsEnergy = (MBinning*) pList->FindObject("binsEnergy");
112 if (!fBinsEnergy)
113 {
114 *fLog << err << "Coarse Energy binning not found... Aborting."
115 << endl;
116 return kFALSE;
117 }
118 }
119
120 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
121 if (fMcEvt)
122 {
123 *fLog << inf << "MMcEvt found... Filling MHMcUnfoldCoeff.fHistSel..." << endl;
124
125 fEnergyEst = (MEnergyEst*)pList->FindObject("MEnergyEst");
126 if(!fEnergyEst)
127 {
128 *fLog << err << "MEnergyEst not found... aborting..." << endl;
129 return kFALSE;
130 }
131
132 return kTRUE;
133 }
134
135 *fLog << inf << "MMcEvt not found... looking for MMcEvtBasic..." << endl;
136
137 fMcEvtBasic = (MMcEvtBasic*) pList->FindObject("MMcEvtBasic");
138
139 if (fMcEvtBasic)
140 {
141 fUnfoldCoeff->SetCoarseBinnings(*fBinsEnergy,*fBinsTheta);
142
143 *fLog << inf << "MMcEvtBasic found... Filling MHMcUnfoldCoeff.fHistAll..." << endl;
144 return kTRUE;
145 }
146
147 *fLog << err << "MMcEvtBasic not found. Aborting..." << endl;
148
149 return kFALSE;
150}
151
152// -------------------------------------------------------------------------
153//
154// Depending on whether fMcEvt or fMcEvtBasic are available, fill
155// MHMcUnfoldCoeff::fHistAll, else, if fMcEvt is available, fill
156// MHMcUnfoldCoeff::fHistMcE and MHMcUnfoldCoeff::fHistEstE
157//
158Int_t MMcUnfoldCoeffCalc::Process()
159{
160 Double_t mcenergy = fMcEvt ? fMcEvt->GetEnergy() : fMcEvtBasic->GetEnergy();
161 Double_t estenergy = fEnergyEst ? fEnergyEst->GetEnergy() : 0;
162 Double_t theta = fMcEvt? fMcEvt->GetTelescopeTheta()*TMath::RadToDeg() : 0; // in deg
163
164 if (fMcEvt)
165 fUnfoldCoeff->FillSel(mcenergy, estenergy, theta);
166 else
167 fUnfoldCoeff->FillAll(mcenergy);
168
169 return kTRUE;
170}
171
172// -------------------------------------------------------------------------
173//
174// Postprocess.
175// If both fHistMcE and fHistEstE are already filled, calculate
176// the coefficients.
177// Else it means we still have to run one more loop, and then the
178// weights are computed
179//
180Int_t MMcUnfoldCoeffCalc::PostProcess()
181{
182 if(((TH2D*)fUnfoldCoeff->GetHistMcE())->GetEntries() > 0 &&
183 ((TH2D*)fUnfoldCoeff->GetHistEstE())->GetEntries() > 0)
184 {
185 fUnfoldCoeff->ComputeCoefficients();
186 return kTRUE;
187 }
188
189 if(((TH1D*)fUnfoldCoeff->GetHistAll())->GetEntries() <= 0)
190 {
191 *fLog << err << "Postprocess reached with no histogram computed. Aborting" << endl;
192 return kFALSE;
193 }
194
195
196 if(fSpectrum)
197 fUnfoldCoeff->ComputeWeights(fSpectrum);
198 else
199 {
200 // default spectrum in case no other one is provided
201 TF1 spectrum("func","pow(x,-1.6)",2.,20000.);
202 fUnfoldCoeff->ComputeWeights(&spectrum);
203 }
204
205 return kTRUE;
206}
Note: See TracBrowser for help on using the repository browser.