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

Last change on this file since 2296 was 2282, 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 fLengthUp = lup;
161 k0 += ncutpar;
162
163 TArrayD wup(ncutpar, par + k0);
164 fWidthUp = wup;
165 k0 += ncutpar;
166
167 TArrayD dup(ncutpar, par + k0);
168 fDistUp = dup;
169 k0 += ncutpar;
170
171 TArrayD llo(ncutpar, par + k0);
172 fLengthLo = llo;
173 k0 += ncutpar;
174
175 TArrayD wlo(ncutpar, par + k0);
176 fWidthLo = wlo;
177 k0 += ncutpar;
178
179 TArrayD dlo(ncutpar, par + k0);
180 fDistLo = dlo;
181 k0 += ncutpar;
182
183 TArrayD aup(ncutpar, par + k0);
184 fAsymUp = aup;
185 k0 += ncutpar;
186
187 TArrayD alo(ncutpar, par + k0);
188 fAsymLo = alo;
189 k0 += ncutpar;
190
191 TArrayD alphaup(ncutpar, par + k0);
192 fAlphaUp = 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 fMatrix = NULL;
281}
282
283// --------------------------------------------------------------------------
284//
285Int_t MCT1SupercutsCalc::PreProcess(MParList *pList)
286{
287 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
288 if (!cam)
289 {
290 *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
291 return kFALSE;
292 }
293
294 fMm2Deg = cam->GetConvMm2Deg();
295
296 fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
297 if (!fHadronness)
298 {
299 *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
300 return kFALSE;
301 }
302
303 if (fMatrix)
304 return kTRUE;
305
306 fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
307 if (!fHil)
308 {
309 *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
310 return kFALSE;
311 }
312
313 fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
314 if (!fHilSrc)
315 {
316 *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
317 return kFALSE;
318 }
319
320 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
321 if (!fMcEvt)
322 {
323 *fLog << err << "MMcEvt not found... aborting." << endl;
324 return kFALSE;
325 }
326
327 return kTRUE;
328}
329
330// --------------------------------------------------------------------------
331//
332// Calculation of upper and lower limits
333//
334Double_t MCT1SupercutsCalc::CtsMCut(
335#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
336const
337#endif
338 TArrayD &a, Double_t ls, Double_t ct,
339 Double_t ls2, Double_t dd2)
340{
341 // define cut-function
342 //
343 // dNOMLOGSIZE = 4.1 (=log(60.0)
344 // dNOMCOSZA = 1.0
345 //
346 // a: array of cut parameters
347 // ls: log(SIZE) - dNOMLOGSIZE
348 // ls2: ls^2
349 // ct: Cos(ZA.) - dNOMCOSZA
350 // dd2: DIST^2
351 const Double_t limit =
352 a[0] + a[1] * dd2 + a[2] * ct +
353 ls * (a[3] + a[4] * dd2 + a[5] * ct) +
354 ls2 * (a[6] + a[7] * dd2);
355
356 //*fLog << "MCT1SupercutsCalc::CtsMCut; *a = "
357 // << *a << ", " << *(a+1) << ", " << *(a+2) << ", "
358 // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", "
359 // << *(a+6) << ", " << *(a+7) << endl;
360
361 //*fLog << "MCT1SupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
362 // << ", " << ls2 << ", " << ct << ", " << dd2 << ", "
363 // << limit << endl;
364
365 return limit;
366}
367
368// --------------------------------------------------------------------------
369//
370// Returns the mapped value from the Matrix
371//
372Double_t MCT1SupercutsCalc::GetVal(Int_t i) const
373{
374 return (*fMatrix)[fMap[i]];
375}
376
377// --------------------------------------------------------------------------
378//
379// You can use this function if you want to use a MHMatrix instead of the
380// given containers. This function adds all necessary columns to the
381// given matrix. Afterward you should fill the matrix with the corresponding
382// data (eg from a file by using MHMatrix::Fill). If you now loop
383// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
384// will take the values from the matrix instead of the containers.
385//
386void MCT1SupercutsCalc::InitMapping(MHMatrix *mat)
387{
388 if (fMatrix)
389 return;
390
391 fMatrix = mat;
392
393 fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta");
394 fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
395 fMap[2] = fMatrix->AddColumn("MHillas.fLength");
396 fMap[3] = fMatrix->AddColumn("MHillas.fSize");
397 fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
398 fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
399 fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
400 fMap[7] = fMatrix->AddColumn("MHillasSrc.fAlpha");
401}
402
403// ---------------------------------------------------------------------------
404//
405// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich)
406// optimized for mkn 421 2001 data
407//
408// set hadronness to 0.25 if cuts are fullfilled
409// 0.75 otherwise
410//
411Int_t MCT1SupercutsCalc::Process()
412{
413 const Double_t kNomLogSize = 4.1;
414 const Double_t kNomCosZA = 1.0;
415
416 const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
417 const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth();
418 const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
419 const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize();
420 const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX();
421 const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY();
422 const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist();
423
424 const Double_t newdist = dist0 * fMm2Deg;
425
426 const Double_t dist2 = meanx*meanx + meany*meany;
427 const Double_t dist = sqrt(dist2) * fMm2Deg;
428 const Double_t dd2 = dist*dist;
429
430
431 const Double_t dmls = log(size) - kNomLogSize;
432 const Double_t dmls2 = dmls * dmls;
433
434 const Double_t dmcza = cos(theta) - kNomCosZA;
435
436 const Double_t length = length0 * fMm2Deg;
437 const Double_t width = width0 * fMm2Deg;
438
439 if (newdist < 1.05 &&
440 newdist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) &&
441 newdist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) &&
442 dist < 1.05 &&
443 length < CtsMCut (fLengthUp, dmls, dmcza, dmls2, dd2) &&
444 length > CtsMCut (fLengthLo, dmls, dmcza, dmls2, dd2) &&
445 width < CtsMCut (fWidthUp, dmls, dmcza, dmls2, dd2) &&
446 width > CtsMCut (fWidthLo, dmls, dmcza, dmls2, dd2) &&
447 //asym < CtsMCut (asymup, dmls, dmcza, dmls2, dd2) &&
448 //asym > CtsMCut (asymlow, dmls, dmcza, dmls2, dd2) &&
449 dist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) &&
450 dist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) )
451 fHadronness->SetHadronness(0.25);
452 else
453 fHadronness->SetHadronness(0.75);
454
455 fHadronness->SetReadyToSave();
456
457 return kTRUE;
458}
459//==========================================================================
Note: See TracBrowser for help on using the repository browser.