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

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