source: branches/MarsMoreSimulationTruth/manalysisct1/MCT1SupercutsCalc.cc

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