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

Last change on this file since 3394 was 3374, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 14.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 << "MRawEvtData not found... aborting." << endl;
167 return kFALSE;
168 }
169
170 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
171 if (!runheader)
172 *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl;
173 else
174 if (runheader->IsMonteCarloRun())
175 return kTRUE;
176
177 fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
178 if (!fCam)
179 return kFALSE;
180
181 fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
182 if (!fPINDiode)
183 return kFALSE;
184
185 fCam->SetPINDiode(fPINDiode);
186
187 fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
188 if (!fBlindPixel)
189 return kFALSE;
190
191 fCam->SetBlindPixel(fBlindPixel);
192
193 fEvtTime = (MTime*)pList->FindObject("MTime");
194
195 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
196 if (!fPedestals)
197 {
198 *fLog << err << "MPedestalCam not found... aborting" << endl;
199 return kFALSE;
200 }
201
202 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
203 if (!fSignals)
204 {
205 *fLog << err << "MExtractedSignalCam not found... aborting" << endl;
206 return kFALSE;
207 }
208
209 return kTRUE;
210}
211
212
213// --------------------------------------------------------------------------
214//
215// The ReInit searches for the following input containers:
216// - MRawRunHeader
217//
218Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
219{
220 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
221 if (!fRunHeader)
222 {
223 *fLog << err << "MRawRunHeader not found... aborting." << endl;
224 return kFALSE;
225 }
226
227 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
228 if (!fGeom)
229 {
230 *fLog << err << "No MGeomCam found... aborting." << endl;
231 return kFALSE;
232 }
233
234 fCam->SetGeomCam(fGeom);
235
236 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
237 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
238 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
239
240 UInt_t npixels = fGeom->GetNumPixels();
241
242 for (UInt_t i=0; i<npixels; i++)
243 {
244 MCalibrationChargePix &pix = (*fCam)[i];
245 pix.DefinePixId(i);
246
247 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
248 fSignals->GetLastUsedSliceHiGain());
249 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
250 fSignals->GetLastUsedSliceLoGain());
251 }
252
253 if (fExcludedPixelsFile.IsNull())
254 return kTRUE;
255
256 //
257 // Look for file to exclude pixels from analysis
258 //
259 gSystem->ExpandPathName(fExcludedPixelsFile);
260
261 //
262 // Initialize reading the file
263 //
264 ifstream in(fExcludedPixelsFile);
265 if (!in)
266 {
267 *fLog << warn << "Cannot open file '" << fExcludedPixelsFile << "'" << endl;
268 return kTRUE;
269 }
270
271 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile << "'" << endl;
272
273 //
274 // Read the file and count the number of entries
275 //
276 UInt_t pixel = 0;
277
278 while (++fNumExcludedPixels)
279 {
280 in >> pixel;
281 if (!in)
282 break;
283
284 //
285 // Check for out of range
286 //
287 if (pixel > npixels)
288 {
289 *fLog << warn << "WARNING: To be excluded pixel: " << pixel << " is out of range " << endl;
290 continue;
291 }
292 //
293 // Exclude pixel
294 //
295 MCalibrationChargePix &pix = (*fCam)[pixel];
296 pix.SetExcluded();
297
298 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
299
300 }
301
302 if (--fNumExcludedPixels == 0)
303 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile << "'" << " is empty " << endl;
304 else
305 fCam->SetNumPixelsExcluded(fNumExcludedPixels);
306
307 return kTRUE;
308}
309
310
311// --------------------------------------------------------------------------
312//
313// Calculate the integral of the FADC time slices and store them as a new
314// pixel in the MCerPhotEvt container.
315//
316Int_t MCalibrationChargeCalc::Process()
317{
318 return kTRUE;
319}
320
321Int_t MCalibrationChargeCalc::PostProcess()
322{
323
324 //
325 // loop over the pedestal events and check if we have calibration
326 //
327 Int_t nvalid = 0;
328 Float_t avinnerped = 0;
329 Float_t avinnerprms = 0;
330 Float_t avinnernum = 0;
331 Float_t avouterped = 0;
332 Float_t avouterprms = 0;
333 Float_t avouternum = 0;
334
335 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
336 {
337
338 MCalibrationChargePix &pix = (*fCam)[pixid];
339
340 //
341 // get the pedestals
342 //
343 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
344 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
345 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
346
347 if (ped == -1.)
348 pix.SetExcluded();
349
350 //
351 // Check if the pixel has been excluded from the fits
352 //
353 if (pix.IsExcluded())
354 continue;
355
356 if (fGeom->GetPixRatio(pixid) == 1.)
357 {
358 avinnerped += ped;
359 avinnerprms += prms;
360 avinnernum++;
361 }
362 else
363 {
364 avouterped += ped;
365 avouterprms += prms;
366 avouternum++;
367 }
368 //
369 // set them in the calibration camera
370 //
371 if (pix.IsHiGainSaturation())
372 {
373 pix.SetPedestal(ped * fNumLoGainSamples,
374 prms * TMath::Sqrt((Float_t)fNumLoGainSamples),
375 prms * fNumLoGainSamples / num);
376 pix.SetNumLoGainSamples((Float_t)fNumLoGainSamples);
377 pix.ApplyLoGainConversion();
378 }
379 else
380 {
381 pix.SetPedestal(ped * fNumHiGainSamples,
382 prms * TMath::Sqrt((Float_t)fNumHiGainSamples),
383 prms * fNumHiGainSamples / num);
384 }
385
386 if (!pix.CheckChargeValidity() || !pix.CheckTimeValidity())
387 continue;
388
389 nvalid++;
390
391 if (!pix.CalcReducedSigma())
392 continue;
393
394 pix.CalcFFactorMethod();
395
396 }
397
398
399
400 //
401 // The Michele check ...
402 //
403 if (nvalid == 0)
404 {
405 *fLog << err << GetDescriptor() << ": Dear Michele! All pixels have non-valid calibration. "
406 << "Did you forget to fill the histograms (filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
407 return kFALSE;
408 }
409
410 MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
411 MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
412 //
413 // set the pedestans in the calibration camera
414 //
415 if (avinnerpix->IsHiGainSaturation())
416 {
417 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumLoGainSamples,
418 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumLoGainSamples),
419 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumLoGainSamples/avinnernum));
420 avinnerpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
421 avinnerpix->ApplyLoGainConversion();
422 }
423 else
424 {
425 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumHiGainSamples,
426 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumHiGainSamples),
427 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumHiGainSamples/avinnernum));
428 }
429
430 if (avouterpix->IsHiGainSaturation())
431 {
432 avouterpix->SetPedestal(avouterped/avouternum * fNumLoGainSamples,
433 avouterprms/avouternum * TMath::Sqrt((Float_t)fNumLoGainSamples),
434 avouterprms/avouternum * TMath::Sqrt((Float_t)fNumLoGainSamples/avouternum));
435 avouterpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
436 avouterpix->ApplyLoGainConversion();
437 }
438 else
439 {
440 avouterpix->SetPedestal(avouterped/avouternum * fNumHiGainSamples,
441 avouterprms/avouternum * TMath::Sqrt((Float_t)fNumHiGainSamples),
442 avouterprms/avouternum * TMath::Sqrt((Float_t)fNumHiGainSamples/avouternum));
443 }
444
445 if (!avinnerpix->CheckChargeValidity() || !avinnerpix->CheckTimeValidity())
446 if (!avinnerpix->CalcReducedSigma())
447 avinnerpix->CalcFFactorMethod();
448
449 if (!avouterpix->CheckChargeValidity() || !avouterpix->CheckTimeValidity())
450 if (!avouterpix->CalcReducedSigma())
451 avouterpix->CalcFFactorMethod();
452
453
454 if (!fBlindPixel->CheckChargeFitValidity())
455 {
456 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
457 fCam->SetBlindPixelMethodValid(kFALSE);
458 }
459 else
460 {
461 if (!fBlindPixel->CalcFluxInsidePlexiglass())
462 {
463 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
464 fCam->SetBlindPixelMethodValid(kFALSE);
465 }
466 else
467 {
468 fCam->SetBlindPixelMethodValid(kTRUE);
469 fCam->ApplyBlindPixelCalibration();
470 }
471 }
472
473 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
474 {
475 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
476 fCam->SetPINDiodeMethodValid(kFALSE);
477 }
478 else
479 {
480 if (!fPINDiode->CalcFluxOutsidePlexiglass())
481 {
482 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
483 fCam->SetPINDiodeMethodValid(kFALSE);
484 }
485 else
486 {
487 fCam->SetPINDiodeMethodValid(kTRUE);
488 fCam->ApplyPINDiodeCalibration();
489 }
490 }
491 fCam->SetReadyToSave();
492
493 return kTRUE;
494}
Note: See TracBrowser for help on using the repository browser.