source: trunk/Mars/manalysis/MEnergyEstParam.cc@ 13736

Last change on this file since 13736 was 8020, checked in by tbretz, 19 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): Thomas Bretz 9/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26// //
27// MEnergyEstParam //
28// //
29// Task to estimate the energy using a parametrization. //
30// Make sure, that all source dependant containers are based on the same //
31// set of basic hillas parameters //
32// //
33/////////////////////////////////////////////////////////////////////////////
34#include "MEnergyEstParam.h"
35
36#include "MParList.h"
37
38#include "MMcEvt.hxx"
39#include "MHMatrix.h"
40#include "MHillasSrc.h"
41#include "MEnergyEst.h"
42
43#include "MLog.h"
44#include "MLogManip.h"
45
46ClassImp(MEnergyEstParam);
47
48using namespace std;
49
50// --------------------------------------------------------------------------
51//
52// Initialize the coefficients with (nonsense) values
53//
54void MEnergyEstParam::InitCoefficients()
55{
56 fA.Set(5);
57 fA[0] = -2539; // [cm]
58 fA[1] = 900; // [cm]
59 fA[2] = 17.5; // [cm]
60 fA[3] = 0.006;
61 fA[4] = 58.3;
62
63 /*
64 fA[0] = 39.667402; // [cm]
65 fA[1] = 143.081912; // [cm]
66 fA[2] = -395.501677; // [cm]
67 fA[3] = 0.006193;
68 fA[4] = 0;
69 */
70
71 fB.Set(7);
72 fB[0] = -8.69; // [GeV]
73 fB[1] = -0.069; // [GeV]
74 fB[2] = 0.000655; // [GeV]
75 fB[3] = 0.0326; // [GeV]
76 fB[4] = 0.000225; // [GeV]
77 fB[5] = 4.13e-8; // [GeV]
78 fB[6] = 0.05;
79
80 /*
81 fB[0] = 45.037867; // [GeV]
82 fB[1] = 0.087222; // [GeV]
83 fB[2] = -0.208065; // [GeV/cm]
84 fB[3] = 78.164942; // [GeV]
85 fB[4] = -159.361283; // [GeV]
86 fB[5] = -0.130455; // [GeV/cm]
87 fB[6] = 0.051139;
88 */
89}
90
91// --------------------------------------------------------------------------
92//
93// Default constructor. Give the name of the parameter container (MHillas)
94// storing width, length and size (Default="MHillas").
95// For the Zenith Angle MMcEvt.fTheta is used.
96//
97MEnergyEstParam::MEnergyEstParam(const char *hillas, const char *name, const char *title)
98 : fMatrix(NULL)
99{
100 fName = name ? name : "MEnergyEstParam";
101 fTitle = title ? title : "Task to estimate the energy";
102
103 fHillasName = hillas;
104
105 fPairs = new TList;
106 fEnergy = new TList;
107 fHillasSrc = new TList;
108
109 fPairs->SetOwner();
110
111 InitCoefficients();
112
113 AddToBranchList("MMcEvt.fTelescopeTheta");
114 AddToBranchList(fHillasName+".fSize");
115 AddToBranchList(fHillasName+".fWidth");
116 AddToBranchList(fHillasName+".fLength");
117}
118
119
120// --------------------------------------------------------------------------
121//
122// Destructor.
123//
124MEnergyEstParam::~MEnergyEstParam()
125{
126 delete fPairs;
127 delete fEnergy;
128 delete fHillasSrc;
129}
130
131// --------------------------------------------------------------------------
132//
133// Check for all necessary parameter containers.
134//
135Int_t MEnergyEstParam::PreProcess(MParList *plist)
136{
137 if (!fMatrix)
138 {
139 fHillas = (MHillas*)plist->FindObject(fHillasName, "MHillas");
140 if (!fHillas)
141 {
142 *fLog << err << dbginf << fHillasName << " [MHillas] not found... aborting." << endl;
143 return kFALSE;
144 }
145
146 fMc = (MMcEvt*)plist->FindObject("MMcEvt");
147 if (!fMc)
148 {
149 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
150 return kFALSE;
151 }
152 }
153
154 TObject *obj = NULL;
155 TIter Next(fPairs);
156
157 while ((obj=Next()))
158 {
159 TObject *o = plist->FindCreateObj(obj->GetTitle(), "MEnergyEst");
160 if (!o)
161 return kFALSE;
162
163 if (!fEnergy->FindObject(obj->GetTitle()))
164 fEnergy->Add(o);
165
166 if (fMatrix)
167 continue;
168
169 o = plist->FindObject(obj->GetName(), "MHillasSrc");
170 if (!o)
171 {
172 *fLog << err << dbginf << obj->GetName() << " not found... aborting." << endl;
173 return kFALSE;
174 }
175 fHillasSrc->Add(o);
176 }
177
178 return kTRUE;
179}
180
181// --------------------------------------------------------------------------
182//
183// Set the five coefficients for the estimation of the impact parameter.
184// Argument must ba a TArrayD of size 5.
185//
186void MEnergyEstParam::SetCoeffA(const TArrayD &arr)
187{
188 if (arr.GetSize() != fA.GetSize())
189 {
190 *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl;
191 return;
192 }
193
194 fA = arr;
195}
196
197// --------------------------------------------------------------------------
198//
199// Set the seven coefficients for the estimation of the energy.
200// Argument must ba a TArrayD of size 7.
201//
202void MEnergyEstParam::SetCoeffB(const TArrayD &arr)
203{
204 if (arr.GetSize() != fB.GetSize())
205 {
206 *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl;
207 return;
208 }
209
210 fB = arr;
211}
212
213// --------------------------------------------------------------------------
214//
215// Set the twelve coefficients for the estimation of impact and energy.
216// Argument must ba a TArrayD of size 12.
217//
218void MEnergyEstParam::SetCoeff(const TArrayD &arr)
219{
220 if (arr.GetSize() != fA.GetSize()+fB.GetSize())
221 {
222 *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl;
223 return;
224 }
225
226 fA = TArrayD(fA.GetSize(), arr.GetArray());
227 fB = TArrayD(fB.GetSize(), arr.GetArray() + fA.GetSize());
228}
229
230// --------------------------------------------------------------------------
231//
232// Returns the mapped value from the Matrix
233//
234Double_t MEnergyEstParam::GetVal(Int_t i) const
235{
236 return (*fMatrix)[fMap[i]];
237}
238
239// --------------------------------------------------------------------------
240//
241// You can use this function if you want to use a MHMatrix instead of the
242// given containers. This function adds all necessary columns to the
243// given matrix. Afterward you should fill the matrix with the corresponding
244// data (eg from a file by using MHMatrix::Fill). If you now loop
245// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
246// will take the values from the matrix instead of the containers.
247//
248void MEnergyEstParam::InitMapping(MHMatrix *mat)
249{
250 if (fMatrix)
251 return;
252
253 fMatrix = mat;
254
255 fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta");
256 fMap[1] = fMatrix->AddColumn(fHillasName+".fWidth");
257 fMap[2] = fMatrix->AddColumn(fHillasName+".fLength");
258 fMap[3] = fMatrix->AddColumn(fHillasName+".fSize");
259
260 Int_t col = 4;
261 TIter Next(fPairs);
262 TObject *o = NULL;
263 while ((o=Next()))
264 fMap[col++] = fMatrix->AddColumn(TString(o->GetName())+".fDist");
265}
266
267void MEnergyEstParam::StopMapping()
268{
269 fMatrix = NULL;
270 fPairs->Clear();
271 fHillasSrc->Clear();
272 fEnergy->Clear();
273}
274
275// --------------------------------------------------------------------------
276//
277// Add a pair of input/output containers.
278// eg. Add("MHillasSrc", "MEnergyEst");
279// Useful if you want to estimate the stuff for the source and antisource
280//
281void MEnergyEstParam::Add(const TString hillas, const TString energy)
282{
283 fPairs->Add(new TNamed(hillas, energy));
284
285 AddToBranchList(hillas+".fDist");
286}
287
288void MEnergyEstParam::Print(Option_t *opt) const
289{
290 for (int i=0; i<fA.GetSize(); i++)
291 *fLog << all << "fA[" << i << "]=" << const_cast<TArrayD&>(fA)[i] << endl;
292 for (int i=0; i<fB.GetSize(); i++)
293 *fLog << all << "fB[" << i << "]=" << const_cast<TArrayD&>(fB)[i] << endl;
294}
295
296// --------------------------------------------------------------------------
297//
298// Estimates the impact parameter and energy using a parameterization
299// (see code)
300//
301Int_t MEnergyEstParam::Process()
302{
303 const Double_t theta = fMatrix ? GetVal(0) : fMc->GetTelescopeTheta();
304 const Double_t width = fMatrix ? GetVal(1) : fHillas->GetWidth();
305 const Double_t length = fMatrix ? GetVal(2) : fHillas->GetLength();
306 const Double_t size = fMatrix ? GetVal(3) : fHillas->GetSize();
307
308 const Double_t k = 1./cos(theta);
309 const Double_t k2 = k*k;
310
311 const Double_t i0 = k * (fA[3]*k+1)/(fA[3]+1);
312 const Double_t i1 = fA[0] + fA[2]*width;
313
314
315 const Double_t e0 = k2 * (fB[6]*k2+1)/(fB[6]+1);
316
317 /* MY PARAM */
318 //const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*width;
319 //const Double_t e2 = fB[2] + fB[5]*size*width;
320
321 /* MARCOS */
322 /*
323 const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*length;
324 const Double_t e2 = fB[2] + fB[5]*length;
325 */
326
327 /* DANIEL */
328 const Double_t e1 = fB[0] + fB[1]*size/(length*width) +
329 fB[3]*length + fB[4]*size/width;
330
331 const Double_t e2 = fB[2] + fB[5]*length;
332
333 TIter NextH(fHillasSrc);
334 TIter NextE(fEnergy);
335
336 Int_t col = 4;
337
338 while (1)
339 {
340 MEnergyEst *est = (MEnergyEst*)NextE();
341 if (!est)
342 break;
343
344 const Double_t dist = fMatrix ? GetVal(col++) : ((MHillasSrc*)NextH())->GetDist();
345
346 /* DANIEL - MARCOS */
347 const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm]
348 Double_t er = e0 * (e1 + e2*ir); // [GeV]
349
350 /* THOMAS BRETZ */
351 // if (width==0) return kCONTINUE;
352 // const Double_t ir = i0 * (i1 + dist*(fA[1]/width + fA[4]/log10(size))); // [cm]
353 //Double_t er = e0 * (e1 + e2*ir); // [GeV]
354
355 /* MKA */
356 //Double_t er = e0 * (fB[0] + fB[1]*size/width + fB[2]*ir /*+ d*leakage*/);
357
358 if (width==0)
359 return kCONTINUE;
360
361 // !Finitite includes IsNaN
362 if (!TMath::Finite(ir))
363 *fLog << all << theta << " " << width << " " << length << " " << size << " " << dist << endl;
364 if (!TMath::Finite(er))
365 {
366 *fLog << all << theta << " " << width << " " << length << " " << size << " " << dist << endl;
367 er = 0;
368 }
369
370 est->SetEnergy(er);
371 est->SetImpact(ir);
372 est->SetReadyToSave();
373 }
374
375 return kTRUE;
376}
377
Note: See TracBrowser for help on using the repository browser.