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

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