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 | ! David Paneque, 02/2004 <mailto:dpaneque@mppmu.mpg.de>
|
---|
20 | !
|
---|
21 | ! Copyright: MAGIC Software Development, 2000-2004
|
---|
22 | !
|
---|
23 | !
|
---|
24 | \* ======================================================================== */
|
---|
25 |
|
---|
26 | /////////////////////////////////////////////////////////////////////////////
|
---|
27 | // //
|
---|
28 | // MSupercutsCalcONOFF //
|
---|
29 | // //
|
---|
30 | // this class calculates the hadronness for the supercuts //
|
---|
31 | // the parameters of the supercuts are taken //
|
---|
32 | // from the container MSupercuts //
|
---|
33 | // //
|
---|
34 | // //
|
---|
35 | /////////////////////////////////////////////////////////////////////////////
|
---|
36 | #include "MSupercutsCalcONOFF.h"
|
---|
37 |
|
---|
38 | #include <math.h>
|
---|
39 | #include <fstream>
|
---|
40 |
|
---|
41 | #include "TFile.h"
|
---|
42 | #include "TArrayD.h"
|
---|
43 |
|
---|
44 | #include "MParList.h"
|
---|
45 | #include "MHillasExt.h"
|
---|
46 | #include "MHillasSrc.h"
|
---|
47 | #include "MNewImagePar.h"
|
---|
48 | #include "MPointingPos.h"
|
---|
49 | #include "MCerPhotEvt.h"
|
---|
50 | #include "MGeomCam.h"
|
---|
51 | #include "MHadronness.h"
|
---|
52 | #include "MTSupercutsApplied.h"
|
---|
53 | #include "MHMatrix.h"
|
---|
54 | #include "MSupercuts.h"
|
---|
55 |
|
---|
56 | #include "MLog.h"
|
---|
57 | #include "MLogManip.h"
|
---|
58 |
|
---|
59 | ClassImp(MSupercutsCalcONOFF);
|
---|
60 |
|
---|
61 | using namespace std;
|
---|
62 |
|
---|
63 |
|
---|
64 | // --------------------------------------------------------------------------
|
---|
65 | //
|
---|
66 | // constructor
|
---|
67 | //
|
---|
68 |
|
---|
69 | MSupercutsCalcONOFF::MSupercutsCalcONOFF(const char *hilname,
|
---|
70 | const char *hilsrcname,
|
---|
71 | const char *name, const char *title)
|
---|
72 | : fHadronnessName("MHadronness"), fSupercutsAppliedName("MSupercutsApplied"),
|
---|
73 | fHilName(hilname), fHilSrcName(hilsrcname),
|
---|
74 | fHilExtName("MHillasExt"), fNewParName("MNewImagePar"),
|
---|
75 | fSuperName("MSupercuts")
|
---|
76 | {
|
---|
77 | fName = name ? name : "MSupercutsCalcONOFF";
|
---|
78 | fTitle = title ? title : "Class to evaluate the Supercuts";
|
---|
79 |
|
---|
80 | fMatrix = NULL;
|
---|
81 |
|
---|
82 | fStoreAppliedSupercuts = kFALSE; // by default, applied SC parameters are not stored
|
---|
83 |
|
---|
84 | fNotUseTheta = kFALSE; // by default, theta info is used in the computation of the cuts
|
---|
85 |
|
---|
86 | fUseStaticCuts = kFALSE; // by default, dynamical cuts are used
|
---|
87 |
|
---|
88 |
|
---|
89 | // Usage of DIST parameter in the parameterization of the cuts.
|
---|
90 | // For the time being is in the constructor. If finally it comes out
|
---|
91 | // that it is important to disable the DIST parameter from the
|
---|
92 | // cut parameterization I will make a function that access this
|
---|
93 | // data variable from outside the class.
|
---|
94 | fUseDist = kTRUE;
|
---|
95 |
|
---|
96 |
|
---|
97 |
|
---|
98 | // OFFSETS FOR THE DYNAMICAL CUTS
|
---|
99 | // Values of Size (photons), Dist (degrees) and Theta (degrees)
|
---|
100 | // for which the value of the dynamical cut is equal to the
|
---|
101 | // non dependent parameter; i.e. to the static cut.
|
---|
102 | // By adjusting these offsets the user can set the value in size,
|
---|
103 | // dist and theta for which the dynamical cuts will be given by
|
---|
104 | // the term that DOES not depend on size, dist and theta.
|
---|
105 |
|
---|
106 | // For the time being, these quantities are set in the constructor.
|
---|
107 | // In future, if they show to be useful, I will make them available
|
---|
108 | // as external variables.
|
---|
109 |
|
---|
110 | fSizeOffset = 3000; // photons
|
---|
111 | fDistOffset = 0.9; // degrees
|
---|
112 | fThetaOffset = 20; // degrees NOT USED FOR THE TIME BEING
|
---|
113 |
|
---|
114 |
|
---|
115 |
|
---|
116 |
|
---|
117 |
|
---|
118 |
|
---|
119 | // Variables defining upper limits for some of the hillas params.
|
---|
120 | // The default values are set to conservative
|
---|
121 | // values so that no gamma showers are removed
|
---|
122 |
|
---|
123 | // If a cut value computed by function MSupercutsCalcONOFF::CtsMCut,
|
---|
124 | // (i.e. widthlow) exceeds (larger or smaller
|
---|
125 | // depending on wether is cutUP or cutLOW) one of these values
|
---|
126 | // (i.e. fWidthLowerLimit), such cut value is replaced by the
|
---|
127 | // the corresponding Limit (i.e. widthlow = fWidthLowerLimit)
|
---|
128 |
|
---|
129 |
|
---|
130 |
|
---|
131 | fDistUpperLimit = 1.4; // in deg
|
---|
132 | fLengthUpperLimit = 0.6;
|
---|
133 | fWidthUpperLimit = 0.4;
|
---|
134 |
|
---|
135 | fLeakage1UpperLimit = 0.25;
|
---|
136 |
|
---|
137 | fLengthOverWidthUpperLimit = 1.0;
|
---|
138 |
|
---|
139 | fDistLowerLimit = 0.1;
|
---|
140 | fLengthLowerLimit = 0.09;
|
---|
141 | fWidthLowerLimit = 0.05;
|
---|
142 |
|
---|
143 |
|
---|
144 |
|
---|
145 |
|
---|
146 |
|
---|
147 |
|
---|
148 | }
|
---|
149 |
|
---|
150 | // --------------------------------------------------------------------------
|
---|
151 | //
|
---|
152 |
|
---|
153 |
|
---|
154 | void MSupercutsCalcONOFF::SetStoreAppliedSupercuts(Bool_t b)
|
---|
155 | {
|
---|
156 | fStoreAppliedSupercuts = b;
|
---|
157 | if (fStoreAppliedSupercuts)
|
---|
158 | {
|
---|
159 | *fLog << "Supercuts applied to all the individual events will be stored "
|
---|
160 | << "in a container of the class MSupercutsApplied." << endl;
|
---|
161 | }
|
---|
162 |
|
---|
163 | }
|
---|
164 |
|
---|
165 |
|
---|
166 | void MSupercutsCalcONOFF::SetVariableNotUseTheta(Bool_t b)
|
---|
167 | {
|
---|
168 | fNotUseTheta = b;
|
---|
169 | if (fNotUseTheta)
|
---|
170 | {
|
---|
171 | *fLog << "Theta variable is NOT used in the computation of the the "
|
---|
172 | << "dynamical cuts" << endl;
|
---|
173 | }
|
---|
174 | else
|
---|
175 | {
|
---|
176 | *fLog << "Theta variable is used in the computation of the the "
|
---|
177 | << "dynamical cuts" << endl;
|
---|
178 | }
|
---|
179 |
|
---|
180 | }
|
---|
181 |
|
---|
182 |
|
---|
183 | void MSupercutsCalcONOFF::SetVariableUseStaticCuts(Bool_t b)
|
---|
184 | {
|
---|
185 | fUseStaticCuts = b;
|
---|
186 | if (fUseStaticCuts)
|
---|
187 | {
|
---|
188 | *fLog << "Supercuts DO NOT take into account Theta, Size and Dist dependence; "
|
---|
189 | << "i.e, they are STATIC CUTS." << endl;
|
---|
190 | }
|
---|
191 | else
|
---|
192 | {
|
---|
193 | *fLog << "Supercuts are computed taking into account Theta, size and Dist dependence."
|
---|
194 | << "i.e., they are DYNAMICAL CUTS" << endl;
|
---|
195 |
|
---|
196 | }
|
---|
197 |
|
---|
198 |
|
---|
199 |
|
---|
200 | }
|
---|
201 |
|
---|
202 |
|
---|
203 | void MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits(Double_t distup,
|
---|
204 | Double_t lengthup,
|
---|
205 | Double_t widthup,
|
---|
206 | Double_t distlow,
|
---|
207 | Double_t lengthlow,
|
---|
208 | Double_t widthlow)
|
---|
209 | {
|
---|
210 |
|
---|
211 | fDistUpperLimit = distup;
|
---|
212 | fLengthUpperLimit = lengthup;
|
---|
213 | fWidthUpperLimit = widthup;
|
---|
214 |
|
---|
215 |
|
---|
216 | fDistLowerLimit = distlow;
|
---|
217 | fLengthLowerLimit = lengthlow;
|
---|
218 | fWidthLowerLimit = widthlow;
|
---|
219 |
|
---|
220 |
|
---|
221 | *fLog << "MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits" << endl
|
---|
222 | << "Upper limits to Hillas parameters Dist, Length and Width "
|
---|
223 | << "are set respectively to : "
|
---|
224 | << fDistUpperLimit << ", "
|
---|
225 | << fLengthUpperLimit << ", "
|
---|
226 | << fWidthUpperLimit << " (in degrees!!)" << endl
|
---|
227 |
|
---|
228 | << "Lower limits to Hillas parameters Dist, Length and Width "
|
---|
229 | << "are set respectively to : "
|
---|
230 | << fDistLowerLimit << ", "
|
---|
231 | << fLengthLowerLimit << ", "
|
---|
232 | << fWidthLowerLimit << " (in degrees!!)" << endl
|
---|
233 |
|
---|
234 | << "The Hadronnes of those images exceeding one of these limits will "
|
---|
235 | << "be set to 0.75." << endl << endl;
|
---|
236 |
|
---|
237 |
|
---|
238 | }
|
---|
239 |
|
---|
240 |
|
---|
241 |
|
---|
242 | Int_t MSupercutsCalcONOFF::PreProcess(MParList *pList)
|
---|
243 | {
|
---|
244 | MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
|
---|
245 | if (!cam)
|
---|
246 | {
|
---|
247 | *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
|
---|
248 | return kFALSE;
|
---|
249 | }
|
---|
250 |
|
---|
251 | fMm2Deg = cam->GetConvMm2Deg();
|
---|
252 |
|
---|
253 | fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
|
---|
254 | if (!fHadronness)
|
---|
255 | {
|
---|
256 | *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
|
---|
257 | return kFALSE;
|
---|
258 | }
|
---|
259 |
|
---|
260 | fSuper = (MSupercuts*)pList->FindObject(fSuperName, "MSupercuts");
|
---|
261 | if (!fSuper)
|
---|
262 | {
|
---|
263 | *fLog << err << fSuperName << " [MSupercuts] not found... aborting." << endl;
|
---|
264 | return kFALSE;
|
---|
265 | }
|
---|
266 |
|
---|
267 | if (fStoreAppliedSupercuts)
|
---|
268 | {
|
---|
269 |
|
---|
270 | fSupercutsApplied = (MTSupercutsApplied*)pList->FindObject(fSupercutsAppliedName,
|
---|
271 | "MTSupercutsApplied");
|
---|
272 | if(!fSupercutsApplied)
|
---|
273 | {
|
---|
274 | *fLog << err << fSupercutsAppliedName
|
---|
275 | << " [MTSupercutsApplied] not found... aborting." << endl;
|
---|
276 | return kFALSE;
|
---|
277 | }
|
---|
278 |
|
---|
279 |
|
---|
280 | }
|
---|
281 |
|
---|
282 | if (fMatrix)
|
---|
283 | return kTRUE;
|
---|
284 |
|
---|
285 | //-----------------------------------------------------------
|
---|
286 | fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
|
---|
287 | if (!fHil)
|
---|
288 | {
|
---|
289 | *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
|
---|
290 | return kFALSE;
|
---|
291 | }
|
---|
292 |
|
---|
293 | fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt");
|
---|
294 | if (!fHilExt)
|
---|
295 | {
|
---|
296 | *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl;
|
---|
297 | return kFALSE;
|
---|
298 | }
|
---|
299 |
|
---|
300 | fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
|
---|
301 | if (!fHilSrc)
|
---|
302 | {
|
---|
303 | *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
|
---|
304 | return kFALSE;
|
---|
305 | }
|
---|
306 |
|
---|
307 | fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar");
|
---|
308 | if (!fNewPar)
|
---|
309 | {
|
---|
310 | *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl;
|
---|
311 | return kFALSE;
|
---|
312 | }
|
---|
313 |
|
---|
314 | fPointPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
|
---|
315 | if (!fPointPos)
|
---|
316 | {
|
---|
317 | *fLog << err << "MSupercutsCalcONOFF::PreProcess; MPointingPos not found... aborting." << endl;
|
---|
318 | return kFALSE;
|
---|
319 | }
|
---|
320 |
|
---|
321 |
|
---|
322 |
|
---|
323 | return kTRUE;
|
---|
324 | }
|
---|
325 |
|
---|
326 | // --------------------------------------------------------------------------
|
---|
327 | //
|
---|
328 | // Calculation of upper and lower limits
|
---|
329 | //
|
---|
330 | Double_t MSupercutsCalcONOFF::CtsMCut(const Double_t* a, Double_t ls, Double_t ct,
|
---|
331 | Double_t ls2, Double_t dd2) const
|
---|
332 | {
|
---|
333 | // define cut-function
|
---|
334 | //
|
---|
335 | //
|
---|
336 | // dNOMCOSZA = 1.0
|
---|
337 | //
|
---|
338 | // a: array of cut parameters
|
---|
339 | // ls: log(SIZE) - log(fSizeOffset)
|
---|
340 | // ls2: ls^2
|
---|
341 | // ct: Cos(ZA.) - Cos(fThetaOffset)
|
---|
342 | // dd2: (DIST- fDistOffset)^2
|
---|
343 |
|
---|
344 | Double_t limit;
|
---|
345 |
|
---|
346 | if(fUseStaticCuts)
|
---|
347 | { // static cuts are used
|
---|
348 | limit = a[0];
|
---|
349 | }
|
---|
350 | else
|
---|
351 | {
|
---|
352 | if (fNotUseTheta)
|
---|
353 | { // Theta info is NOT used in the computation of the dynamical cuts
|
---|
354 | // For the time being dist info will not be used
|
---|
355 |
|
---|
356 | if(fUseDist)
|
---|
357 | {
|
---|
358 |
|
---|
359 |
|
---|
360 | limit = a[0] + a[1] * dd2 +
|
---|
361 | ls * (a[3] + a[4] * dd2) +
|
---|
362 | ls2 * (a[6] + a[7] * dd2);
|
---|
363 |
|
---|
364 |
|
---|
365 | }
|
---|
366 | else
|
---|
367 | {
|
---|
368 | limit = a[0] + ls * (a[3])+ ls2 * (a[6]);
|
---|
369 | }
|
---|
370 |
|
---|
371 |
|
---|
372 | }
|
---|
373 | else
|
---|
374 | {// Theta info IS used in the computation of the dynamical cuts
|
---|
375 |
|
---|
376 | limit =
|
---|
377 | a[0] + a[1] * dd2 + a[2] * ct +
|
---|
378 | ls * (a[3] + a[4] * dd2 + a[5] * ct) +
|
---|
379 | ls2 * (a[6] + a[7] * dd2);
|
---|
380 | }
|
---|
381 |
|
---|
382 | }
|
---|
383 |
|
---|
384 |
|
---|
385 |
|
---|
386 | //*fLog << "MSupercutsCalcONOFF::CtsMCut; *a = "
|
---|
387 | // << *a << ", " << *(a+1) << ", " << *(a+2) << ", "
|
---|
388 | // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", "
|
---|
389 | // << *(a+6) << ", " << *(a+7) << endl;
|
---|
390 |
|
---|
391 | //*fLog << "MSupercutsCalcONOFF::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
|
---|
392 | // << ", " << ls2 << ", " << ct << ", " << dd2 << ", "
|
---|
393 | // << limit << endl;
|
---|
394 |
|
---|
395 | return limit;
|
---|
396 | }
|
---|
397 |
|
---|
398 | // --------------------------------------------------------------------------
|
---|
399 | //
|
---|
400 | // Returns the mapped value from the Matrix
|
---|
401 | //
|
---|
402 | Double_t MSupercutsCalcONOFF::GetVal(Int_t i) const
|
---|
403 | {
|
---|
404 |
|
---|
405 | Double_t val = (*fMatrix)[fMap[i]];
|
---|
406 |
|
---|
407 |
|
---|
408 | //*fLog << "MSupercutsCalcONOFF::GetVal; i, fMatrix, fMap, val = "
|
---|
409 | // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl;
|
---|
410 |
|
---|
411 |
|
---|
412 | return val;
|
---|
413 | }
|
---|
414 |
|
---|
415 | // --------------------------------------------------------------------------
|
---|
416 | //
|
---|
417 | // You can use this function if you want to use a MHMatrix instead of the
|
---|
418 | // given containers. This function adds all necessary columns to the
|
---|
419 | // given matrix. Afterward you should fill the matrix with the corresponding
|
---|
420 | // data (eg from a file by using MHMatrix::Fill). If you now loop
|
---|
421 | // through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
|
---|
422 | // will take the values from the matrix instead of the containers.
|
---|
423 | //
|
---|
424 |
|
---|
425 |
|
---|
426 | void MSupercutsCalcONOFF::InitMapping(MHMatrix *mat)
|
---|
427 | {
|
---|
428 | if (fMatrix)
|
---|
429 | return;
|
---|
430 |
|
---|
431 | fMatrix = mat;
|
---|
432 |
|
---|
433 | fMap[0] = fMatrix->AddColumn("MPointingPos.fZd"); //deg
|
---|
434 | fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
|
---|
435 | fMap[2] = fMatrix->AddColumn("MHillas.fLength");
|
---|
436 | fMap[3] = fMatrix->AddColumn("MHillas.fSize");
|
---|
437 | fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
|
---|
438 | fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
|
---|
439 | fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
|
---|
440 | fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)");
|
---|
441 | fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
|
---|
442 | fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc");
|
---|
443 | fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1");
|
---|
444 | }
|
---|
445 |
|
---|
446 |
|
---|
447 | // ---------------------------------------------------------------------------
|
---|
448 | //
|
---|
449 | // Evaluate dynamical supercuts
|
---|
450 | //
|
---|
451 | // set hadronness to 0.25 if cuts are fullfilled
|
---|
452 | // 0.75 otherwise
|
---|
453 | //
|
---|
454 | Int_t MSupercutsCalcONOFF::Process()
|
---|
455 | {
|
---|
456 | // const Double_t kNomLogSize = 4.1;
|
---|
457 | const Double_t kNomCosZA = 1.0;
|
---|
458 |
|
---|
459 | //const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
|
---|
460 | // For the time being I do it in a way that theta must be
|
---|
461 | // the value used in the row 0 of fMatrix. That means that if there is
|
---|
462 | // no matrix, the program will complain (and crash). And tehn I will know
|
---|
463 | // that value 0 (supposed to be fThetaOrig.Val) could not be taken.
|
---|
464 | const Double_t theta = GetVal(0);
|
---|
465 | const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth();
|
---|
466 | const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
|
---|
467 | const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize();
|
---|
468 | const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX();
|
---|
469 | const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY();
|
---|
470 | const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist();
|
---|
471 |
|
---|
472 | const Double_t alpha = fMatrix ? GetVal(7) : fHilSrc->GetAlpha();
|
---|
473 |
|
---|
474 | Double_t help;
|
---|
475 | if (!fMatrix)
|
---|
476 | help = TMath::Sign(fHilExt->GetM3Long(),
|
---|
477 | fHilSrc->GetCosDeltaAlpha());
|
---|
478 | const Double_t asym0 = fMatrix ? GetVal(8) : help;
|
---|
479 | const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc();
|
---|
480 | const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1();
|
---|
481 |
|
---|
482 | const Double_t newdist = dist0 * fMm2Deg;
|
---|
483 |
|
---|
484 | const Double_t dist2 = meanx*meanx + meany*meany;
|
---|
485 | const Double_t dist = sqrt(dist2) * fMm2Deg;
|
---|
486 |
|
---|
487 | // const Double_t dd2 = dist*dist;
|
---|
488 |
|
---|
489 |
|
---|
490 | // in the parameterization of the cuts
|
---|
491 | // the dist parameter used is the one computed from the source position,
|
---|
492 | // and not the one computed from the camera center.
|
---|
493 |
|
---|
494 | // Actually the value used in the parameterization is newdist^2,
|
---|
495 | // minus the offset_in_dist^2
|
---|
496 |
|
---|
497 | const Double_t dd2 = newdist*newdist - fDistOffset*fDistOffset;
|
---|
498 |
|
---|
499 |
|
---|
500 |
|
---|
501 | const Double_t dmls = log(size) - log(fSizeOffset);
|
---|
502 | const Double_t dmls2 = dmls * dmls;
|
---|
503 |
|
---|
504 | const Double_t dmcza = cos(theta) - kNomCosZA;
|
---|
505 |
|
---|
506 | const Double_t length = length0 * fMm2Deg;
|
---|
507 | const Double_t width = width0 * fMm2Deg;
|
---|
508 | const Double_t asym = asym0 * fMm2Deg;
|
---|
509 |
|
---|
510 |
|
---|
511 |
|
---|
512 | // computation of the cut limits
|
---|
513 |
|
---|
514 | Double_t lengthup = CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2);
|
---|
515 | Double_t lengthlow = CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2);
|
---|
516 |
|
---|
517 | Double_t widthup = CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2);
|
---|
518 | Double_t widthlow = CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2);
|
---|
519 |
|
---|
520 | Double_t distup = CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2);
|
---|
521 | Double_t distlow = CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2);
|
---|
522 |
|
---|
523 | Double_t asymup = CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2);
|
---|
524 | Double_t asymlow = CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2);
|
---|
525 |
|
---|
526 | Double_t concup = CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2);
|
---|
527 | Double_t conclow = CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2);
|
---|
528 |
|
---|
529 | Double_t leakageup = CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2);
|
---|
530 | Double_t leakagelow = CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2);
|
---|
531 |
|
---|
532 |
|
---|
533 | Double_t lengthoverwidth = length/width;
|
---|
534 |
|
---|
535 | Double_t hadronness = 0.0;
|
---|
536 |
|
---|
537 |
|
---|
538 | // If the cut values computed before for Dist, Length and Width exceed the upper limits set
|
---|
539 | // by the user through function MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits,
|
---|
540 | // (or the default values defined in constructor); such cut values are replaced by the
|
---|
541 | // the limit values
|
---|
542 |
|
---|
543 |
|
---|
544 |
|
---|
545 | if (distup > fDistUpperLimit)
|
---|
546 | {
|
---|
547 | distup = fDistUpperLimit;
|
---|
548 | }
|
---|
549 |
|
---|
550 |
|
---|
551 | if (lengthup > fLengthUpperLimit)
|
---|
552 | {
|
---|
553 | lengthup = fLengthUpperLimit;
|
---|
554 | }
|
---|
555 |
|
---|
556 | if (widthup > fWidthUpperLimit)
|
---|
557 | {
|
---|
558 | widthup = fWidthUpperLimit;
|
---|
559 | }
|
---|
560 |
|
---|
561 |
|
---|
562 |
|
---|
563 | if (distlow < fDistLowerLimit)
|
---|
564 | {
|
---|
565 | distlow = fDistLowerLimit;
|
---|
566 | }
|
---|
567 |
|
---|
568 | if (lengthlow < fLengthLowerLimit)
|
---|
569 | {
|
---|
570 | lengthlow = fLengthLowerLimit;
|
---|
571 | }
|
---|
572 |
|
---|
573 | if (widthlow < fWidthLowerLimit)
|
---|
574 | {
|
---|
575 | widthlow = fWidthLowerLimit;
|
---|
576 | }
|
---|
577 |
|
---|
578 |
|
---|
579 | if (leakageup > fLeakage1UpperLimit)
|
---|
580 | {
|
---|
581 | leakageup = fLeakage1UpperLimit;
|
---|
582 | }
|
---|
583 |
|
---|
584 |
|
---|
585 |
|
---|
586 |
|
---|
587 |
|
---|
588 | // Upper cut in leakage 1 also implemented
|
---|
589 |
|
---|
590 |
|
---|
591 |
|
---|
592 |
|
---|
593 | if (//
|
---|
594 | newdist > distup ||
|
---|
595 | newdist < distlow ||
|
---|
596 | length > lengthup ||
|
---|
597 | length < lengthlow ||
|
---|
598 | width > widthup ||
|
---|
599 | width < widthlow ||
|
---|
600 | leakage > leakageup ||
|
---|
601 | leakage < leakagelow ||
|
---|
602 | lengthoverwidth < fLengthOverWidthUpperLimit
|
---|
603 | //
|
---|
604 | //asym < asymup &&
|
---|
605 | //asym > asymlow &&
|
---|
606 |
|
---|
607 | //dist < distup &&
|
---|
608 | //dist > distlow &&
|
---|
609 |
|
---|
610 | //conc < concup &&
|
---|
611 | //conc > conclow &&
|
---|
612 | //
|
---|
613 | )
|
---|
614 |
|
---|
615 | {hadronness = 0.75;}
|
---|
616 | else
|
---|
617 | {
|
---|
618 | hadronness = 0.25;
|
---|
619 | }
|
---|
620 |
|
---|
621 |
|
---|
622 | fHadronness->SetHadronness(hadronness);
|
---|
623 | fHadronness->SetReadyToSave();
|
---|
624 |
|
---|
625 |
|
---|
626 | if(fStoreAppliedSupercuts)
|
---|
627 | {
|
---|
628 |
|
---|
629 | // TArrayD vector with the shower parameters (matching the ones
|
---|
630 | // specified in function MTSupercutsApplied::CreateTreeBranches)
|
---|
631 | // is created and filled with this event features
|
---|
632 |
|
---|
633 | // Eventually this definition might be done from outside this
|
---|
634 | // class, allowing for adding removing parameters without
|
---|
635 | // recompiling mars. yet for the time being, let's keep it simple.
|
---|
636 |
|
---|
637 | TArrayD ShowerParams(10);
|
---|
638 | ShowerParams[0] = length;
|
---|
639 | ShowerParams[1] = width;
|
---|
640 | ShowerParams[2] = dist;
|
---|
641 | ShowerParams[3] = newdist;
|
---|
642 | ShowerParams[4] = asym;
|
---|
643 | ShowerParams[5] = conc;
|
---|
644 | ShowerParams[6] = leakage;
|
---|
645 | ShowerParams[7] = theta;
|
---|
646 | ShowerParams[8] = size;
|
---|
647 | ShowerParams[9] = alpha;
|
---|
648 |
|
---|
649 |
|
---|
650 | // TArrayD vector with the cut parameters (matching the ones
|
---|
651 | // specified in function MTSupercutsApplied::CreateTreeBranches)
|
---|
652 | // is created and filled with this event features
|
---|
653 |
|
---|
654 | // Eventually this definition might be done from outside this
|
---|
655 | // class, allowing for adding removing parameters without
|
---|
656 | // recompiling mars. yet for the time being, let's keep it simple.
|
---|
657 |
|
---|
658 |
|
---|
659 | TArrayD SuperCutParams(13);
|
---|
660 | SuperCutParams[0] = lengthup;
|
---|
661 | SuperCutParams[1] = lengthlow;
|
---|
662 | SuperCutParams[2] = widthup;
|
---|
663 | SuperCutParams[3] = widthlow;
|
---|
664 | SuperCutParams[4] = distup;
|
---|
665 | SuperCutParams[5] = distlow;
|
---|
666 | SuperCutParams[6] = asymup;
|
---|
667 | SuperCutParams[7] = asymlow;
|
---|
668 | SuperCutParams[8] = concup;
|
---|
669 | SuperCutParams[9] = conclow;
|
---|
670 | SuperCutParams[10] = leakageup;
|
---|
671 | SuperCutParams[11] = leakagelow;
|
---|
672 | SuperCutParams[12] = hadronness;
|
---|
673 |
|
---|
674 | // SC parameters applied to this event, as well as the
|
---|
675 | // shower parameters, are stored in the branches of the
|
---|
676 | // TTree object of a MTSupercutsApplied object
|
---|
677 |
|
---|
678 |
|
---|
679 | if (!fSupercutsApplied ->FillTreeBranches(SuperCutParams, ShowerParams))
|
---|
680 | {
|
---|
681 | *fLog << "MSupercutsCalcONOFF::Process()" << endl
|
---|
682 | << "Supercuts applied could not be stored in tree..."
|
---|
683 | << endl;
|
---|
684 |
|
---|
685 | return kFALSE;
|
---|
686 | }
|
---|
687 |
|
---|
688 |
|
---|
689 | }
|
---|
690 |
|
---|
691 |
|
---|
692 |
|
---|
693 | return kTRUE;
|
---|
694 |
|
---|
695 |
|
---|
696 | // OLD STUFF
|
---|
697 |
|
---|
698 | /*
|
---|
699 | *fLog << "newdist, length, width, asym, dist, conc, leakage = "
|
---|
700 | << newdist << ", " << length << ", " << width << ", "
|
---|
701 | << asym << ", " << dist << ", " << conc << ", " << leakage
|
---|
702 | << endl;
|
---|
703 |
|
---|
704 | *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = "
|
---|
705 | << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
706 | << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
707 |
|
---|
708 | << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
709 | << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
710 |
|
---|
711 | << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
712 | << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
713 |
|
---|
714 | << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
715 | << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
716 |
|
---|
717 | << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
718 | << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
719 |
|
---|
720 | << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
721 | << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
722 |
|
---|
723 | << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
724 | << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", "
|
---|
725 | << endl;
|
---|
726 | */
|
---|
727 |
|
---|
728 | /*
|
---|
729 | if (
|
---|
730 | //dist < 1.05 &&
|
---|
731 | //newdist < 1.05 &&
|
---|
732 |
|
---|
733 | newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) &&
|
---|
734 | newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) &&
|
---|
735 |
|
---|
736 | length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) &&
|
---|
737 | length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) &&
|
---|
738 |
|
---|
739 | width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) &&
|
---|
740 | width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) &&
|
---|
741 |
|
---|
742 | asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) &&
|
---|
743 | asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) &&
|
---|
744 |
|
---|
745 | dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) &&
|
---|
746 | dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) &&
|
---|
747 |
|
---|
748 | conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) &&
|
---|
749 | conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) &&
|
---|
750 |
|
---|
751 | leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) &&
|
---|
752 | leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) )
|
---|
753 |
|
---|
754 | fHadronness->SetHadronness(0.25);
|
---|
755 | else
|
---|
756 | fHadronness->SetHadronness(0.75);
|
---|
757 |
|
---|
758 |
|
---|
759 |
|
---|
760 | // *fLog << "SChadroness = " << fHadronness->GetHadronness() << endl;
|
---|
761 |
|
---|
762 | fHadronness->SetReadyToSave();
|
---|
763 | */
|
---|
764 |
|
---|
765 | // END OF OLD STUFF
|
---|
766 | }
|
---|
767 |
|
---|
768 | /*
|
---|
769 | Bool_t MSupercutsCalcONOFF::StoreSupercutsAppliedToThisEvent(TArrayD CutParams,
|
---|
770 | TArrayD ShowerParams)
|
---|
771 | {
|
---|
772 | // SC parameters applied to this event are stored in
|
---|
773 | // TTree object of MTSupercutsApplied object
|
---|
774 |
|
---|
775 |
|
---|
776 | if (!fSupercutsApplied ->FillTreeBranches(CutParams, ShowerParams))
|
---|
777 | {
|
---|
778 | *fLog << "MSupercutsCalcONOFF::StoreSupercutsAppliedToThisEvent" << endl
|
---|
779 | << "Supercuts applied could not be stored in tree..."
|
---|
780 | << endl;
|
---|
781 |
|
---|
782 | return kFALSE;
|
---|
783 | }
|
---|
784 |
|
---|
785 |
|
---|
786 | return kTRUE;
|
---|
787 | }
|
---|
788 |
|
---|
789 | */
|
---|
790 |
|
---|
791 |
|
---|
792 |
|
---|
793 |
|
---|
794 |
|
---|
795 |
|
---|
796 |
|
---|
797 |
|
---|