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

Last change on this file since 3350 was 3350, checked in by gaug, 21 years ago
*** empty log message ***
File size: 14.8 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 Float_t avped = 0;
349 Float_t avprms = 0;
350 Float_t avnum = 0;
351 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
352 {
353
354 MCalibrationChargePix &pix = (*fCam)[pixid];
355
356 //
357 // Check if the pixel has been excluded from the fits
358 //
359 if (pix.IsExcluded())
360 continue;
361
362 //
363 // get the pedestals
364 //
365 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
366 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
367 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
368
369 avped += ped;
370 avprms += prms;
371 avnum += num;
372 //
373 // set them in the calibration camera
374 //
375 if (pix.IsHiGainSaturation())
376 {
377 pix.SetPedestal(ped * fNumLoGainSamples,
378 prms * TMath::Sqrt((Float_t)fNumLoGainSamples),
379 prms * fNumLoGainSamples / num);
380 pix.SetNumLoGainSamples((Float_t)fNumLoGainSamples);
381 pix.ApplyLoGainConversion();
382 }
383 else
384 {
385 pix.SetPedestal(ped * fNumHiGainSamples,
386 prms * TMath::Sqrt((Float_t)fNumHiGainSamples),
387 prms * fNumHiGainSamples / num);
388 }
389
390 if (!pix.CheckChargeValidity() || !pix.CheckTimeValidity())
391 continue;
392
393 nvalid++;
394
395 if (!pix.CalcReducedSigma())
396 continue;
397
398 pix.CalcFFactorMethod();
399
400 }
401
402
403
404 //
405 // The Michele check ...
406 //
407 if (nvalid == 0)
408 {
409 *fLog << err << GetDescriptor() << ": Dear Michele! All pixels have non-valid calibration. "
410 << "Did you forget to fill the histograms (filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
411 return kFALSE;
412 }
413
414 MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
415 MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
416 //
417 // set the pedestans in the calibration camera
418 //
419 if (avinnerpix->IsHiGainSaturation())
420 {
421 avinnerpix->SetPedestal(avped * fNumLoGainSamples,
422 avprms * TMath::Sqrt((Float_t)fNumLoGainSamples),
423 avprms * fNumLoGainSamples / avnum);
424 avinnerpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
425 avinnerpix->ApplyLoGainConversion();
426 }
427 else
428 {
429 avinnerpix->SetPedestal(avped * fNumHiGainSamples,
430 avprms * TMath::Sqrt((Float_t)fNumHiGainSamples),
431 avprms * fNumHiGainSamples / avnum);
432 }
433
434 if (avouterpix->IsHiGainSaturation())
435 {
436 avouterpix->SetPedestal(avped * fNumLoGainSamples,
437 avprms * TMath::Sqrt((Float_t)fNumLoGainSamples),
438 avprms * fNumLoGainSamples / avnum);
439 avouterpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
440 avouterpix->ApplyLoGainConversion();
441 }
442 else
443 {
444 avouterpix->SetPedestal(avped * fNumHiGainSamples,
445 avprms * TMath::Sqrt((Float_t)fNumHiGainSamples),
446 avprms * fNumHiGainSamples / avnum);
447 }
448
449 if (!avinnerpix->CheckChargeValidity() || !avinnerpix->CheckTimeValidity())
450 if (!avinnerpix->CalcReducedSigma())
451 avinnerpix->CalcFFactorMethod();
452
453 if (!avouterpix->CheckChargeValidity() || !avouterpix->CheckTimeValidity())
454 if (!avouterpix->CalcReducedSigma())
455 avouterpix->CalcFFactorMethod();
456
457
458 if (!fBlindPixel->CheckChargeFitValidity())
459 {
460 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
461 fCam->SetBlindPixelMethodValid(kFALSE);
462 }
463 else
464 {
465 if (!fBlindPixel->CalcFluxInsidePlexiglass())
466 {
467 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
468 fCam->SetBlindPixelMethodValid(kFALSE);
469 }
470 else
471 {
472 fCam->SetBlindPixelMethodValid(kTRUE);
473 fCam->ApplyBlindPixelCalibration();
474 }
475 }
476
477 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
478 {
479 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
480 fCam->SetPINDiodeMethodValid(kFALSE);
481 }
482 else
483 {
484 if (!fPINDiode->CalcFluxOutsidePlexiglass())
485 {
486 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
487 fCam->SetPINDiodeMethodValid(kFALSE);
488 }
489 else
490 {
491 fCam->SetPINDiodeMethodValid(kTRUE);
492 fCam->ApplyPINDiodeCalibration();
493 }
494 }
495 fCam->SetReadyToSave();
496
497 return kTRUE;
498}
499
500
501
502
503
504
505
Note: See TracBrowser for help on using the repository browser.