source: trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc@ 9227

Last change on this file since 9227 was 8890, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 21.9 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 = log10(MNewImagePar.fLeakage1+1)
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]*leak +
54// (lgsize>c[10])*c[11]*(lgsize-c[10])^2;
55// p = xi*(1-width/length);
56// sign1 = m3long-c[5]
57// sign2 = (dist-c[7])*c[6]-slope
58// disp = sign1<0 ||sign2<0 ? -disp : disp
59// thetasq = disp^2 + dist^2 - 2*disp*dist*cos(alpha)
60//
61// And the values with respect to the antisource position respectively.
62//
63//
64// c[1]:
65// thetasq < c[1]*c[1]
66//
67// AreaCut:
68// c[2], c[3], c[4]:
69// A = c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]))
70// area < A
71//
72// HadronnessCut:
73// c[13], c[14]:
74// had <= c[13]
75// 10^lgsize >= c[14]
76//
77//
78// Options:
79// --------
80//
81// HadronnessCut:
82// - none/nocut: No area cut
83// - area: Area cut <default>
84// - hadronness: Cut in size and hadronness (c[10], c[11])
85// - all: area + hadronness
86//
87// ThetaCut:
88// - none/nocut: no theta cut <default>
89// - on: cut in theta only
90// - off: anti-cut in anti-theta only
91// - wobble: same as on|off (both)
92//
93// CalcDisp:
94// - yes: Disp is calculated as defined above <default>
95// - no: abs(Disp.fVal) from the parameter list is used instead
96//
97// CalcGhostbuster:
98// - yes: The ghostbuster is calculated as defined above <default>
99// - no: Ghostbuster.fVal<c[12] is used as ghostbusting condition
100//
101//
102// Class Version 2:
103// ----------------
104// + Bool_t fCalcDisp; // Should we use Disp from the parameterlist?
105// + Bool_t fCalcGhostbuster; // Should we use Ghostbuster from the parameterlist?
106//
107//
108// Input Container:
109// ------
110// MGeomCam
111// MHillas
112// MHillasExt
113// MNewImagePar
114// MHillasSrc
115// [MHillasSrcAnti [MHillasSrc]]
116// [Disp [MParameterD]]
117// [Ghostbuster [MParameterD]]
118//
119// Output:
120// -------
121// ThetaSquared [MParameterD] // stores thetasq
122// Disp [MParameterD] // stores -p*sign(MHillasSrc.fCosDeltaAlpha)
123// ThetaSquaredCut [MParameterD] // stores c[1]*c[1]
124//
125/////////////////////////////////////////////////////////////////////////////
126#include "MFMagicCuts.h"
127
128#include <fstream>
129#include <errno.h>
130
131#include "MHMatrix.h"
132
133#include "MLog.h"
134#include "MLogManip.h"
135
136#include "MMath.h"
137
138#include "MParList.h"
139
140#include "MGeomCam.h"
141#include "MHillas.h"
142#include "MHillasSrc.h"
143#include "MHillasExt.h"
144#include "MNewImagePar.h"
145#include "MParameters.h"
146
147ClassImp(MFMagicCuts);
148
149using namespace std;
150
151// --------------------------------------------------------------------------
152//
153// constructor
154//
155MFMagicCuts::MFMagicCuts(const char *name, const char *title)
156 : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0),
157 fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(20),
158 fThetaCut(kNone), fHadronnessCut(kArea), fCalcDisp(kTRUE),
159 fCalcGhostbuster(kTRUE)
160{
161 fName = name ? name : "MFMagicCuts";
162 fTitle = title ? title : "Class to evaluate the MagicCuts";
163
164 fMatrix = NULL;
165
166 AddToBranchList("MHillas.fWidth");
167 AddToBranchList("MHillas.fLength");
168 AddToBranchList("MHillas.fSize");
169 AddToBranchList("MHillasSrc.fAlpha");
170 AddToBranchList("MHillasSrc.fDist");
171 AddToBranchList("MHillasSrc.fCosDeltaAlpha");
172 AddToBranchList("MHillasSrcAnti.fAlpha");
173 AddToBranchList("MHillasSrcAnti.fDist");
174 AddToBranchList("MHillasSrcAnti.fCosDeltaAlpha");
175 AddToBranchList("MHillasExt.fM3Long");
176 AddToBranchList("MHillasExt.fSlopeLong");
177 AddToBranchList("MNewImagePar.fLeakage1");
178 AddToBranchList("Hadronness.fVal");
179 AddToBranchList("Disp.fVal");
180 AddToBranchList("Ghostbuster.fVal");
181}
182
183// --------------------------------------------------------------------------
184//
185// The theta cut value GetThetaSqCut() is propageted to the parameter list
186// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else.
187//
188Int_t MFMagicCuts::PreProcess(MParList *pList)
189{
190 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
191 if (!cam)
192 {
193 *fLog << err << "MGeomCam not found... aborting." << endl;
194 return kFALSE;
195 }
196
197 // Geometry in SPONDE?
198
199 fMm2Deg = cam->GetConvMm2Deg();
200
201 fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
202 if (!fThetaSq)
203 return kFALSE;
204
205 if (!fCalcDisp)
206 fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD");
207 else
208 fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
209 if (!fDisp)
210 {
211 *fLog << err << "Disp [MParameterD] not found... aborting." << endl;
212 return kFALSE;
213 }
214
215 if (!fCalcGhostbuster)
216 fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD");
217 else
218 fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster");
219 if (!fGhostbuster)
220 {
221 *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl;
222 return kFALSE;
223 }
224
225 // propagate Theta cut to the parameter list
226 // so that it can be used somewhere else
227 MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
228 if (!thetacut)
229 return kFALSE;
230 thetacut->SetVal(GetThetaSqCut());
231
232 Print();
233
234 if (fMatrix)
235 return kTRUE;
236
237 //-----------------------------------------------------------
238
239 fHil = (MHillas*)pList->FindObject("MHillas");
240 if (!fHil)
241 {
242 *fLog << err << "MHillas not found... aborting." << endl;
243 return kFALSE;
244 }
245
246 fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
247 if (!fHilExt)
248 {
249 *fLog << err << "MHillasExt not found... aborting." << endl;
250 return kFALSE;
251 }
252
253 fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
254 if (!fNewImgPar)
255 {
256 *fLog << err << "MNewImagePar not found... aborting." << endl;
257 return kFALSE;
258 }
259
260 fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
261 if (!fHilSrc)
262 {
263 *fLog << err << "MHillasSrc not found... aborting." << endl;
264 return kFALSE;
265 }
266
267 if (fThetaCut&kOff)
268 {
269 fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
270 if (!fHilAnti)
271 {
272 *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
273 return kFALSE;
274 }
275 }
276
277 if (fHadronnessCut&kHadronness)
278 {
279 fHadronness = (MParameterD*)pList->FindObject("Hadronness", "MParameterD");
280 if (!fHadronness)
281 {
282 *fLog << warn << "Hadronness [MParameterD] not found... aborting." << endl;
283 return kFALSE;
284 }
285 }
286
287 return kTRUE;
288}
289
290// --------------------------------------------------------------------------
291//
292// Returns the mapped value from the Matrix
293//
294Double_t MFMagicCuts::GetVal(Int_t i) const
295{
296 return (*fMatrix)[fMap[i]];
297}
298
299// --------------------------------------------------------------------------
300//
301// You can use this function if you want to use a MHMatrix instead of the
302// given containers. This function adds all necessary columns to the
303// given matrix. Afterward you should fill the matrix with the corresponding
304// data (eg from a file by using MHMatrix::Fill). If you now loop
305// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
306// will take the values from the matrix instead of the containers.
307//
308void MFMagicCuts::InitMapping(MHMatrix *mat)
309{
310 if (fMatrix)
311 return;
312
313 fMatrix = mat;
314
315// fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
316// fMap[kEWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
317
318 fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
319 fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)");
320 fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
321
322 fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");
323
324 fMap[kEAlpha] = fMatrix->AddColumn("MHillasSrc.fAlpha");
325 fMap[kEDist] = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
326 fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
327
328 fMap[kESign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
329
330 fMap[kESlope] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
331
332 if (!fCalcDisp)
333 fMap[kEDisp] = fMatrix->AddColumn("abs(Disp.fVal)");
334
335 if (!fCalcGhostbuster)
336 fMap[kEGhostbuster] = fMatrix->AddColumn("Ghostbuster.fVal");
337
338 if (fThetaCut&kOff)
339 {
340 fMap[kEAlphaAnti] = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
341 fMap[kEDistAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
342 fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
343 fMap[kESlopeAnti] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
344 }
345
346 if (fHadronnessCut&kHadronness)
347 fMap[kEHadronness] = fMatrix->AddColumn("Hadronness.fVal");
348}
349
350// --------------------------------------------------------------------------
351//
352// Returns theta squared based on the formula:
353// p := c*(1-width/length)
354// return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
355//
356// If par!=NULL p is stored in this MParameterD
357//
358Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t d, Double_t a) const
359{
360 // wl = width/l [1]
361 // d = dist [deg]
362 // a = alpha [deg]
363 return d*d + p*p - 2*d*p*TMath::Cos(a*TMath::DegToRad());
364}
365
366// --------------------------------------------------------------------------
367//
368// Returns squared cut value in theta: fVariables[1]^2
369//
370Double_t MFMagicCuts::GetThetaSqCut() const
371{
372 return fVariables[1]*fVariables[1];
373}
374
375// --------------------------------------------------------------------------
376//
377// Return abs(Disp.fVal) if disp calculation is switched off.
378// Otherwise return (c0+c6*leak^c7)*(1-width/length)
379//
380Double_t MFMagicCuts::GetDisp(Double_t slope, Double_t lgsize) const
381{
382 if (!fCalcDisp)
383 return fMatrix ? GetVal(kEDisp) : TMath::Abs(fDisp->GetVal());
384
385 // Get some variables
386 const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
387 const Double_t leak = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
388
389 // For simplicity
390 const Double_t *c = fVariables.GetArray();
391
392 // As rule for root or MDataPhrase:
393 // ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0]
394 //
395 Double_t xi = c[0] + c[8]*slope + c[9]*leak;
396 if (lgsize>c[10])
397 xi += (lgsize-c[10])*(lgsize-c[10])*c[11];
398
399 const Double_t disp = xi*(1-wdivl);
400
401 return disp;
402}
403
404Bool_t MFMagicCuts::IsGhost(Double_t m3long, Double_t slope, Double_t dist) const
405{
406 // For simplicity
407 const Double_t *c = fVariables.GetArray();
408
409 if (!fCalcGhostbuster)
410 return (fMatrix ? GetVal(kEGhostbuster) : fGhostbuster->GetVal())<c[12];
411
412 // Do Ghostbusting with slope and m3l
413 const Bool_t sign1 = m3long < c[5];
414 const Bool_t sign2 = slope > (dist-c[7])*c[6];
415
416 return sign1 || sign2;
417}
418
419// ---------------------------------------------------------------------------
420//
421// Evaluate dynamical magic-cuts
422//
423Int_t MFMagicCuts::Process()
424{
425 // For simplicity
426 const Double_t *c = fVariables.GetArray();
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.