source: branches/MarsMoreSimulationTruth/mfilter/MFMagicCuts.cc

Last change on this file was 18157, checked in by Daniela Dorner, 10 years ago
fixed typo
File size: 22.1 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): Thomas Bretz, 03/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2007
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MFMagicCuts
28//
29// Predefinitions:
30// ---------------
31// width/length = MHillas.fWidth/MHillas.fLength
32// area = MHillas.GetArea*fMm2Deg*fMm2Deg
33// lgsize = log10(MHillas.fSize)
34//
35// leakage1 = MNewImagePar.fLeakage1
36//
37// alpha = MHillasSrc.fAlpha
38// dist = MHillasSrc.fDist*fMm2Deg
39// m3long = MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*fMm2Deg
40// slope = MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/fMm2Deg
41//
42// antialpha = MHillasSrcAnti.fAlpha
43// antidist = MHillasSrcAnti.fDist*fMm2Deg
44// antim3long = MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMm2Deg
45// antislope = MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/fMm2Deg
46//
47// had = Hadronness.fVal
48//
49// Coefficients/Cuts:
50// ------------------
51//
52// c[0], c[5], c[6], c[7], c[8], c[9]:
53// xi = c[0] + c[8]*slope + c[9]*(1-1/(1+c[10]*leak))
54// p = xi*(1-width/length);
55// sign1 = m3long-c[5]
56// sign2 = (dist-c[7])*c[6]-slope
57// disp = sign1<0 ||sign2<0 ? -disp : disp
58// thetasq = disp^2 + dist^2 - 2*disp*dist*cos(alpha)
59//
60// And the values with respect to the antisource position respectively.
61//
62//
63// c[1]:
64// thetasq < c[1]*c[1]
65//
66// AreaCut:
67// c[2], c[3], c[4]:
68// A = c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]))
69// area < A
70//
71// HadronnessCut:
72// c[13], c[14]:
73// had <= c[13]
74// 10^lgsize >= c[14]
75//
76//
77// Options:
78// --------
79//
80// HadronnessCut:
81// - none/nocut: No area cut
82// - area: Area cut <default>
83// - hadronness: Cut in size and hadronness (c[10], c[11])
84// - all: area + hadronness
85//
86// ThetaCut:
87// - none/nocut: no theta cut <default>
88// - on: cut in theta only
89// - off: anti-cut in anti-theta only
90// - wobble: same as on|off (both)
91//
92// CalcDisp:
93// - yes: Disp is calculated as defined above <default>
94// - no: abs(Disp.fVal) from the parameter list is used instead
95//
96// CalcGhostbuster:
97// - yes: The ghostbuster is calculated as defined above <default>
98// - no: Ghostbuster.fVal<c[12] is used as ghostbusting condition
99//
100//
101// Class Version 2:
102// ----------------
103// + Bool_t fCalcDisp; // Should we use Disp from the parameterlist?
104// + Bool_t fCalcGhostbuster; // Should we use Ghostbuster from the parameterlist?
105//
106//
107// Input Container:
108// ------
109// MGeomCam
110// MHillas
111// MHillasExt
112// MNewImagePar
113// MHillasSrc
114// [MHillasSrcAnti [MHillasSrc]]
115// [Disp [MParameterD]]
116// [Ghostbuster [MParameterD]]
117//
118// Output:
119// -------
120// ThetaSquared [MParameterD] // stores thetasq
121// Disp [MParameterD] // stores -p*sign(MHillasSrc.fCosDeltaAlpha)
122// ThetaSquaredCut [MParameterD] // stores c[1]*c[1]
123//
124/////////////////////////////////////////////////////////////////////////////
125#include "MFMagicCuts.h"
126
127#include <fstream>
128#include <errno.h>
129
130#include "MHMatrix.h"
131
132#include "MLog.h"
133#include "MLogManip.h"
134
135#include "MMath.h"
136
137#include "MParList.h"
138
139#include "MGeomCam.h"
140#include "MHillas.h"
141#include "MHillasSrc.h"
142#include "MHillasExt.h"
143#include "MNewImagePar.h"
144#include "MParameters.h"
145
146ClassImp(MFMagicCuts);
147
148using namespace std;
149
150// --------------------------------------------------------------------------
151//
152// constructor
153//
154MFMagicCuts::MFMagicCuts(const char *name, const char *title)
155 : fGeom(0), fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0),
156 fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(20),
157 fThetaCut(kNone), fHadronnessCut(kArea), fCalcDisp(kTRUE),
158 fCalcGhostbuster(kTRUE)
159{
160 fName = name ? name : "MFMagicCuts";
161 fTitle = title ? title : "Class to evaluate the MagicCuts";
162
163 fMatrix = NULL;
164
165 AddToBranchList("MHillas.fWidth");
166 AddToBranchList("MHillas.fLength");
167 AddToBranchList("MHillas.fSize");
168 AddToBranchList("MHillasSrc.fAlpha");
169 AddToBranchList("MHillasSrc.fDist");
170 AddToBranchList("MHillasSrc.fCosDeltaAlpha");
171 AddToBranchList("MHillasSrcAnti.fAlpha");
172 AddToBranchList("MHillasSrcAnti.fDist");
173 AddToBranchList("MHillasSrcAnti.fCosDeltaAlpha");
174 AddToBranchList("MHillasExt.fM3Long");
175 AddToBranchList("MHillasExt.fSlopeLong");
176 AddToBranchList("MNewImagePar.fLeakage1");
177 AddToBranchList("Hadronness.fVal");
178 AddToBranchList("Disp.fVal");
179 AddToBranchList("Ghostbuster.fVal");
180}
181
182// --------------------------------------------------------------------------
183//
184// The theta cut value GetThetaSqCut() is propageted to the parameter list
185// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else.
186//
187Int_t MFMagicCuts::PreProcess(MParList *pList)
188{
189 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
190 if (!fGeom)
191 {
192 *fLog << err << "MGeomCam not found... aborting." << endl;
193 return kFALSE;
194 }
195
196 fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
197 if (!fThetaSq)
198 return kFALSE;
199
200 if (!fCalcDisp)
201 fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD");
202 else
203 fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
204 if (!fDisp)
205 {
206 *fLog << err << "Disp [MParameterD] not found... aborting." << endl;
207 return kFALSE;
208 }
209
210 if (!fCalcGhostbuster)
211 fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD");
212 else
213 fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster");
214 if (!fGhostbuster)
215 {
216 *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl;
217 return kFALSE;
218 }
219
220 // propagate Theta cut to the parameter list
221 // so that it can be used somewhere else
222 MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
223 if (!thetacut)
224 return kFALSE;
225 thetacut->SetVal(GetThetaSqCut());
226
227 Print();
228
229 if (fMatrix)
230 return kTRUE;
231
232 //-----------------------------------------------------------
233
234 fHil = (MHillas*)pList->FindObject("MHillas");
235 if (!fHil)
236 {
237 *fLog << err << "MHillas not found... aborting." << endl;
238 return kFALSE;
239 }
240
241 fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
242 if (!fHilExt)
243 {
244 *fLog << err << "MHillasExt not found... aborting." << endl;
245 return kFALSE;
246 }
247
248 fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
249 if (!fNewImgPar)
250 {
251 *fLog << err << "MNewImagePar not found... aborting." << endl;
252 return kFALSE;
253 }
254
255 fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
256 if (!fHilSrc)
257 {
258 *fLog << err << "MHillasSrc not found... aborting." << endl;
259 return kFALSE;
260 }
261
262 if (fThetaCut&kOff)
263 {
264 fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
265 if (!fHilAnti)
266 {
267 *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
268 return kFALSE;
269 }
270 }
271
272 if (fHadronnessCut&kHadronness)
273 {
274 fHadronness = (MParameterD*)pList->FindObject("Hadronness", "MParameterD");
275 if (!fHadronness)
276 {
277 *fLog << warn << "Hadronness [MParameterD] not found... aborting." << endl;
278 return kFALSE;
279 }
280 }
281
282 return kTRUE;
283}
284
285// --------------------------------------------------------------------------
286//
287// Returns the mapped value from the Matrix
288//
289Double_t MFMagicCuts::GetVal(Int_t i) const
290{
291 return (*fMatrix)[fMap[i]];
292}
293
294// --------------------------------------------------------------------------
295//
296// You can use this function if you want to use a MHMatrix instead of the
297// given containers. This function adds all necessary columns to the
298// given matrix. Afterward you should fill the matrix with the corresponding
299// data (eg from a file by using MHMatrix::Fill). If you now loop
300// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
301// will take the values from the matrix instead of the containers.
302//
303void MFMagicCuts::InitMapping(MHMatrix *mat)
304{
305 if (fMatrix)
306 return;
307
308 fMatrix = mat;
309
310// fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
311// fMap[kEWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
312
313 fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
314 fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)");
315 fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
316
317 //fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");
318 fMap[kELeakage] = fMatrix->AddColumn("MNewImagePar.fLeakage1");
319
320 fMap[kEAlpha] = fMatrix->AddColumn("MHillasSrc.fAlpha");
321 fMap[kEDist] = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
322 fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
323
324 fMap[kESign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
325
326 fMap[kESlope] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
327
328 if (!fCalcDisp)
329 fMap[kEDisp] = fMatrix->AddColumn("abs(Disp.fVal)");
330
331 if (!fCalcGhostbuster)
332 fMap[kEGhostbuster] = fMatrix->AddColumn("Ghostbuster.fVal");
333
334 if (fThetaCut&kOff)
335 {
336 fMap[kEAlphaAnti] = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
337 fMap[kEDistAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
338 fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
339 fMap[kESlopeAnti] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
340 }
341
342 if (fHadronnessCut&kHadronness)
343 fMap[kEHadronness] = fMatrix->AddColumn("Hadronness.fVal");
344}
345
346// --------------------------------------------------------------------------
347//
348// Returns theta squared based on the formula:
349// p := c*(1-width/length)
350// return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
351//
352// If par!=NULL p is stored in this MParameterD
353//
354Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t d, Double_t a) const
355{
356 // wl = width/l [1]
357 // d = dist [deg]
358 // a = alpha [deg]
359 return d*d + p*p - 2*d*p*TMath::Cos(a*TMath::DegToRad());
360}
361
362// --------------------------------------------------------------------------
363//
364// Returns squared cut value in theta: fVariables[1]^2
365//
366Double_t MFMagicCuts::GetThetaSqCut() const
367{
368 return fVariables[1]*fVariables[1];
369}
370
371// --------------------------------------------------------------------------
372//
373// Return abs(Disp.fVal) if disp calculation is switched off.
374// Otherwise return (c0+c6*leak^c7)*(1-width/length)
375//
376Double_t MFMagicCuts::GetDisp(Double_t slope, Double_t lgsize) const
377{
378 if (!fCalcDisp)
379 return fMatrix ? GetVal(kEDisp) : TMath::Abs(fDisp->GetVal());
380
381 // Get some variables
382 const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
383 //const Double_t leak = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
384 const Double_t leak = fMatrix ? GetVal(kELeakage) : fNewImgPar->GetLeakage1();
385
386 // For simplicity
387 const Double_t *c = fVariables.GetArray();
388
389 // As rule for root or MDataPhrase:
390 // ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0]
391 //char *r = "M[0]*([0] + [1]*M[1] + [2]*(1-1/(1+[3]*M[2])))";
392 //
393 Double_t xi = c[0] + c[8]*slope + c[9]*(1-1/(1+c[10]*leak));
394// if (lgsize>c[10])
395// xi += (lgsize-c[10])*(lgsize-c[10])*c[11];
396
397 const Double_t disp = xi*(1-wdivl);
398
399 return disp;
400}
401
402Bool_t MFMagicCuts::IsGhost(Double_t m3long, Double_t slope, Double_t dist) const
403{
404 // For simplicity
405 const Double_t *c = fVariables.GetArray();
406
407 if (!fCalcGhostbuster)
408 return (fMatrix ? GetVal(kEGhostbuster) : fGhostbuster->GetVal())<c[12];
409
410 // Do Ghostbusting with slope and m3l
411 const Bool_t sign1 = m3long < c[5];
412 const Bool_t sign2 = slope > (dist-c[7])*c[6];
413
414 return sign1 || sign2;
415}
416
417// ---------------------------------------------------------------------------
418//
419// Evaluate dynamical magic-cuts
420//
421Int_t MFMagicCuts::Process()
422{
423 // For simplicity
424 const Double_t *c = fVariables.GetArray();
425
426 const Float_t fMm2Deg = fGeom->GetConvMm2Deg();
427
428 // Default if we return before the end
429 fResult = kFALSE;
430
431 // Value for Theta cut (Disp parametrization)
432 const Double_t cut = GetThetaSqCut();
433
434 // ------------------- Most efficient cut -----------------------
435 // ---------------------- Theta cut ON --------------------------
436 const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg;
437 const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha();
438 const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
439 const Double_t slope = fMatrix ? GetVal(kESlope) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
440 const Double_t sign = fMatrix ? GetVal(kESign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
441 const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize());
442
443 // Calculate disp including sign
444 const Double_t disp = GetDisp(slope, lgsize);
445 const Double_t ghost = IsGhost(m3long, slope, dist);
446
447 const Double_t p = ghost ? -disp : disp;
448
449 // Align sign of disp along shower axis like m3long
450 fDisp->SetVal(-p*sign);
451
452 // FIXME: Allow to use ThetaSq from INPUT!
453
454 // Calculate corresponding Theta^2
455 const Double_t thetasq = GetThetaSq(p, dist, alpha);
456 fThetaSq->SetVal(thetasq);
457
458 // Perform theta cut
459 if (fThetaCut&kOn && thetasq >= cut)
460 return kTRUE;
461
462 // --------------------- Efficient cut -------------------------
463 // ---------------------- Area/M3l cut --------------------------
464 if (fHadronnessCut&kArea)
465 {
466 const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
467 const Double_t A = (c[2]*c[4]>0 ? (lgsize>c[3]) : (lgsize<c[3])) ?
468 c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
469 //const Double_t A = lgsize>c[3] ? c[2] : c[2] - c[4]/c[2]*(lgsize-c[3])*(lgsize-c[3]));
470 if (area>=A)
471 return kTRUE;
472 }
473
474 if (fHadronnessCut&kAreaLin)
475 {
476 const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
477 const Double_t A = c[2] + c[4]*(lgsize-c[3]);
478 if (area>=A)
479 return kTRUE;
480 }
481
482 // ---------------------------------------------------------------
483 // ---------------------------------------------------------------
484
485 if (fHadronnessCut&kHadronness)
486 {
487 const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
488 if (had>c[13])
489 return kTRUE;
490
491 if (TMath::Power(10, lgsize)<c[14])
492 return kTRUE;
493 }
494
495 // ------------------- Least efficient cut -----------------------
496 // ---------------------- Theta cut OFF --------------------------
497
498 if (fThetaCut&kOff)
499 {
500 const Double_t distanti = fMatrix ? GetVal(kEDistAnti) : fHilAnti->GetDist()*fMm2Deg;
501 const Double_t alphaanti = fMatrix ? GetVal(kEAlphaAnti) : fHilAnti->GetAlpha();
502 const Double_t m3lanti = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
503 const Double_t slopeanti = fMatrix ? GetVal(kESlopeAnti) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
504
505 const Double_t dispanti = GetDisp(slopeanti, lgsize);
506 const Double_t ghostanti = IsGhost(m3lanti, slopeanti, lgsize);
507
508 const Double_t panti = ghostanti ? -dispanti : dispanti;
509
510 // Align disp along source direction depending on third moment
511 const Double_t thetasqanti = GetThetaSq(panti, distanti, alphaanti);
512
513 if (thetasqanti<cut)
514 return kTRUE;
515
516 // This cut throws away gamma-like events which would be assigned to the
517 // anti-source. It is not necessary but it offers the possibility to
518 // fill the MHDisp plot in one run (Having a hole instead of eg. Crab
519 // the data can be rotated and subtracted)
520 }
521
522 fResult = kTRUE;
523
524 return kTRUE;
525}
526
527void MFMagicCuts::SetVariables(const TArrayD &arr)
528{
529 fVariables = arr;
530}
531
532TString MFMagicCuts::GetParam(Int_t i) const
533{
534 if (i>=fVariables.GetSize())
535 return "";
536
537 return Form("Param[%2d]=%+f", i, fVariables[i]);
538}
539
540void MFMagicCuts::Print(Option_t *o) const
541{
542 *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;
543
544 *fLog << "Using Theta cut: ";
545 switch (fThetaCut)
546 {
547 case kNone:
548 *fLog << "none" << endl;
549 break;
550 case kOn:
551 *fLog << "on" << endl;
552 break;
553 case kOff:
554 *fLog << "off" << endl;
555 break;
556 case kWobble:
557 *fLog << "on and off (wobble)" << endl;
558 break;
559 }
560 *fLog << "Using hadronness cut: ";
561 switch (fHadronnessCut)
562 {
563 case kNoCut:
564 *fLog << "none" << endl;
565 break;
566 case kArea:
567 *fLog << "area" << endl;
568 break;
569 case kHadronness:
570 *fLog << "hadronness" << endl;
571 break;
572 case kAreaLin:
573 *fLog << "arealin" << endl;
574 break;
575 case kAll:
576 *fLog << "all" << endl;
577 break;
578 }
579 if (fCalcDisp)
580 *fLog << "Disp is calculated from c0, c8-c11." << endl;
581 else
582 *fLog << "Disp.fVal from the parameter list is used as disp." << endl;
583 if (fCalcGhostbuster)
584 *fLog << "Ghostbusting is calculated from c5-c7." << endl;
585 else
586 *fLog << "Ghostbuster.fVal from the parameter list is used for ghostbusting." << endl;
587
588 *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
589
590 const Int_t n = fVariables.GetSize();
591 for (int i=0; i<n; i+=3)
592 {
593 *fLog << setw(25) << GetParam(i);
594 *fLog << setw(25) << GetParam(i+1);
595 *fLog << setw(25) << GetParam(i+2);
596 *fLog << endl;
597 }
598 *fLog << setiosflags(ios::right);
599}
600
601Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
602{
603 ifstream fin(fname);
604 if (!fin)
605 {
606 *fLog << err << "Cannot open file " << fname << ": ";
607 *fLog << strerror(errno) << endl;
608 return kFALSE;
609 }
610
611 for (int i=0; i<fVariables.GetSize(); i++)
612 {
613 fin >> fVariables[i];
614 if (!fin)
615 {
616 *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
617 return kERROR;
618 }
619 }
620 return kTRUE;
621}
622
623Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
624{
625 ofstream fout(fname);
626 if (!fout)
627 {
628 *fLog << err << "Cannot open file " << fname << ": ";
629 *fLog << strerror(errno) << endl;
630 return kFALSE;
631 }
632
633 for (int i=0; i<fVariables.GetSize(); i++)
634 fout << setw(11) << fVariables[i] << endl;
635
636 return kTRUE;
637}
638
639Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
640{
641 if (MFilter::ReadEnv(env, prefix, print)==kERROR)
642 return kERROR;
643
644
645 Bool_t rc = kFALSE;
646 if (IsEnvDefined(env, prefix, "ThetaCut", print))
647 {
648 TString str = GetEnvValue(env, prefix, "ThetaCut", "");
649 str.ToLower();
650 str=str.Strip(TString::kBoth);
651
652 if (str==(TString)"nocut" || str==(TString)"none")
653 fThetaCut = kNone;
654 if (str==(TString)"on")
655 fThetaCut = kOn;
656 if (str==(TString)"off")
657 fThetaCut = kOff;
658 if (str==(TString)"wobble")
659 fThetaCut = kWobble;
660 rc = kTRUE;
661 }
662
663 if (IsEnvDefined(env, prefix, "HadronnessCut", print))
664 {
665 TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
666 str.ToLower();
667 str=str.Strip(TString::kBoth);
668
669 if (str==(TString)"nocut" || str==(TString)"none")
670 fHadronnessCut = kNoCut;
671 if (str==(TString)"area")
672 fHadronnessCut = kArea;
673 if (str==(TString)"hadronness")
674 fHadronnessCut = kHadronness;
675 if (str==(TString)"all")
676 fHadronnessCut = kAll;
677
678 rc = kTRUE;
679 }
680
681 if (IsEnvDefined(env, prefix, "CalcDisp", print))
682 {
683 fCalcDisp = GetEnvValue(env, prefix, "CalcDisp", fCalcDisp);
684 rc = kTRUE;
685 }
686
687 if (IsEnvDefined(env, prefix, "CalcGhostbuster", print))
688 {
689 fCalcGhostbuster = GetEnvValue(env, prefix, "CalcGhostbuster", fCalcGhostbuster);
690 rc = kTRUE;
691 }
692
693 if (IsEnvDefined(env, prefix, "File", print))
694 {
695 const TString fname = GetEnvValue(env, prefix, "File", "");
696 if (!CoefficentsRead(fname))
697 return kERROR;
698
699 return kTRUE;
700 }
701
702 for (int i=0; i<fVariables.GetSize(); i++)
703 {
704 if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
705 {
706 // It is important that the last argument is a floating point
707 fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
708 rc = kTRUE;
709 }
710 }
711 return rc;
712}
Note: See TracBrowser for help on using the repository browser.