source: trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc@ 2273

Last change on this file since 2273 was 2272, checked in by wittek, 22 years ago
*** empty log message ***
File size: 12.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): 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/////////////////////////////////////////////////////////////////////////////
30#include "MCT1SupercutsCalc.h"
31
32#include <math.h>
33#include <fstream>
34
35#include "MParList.h"
36#include "MHillasExt.h"
37#include "MHillasSrc.h"
38#include "MMcEvt.hxx"
39#include "MCerPhotEvt.h"
40#include "MGeomCam.h"
41#include "MHadronness.h"
42#include "MHMatrix.h"
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47ClassImp(MCT1SupercutsCalc);
48
49using namespace std;
50
51void MCT1SupercutsCalc::InitParams()
52{
53 fLengthUp.Set(8);
54 fLengthLo.Set(8);
55 fWidthUp.Set(8);
56 fWidthLo.Set(8);
57 fDistUp.Set(8);
58 fDistLo.Set(8);
59 fAsymUp.Set(8);
60 fAsymLo.Set(8);
61 fAlphaUp.Set(8);
62
63 //---------------------------------
64 // cut parameters
65 fLengthUp[0] = 0.315585;
66 fLengthUp[1] = 0.001455;
67 fLengthUp[2] = 0.203198;
68 fLengthUp[3] = 0.005532;
69 fLengthUp[4] =-0.001670;
70 fLengthUp[5] =-0.020362;
71 fLengthUp[6] = 0.007388;
72 fLengthUp[7] =-0.013463;
73
74 fWidthUp[0] = 0.145412;
75 fWidthUp[1] =-0.001771;
76 fWidthUp[2] = 0.054462;
77 fWidthUp[3] = 0.022280;
78 fWidthUp[4] =-0.009893;
79 fWidthUp[5] = 0.056353;
80 fWidthUp[6] = 0.020711;
81 fWidthUp[7] =-0.016703;
82
83 fDistUp[0] = 1.787943;
84 fDistUp[1] = 0.;
85 fDistUp[2] = 2.942310;
86 fDistUp[3] = 0.199815;
87 fDistUp[4] = 0.;
88 fDistUp[5] = 0.249909;
89 fDistUp[6] = 0.189697;
90 fDistUp[7] = 0.;
91
92 fLengthLo[0] = 0.151530;
93 fLengthLo[1] = 0.028323;
94 fLengthLo[2] = 0.510707;
95 fLengthLo[3] = 0.053089;
96 fLengthLo[4] = 0.013708;
97 fLengthLo[5] = 2.357993;
98 fLengthLo[6] = 0.000080;
99 fLengthLo[7] =-0.007157;
100
101 fWidthLo[0] = 0.089187;
102 fWidthLo[1] =-0.006430;
103 fWidthLo[2] = 0.074442;
104 fWidthLo[3] = 0.003738;
105 fWidthLo[4] =-0.004256;
106 fWidthLo[5] =-0.014101;
107 fWidthLo[6] = 0.006126;
108 fWidthLo[7] =-0.002849;
109
110 fDistLo[0] = 0.589406;
111 fDistLo[1] = 0.;
112 fDistLo[2] =-0.083964;
113 fDistLo[3] =-0.007975;
114 fDistLo[4] = 0.;
115 fDistLo[5] = 0.045374;
116 fDistLo[6] =-0.001750;
117 fDistLo[7] = 0.;
118
119 fAsymUp[0] = 0.061267;
120 fAsymUp[1] = 0.014462;
121 fAsymUp[2] = 0.014327;
122 fAsymUp[3] = 0.014540;
123 fAsymUp[4] = 0.013391;
124 fAsymUp[5] = 0.012319;
125 fAsymUp[6] = 0.010444;
126 fAsymUp[7] = 0.008328;
127
128 fAsymLo[0] =-0.012055;
129 fAsymLo[1] = 0.009157;
130 fAsymLo[2] = 0.005441;
131 fAsymLo[3] = 0.000399;
132 fAsymLo[4] = 0.001433;
133 fAsymLo[5] =-0.002050;
134 fAsymLo[6] =-0.000104;
135 fAsymLo[7] =-0.001188;
136
137 fAlphaUp[0] = 13.123440;
138 fAlphaUp[1] = 0.;
139 fAlphaUp[2] = 0.;
140 fAlphaUp[3] = 0.;
141 fAlphaUp[4] = 0.;
142 fAlphaUp[5] = 0.;
143 fAlphaUp[6] = 0.;
144 fAlphaUp[7] = 0.;
145 //---------------------------------
146}
147
148// --------------------------------------------------------------------------
149//
150// Set the parameter values from vector 'par'
151//
152// Attention : it is assumed that there are (9*ncutpar) values
153//
154void MCT1SupercutsCalc::SetParams(Double_t *par)
155{
156 UInt_t ncutpar = fLengthUp.GetSize();
157 UInt_t k0 = 0;
158
159 TArrayD lup(ncutpar, par + k0);
160 SetLengthUp(lup);
161 k0 += ncutpar;
162
163 TArrayD wup(ncutpar, par + k0);
164 SetWidthUp(wup);
165 k0 += ncutpar;
166
167 TArrayD dup(ncutpar, par + k0);
168 SetDistUp(dup);
169 k0 += ncutpar;
170
171 TArrayD llo(ncutpar, par + k0);
172 SetLengthLo(llo);
173 k0 += ncutpar;
174
175 TArrayD wlo(ncutpar, par + k0);
176 SetWidthLo(wlo);
177 k0 += ncutpar;
178
179 TArrayD dlo(ncutpar, par + k0);
180 SetDistLo(dlo);
181 k0 += ncutpar;
182
183 TArrayD aup(ncutpar, par + k0);
184 SetAsymUp(aup);
185 k0 += ncutpar;
186
187 TArrayD alo(ncutpar, par + k0);
188 SetAsymLo(alo);
189 k0 += ncutpar;
190
191 TArrayD alphaup(ncutpar, par + k0);
192 SetAlphaUp(alphaup);
193}
194
195// --------------------------------------------------------------------------
196//
197// Get the parameter values
198//
199// Attention : it is assumed that there are (9*ncutpar) values
200//
201void MCT1SupercutsCalc::GetParams(Double_t *par)
202{
203 UInt_t ncutpar = fLengthUp.GetSize();
204 UInt_t k0 = 0;
205
206 for (UInt_t j=0; j<ncutpar; j++)
207 {
208 UInt_t k = k0 + j;
209 par[k] = fLengthUp[j];
210 }
211 k0 += ncutpar;
212
213 for (UInt_t j=0; j<ncutpar; j++)
214 {
215 UInt_t k = k0 + j;
216 par[k] = fWidthUp[j];
217 }
218 k0 += ncutpar;
219
220 for (UInt_t j=0; j<ncutpar; j++)
221 {
222 UInt_t k = k0 + j;
223 par[k] = fDistUp[j];
224 }
225 k0 += ncutpar;
226
227 for (UInt_t j=0; j<ncutpar; j++)
228 {
229 UInt_t k = k0 + j;
230 par[k] = fLengthLo[j];
231 }
232 k0 += ncutpar;
233
234 for (UInt_t j=0; j<ncutpar; j++)
235 {
236 UInt_t k = k0 + j;
237 par[k] = fWidthLo[j];
238 }
239 k0 += ncutpar;
240
241 for (UInt_t j=0; j<ncutpar; j++)
242 {
243 UInt_t k = k0 + j;
244 par[k] = fDistLo[j];
245 }
246 k0 += ncutpar;
247
248 for (UInt_t j=0; j<ncutpar; j++)
249 {
250 UInt_t k = k0 + j;
251 par[k] = fAsymUp[j];
252 }
253 k0 += ncutpar;
254
255 for (UInt_t j=0; j<ncutpar; j++)
256 {
257 UInt_t k = k0 + j;
258 par[k] = fAsymLo[j];
259 }
260 k0 += ncutpar;
261
262 for (UInt_t j=0; j<ncutpar; j++)
263 {
264 UInt_t k = k0 + j;
265 par[k] = fAlphaUp[j];
266 }
267}
268
269// --------------------------------------------------------------------------
270//
271MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname,
272 const char *hilsrcname, const char *name, const char *title)
273: fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname)
274{
275 fName = name ? name : "MCT1SupercutsCalc";
276 fTitle = title ? title : "Class to evaluate the Supercuts";
277
278 InitParams();
279}
280
281// --------------------------------------------------------------------------
282//
283Int_t MCT1SupercutsCalc::PreProcess(MParList *pList)
284{
285 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
286 if (!cam)
287 {
288 *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
289 return kFALSE;
290 }
291
292 fMm2Deg = cam->GetConvMm2Deg();
293
294 fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
295 if (!fHadronness)
296 {
297 *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
298 return kFALSE;
299 }
300
301 if (fMatrix)
302 return kTRUE;
303
304 fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
305 if (!fHil)
306 {
307 *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
308 return kFALSE;
309 }
310
311 fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
312 if (!fHilSrc)
313 {
314 *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
315 return kFALSE;
316 }
317
318 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
319 if (!fMcEvt)
320 {
321 *fLog << err << "MMcEvt not found... aborting." << endl;
322 return kFALSE;
323 }
324
325 return kTRUE;
326}
327
328// --------------------------------------------------------------------------
329//
330// Calculation of upper and lower limits
331//
332Double_t MCT1SupercutsCalc::CtsMCut(
333#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
334const
335#endif
336 TArrayD &a, Double_t ls, Double_t ct,
337 Double_t ls2, Double_t dd2)
338{
339 // define cut-function
340 //
341 // dNOMLOGSIZE = 4.1 (=log(60.0)
342 // dNOMCOSZA = 1.0
343 //
344 // a: array of cut parameters
345 // ls: log(SIZE) - dNOMLOGSIZE
346 // ls2: ls^2
347 // ct: Cos(ZA.) - dNOMCOSZA
348 // dd2: DIST^2
349 const Double_t limit =
350 a[0] + a[1] * dd2 + a[2] * ct +
351 ls * (a[3] + a[4] * dd2 + a[5] * ct) +
352 ls2 * (a[6] + a[7] * dd2);
353
354 //*fLog << "MCT1SupercutsCalc::CtsMCut; *a = "
355 // << *a << ", " << *(a+1) << ", " << *(a+2) << ", "
356 // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", "
357 // << *(a+6) << ", " << *(a+7) << endl;
358
359 //*fLog << "MCT1SupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
360 // << ", " << ls2 << ", " << ct << ", " << dd2 << ", "
361 // << limit << endl;
362
363 return limit;
364}
365
366// --------------------------------------------------------------------------
367//
368// Returns the mapped value from the Matrix
369//
370Double_t MCT1SupercutsCalc::GetVal(Int_t i) const
371{
372 return (*fMatrix)[fMap[i]];
373}
374
375// --------------------------------------------------------------------------
376//
377// You can use this function if you want to use a MHMatrix instead of the
378// given containers. This function adds all necessary columns to the
379// given matrix. Afterward you should fill the matrix with the corresponding
380// data (eg from a file by using MHMatrix::Fill). If you now loop
381// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
382// will take the values from the matrix instead of the containers.
383//
384void MCT1SupercutsCalc::InitMapping(MHMatrix *mat)
385{
386 if (fMatrix)
387 return;
388
389 fMatrix = mat;
390
391 fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta");
392 fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
393 fMap[2] = fMatrix->AddColumn("MHillas.fLength");
394 fMap[3] = fMatrix->AddColumn("MHillas.fSize");
395 fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
396 fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
397 fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
398 fMap[7] = fMatrix->AddColumn("MHillasSrc.fAlpha");
399}
400
401// ---------------------------------------------------------------------------
402//
403// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich)
404// optimized for mkn 421 2001 data
405//
406// set hadronness to 0.25 if cuts are fullfilled
407// 0.75 otherwise
408//
409Int_t MCT1SupercutsCalc::Process()
410{
411 const Double_t kNomLogSize = 4.1;
412 const Double_t kNomCosZA = 1.0;
413
414 const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
415 const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth();
416 const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
417 const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize();
418 const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX();
419 const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY();
420 const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist();
421
422 const Double_t newdist = dist0 * fMm2Deg;
423
424 const Double_t dist2 = meanx*meanx + meany*meany;
425 const Double_t dd2 = dist2*fMm2Deg;
426 const Double_t dist = sqrt(dist2) * fMm2Deg;
427
428 const Double_t dmls = log(size) - kNomLogSize;
429 const Double_t dmls2 = dmls * dmls;
430
431 const Double_t dmcza = cos(theta) - kNomCosZA;
432
433 const Double_t length = length0 * fMm2Deg;
434 const Double_t width = width0 * fMm2Deg;
435
436 if (newdist < 1.05 &&
437 newdist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) &&
438 newdist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) &&
439 dist < 1.05 &&
440 length < CtsMCut (fLengthUp, dmls, dmcza, dmls2, dd2) &&
441 length > CtsMCut (fLengthLo, dmls, dmcza, dmls2, dd2) &&
442 width < CtsMCut (fWidthUp, dmls, dmcza, dmls2, dd2) &&
443 width > CtsMCut (fWidthLo, dmls, dmcza, dmls2, dd2) &&
444 //asym < CtsMCut (asymup, dmls, dmcza, dmls2, dd2) &&
445 //asym > CtsMCut (asymlow, dmls, dmcza, dmls2, dd2) &&
446 dist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) &&
447 dist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) )
448 fHadronness->SetHadronness(0.25);
449 else
450 fHadronness->SetHadronness(0.75);
451
452 fHadronness->SetReadyToSave();
453
454 return kTRUE;
455}
456//==========================================================================
Note: See TracBrowser for help on using the repository browser.