source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc@ 3324

Last change on this file since 3324 was 3312, checked in by gaug, 21 years ago
*** empty log message ***
File size: 13.2 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): Markus Gaug 02/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MCalibrationChargeCalc
28//
29// Task to calculate the calibration conversion factors from the FADC
30// time slices. The integrated time slices have to be delivered by an
31// MExtractedSignalCam. The pedestals by an MPedestalCam.
32//
33// The output container MCalibrationCam holds one entry of type MCalibrationChargePix
34// for every pixel. It is filled in the following way:
35//
36// ProProcess: Search for MPedestalCam, MExtractedSignalCam and MExtractedSignalBlindPixel
37// Initialize MCalibrationCam
38// Initialize pulser light wavelength
39//
40// ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates
41// memory in a TClonesArray of type MCalibrationChargePix
42// Initialize number of used FADC slices
43// Optionally exclude pixels from calibration
44//
45// Process: Every MCalibrationChargePix holds a histogram class,
46// MHCalibrationPixel which itself hold histograms of type:
47// HCharge(npix) (distribution of summed FADC time slice
48// entries)
49// HTime(npix) (distribution of position of maximum)
50// HChargevsN(npix) (distribution of charges vs. event number.
51//
52// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian
53// All histograms HTime(npix) are fitted to a Gaussian
54// The histogram HBlindPixelCharge (blind pixel) is fitted to
55// a single PhE fit
56//
57// The histograms of the PIN Diode are fitted to Gaussians
58//
59// Fits can be excluded via the commands:
60// MalibrationCam::SkipBlindPixelFits() (skip all blind
61// pixel fits)
62//
63// Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
64// can be skipped with the command:
65// MalibrationCam::SkipHiLoGainCalibration()
66//
67// Input Containers:
68// MRawEvtData
69// MPedestalCam
70// MExtractedSignalCam
71// MExtractedSignalBlindPixel
72//
73// Output Containers:
74// MCalibrationCam
75//
76//////////////////////////////////////////////////////////////////////////////
77#include "MCalibrationChargeCalc.h"
78
79// FXIME: Usage of fstream is a preliminary workaround!
80#include <fstream>
81
82#include <TSystem.h>
83#include <TH1.h>
84
85#include "MLog.h"
86#include "MLogManip.h"
87
88#include "MParList.h"
89
90#include "MGeomCam.h"
91#include "MRawRunHeader.h"
92#include "MRawEvtPixelIter.h"
93
94#include "MPedestalCam.h"
95#include "MPedestalPix.h"
96
97#include "MCalibrationChargeCam.h"
98#include "MCalibrationChargePix.h"
99#include "MCalibrationChargePINDiode.h"
100#include "MCalibrationChargeBlindPix.h"
101
102#include "MExtractedSignalCam.h"
103#include "MExtractedSignalPix.h"
104
105
106ClassImp(MCalibrationChargeCalc);
107
108using namespace std;
109
110// --------------------------------------------------------------------------
111//
112// Default constructor.
113//
114MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
115 : fPedestals(NULL), fCam(NULL),
116 fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL),
117 fSignals(NULL), fPINDiode(NULL), fBlindPixel(NULL)
118{
119
120 fName = name ? name : "MCalibrationChargeCalc";
121 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
122
123 AddToBranchList("MRawEvtData.fHiGainPixId");
124 AddToBranchList("MRawEvtData.fLoGainPixId");
125 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
126 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
127
128 Clear();
129}
130
131void MCalibrationChargeCalc::Clear(const Option_t *o)
132{
133
134 SETBIT(fFlags, kUseQualityChecks);
135 SETBIT(fFlags, kHiLoGainCalibration);
136
137 fNumHiGainSamples = 0;
138 fNumLoGainSamples = 0;
139 fConversionHiLo = 0;
140 fNumExcludedPixels = 0;
141
142}
143
144
145// --------------------------------------------------------------------------
146//
147// The PreProcess searches for the following input containers:
148// - MRawEvtData
149// - MPedestalCam
150//
151// The following output containers are also searched and created if
152// they were not found:
153//
154// - MCalibrationCam
155//
156// The following output containers are only searched, but not created
157//
158// - MTime
159//
160Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
161{
162
163 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
164 if (!fRawEvt)
165 {
166 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
167 return kFALSE;
168 }
169
170 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
171 if (!runheader)
172 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
173 else
174 if (runheader->GetRunType() == kRTMonteCarlo)
175 {
176 return kTRUE;
177 }
178
179 fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
180 if (!fCam)
181 {
182 *fLog << err << dbginf << "MCalibrationChargeCam could not be created ... aborting." << endl;
183 return kFALSE;
184 }
185
186 fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
187 if (!fPINDiode)
188 {
189 *fLog << err << dbginf << "MCalibrationChargePINDiode could not be created ... aborting." << endl;
190 return kFALSE;
191 }
192
193 fCam->SetPINDiode(fPINDiode);
194
195 fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
196 if (!fBlindPixel)
197 {
198 *fLog << err << dbginf << "MCalibrationChargeBlindPix could not be created ... aborting." << endl;
199 return kFALSE;
200 }
201
202 fCam->SetBlindPixel(fBlindPixel);
203
204 fEvtTime = (MTime*)pList->FindObject("MTime");
205
206 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
207 if (!fPedestals)
208 {
209 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
210 return kFALSE;
211 }
212
213
214 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
215 if (!fSignals)
216 {
217 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
218 return kFALSE;
219 }
220
221 return kTRUE;
222}
223
224
225// --------------------------------------------------------------------------
226//
227// The ReInit searches for the following input containers:
228// - MRawRunHeader
229//
230Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
231{
232
233 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
234 if (!fRunHeader)
235 {
236 *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
237 return kFALSE;
238 }
239
240
241 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
242 if (!cam)
243 {
244 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
245 return kFALSE;
246 }
247
248 fCam->SetGeomCam(cam);
249
250 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
251 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
252 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
253
254 UInt_t npixels = cam->GetNumPixels();
255
256 for (UInt_t i=0; i<npixels; i++)
257 {
258
259 MCalibrationChargePix &pix = (*fCam)[i];
260 pix.DefinePixId(i);
261
262 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
263 fSignals->GetLastUsedSliceHiGain());
264 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
265 fSignals->GetLastUsedSliceLoGain());
266
267 }
268
269 //
270 // Look for file to exclude pixels from analysis
271 //
272 if (!fExcludedPixelsFile.IsNull())
273 {
274
275 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
276
277 //
278 // Initialize reading the file
279 //
280 ifstream in(fExcludedPixelsFile.Data(),ios::in);
281
282 if (in)
283 {
284 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
285 //
286 // Read the file and count the number of entries
287 //
288 UInt_t pixel = 0;
289
290 while (++fNumExcludedPixels)
291 {
292
293 in >> pixel;
294
295 if (!in.good())
296 break;
297 //
298 // Check for out of range
299 //
300 if (pixel > npixels)
301 {
302 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
303 << " is out of range " << endl;
304 continue;
305 }
306 //
307 // Exclude pixel
308 //
309 MCalibrationChargePix &pix = (*fCam)[pixel];
310 pix.SetExcluded();
311
312 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
313
314 }
315
316 if (--fNumExcludedPixels == 0)
317 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
318 << "'" << " is empty " << endl;
319 else
320 fCam->SetNumPixelsExcluded(fNumExcludedPixels);
321
322 }
323 else
324 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
325 }
326
327 return kTRUE;
328}
329
330
331// --------------------------------------------------------------------------
332//
333// Calculate the integral of the FADC time slices and store them as a new
334// pixel in the MCerPhotEvt container.
335//
336Int_t MCalibrationChargeCalc::Process()
337{
338 return kTRUE;
339}
340
341Int_t MCalibrationChargeCalc::PostProcess()
342{
343
344 //
345 // loop over the pedestal events and check if we have calibration
346 //
347 Int_t nvalid = 0;
348 for (UInt_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
349 {
350
351 MCalibrationChargePix &pix = (*fCam)[pixid];
352
353 //
354 // Check if the pixel has been excluded from the fits
355 //
356 if (pix.IsExcluded())
357 continue;
358
359 //
360 // get the pedestals
361 //
362 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
363 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
364 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
365
366 //
367 // set them in the calibration camera
368 //
369 if (pix.IsHiGainSaturation())
370 {
371 pix.SetPedestal(ped * fNumLoGainSamples,
372 prms * TMath::Sqrt((Float_t)fNumLoGainSamples),
373 prms * fNumLoGainSamples / num);
374 pix.SetNumLoGainSamples((Float_t)fNumLoGainSamples);
375 pix.ApplyLoGainConversion();
376 }
377 else
378 {
379 pix.SetPedestal(ped * fNumHiGainSamples,
380 prms * TMath::Sqrt((Float_t)fNumHiGainSamples),
381 prms * fNumHiGainSamples / num);
382 }
383
384 if (!pix.CheckChargeValidity() || !pix.CheckTimeValidity())
385 continue;
386
387 nvalid++;
388
389 if (!pix.CalcReducedSigma())
390 continue;
391
392 pix.CalcFFactorMethod();
393
394 }
395
396 //
397 // The Michele check ...
398 //
399 if (nvalid == 0)
400 {
401 *fLog << err << GetDescriptor() << ": Dear Michele! All pixels have non-valid calibration. "
402 << "Did you forget to fill the histograms (filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
403 return kFALSE;
404 }
405
406 if (!fBlindPixel->CheckChargeFitValidity())
407 {
408 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
409 fCam->SetBlindPixelMethodValid(kFALSE);
410 }
411 else
412 {
413 if (!fBlindPixel->CalcFluxInsidePlexiglass())
414 {
415 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
416 fCam->SetBlindPixelMethodValid(kFALSE);
417 }
418 else
419 {
420 fCam->SetBlindPixelMethodValid(kTRUE);
421 fCam->ApplyBlindPixelCalibration();
422 }
423 }
424
425 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
426 {
427 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
428 fCam->SetPINDiodeMethodValid(kFALSE);
429 }
430 else
431 {
432 if (!fPINDiode->CalcFluxOutsidePlexiglass())
433 {
434 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
435 fCam->SetPINDiodeMethodValid(kFALSE);
436 }
437 else
438 {
439 fCam->SetPINDiodeMethodValid(kTRUE);
440 fCam->ApplyPINDiodeCalibration();
441 }
442 }
443 fCam->SetReadyToSave();
444
445 return kTRUE;
446}
447
448
449
450
451
452
453
Note: See TracBrowser for help on using the repository browser.