source: trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc@ 1861

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