source: trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc@ 6977

Last change on this file since 6977 was 5955, checked in by domingo, 20 years ago
*** empty log message ***
File size: 9.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): Eva Domingo, 12/2004 <mailto:domingo@ifae.es>
19! Wolfgang Wittek, 12/2004 <mailto:wittek@mppmu.mpg.de>
20!
21!
22! Copyright: MAGIC Software Development, 2000-2005
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28// //
29// MFDisp //
30// //
31// This is a class to set cuts to select an event sample //
32// for the Disp optimization //
33// //
34// //
35//////////////////////////////////////////////////////////////////////////////
36
37#include "MFDisp.h"
38
39#include "MGeomCam.h"
40#include "MHillas.h"
41#include "MHillasSrc.h"
42#include "MImagePar.h"
43#include "MNewImagePar.h"
44
45#include "MParList.h"
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50
51ClassImp(MFDisp);
52
53using namespace std;
54
55// --------------------------------------------------------------------------
56//
57// Default constructor.
58//
59MFDisp::MFDisp(const char *name, const char *title)
60{
61 fName = name ? name : "MFDisp";
62 fTitle = title ? title : "Cuts for Disp optimization";
63
64 // default values of cuts
65 SetCuts(0, 1000, 0, 1000, 0, 1000,
66 0.0, 1.e10, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
67}
68
69
70// --------------------------------------------------------------------------
71//
72// Set the values for the cuts
73//
74void MFDisp::SetCuts(Int_t islandsmin, Int_t islandsmax,
75 Int_t usedpixelsmin, Int_t usedpixelsmax,
76 Int_t corepixelsmin, Int_t corepixelsmax,
77 Float_t sizemin, Float_t sizemax,
78 Float_t leakage1min, Float_t leakage1max,
79 Float_t leakage2min, Float_t leakage2max,
80 Float_t lengthmin, Float_t widthmin)
81{
82 fIslandsMin = islandsmin;
83 fIslandsMax = islandsmax;
84
85 fUsedPixelsMin = usedpixelsmin;
86 fUsedPixelsMax = usedpixelsmax;
87
88 fCorePixelsMin = corepixelsmin;
89 fCorePixelsMax = corepixelsmax;
90
91 fSizeMin = sizemin;
92 fSizeMax = sizemax;
93
94 fLeakage1Min = leakage1min;
95 fLeakage1Max = leakage1max;
96
97 fLeakage2Min = leakage2min;
98 fLeakage2Max = leakage2max;
99
100 fLengthMin = lengthmin;
101 fWidthMin = widthmin;
102}
103
104
105// --------------------------------------------------------------------------
106//
107Int_t MFDisp::PreProcess(MParList *pList)
108{
109 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
110 if (!cam)
111 {
112 *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
113 return kFALSE;
114 }
115
116 fMm2Deg = cam->GetConvMm2Deg();
117
118
119 fHil = (MHillas*)pList->FindObject("MHillas");
120 if (!fHil)
121 {
122 *fLog << err << "MHillas not found... aborting." << endl;
123 return kFALSE;
124 }
125
126 fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
127 if (!fHilSrc)
128 {
129 *fLog << err << "MHillasSrc not found... aborting." << endl;
130 return kFALSE;
131 }
132
133 fImgPar = (MImagePar*)pList->FindObject("MImagePar");
134 if (!fImgPar)
135 {
136 *fLog << err << "MImagePar not found... aborting." << endl;
137 return kFALSE;
138 }
139
140 fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
141 if (!fNewImgPar)
142 {
143 *fLog << err << "MNewImagePar not found... aborting." << endl;
144 return kFALSE;
145 }
146
147
148 memset(fCut, 0, sizeof(fCut));
149
150
151 //--------------------
152 *fLog << inf << "MFDisp cut values :" << endl;
153
154 *fLog << inf << " fIslandsMin, fIslandsMax = "
155 << fIslandsMin << ", " << fIslandsMax << endl;
156
157 *fLog << inf << " fUsedPixelsMin, fUsedPixelsMax = "
158 << fUsedPixelsMin << ", " << fUsedPixelsMax << endl;
159
160 *fLog << inf << " fCorePixelsMin, fCorePixelsMax = "
161 << fCorePixelsMin << ", " << fCorePixelsMax << endl;
162
163 *fLog << inf << " fSizeMin, fSizeMax = "
164 << fSizeMin << ", " << fSizeMax << endl;
165
166 *fLog << inf << " fLeakage1Min, fLeakage1Max = "
167 << fLeakage1Min << ", " << fLeakage1Max << endl;
168
169 *fLog << inf << " fLeakage2Min, fLeakage2Max = "
170 << fLeakage2Min << ", " << fLeakage2Max << endl;
171
172 *fLog << inf << " fLengthMin, fWidthMin = "
173 << fLengthMin << ", " << fWidthMin << endl;
174 //--------------------
175
176 return kTRUE;
177}
178
179
180// --------------------------------------------------------------------------
181//
182// Sum up the number of events passing each cut
183//
184Bool_t MFDisp::Set(Int_t rc)
185{
186 fResult = kTRUE;
187 fCut[rc]++;
188 return kTRUE;
189}
190
191
192// --------------------------------------------------------------------------
193//
194// Evaluate cuts for Disp optimization
195// if selections are fulfilled set fResult = kTRUE;
196//
197Int_t MFDisp::Process()
198{
199 const Double_t length = fHil->GetLength() * fMm2Deg;
200 const Double_t width = fHil->GetWidth() * fMm2Deg;
201 //const Double_t dist = fHilSrc->GetDist()* fMm2Deg;
202 //const Double_t delta = fHil->GetDelta() * kRad2Deg;
203 const Double_t size = fHil->GetSize();
204
205 const Double_t leakage1 = fNewImgPar->GetLeakage1();
206 const Double_t leakage2 = fNewImgPar->GetLeakage2();
207
208 const Int_t numusedpixels = fNewImgPar->GetNumUsedPixels();
209 const Int_t numcorepixels = fNewImgPar->GetNumCorePixels();
210 const Int_t numislands = fImgPar->GetNumIslands();
211
212 fResult = kFALSE;
213
214 if (numislands<fIslandsMin || numislands>fIslandsMax )
215 return Set(1);
216
217 if (numusedpixels<fUsedPixelsMin || numusedpixels>fUsedPixelsMax )
218 return Set(2);
219
220 if (numcorepixels<fCorePixelsMin || numcorepixels>fCorePixelsMax )
221 return Set(3);
222
223 if (size<fSizeMin || size>fSizeMax)
224 return Set(4);
225
226 if (leakage1<fLeakage1Min || leakage1>fLeakage1Max)
227 return Set(5);
228
229 if (leakage2<fLeakage2Min || leakage2>fLeakage2Max)
230 return Set(6);
231
232 if (length<fLengthMin || width<fWidthMin)
233 return Set(7);
234
235 fCut[0]++;
236 return kTRUE;
237}
238
239
240// --------------------------------------------------------------------------
241//
242// Prints some statistics about the cuts selections
243//
244Int_t MFDisp::PostProcess()
245{
246 if (GetNumExecutions()==0)
247 return kTRUE;
248
249 *fLog << inf << endl;
250 *fLog << GetDescriptor() << " execution statistics:" << endl;
251
252 *fLog << dec << setfill(' ');
253 *fLog << " " << setw(7) << fCut[1] << " (" << setw(3);
254 *fLog << (int)(fCut[1]*100/GetNumExecutions()+0.5) ;
255 *fLog << "%) Evts skipped due to: No.of islands < " << fIslandsMin ;
256 *fLog << " or > " << fIslandsMax << endl;
257
258 *fLog << " " << setw(7) << fCut[2] << " (" << setw(3);
259 *fLog << (int)(fCut[2]*100/GetNumExecutions()+0.5) ;
260 *fLog << "%) Evts skipped due to: No.of used pixels < " << fUsedPixelsMin ;
261 *fLog << " or > " << fUsedPixelsMax << endl;
262
263 *fLog << " " << setw(7) << fCut[3] << " (" << setw(3);
264 *fLog << (int)(fCut[3]*100/GetNumExecutions()+0.5) ;
265 *fLog << "%) Evts skipped due to: No.of core pixels < " << fCorePixelsMin ;
266 *fLog << " or > " << fCorePixelsMax << endl;
267
268 *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ;
269 *fLog << (int)(fCut[4]*100/GetNumExecutions()+0.5) ;
270 *fLog << "%) Evts skipped due to: SIZE < " << fSizeMin;
271 *fLog << " or > " << fSizeMax << endl;
272
273 *fLog << " " << setw(7) << fCut[5] << " (" << setw(3) ;
274 *fLog << (int)(fCut[5]*100/GetNumExecutions()+0.5) ;
275 *fLog << "%) Evts skipped due to: Leakage1 < " << fLeakage1Min;
276 *fLog << " or > " << fLeakage1Max << endl;
277
278 *fLog << " " << setw(7) << fCut[6] << " (" << setw(3) ;
279 *fLog << (int)(fCut[6]*100/GetNumExecutions()+0.5) ;
280 *fLog << "%) Evts skipped due to: Leakage2 < " << fLeakage2Min;
281 *fLog << " or > " << fLeakage2Max << endl;
282
283 *fLog << " " << setw(7) << fCut[7] << " (" << setw(3) ;
284 *fLog << (int)(fCut[7]*100/GetNumExecutions()+0.5) ;
285 *fLog << "%) Evts skipped due to: LENGTH <= " << fLengthMin;
286 *fLog << " or WIDTH <= " << fWidthMin << endl;
287
288 *fLog << " " << fCut[0] << " (" ;
289 *fLog << (int)(fCut[0]*100/GetNumExecutions()+0.5) ;
290 *fLog << "%) Evts survived the Disp cuts!" << endl;
291 *fLog << endl;
292
293 return kTRUE;
294}
295
296
297
298
299
Note: See TracBrowser for help on using the repository browser.