source: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.cc@ 6232

Last change on this file since 6232 was 6232, checked in by mazin, 20 years ago
*** empty log message ***
File size: 25.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! 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
59ClassImp(MSupercutsCalcONOFF);
60
61using namespace std;
62
63
64// --------------------------------------------------------------------------
65//
66// constructor
67//
68
69MSupercutsCalcONOFF::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 // Temp upper limit to get rid of triangular events
136 fLog10ConcUpperLimit = -0.35;
137
138 fLeakage1UpperLimit = 0.25;
139
140 fLengthOverWidthUpperLimit = 1.0;
141
142 fDistLowerLimit = 0.1;
143 fLengthLowerLimit = 0.05; // values at October 21th 2004
144 fWidthLowerLimit = 0.03; // values at October 21th 2004
145
146
147
148
149
150
151}
152
153// --------------------------------------------------------------------------
154//
155
156
157void MSupercutsCalcONOFF::SetStoreAppliedSupercuts(Bool_t b)
158{
159 fStoreAppliedSupercuts = b;
160 if (fStoreAppliedSupercuts)
161 {
162 *fLog << "Supercuts applied to all the individual events will be stored "
163 << "in a container of the class MSupercutsApplied." << endl;
164 }
165
166}
167
168
169void MSupercutsCalcONOFF::SetVariableNotUseTheta(Bool_t b)
170{
171 fNotUseTheta = b;
172 if (fNotUseTheta)
173 {
174 *fLog << "Theta variable is NOT used in the computation of the the "
175 << "dynamical cuts" << endl;
176 }
177 else
178 {
179 *fLog << "Theta variable is used in the computation of the the "
180 << "dynamical cuts" << endl;
181 }
182
183}
184
185
186void MSupercutsCalcONOFF::SetVariableUseStaticCuts(Bool_t b)
187{
188 fUseStaticCuts = b;
189 if (fUseStaticCuts)
190 {
191 *fLog << "Supercuts DO NOT take into account Theta, Size and Dist dependence; "
192 << "i.e, they are STATIC CUTS." << endl;
193 }
194 else
195 {
196 *fLog << "Supercuts are computed taking into account Theta, size and Dist dependence."
197 << "i.e., they are DYNAMICAL CUTS" << endl;
198
199 }
200
201
202
203}
204
205
206void MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits(Double_t distup,
207 Double_t lengthup,
208 Double_t widthup,
209 Double_t distlow,
210 Double_t lengthlow,
211 Double_t widthlow)
212{
213
214 fDistUpperLimit = distup;
215 fLengthUpperLimit = lengthup;
216 fWidthUpperLimit = widthup;
217
218
219 fDistLowerLimit = distlow;
220 fLengthLowerLimit = lengthlow;
221 fWidthLowerLimit = widthlow;
222
223
224 *fLog << "MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits" << endl
225 << "Upper limits to Hillas parameters Dist, Length and Width "
226 << "are set respectively to : "
227 << fDistUpperLimit << ", "
228 << fLengthUpperLimit << ", "
229 << fWidthUpperLimit << " (in degrees!!)" << endl
230
231 << "Lower limits to Hillas parameters Dist, Length and Width "
232 << "are set respectively to : "
233 << fDistLowerLimit << ", "
234 << fLengthLowerLimit << ", "
235 << fWidthLowerLimit << " (in degrees!!)" << endl
236
237 << "The Hadronnes of those images exceeding one of these limits will "
238 << "be set to 0.75." << endl << endl;
239
240
241}
242
243
244// Function implementing a filter for muon induced images
245// (Values from Keichi cuts October 18th 2004)
246
247// Function returns kTRUE if the event does NOT survive the
248// muon filter. So kTRUE means event MUST be rejected
249
250
251Bool_t MSupercutsCalcONOFF::MuonFilter (Double_t Size,
252 Double_t Length)
253{
254 // Parametrization for cut in length.
255 // Event is rejected if
256 // Log10(Length) < k1*(log10(Size)) + k2 (Size<10000photons)
257 // k1 = 0.286; k2 = -1.91 (October 18th 2004)
258
259 // Length < 0.2 deg (Size > 10000)
260
261
262 // Parametrization for cut in Length/Size
263 // Event is rejected if
264 // log10(LenOverSize) < k3*(log10(Size)) + k4 (Size<4000photons)
265 // k3 = -0.780; k4 = -1.71
266
267
268 Double_t LengthUpperLimit = 0.0;
269 Double_t LengthOverSizeUpperLimit = 0.0;
270
271
272 Double_t SizeUpperLimitForCutInLength = 10000; // in photons
273 Double_t SizeUpperLimitForCutInLengthOverSize = 4000; // in photons
274
275 Double_t k1 = 0.286;
276 Double_t k2 = -1.91;
277 Double_t k3 = -0.78;
278 Double_t k4 = -1.71;
279
280
281
282
283
284 if (Size <= SizeUpperLimitForCutInLength)
285 {
286 LengthUpperLimit = k1 * TMath::Log10(Size) + k2;
287 LengthUpperLimit = TMath::Power(10, LengthUpperLimit);
288 }
289 else
290 {
291 LengthUpperLimit = 0.2;
292 }
293
294 // Apply cut in Length
295
296 if (Length < LengthUpperLimit)
297 return kTRUE;
298
299
300
301 if (Size < SizeUpperLimitForCutInLengthOverSize)
302 {
303 LengthOverSizeUpperLimit = k3 * TMath::Log10(Size) + k4;
304 LengthOverSizeUpperLimit = TMath::Power(10, LengthOverSizeUpperLimit);
305
306 if (Length/Size < LengthOverSizeUpperLimit)
307 return kTRUE;
308 }
309
310
311 return kFALSE;
312
313}
314
315
316
317
318
319
320
321
322Int_t MSupercutsCalcONOFF::PreProcess(MParList *pList)
323{
324 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
325 if (!cam)
326 {
327 *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
328 return kFALSE;
329 }
330
331 fMm2Deg = cam->GetConvMm2Deg();
332
333 fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
334 if (!fHadronness)
335 {
336 *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
337 return kFALSE;
338 }
339
340 fSuper = (MSupercuts*)pList->FindObject(fSuperName, "MSupercuts");
341 if (!fSuper)
342 {
343 *fLog << err << fSuperName << " [MSupercuts] not found... aborting." << endl;
344 return kFALSE;
345 }
346
347 if (fStoreAppliedSupercuts)
348 {
349
350 fSupercutsApplied = (MTSupercutsApplied*)pList->FindObject(fSupercutsAppliedName,
351 "MTSupercutsApplied");
352 if(!fSupercutsApplied)
353 {
354 *fLog << err << fSupercutsAppliedName
355 << " [MTSupercutsApplied] not found... aborting." << endl;
356 return kFALSE;
357 }
358
359
360 }
361
362 if (fMatrix)
363 return kTRUE;
364
365 //-----------------------------------------------------------
366 fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
367 if (!fHil)
368 {
369 *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
370 return kFALSE;
371 }
372
373 fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt");
374 if (!fHilExt)
375 {
376 *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl;
377 return kFALSE;
378 }
379
380 fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
381 if (!fHilSrc)
382 {
383 *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
384 return kFALSE;
385 }
386
387 fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar");
388 if (!fNewPar)
389 {
390 *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl;
391 return kFALSE;
392 }
393
394 fPointPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
395 if (!fPointPos)
396 {
397 *fLog << err << "MSupercutsCalcONOFF::PreProcess; MPointingPos not found... aborting." << endl;
398 return kFALSE;
399 }
400
401
402
403 return kTRUE;
404}
405
406// --------------------------------------------------------------------------
407//
408// Calculation of upper and lower limits
409//
410Double_t MSupercutsCalcONOFF::CtsMCut(const Double_t* a, Double_t ls, Double_t ct,
411 Double_t ls2, Double_t dd2) const
412{
413 // define cut-function
414 //
415 //
416 // dNOMCOSZA = 1.0
417 //
418 // a: array of cut parameters
419 // ls: log(SIZE) - log(fSizeOffset)
420 // ls2: ls^2
421 // ct: Cos(ZA.) - Cos(fThetaOffset)
422 // dd2: DIST^2 - fDistOffset^2
423
424 // lconc: log10CONC
425
426
427 Double_t limit;
428
429 if(fUseStaticCuts)
430 { // static cuts are used
431 limit = a[0];
432 }
433 else
434 {
435 if (fNotUseTheta)
436 { // Theta info is NOT used in the computation of the dynamical cuts
437 // For the time being dist info will not be used
438
439 if(fUseDist)
440 {
441
442
443 limit = a[0] + a[1] * dd2 +
444 ls * (a[3] + a[4] * dd2) +
445 ls2 * (a[6] + a[7] * dd2);
446
447
448 }
449 else
450 {
451 limit = a[0] + ls * (a[3])+ ls2 * (a[6]);
452 }
453
454
455 }
456 else
457 {// Theta info IS used in the computation of the dynamical cuts
458
459 limit =
460 a[0] + a[1] * dd2 + a[2] * ct +
461 ls * (a[3] + a[4] * dd2 + a[5] * ct) +
462 ls2 * (a[6] + a[7] * dd2);
463 }
464
465 }
466
467
468
469 //*fLog << "MSupercutsCalcONOFF::CtsMCut; *a = "
470 // << *a << ", " << *(a+1) << ", " << *(a+2) << ", "
471 // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", "
472 // << *(a+6) << ", " << *(a+7) << endl;
473
474 //*fLog << "MSupercutsCalcONOFF::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
475 // << ", " << ls2 << ", " << ct << ", " << dd2 << ", "
476 // << limit << endl;
477
478 return limit;
479}
480
481// --------------------------------------------------------------------------
482//
483// Returns the mapped value from the Matrix
484//
485Double_t MSupercutsCalcONOFF::GetVal(Int_t i) const
486{
487
488 Double_t val = (*fMatrix)[fMap[i]];
489
490
491 //*fLog << "MSupercutsCalcONOFF::GetVal; i, fMatrix, fMap, val = "
492 // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl;
493
494
495 return val;
496}
497
498// --------------------------------------------------------------------------
499//
500// You can use this function if you want to use a MHMatrix instead of the
501// given containers. This function adds all necessary columns to the
502// given matrix. Afterward you should fill the matrix with the corresponding
503// data (eg from a file by using MHMatrix::Fill). If you now loop
504// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
505// will take the values from the matrix instead of the containers.
506//
507
508
509void MSupercutsCalcONOFF::InitMapping(MHMatrix *mat)
510{
511 if (fMatrix)
512 return;
513
514 fMatrix = mat;
515
516 fMap[0] = fMatrix->AddColumn("MPointingPos.fZd"); //deg
517 fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
518 fMap[2] = fMatrix->AddColumn("MHillas.fLength");
519 fMap[3] = fMatrix->AddColumn("MHillas.fSize");
520 fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
521 fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
522 fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
523 fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)");
524 fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
525 fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc");
526 fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1");
527}
528
529
530// ---------------------------------------------------------------------------
531//
532// Evaluate dynamical supercuts
533//
534// set hadronness to 0.25 if cuts are fullfilled
535// 0.75 otherwise
536//
537Int_t MSupercutsCalcONOFF::Process()
538{
539 // const Double_t kNomLogSize = 4.1;
540 const Double_t kNomCosZA = 1.0;
541
542 //const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
543 // For the time being I do it in a way that theta must be
544 // the value used in the row 0 of fMatrix. That means that if there is
545 // no matrix, the program will complain (and crash). And tehn I will know
546 // that value 0 (supposed to be fThetaOrig.Val) could not be taken.
547 const Double_t theta = GetVal(0);
548 const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth();
549 const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
550 const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize();
551 const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX();
552 const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY();
553 const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist();
554
555 const Double_t alpha = fMatrix ? GetVal(7) : fHilSrc->GetAlpha();
556
557 Double_t help;
558 if (!fMatrix)
559 help = TMath::Sign(fHilExt->GetM3Long(),
560 fHilSrc->GetCosDeltaAlpha());
561 const Double_t asym0 = fMatrix ? GetVal(8) : help;
562 const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc();
563 const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1();
564
565 const Double_t newdist = dist0 * fMm2Deg;
566
567 const Double_t dist2 = meanx*meanx + meany*meany;
568 const Double_t dist = sqrt(dist2) * fMm2Deg;
569
570 // const Double_t dd2 = dist*dist;
571
572
573
574 // The parametrization of upper cut in CONC is just provisional.
575 // The aim is to get rid of triangular events.
576 // The parametrization is
577 // Log10(Conc)Uppercut = m*(log(Size) - log(SizeOffset)) + b
578
579
580 Double_t log10conc = TMath::Log10(conc);
581
582
583 // in the parameterization of the cuts
584 // the dist parameter used is the one computed from the source position,
585 // and not the one computed from the camera center.
586
587 // Actually the value used in the parameterization is newdist^2,
588 // minus the offset_in_dist^2
589
590 const Double_t dd2 = newdist*newdist - fDistOffset*fDistOffset;
591
592
593
594 const Double_t dmls = log(size) - log(fSizeOffset);
595 const Double_t dmls2 = dmls * dmls;
596
597 const Double_t dmcza = cos(theta) - kNomCosZA;
598
599 const Double_t length = length0 * fMm2Deg;
600 const Double_t width = width0 * fMm2Deg;
601 const Double_t asym = asym0 * fMm2Deg;
602
603
604 // MuuonLikeEvent is a variable which is set
605 // to kFALSE/kTRUE by the filter cuts implemented
606 // in MuonFilter function. If kTRUE, the event
607 // will be rejected.
608
609// DM:
610 //Bool_t MuonLikeEvent = MuonFilter (size, length);
611
612
613 // computation of the cut limits
614
615 Double_t lengthup = CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2);
616 Double_t lengthlow = CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2);
617
618 Double_t widthup = CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2);
619 Double_t widthlow = CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2);
620
621 Double_t distup = CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2);
622 Double_t distlow = CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2);
623
624 Double_t asymup = CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2);
625 Double_t asymlow = CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2);
626
627 Double_t log10concup = CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2);
628 Double_t log10conclow = CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2);
629
630 Double_t leakageup = CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2);
631 Double_t leakagelow = CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2);
632
633
634 Double_t lengthoverwidth = length/width;
635
636
637 // Cut in CONC is parametrized
638
639 Double_t hadronness = 0.0;
640
641
642 // If the cut values computed before for Dist,
643 // Length and Width exceed the upper limits set
644 // by the user through function
645 // MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits,
646 // (or the default values defined in constructor); such cut values
647 // are replaced by the
648 // the limit values
649
650
651
652 if (distup > fDistUpperLimit)
653 {
654 distup = fDistUpperLimit;
655 }
656
657
658 if (lengthup > fLengthUpperLimit)
659 {
660 lengthup = fLengthUpperLimit;
661 }
662
663 if (widthup > fWidthUpperLimit)
664 {
665 widthup = fWidthUpperLimit;
666 }
667
668
669
670 if (distlow < fDistLowerLimit)
671 {
672 distlow = fDistLowerLimit;
673 }
674
675 if (lengthlow < fLengthLowerLimit)
676 {
677 lengthlow = fLengthLowerLimit;
678 }
679
680 if (widthlow < fWidthLowerLimit)
681 {
682 widthlow = fWidthLowerLimit;
683 }
684
685 if (log10concup > fLog10ConcUpperLimit)
686 {
687 log10concup = fLog10ConcUpperLimit;
688 }
689
690
691 if (leakageup > fLeakage1UpperLimit)
692 {
693 leakageup = fLeakage1UpperLimit;
694 }
695
696
697
698
699 // Upper cut in leakage 1 also implemented
700
701 if (//
702 newdist > distup ||
703 newdist < distlow ||
704 length > lengthup ||
705 length < lengthlow ||
706 width > widthup ||
707 width < widthlow ||
708 leakage > leakageup ||
709 leakage < leakagelow ||
710 lengthoverwidth < fLengthOverWidthUpperLimit ||
711 //
712 //asym < asymup &&
713 //asym > asymlow &&
714
715 //dist < distup &&
716 //dist > distlow &&
717
718 log10conc > log10concup
719 // log10conc < log10conclow
720 //
721 //log10conc > log10concup ||
722 // DM MuonLikeEvent // if KTRUE, event is rejected
723 )
724
725 {hadronness = 0.75;}
726 else
727 {
728 hadronness = 0.25;
729 }
730
731
732 fHadronness->SetHadronness(hadronness);
733 fHadronness->SetReadyToSave();
734
735
736 if(fStoreAppliedSupercuts)
737 {
738
739 // TArrayD vector with the shower parameters (matching the ones
740 // specified in function MTSupercutsApplied::CreateTreeBranches)
741 // is created and filled with this event features
742
743 // Eventually this definition might be done from outside this
744 // class, allowing for adding removing parameters without
745 // recompiling mars. yet for the time being, let's keep it simple.
746
747 TArrayD ShowerParams(10);
748 ShowerParams[0] = length;
749 ShowerParams[1] = width;
750 ShowerParams[2] = dist;
751 ShowerParams[3] = newdist;
752 ShowerParams[4] = asym;
753 ShowerParams[5] = log10conc;
754 ShowerParams[6] = leakage;
755 ShowerParams[7] = theta;
756 ShowerParams[8] = size;
757 ShowerParams[9] = alpha;
758
759
760 // TArrayD vector with the cut parameters (matching the ones
761 // specified in function MTSupercutsApplied::CreateTreeBranches)
762 // is created and filled with this event features
763
764 // Eventually this definition might be done from outside this
765 // class, allowing for adding removing parameters without
766 // recompiling mars. yet for the time being, let's keep it simple.
767
768
769 TArrayD SuperCutParams(13);
770 SuperCutParams[0] = lengthup;
771 SuperCutParams[1] = lengthlow;
772 SuperCutParams[2] = widthup;
773 SuperCutParams[3] = widthlow;
774 SuperCutParams[4] = distup;
775 SuperCutParams[5] = distlow;
776 SuperCutParams[6] = asymup;
777 SuperCutParams[7] = asymlow;
778 SuperCutParams[8] = log10concup;
779 SuperCutParams[9] = log10conclow;
780 SuperCutParams[10] = leakageup;
781 SuperCutParams[11] = leakagelow;
782 SuperCutParams[12] = hadronness;
783
784 // SC parameters applied to this event, as well as the
785 // shower parameters, are stored in the branches of the
786 // TTree object of a MTSupercutsApplied object
787
788
789 if (!fSupercutsApplied ->FillTreeBranches(SuperCutParams, ShowerParams))
790 {
791 *fLog << "MSupercutsCalcONOFF::Process()" << endl
792 << "Supercuts applied could not be stored in tree..."
793 << endl;
794
795 return kFALSE;
796 }
797
798
799 }
800
801
802
803 return kTRUE;
804
805
806 // OLD STUFF
807
808 /*
809 *fLog << "newdist, length, width, asym, dist, conc, leakage = "
810 << newdist << ", " << length << ", " << width << ", "
811 << asym << ", " << dist << ", " << conc << ", " << leakage
812 << endl;
813
814 *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = "
815 << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", "
816 << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", "
817
818 << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", "
819 << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", "
820
821 << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", "
822 << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", "
823
824 << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", "
825 << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", "
826
827 << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", "
828 << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", "
829
830 << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", "
831 << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", "
832
833 << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", "
834 << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", "
835 << endl;
836 */
837
838 /*
839 if (
840 //dist < 1.05 &&
841 //newdist < 1.05 &&
842
843 newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) &&
844 newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) &&
845
846 length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) &&
847 length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) &&
848
849 width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) &&
850 width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) &&
851
852 asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) &&
853 asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) &&
854
855 dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) &&
856 dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) &&
857
858 conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) &&
859 conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) &&
860
861 leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) &&
862 leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) )
863
864 fHadronness->SetHadronness(0.25);
865 else
866 fHadronness->SetHadronness(0.75);
867
868
869
870 // *fLog << "SChadroness = " << fHadronness->GetHadronness() << endl;
871
872 fHadronness->SetReadyToSave();
873 */
874
875 // END OF OLD STUFF
876}
877
878/*
879Bool_t MSupercutsCalcONOFF::StoreSupercutsAppliedToThisEvent(TArrayD CutParams,
880 TArrayD ShowerParams)
881{
882 // SC parameters applied to this event are stored in
883 // TTree object of MTSupercutsApplied object
884
885
886 if (!fSupercutsApplied ->FillTreeBranches(CutParams, ShowerParams))
887 {
888 *fLog << "MSupercutsCalcONOFF::StoreSupercutsAppliedToThisEvent" << endl
889 << "Supercuts applied could not be stored in tree..."
890 << endl;
891
892 return kFALSE;
893 }
894
895
896 return kTRUE;
897}
898
899*/
900
901
902
903
904
905
906
907
Note: See TracBrowser for help on using the repository browser.