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*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 |
|
---|
147 | ClassImp(MFMagicCuts);
|
---|
148 |
|
---|
149 | using namespace std;
|
---|
150 |
|
---|
151 | // --------------------------------------------------------------------------
|
---|
152 | //
|
---|
153 | // constructor
|
---|
154 | //
|
---|
155 | MFMagicCuts::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 | //
|
---|
188 | Int_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 | fMm2Deg = cam->GetConvMm2Deg();
|
---|
198 |
|
---|
199 | fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
|
---|
200 | if (!fThetaSq)
|
---|
201 | return kFALSE;
|
---|
202 |
|
---|
203 | if (!fCalcDisp)
|
---|
204 | fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD");
|
---|
205 | else
|
---|
206 | fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
|
---|
207 | if (!fDisp)
|
---|
208 | {
|
---|
209 | *fLog << err << "Disp [MParameterD] not found... aborting." << endl;
|
---|
210 | return kFALSE;
|
---|
211 | }
|
---|
212 |
|
---|
213 | if (!fCalcGhostbuster)
|
---|
214 | fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD");
|
---|
215 | else
|
---|
216 | fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster");
|
---|
217 | if (!fGhostbuster)
|
---|
218 | {
|
---|
219 | *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl;
|
---|
220 | return kFALSE;
|
---|
221 | }
|
---|
222 |
|
---|
223 | // propagate Theta cut to the parameter list
|
---|
224 | // so that it can be used somewhere else
|
---|
225 | MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
|
---|
226 | if (!thetacut)
|
---|
227 | return kFALSE;
|
---|
228 | thetacut->SetVal(GetThetaSqCut());
|
---|
229 |
|
---|
230 | Print();
|
---|
231 |
|
---|
232 | if (fMatrix)
|
---|
233 | return kTRUE;
|
---|
234 |
|
---|
235 | //-----------------------------------------------------------
|
---|
236 |
|
---|
237 | fHil = (MHillas*)pList->FindObject("MHillas");
|
---|
238 | if (!fHil)
|
---|
239 | {
|
---|
240 | *fLog << err << "MHillas not found... aborting." << endl;
|
---|
241 | return kFALSE;
|
---|
242 | }
|
---|
243 |
|
---|
244 | fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
|
---|
245 | if (!fHilExt)
|
---|
246 | {
|
---|
247 | *fLog << err << "MHillasExt not found... aborting." << endl;
|
---|
248 | return kFALSE;
|
---|
249 | }
|
---|
250 |
|
---|
251 | fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
|
---|
252 | if (!fNewImgPar)
|
---|
253 | {
|
---|
254 | *fLog << err << "MNewImagePar not found... aborting." << endl;
|
---|
255 | return kFALSE;
|
---|
256 | }
|
---|
257 |
|
---|
258 | fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
|
---|
259 | if (!fHilSrc)
|
---|
260 | {
|
---|
261 | *fLog << err << "MHillasSrc not found... aborting." << endl;
|
---|
262 | return kFALSE;
|
---|
263 | }
|
---|
264 |
|
---|
265 | if (fThetaCut&kOff)
|
---|
266 | {
|
---|
267 | fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
|
---|
268 | if (!fHilAnti)
|
---|
269 | {
|
---|
270 | *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
|
---|
271 | return kFALSE;
|
---|
272 | }
|
---|
273 | }
|
---|
274 |
|
---|
275 | if (fHadronnessCut&kHadronness)
|
---|
276 | {
|
---|
277 | fHadronness = (MParameterD*)pList->FindObject("Hadronness", "MParameterD");
|
---|
278 | if (!fHadronness)
|
---|
279 | {
|
---|
280 | *fLog << warn << "Hadronness [MParameterD] not found... aborting." << endl;
|
---|
281 | return kFALSE;
|
---|
282 | }
|
---|
283 | }
|
---|
284 |
|
---|
285 | return kTRUE;
|
---|
286 | }
|
---|
287 |
|
---|
288 | // --------------------------------------------------------------------------
|
---|
289 | //
|
---|
290 | // Returns the mapped value from the Matrix
|
---|
291 | //
|
---|
292 | Double_t MFMagicCuts::GetVal(Int_t i) const
|
---|
293 | {
|
---|
294 | return (*fMatrix)[fMap[i]];
|
---|
295 | }
|
---|
296 |
|
---|
297 | // --------------------------------------------------------------------------
|
---|
298 | //
|
---|
299 | // You can use this function if you want to use a MHMatrix instead of the
|
---|
300 | // given containers. This function adds all necessary columns to the
|
---|
301 | // given matrix. Afterward you should fill the matrix with the corresponding
|
---|
302 | // data (eg from a file by using MHMatrix::Fill). If you now loop
|
---|
303 | // through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
|
---|
304 | // will take the values from the matrix instead of the containers.
|
---|
305 | //
|
---|
306 | void MFMagicCuts::InitMapping(MHMatrix *mat)
|
---|
307 | {
|
---|
308 | if (fMatrix)
|
---|
309 | return;
|
---|
310 |
|
---|
311 | fMatrix = mat;
|
---|
312 |
|
---|
313 | // fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
|
---|
314 | // fMap[kEWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
|
---|
315 |
|
---|
316 | fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
|
---|
317 | fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)");
|
---|
318 | fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
|
---|
319 |
|
---|
320 | fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");
|
---|
321 |
|
---|
322 | fMap[kEAlpha] = fMatrix->AddColumn("MHillasSrc.fAlpha");
|
---|
323 | fMap[kEDist] = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
|
---|
324 | fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
|
---|
325 |
|
---|
326 | fMap[kESign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
|
---|
327 |
|
---|
328 | fMap[kESlope] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
|
---|
329 |
|
---|
330 | if (!fCalcDisp)
|
---|
331 | fMap[kEDisp] = fMatrix->AddColumn("abs(Disp.fVal)");
|
---|
332 |
|
---|
333 | if (!fCalcGhostbuster)
|
---|
334 | fMap[kEGhostbuster] = fMatrix->AddColumn("Ghostbuster.fVal");
|
---|
335 |
|
---|
336 | if (fThetaCut&kOff)
|
---|
337 | {
|
---|
338 | fMap[kEAlphaAnti] = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
|
---|
339 | fMap[kEDistAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
|
---|
340 | fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
|
---|
341 | fMap[kESlopeAnti] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
|
---|
342 | }
|
---|
343 |
|
---|
344 | if (fHadronnessCut&kHadronness)
|
---|
345 | fMap[kEHadronness] = fMatrix->AddColumn("Hadronness.fVal");
|
---|
346 | }
|
---|
347 |
|
---|
348 | // --------------------------------------------------------------------------
|
---|
349 | //
|
---|
350 | // Returns theta squared based on the formula:
|
---|
351 | // p := c*(1-width/length)
|
---|
352 | // return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
|
---|
353 | //
|
---|
354 | // If par!=NULL p is stored in this MParameterD
|
---|
355 | //
|
---|
356 | Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t d, Double_t a) const
|
---|
357 | {
|
---|
358 | // wl = width/l [1]
|
---|
359 | // d = dist [deg]
|
---|
360 | // a = alpha [deg]
|
---|
361 | return d*d + p*p - 2*d*p*TMath::Cos(a*TMath::DegToRad());
|
---|
362 | }
|
---|
363 |
|
---|
364 | // --------------------------------------------------------------------------
|
---|
365 | //
|
---|
366 | // Returns squared cut value in theta: fVariables[1]^2
|
---|
367 | //
|
---|
368 | Double_t MFMagicCuts::GetThetaSqCut() const
|
---|
369 | {
|
---|
370 | return fVariables[1]*fVariables[1];
|
---|
371 | }
|
---|
372 |
|
---|
373 | // --------------------------------------------------------------------------
|
---|
374 | //
|
---|
375 | // Return abs(Disp.fVal) if disp calculation is switched off.
|
---|
376 | // Otherwise return (c0+c6*leak^c7)*(1-width/length)
|
---|
377 | //
|
---|
378 | Double_t MFMagicCuts::GetDisp(Double_t slope, Double_t lgsize) const
|
---|
379 | {
|
---|
380 | if (!fCalcDisp)
|
---|
381 | return fMatrix ? GetVal(kEDisp) : TMath::Abs(fDisp->GetVal());
|
---|
382 |
|
---|
383 | // Get some variables
|
---|
384 | const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
|
---|
385 | const Double_t leak = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
|
---|
386 |
|
---|
387 | // For simplicity
|
---|
388 | const Double_t *c = fVariables.GetArray();
|
---|
389 |
|
---|
390 | // As rule for root or MDataPhrase:
|
---|
391 | // ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0]
|
---|
392 | //
|
---|
393 | Double_t xi = c[0] + c[8]*slope + c[9]*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 |
|
---|
402 | Bool_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 | //
|
---|
421 | Int_t MFMagicCuts::Process()
|
---|
422 | {
|
---|
423 | // For simplicity
|
---|
424 | const Double_t *c = fVariables.GetArray();
|
---|
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 | // Calculate corresponding Theta^2
|
---|
451 | const Double_t thetasq = GetThetaSq(p, dist, alpha);
|
---|
452 | fThetaSq->SetVal(thetasq);
|
---|
453 |
|
---|
454 | // Perform theta cut
|
---|
455 | if (fThetaCut&kOn && thetasq >= cut)
|
---|
456 | return kTRUE;
|
---|
457 |
|
---|
458 | // --------------------- Efficient cut -------------------------
|
---|
459 | // ---------------------- Area/M3l cut --------------------------
|
---|
460 | if (fHadronnessCut&kArea)
|
---|
461 | {
|
---|
462 | const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
|
---|
463 | const Double_t A = lgsize>c[3] ? c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
|
---|
464 | if (area>=A)
|
---|
465 | return kTRUE;
|
---|
466 | }
|
---|
467 |
|
---|
468 | // ---------------------------------------------------------------
|
---|
469 | // ---------------------------------------------------------------
|
---|
470 |
|
---|
471 | if (fHadronnessCut&kHadronness)
|
---|
472 | {
|
---|
473 | const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
|
---|
474 | if (had>c[13])
|
---|
475 | return kTRUE;
|
---|
476 |
|
---|
477 | if (TMath::Power(10, lgsize)<c[14])
|
---|
478 | return kTRUE;
|
---|
479 | }
|
---|
480 |
|
---|
481 | // ------------------- Least efficient cut -----------------------
|
---|
482 | // ---------------------- Theta cut OFF --------------------------
|
---|
483 |
|
---|
484 | if (fThetaCut&kOff)
|
---|
485 | {
|
---|
486 | const Double_t distanti = fMatrix ? GetVal(kEDistAnti) : fHilAnti->GetDist()*fMm2Deg;
|
---|
487 | const Double_t alphaanti = fMatrix ? GetVal(kEAlphaAnti) : fHilAnti->GetAlpha();
|
---|
488 | const Double_t m3lanti = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
|
---|
489 | const Double_t slopeanti = fMatrix ? GetVal(kESlopeAnti) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
|
---|
490 |
|
---|
491 | const Double_t dispanti = GetDisp(slopeanti, lgsize);
|
---|
492 | const Double_t ghostanti = IsGhost(m3lanti, slopeanti, lgsize);
|
---|
493 |
|
---|
494 | const Double_t panti = ghostanti ? -dispanti : dispanti;
|
---|
495 |
|
---|
496 | // Align disp along source direction depending on third moment
|
---|
497 | const Double_t thetasqanti = GetThetaSq(panti, distanti, alphaanti);
|
---|
498 |
|
---|
499 | if (thetasqanti<cut)
|
---|
500 | return kTRUE;
|
---|
501 |
|
---|
502 | // This cut throws away gamma-like events which would be assigned to the
|
---|
503 | // anti-source. It is not necessary but it offers the possibility to
|
---|
504 | // fill the MHDisp plot in one run (Having a hole instead of eg. Crab
|
---|
505 | // the data can be rotated and subtracted)
|
---|
506 | }
|
---|
507 |
|
---|
508 | fResult = kTRUE;
|
---|
509 |
|
---|
510 | return kTRUE;
|
---|
511 | }
|
---|
512 |
|
---|
513 | void MFMagicCuts::SetVariables(const TArrayD &arr)
|
---|
514 | {
|
---|
515 | fVariables = arr;
|
---|
516 | }
|
---|
517 |
|
---|
518 | TString MFMagicCuts::GetParam(Int_t i) const
|
---|
519 | {
|
---|
520 | if (i>=fVariables.GetSize())
|
---|
521 | return "";
|
---|
522 |
|
---|
523 | return Form("Param[%2d]=%+f", i, fVariables[i]);
|
---|
524 | }
|
---|
525 |
|
---|
526 | void MFMagicCuts::Print(Option_t *o) const
|
---|
527 | {
|
---|
528 | *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;
|
---|
529 |
|
---|
530 | *fLog << "Using Theta cut: ";
|
---|
531 | switch (fThetaCut)
|
---|
532 | {
|
---|
533 | case kNone:
|
---|
534 | *fLog << "none" << endl;
|
---|
535 | break;
|
---|
536 | case kOn:
|
---|
537 | *fLog << "on" << endl;
|
---|
538 | break;
|
---|
539 | case kOff:
|
---|
540 | *fLog << "off" << endl;
|
---|
541 | break;
|
---|
542 | case kWobble:
|
---|
543 | *fLog << "on and off (wobble)" << endl;
|
---|
544 | break;
|
---|
545 | }
|
---|
546 | *fLog << "Using hadronness cut: ";
|
---|
547 | switch (fHadronnessCut)
|
---|
548 | {
|
---|
549 | case kNoCut:
|
---|
550 | *fLog << "none" << endl;
|
---|
551 | break;
|
---|
552 | case kArea:
|
---|
553 | *fLog << "area" << endl;
|
---|
554 | break;
|
---|
555 | case kHadronness:
|
---|
556 | *fLog << "hadronness" << endl;
|
---|
557 | break;
|
---|
558 | case kAll:
|
---|
559 | *fLog << "all" << endl;
|
---|
560 | break;
|
---|
561 | }
|
---|
562 | if (fCalcDisp)
|
---|
563 | *fLog << "Disp is calculated from c0, c7, c8." << endl;
|
---|
564 | else
|
---|
565 | *fLog << "Disp.fVal from the parameter list is used as disp." << endl;
|
---|
566 | if (fCalcGhostbuster)
|
---|
567 | *fLog << "Ghostbusting is calculated from c5, c6, c7." << endl;
|
---|
568 | else
|
---|
569 | *fLog << "Ghostbuster.fVal from the parameter list is used for ghostbusting." << endl;
|
---|
570 |
|
---|
571 | *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
|
---|
572 |
|
---|
573 | const Int_t n = fVariables.GetSize();
|
---|
574 | for (int i=0; i<n; i+=3)
|
---|
575 | {
|
---|
576 | *fLog << setw(25) << GetParam(i);
|
---|
577 | *fLog << setw(25) << GetParam(i+1);
|
---|
578 | *fLog << setw(25) << GetParam(i+2);
|
---|
579 | *fLog << endl;
|
---|
580 | }
|
---|
581 | *fLog << setiosflags(ios::right);
|
---|
582 | }
|
---|
583 |
|
---|
584 | Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
|
---|
585 | {
|
---|
586 | ifstream fin(fname);
|
---|
587 | if (!fin)
|
---|
588 | {
|
---|
589 | *fLog << err << "Cannot open file " << fname << ": ";
|
---|
590 | *fLog << strerror(errno) << endl;
|
---|
591 | return kFALSE;
|
---|
592 | }
|
---|
593 |
|
---|
594 | for (int i=0; i<fVariables.GetSize(); i++)
|
---|
595 | {
|
---|
596 | fin >> fVariables[i];
|
---|
597 | if (!fin)
|
---|
598 | {
|
---|
599 | *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
|
---|
600 | return kERROR;
|
---|
601 | }
|
---|
602 | }
|
---|
603 | return kTRUE;
|
---|
604 | }
|
---|
605 |
|
---|
606 | Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
|
---|
607 | {
|
---|
608 | ofstream fout(fname);
|
---|
609 | if (!fout)
|
---|
610 | {
|
---|
611 | *fLog << err << "Cannot open file " << fname << ": ";
|
---|
612 | *fLog << strerror(errno) << endl;
|
---|
613 | return kFALSE;
|
---|
614 | }
|
---|
615 |
|
---|
616 | for (int i=0; i<fVariables.GetSize(); i++)
|
---|
617 | fout << setw(11) << fVariables[i] << endl;
|
---|
618 |
|
---|
619 | return kTRUE;
|
---|
620 | }
|
---|
621 |
|
---|
622 | Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
---|
623 | {
|
---|
624 | if (MFilter::ReadEnv(env, prefix, print)==kERROR)
|
---|
625 | return kERROR;
|
---|
626 |
|
---|
627 |
|
---|
628 | Bool_t rc = kFALSE;
|
---|
629 | if (IsEnvDefined(env, prefix, "ThetaCut", print))
|
---|
630 | {
|
---|
631 | TString str = GetEnvValue(env, prefix, "ThetaCut", "");
|
---|
632 | str.ToLower();
|
---|
633 | str=str.Strip(TString::kBoth);
|
---|
634 |
|
---|
635 | if (str==(TString)"nocut" || str==(TString)"none")
|
---|
636 | fThetaCut = kNone;
|
---|
637 | if (str==(TString)"on")
|
---|
638 | fThetaCut = kOn;
|
---|
639 | if (str==(TString)"off")
|
---|
640 | fThetaCut = kOff;
|
---|
641 | if (str==(TString)"wobble")
|
---|
642 | fThetaCut = kWobble;
|
---|
643 | rc = kTRUE;
|
---|
644 | }
|
---|
645 |
|
---|
646 | if (IsEnvDefined(env, prefix, "HadronnessCut", print))
|
---|
647 | {
|
---|
648 | TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
|
---|
649 | str.ToLower();
|
---|
650 | str=str.Strip(TString::kBoth);
|
---|
651 |
|
---|
652 | if (str==(TString)"nocut" || str==(TString)"none")
|
---|
653 | fHadronnessCut = kNoCut;
|
---|
654 | if (str==(TString)"area")
|
---|
655 | fHadronnessCut = kArea;
|
---|
656 | if (str==(TString)"hadronness")
|
---|
657 | fHadronnessCut = kHadronness;
|
---|
658 | if (str==(TString)"all")
|
---|
659 | fHadronnessCut = kAll;
|
---|
660 |
|
---|
661 | rc = kTRUE;
|
---|
662 | }
|
---|
663 |
|
---|
664 | if (IsEnvDefined(env, prefix, "CalcDisp", print))
|
---|
665 | {
|
---|
666 | fCalcDisp = GetEnvValue(env, prefix, "CalcDisp", fCalcDisp);
|
---|
667 | rc = kTRUE;
|
---|
668 | }
|
---|
669 |
|
---|
670 | if (IsEnvDefined(env, prefix, "CalcGhostbuster", print))
|
---|
671 | {
|
---|
672 | fCalcGhostbuster = GetEnvValue(env, prefix, "CalcGhostbuster", fCalcGhostbuster);
|
---|
673 | rc = kTRUE;
|
---|
674 | }
|
---|
675 |
|
---|
676 | if (IsEnvDefined(env, prefix, "File", print))
|
---|
677 | {
|
---|
678 | const TString fname = GetEnvValue(env, prefix, "File", "");
|
---|
679 | if (!CoefficentsRead(fname))
|
---|
680 | return kERROR;
|
---|
681 |
|
---|
682 | return kTRUE;
|
---|
683 | }
|
---|
684 |
|
---|
685 | for (int i=0; i<fVariables.GetSize(); i++)
|
---|
686 | {
|
---|
687 | if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
|
---|
688 | {
|
---|
689 | // It is important that the last argument is a floating point
|
---|
690 | fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
|
---|
691 | rc = kTRUE;
|
---|
692 | }
|
---|
693 | }
|
---|
694 | return rc;
|
---|
695 | }
|
---|