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

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