source: branches/MarsISDCBranchBasedOn17887/mfilter/MFMagicCuts.cc@ 18815

Last change on this file since 18815 was 9339, 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 : fGeom(0), 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 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
191 if (!fGeom)
192 {
193 *fLog << err << "MGeomCam not found... aborting." << endl;
194 return kFALSE;
195 }
196
197 fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
198 if (!fThetaSq)
199 return kFALSE;
200
201 if (!fCalcDisp)
202 fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD");
203 else
204 fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
205 if (!fDisp)
206 {
207 *fLog << err << "Disp [MParameterD] not found... aborting." << endl;
208 return kFALSE;
209 }
210
211 if (!fCalcGhostbuster)
212 fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD");
213 else
214 fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster");
215 if (!fGhostbuster)
216 {
217 *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl;
218 return kFALSE;
219 }
220
221 // propagate Theta cut to the parameter list
222 // so that it can be used somewhere else
223 MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
224 if (!thetacut)
225 return kFALSE;
226 thetacut->SetVal(GetThetaSqCut());
227
228 Print();
229
230 if (fMatrix)
231 return kTRUE;
232
233 //-----------------------------------------------------------
234
235 fHil = (MHillas*)pList->FindObject("MHillas");
236 if (!fHil)
237 {
238 *fLog << err << "MHillas not found... aborting." << endl;
239 return kFALSE;
240 }
241
242 fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
243 if (!fHilExt)
244 {
245 *fLog << err << "MHillasExt not found... aborting." << endl;
246 return kFALSE;
247 }
248
249 fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
250 if (!fNewImgPar)
251 {
252 *fLog << err << "MNewImagePar not found... aborting." << endl;
253 return kFALSE;
254 }
255
256 fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
257 if (!fHilSrc)
258 {
259 *fLog << err << "MHillasSrc not found... aborting." << endl;
260 return kFALSE;
261 }
262
263 if (fThetaCut&kOff)
264 {
265 fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
266 if (!fHilAnti)
267 {
268 *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
269 return kFALSE;
270 }
271 }
272
273 if (fHadronnessCut&kHadronness)
274 {
275 fHadronness = (MParameterD*)pList->FindObject("Hadronness", "MParameterD");
276 if (!fHadronness)
277 {
278 *fLog << warn << "Hadronness [MParameterD] not found... aborting." << endl;
279 return kFALSE;
280 }
281 }
282
283 return kTRUE;
284}
285
286// --------------------------------------------------------------------------
287//
288// Returns the mapped value from the Matrix
289//
290Double_t MFMagicCuts::GetVal(Int_t i) const
291{
292 return (*fMatrix)[fMap[i]];
293}
294
295// --------------------------------------------------------------------------
296//
297// You can use this function if you want to use a MHMatrix instead of the
298// given containers. This function adds all necessary columns to the
299// given matrix. Afterward you should fill the matrix with the corresponding
300// data (eg from a file by using MHMatrix::Fill). If you now loop
301// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
302// will take the values from the matrix instead of the containers.
303//
304void MFMagicCuts::InitMapping(MHMatrix *mat)
305{
306 if (fMatrix)
307 return;
308
309 fMatrix = mat;
310
311// fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
312// fMap[kEWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
313
314 fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
315 fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)");
316 fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
317
318 fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");
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
385 // For simplicity
386 const Double_t *c = fVariables.GetArray();
387
388 // As rule for root or MDataPhrase:
389 // ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0]
390 //
391 Double_t xi = c[0] + c[8]*slope + c[9]*leak;
392 if (lgsize>c[10])
393 xi += (lgsize-c[10])*(lgsize-c[10])*c[11];
394
395 const Double_t disp = xi*(1-wdivl);
396
397 return disp;
398}
399
400Bool_t MFMagicCuts::IsGhost(Double_t m3long, Double_t slope, Double_t dist) const
401{
402 // For simplicity
403 const Double_t *c = fVariables.GetArray();
404
405 if (!fCalcGhostbuster)
406 return (fMatrix ? GetVal(kEGhostbuster) : fGhostbuster->GetVal())<c[12];
407
408 // Do Ghostbusting with slope and m3l
409 const Bool_t sign1 = m3long < c[5];
410 const Bool_t sign2 = slope > (dist-c[7])*c[6];
411
412 return sign1 || sign2;
413}
414
415// ---------------------------------------------------------------------------
416//
417// Evaluate dynamical magic-cuts
418//
419Int_t MFMagicCuts::Process()
420{
421 // For simplicity
422 const Double_t *c = fVariables.GetArray();
423
424 const Float_t fMm2Deg = fGeom->GetConvMm2Deg();
425
426 // Default if we return before the end
427 fResult = kFALSE;
428
429 // Value for Theta cut (Disp parametrization)
430 const Double_t cut = GetThetaSqCut();
431
432 // ------------------- Most efficient cut -----------------------
433 // ---------------------- Theta cut ON --------------------------
434 const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg;
435 const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha();
436 const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
437 const Double_t slope = fMatrix ? GetVal(kESlope) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
438 const Double_t sign = fMatrix ? GetVal(kESign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
439 const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize());
440
441 // Calculate disp including sign
442 const Double_t disp = GetDisp(slope, lgsize);
443 const Double_t ghost = IsGhost(m3long, slope, dist);
444
445 const Double_t p = ghost ? -disp : disp;
446
447 // Align sign of disp along shower axis like m3long
448 fDisp->SetVal(-p*sign);
449
450 // FIXME: Allow to use ThetaSq from INPUT!
451
452 // Calculate corresponding Theta^2
453 const Double_t thetasq = GetThetaSq(p, dist, alpha);
454 fThetaSq->SetVal(thetasq);
455
456 // Perform theta cut
457 if (fThetaCut&kOn && thetasq >= cut)
458 return kTRUE;
459
460 // --------------------- Efficient cut -------------------------
461 // ---------------------- Area/M3l cut --------------------------
462 if (fHadronnessCut&kArea)
463 {
464 const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
465 const Double_t A = (c[2]*c[4]>0 ? (lgsize>c[3]) : (lgsize<c[3])) ?
466 c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
467 //const Double_t A = lgsize>c[3] ? c[2] : c[2] - c[4]/c[2]*(lgsize-c[3])*(lgsize-c[3]));
468 if (area>=A)
469 return kTRUE;
470 }
471
472 if (fHadronnessCut&kAreaLin)
473 {
474 const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
475 const Double_t A = c[2] + c[4]*(lgsize-c[3]);
476 if (area>=A)
477 return kTRUE;
478 }
479
480 // ---------------------------------------------------------------
481 // ---------------------------------------------------------------
482
483 if (fHadronnessCut&kHadronness)
484 {
485 const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
486 if (had>c[13])
487 return kTRUE;
488
489 if (TMath::Power(10, lgsize)<c[14])
490 return kTRUE;
491 }
492
493 // ------------------- Least efficient cut -----------------------
494 // ---------------------- Theta cut OFF --------------------------
495
496 if (fThetaCut&kOff)
497 {
498 const Double_t distanti = fMatrix ? GetVal(kEDistAnti) : fHilAnti->GetDist()*fMm2Deg;
499 const Double_t alphaanti = fMatrix ? GetVal(kEAlphaAnti) : fHilAnti->GetAlpha();
500 const Double_t m3lanti = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
501 const Double_t slopeanti = fMatrix ? GetVal(kESlopeAnti) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
502
503 const Double_t dispanti = GetDisp(slopeanti, lgsize);
504 const Double_t ghostanti = IsGhost(m3lanti, slopeanti, lgsize);
505
506 const Double_t panti = ghostanti ? -dispanti : dispanti;
507
508 // Align disp along source direction depending on third moment
509 const Double_t thetasqanti = GetThetaSq(panti, distanti, alphaanti);
510
511 if (thetasqanti<cut)
512 return kTRUE;
513
514 // This cut throws away gamma-like events which would be assigned to the
515 // anti-source. It is not necessary but it offers the possibility to
516 // fill the MHDisp plot in one run (Having a hole instead of eg. Crab
517 // the data can be rotated and subtracted)
518 }
519
520 fResult = kTRUE;
521
522 return kTRUE;
523}
524
525void MFMagicCuts::SetVariables(const TArrayD &arr)
526{
527 fVariables = arr;
528}
529
530TString MFMagicCuts::GetParam(Int_t i) const
531{
532 if (i>=fVariables.GetSize())
533 return "";
534
535 return Form("Param[%2d]=%+f", i, fVariables[i]);
536}
537
538void MFMagicCuts::Print(Option_t *o) const
539{
540 *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;
541
542 *fLog << "Using Theta cut: ";
543 switch (fThetaCut)
544 {
545 case kNone:
546 *fLog << "none" << endl;
547 break;
548 case kOn:
549 *fLog << "on" << endl;
550 break;
551 case kOff:
552 *fLog << "off" << endl;
553 break;
554 case kWobble:
555 *fLog << "on and off (wobble)" << endl;
556 break;
557 }
558 *fLog << "Using hadronness cut: ";
559 switch (fHadronnessCut)
560 {
561 case kNoCut:
562 *fLog << "none" << endl;
563 break;
564 case kArea:
565 *fLog << "area" << endl;
566 break;
567 case kHadronness:
568 *fLog << "hadronness" << endl;
569 break;
570 case kAreaLin:
571 *fLog << "arealin" << endl;
572 break;
573 case kAll:
574 *fLog << "all" << endl;
575 break;
576 }
577 if (fCalcDisp)
578 *fLog << "Disp is calculated from c0, c8-c11." << endl;
579 else
580 *fLog << "Disp.fVal from the parameter list is used as disp." << endl;
581 if (fCalcGhostbuster)
582 *fLog << "Ghostbusting is calculated from c5-c7." << endl;
583 else
584 *fLog << "Ghostbuster.fVal from the parameter list is used for ghostbusting." << endl;
585
586 *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
587
588 const Int_t n = fVariables.GetSize();
589 for (int i=0; i<n; i+=3)
590 {
591 *fLog << setw(25) << GetParam(i);
592 *fLog << setw(25) << GetParam(i+1);
593 *fLog << setw(25) << GetParam(i+2);
594 *fLog << endl;
595 }
596 *fLog << setiosflags(ios::right);
597}
598
599Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
600{
601 ifstream fin(fname);
602 if (!fin)
603 {
604 *fLog << err << "Cannot open file " << fname << ": ";
605 *fLog << strerror(errno) << endl;
606 return kFALSE;
607 }
608
609 for (int i=0; i<fVariables.GetSize(); i++)
610 {
611 fin >> fVariables[i];
612 if (!fin)
613 {
614 *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
615 return kERROR;
616 }
617 }
618 return kTRUE;
619}
620
621Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
622{
623 ofstream fout(fname);
624 if (!fout)
625 {
626 *fLog << err << "Cannot open file " << fname << ": ";
627 *fLog << strerror(errno) << endl;
628 return kFALSE;
629 }
630
631 for (int i=0; i<fVariables.GetSize(); i++)
632 fout << setw(11) << fVariables[i] << endl;
633
634 return kTRUE;
635}
636
637Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
638{
639 if (MFilter::ReadEnv(env, prefix, print)==kERROR)
640 return kERROR;
641
642
643 Bool_t rc = kFALSE;
644 if (IsEnvDefined(env, prefix, "ThetaCut", print))
645 {
646 TString str = GetEnvValue(env, prefix, "ThetaCut", "");
647 str.ToLower();
648 str=str.Strip(TString::kBoth);
649
650 if (str==(TString)"nocut" || str==(TString)"none")
651 fThetaCut = kNone;
652 if (str==(TString)"on")
653 fThetaCut = kOn;
654 if (str==(TString)"off")
655 fThetaCut = kOff;
656 if (str==(TString)"wobble")
657 fThetaCut = kWobble;
658 rc = kTRUE;
659 }
660
661 if (IsEnvDefined(env, prefix, "HadronnessCut", print))
662 {
663 TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
664 str.ToLower();
665 str=str.Strip(TString::kBoth);
666
667 if (str==(TString)"nocut" || str==(TString)"none")
668 fHadronnessCut = kNoCut;
669 if (str==(TString)"area")
670 fHadronnessCut = kArea;
671 if (str==(TString)"hadronness")
672 fHadronnessCut = kHadronness;
673 if (str==(TString)"all")
674 fHadronnessCut = kAll;
675
676 rc = kTRUE;
677 }
678
679 if (IsEnvDefined(env, prefix, "CalcDisp", print))
680 {
681 fCalcDisp = GetEnvValue(env, prefix, "CalcDisp", fCalcDisp);
682 rc = kTRUE;
683 }
684
685 if (IsEnvDefined(env, prefix, "CalcGhostbuster", print))
686 {
687 fCalcGhostbuster = GetEnvValue(env, prefix, "CalcGhostbuster", fCalcGhostbuster);
688 rc = kTRUE;
689 }
690
691 if (IsEnvDefined(env, prefix, "File", print))
692 {
693 const TString fname = GetEnvValue(env, prefix, "File", "");
694 if (!CoefficentsRead(fname))
695 return kERROR;
696
697 return kTRUE;
698 }
699
700 for (int i=0; i<fVariables.GetSize(); i++)
701 {
702 if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
703 {
704 // It is important that the last argument is a floating point
705 fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
706 rc = kTRUE;
707 }
708 }
709 return rc;
710}
Note: See TracBrowser for help on using the repository browser.