source: trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc@ 9223

Last change on this file since 9223 was 6855, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 11.4 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): Wolfgang Wittek, 04/2003 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MSupercutsCalc
28//
29// this class calculates the hadronness for the supercuts
30// the parameters of the supercuts are taken
31// from the container MSupercuts
32//
33/////////////////////////////////////////////////////////////////////////////
34#include "MSupercutsCalc.h"
35
36#include <math.h>
37#include <fstream>
38
39#include "TFile.h"
40#include "TArrayD.h"
41
42#include "MParList.h"
43#include "MHillasExt.h"
44#include "MHillasSrc.h"
45#include "MNewImagePar.h"
46#include "MMcEvt.hxx"
47#include "MGeomCam.h"
48#include "MHadronness.h"
49#include "MHMatrix.h"
50#include "MSupercuts.h"
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55ClassImp(MSupercutsCalc);
56
57using namespace std;
58
59
60// --------------------------------------------------------------------------
61//
62// constructor
63//
64
65MSupercutsCalc::MSupercutsCalc(const char *hilname,
66 const char *hilsrcname,
67 const char *name, const char *title)
68 : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname),
69 fHilExtName("MHillasExt"), fNewParName("MNewImagePar"),
70 fSuperName("MSupercuts")
71{
72 fName = name ? name : "MSupercutsCalc";
73 fTitle = title ? title : "Class to evaluate the Supercuts";
74
75 fMatrix = NULL;
76}
77
78// --------------------------------------------------------------------------
79//
80Int_t MSupercutsCalc::PreProcess(MParList *pList)
81{
82 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
83 if (!cam)
84 {
85 *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
86 return kFALSE;
87 }
88
89 fMm2Deg = cam->GetConvMm2Deg();
90
91 fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
92 if (!fHadronness)
93 {
94 *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
95 return kFALSE;
96 }
97
98 fSuper = (MSupercuts*)pList->FindObject(fSuperName, "MSupercuts");
99 if (!fSuper)
100 {
101 *fLog << err << fSuperName << " [MSupercuts] not found... aborting." << endl;
102 return kFALSE;
103 }
104
105 if (fMatrix)
106 return kTRUE;
107
108 //-----------------------------------------------------------
109 fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
110 if (!fHil)
111 {
112 *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
113 return kFALSE;
114 }
115
116 fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt");
117 if (!fHilExt)
118 {
119 *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl;
120 return kFALSE;
121 }
122
123 fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
124 if (!fHilSrc)
125 {
126 *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
127 return kFALSE;
128 }
129
130 fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar");
131 if (!fNewPar)
132 {
133 *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl;
134 return kFALSE;
135 }
136
137 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
138 if (!fMcEvt)
139 {
140 *fLog << err << "MMcEvt not found... aborting." << endl;
141 return kFALSE;
142 }
143
144
145 return kTRUE;
146}
147
148// --------------------------------------------------------------------------
149//
150// Calculation of upper and lower limits
151//
152Double_t MSupercutsCalc::CtsMCut(const Double_t* a, Double_t ls, Double_t ct,
153 Double_t ls2, Double_t dd2) const
154{
155 // define cut-function
156 //
157 // dNOMLOGSIZE = 5.0 (=log(150.0)
158 // dNOMCOSZA = 1.0
159 //
160 // a: array of cut parameters
161 // ls: log(SIZE) - dNOMLOGSIZE
162 // ls2: ls^2
163 // ct: Cos(ZA.) - dNOMCOSZA
164 // dd2: DIST^2
165 const Double_t limit =
166 a[0] + a[1] * dd2 + a[2] * ct +
167 ls * (a[3] + a[4] * dd2 + a[5] * ct) +
168 ls2 * (a[6] + a[7] * dd2);
169
170 //*fLog << "MSupercutsCalc::CtsMCut; *a = "
171 // << *a << ", " << *(a+1) << ", " << *(a+2) << ", "
172 // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", "
173 // << *(a+6) << ", " << *(a+7) << endl;
174
175 //*fLog << "MSupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
176 // << ", " << ls2 << ", " << ct << ", " << dd2 << ", "
177 // << limit << endl;
178
179 return limit;
180}
181
182// --------------------------------------------------------------------------
183//
184// Returns the mapped value from the Matrix
185//
186Double_t MSupercutsCalc::GetVal(Int_t i) const
187{
188
189 Double_t val = (*fMatrix)[fMap[i]];
190
191
192 //*fLog << "MSupercutsCalc::GetVal; i, fMatrix, fMap, val = "
193 // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl;
194
195
196 return val;
197}
198
199// --------------------------------------------------------------------------
200//
201// You can use this function if you want to use a MHMatrix instead of the
202// given containers. This function adds all necessary columns to the
203// given matrix. Afterward you should fill the matrix with the corresponding
204// data (eg from a file by using MHMatrix::Fill). If you now loop
205// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
206// will take the values from the matrix instead of the containers.
207//
208void MSupercutsCalc::InitMapping(MHMatrix *mat)
209{
210 if (fMatrix)
211 return;
212
213 fMatrix = mat;
214
215 fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta");
216 fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
217 fMap[2] = fMatrix->AddColumn("MHillas.fLength");
218 fMap[3] = fMatrix->AddColumn("MHillas.fSize");
219 fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
220 fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
221 fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
222 fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)");
223 fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
224 fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc");
225 fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1");
226}
227
228// ---------------------------------------------------------------------------
229//
230// Evaluate dynamical supercuts
231//
232// set hadronness to 0.25 if cuts are fullfilled
233// 0.75 otherwise
234//
235Int_t MSupercutsCalc::Process()
236{
237 const Double_t kNomLogSize = 4.1;
238 const Double_t kNomCosZA = 1.0;
239
240 const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
241 const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth();
242 const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
243 const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize();
244 const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX();
245 const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY();
246 const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist();
247
248 Double_t help;
249 if (!fMatrix)
250 help = TMath::Sign(fHilExt->GetM3Long(),
251 fHilSrc->GetCosDeltaAlpha());
252 const Double_t asym0 = fMatrix ? GetVal(8) : help;
253 const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc();
254 const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1();
255
256 const Double_t newdist = dist0 * fMm2Deg;
257
258 const Double_t dist2 = meanx*meanx + meany*meany;
259 const Double_t dist = sqrt(dist2) * fMm2Deg;
260 const Double_t dd2 = dist*dist;
261
262
263 const Double_t dmls = log(size) - kNomLogSize;
264 const Double_t dmls2 = dmls * dmls;
265
266 const Double_t dmcza = cos(theta) - kNomCosZA;
267
268 const Double_t length = length0 * fMm2Deg;
269 const Double_t width = width0 * fMm2Deg;
270 const Double_t asym = asym0 * fMm2Deg;
271
272 /*
273 *fLog << "newdist, length, width, asym, dist, conc, leakage = "
274 << newdist << ", " << length << ", " << width << ", "
275 << asym << ", " << dist << ", " << conc << ", " << leakage
276 << endl;
277
278 *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = "
279 << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", "
280 << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", "
281
282 << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", "
283 << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", "
284
285 << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", "
286 << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", "
287
288 << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", "
289 << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", "
290
291 << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", "
292 << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", "
293
294 << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", "
295 << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", "
296
297 << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", "
298 << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", "
299 << endl;
300 */
301
302
303 if (
304 //dist < 1.05 &&
305 //newdist < 1.05 &&
306
307 newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) &&
308 newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) &&
309
310 length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) &&
311 length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) &&
312
313 width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) &&
314 width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) &&
315
316 asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) &&
317 asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) &&
318
319 dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) &&
320 dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) &&
321
322 conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) &&
323 conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) &&
324
325 leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) &&
326 leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) )
327
328 fHadronness->SetHadronness(0.25);
329 else
330 fHadronness->SetHadronness(0.75);
331
332 //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl;
333
334 fHadronness->SetReadyToSave();
335
336 return kTRUE;
337}
338
339
340
341
342
343
344
345
346
Note: See TracBrowser for help on using the repository browser.