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

Last change on this file since 3556 was 3556, checked in by gaug, 21 years ago
*** empty log message ***
File size: 20.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: Initialize MCalibrationCam
37// Initialize pulser light wavelength
38//
39// ReInit: MCalibrationCam::InitSize(NumPixels) is called from MGeomApply (which allocates
40// memory in a TClonesArray of type MCalibrationChargePix)
41// Initializes pointer to MBadPixelsCam
42//
43// Process: Nothing done by this class, histograms are filled by
44// MHCalibrationChargeCam
45//
46// PostProcess: Fit results from MHCalibrationChargeCam are retrieved
47// and used for the calculation of the reduced sigma,
48// the F-Factor method, the blind pixel method (photon flux
49// inside plexiglass) and
50// the PINDiode method (photon flux
51// outside plexiglass)
52//
53// Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
54// can be skipped with the command:
55// MalibrationCam::SkipHiLoGainCalibration()
56//
57// Input Containers:
58// MRawEvtData
59// MPedestalCam
60// MBadPixelsCam
61//
62// Output Containers:
63// MCalibrationCam
64// MCalibrationQECam
65// MBadPixelsCam
66//
67//
68// Preliminary description of the calibration in photons (email from 12/02/04)
69//
70// Why calibrating in photons:
71// ===========================
72//
73// At the Barcelona meeting in 2002, we decided to calibrate the camera in
74// photons. This for the following reasons:
75//
76// * The physical quantity arriving at the camera are photons. This is
77// the direct physical information from the air shower. The photons
78// have a flux and a spectrum.
79//
80// * The photon fluxes depend mostly on the shower energy (with
81// corrections deriving from the observation conditions), while the photon
82// spectra depend mostly on the observation conditions: zenith angle,
83// quality of the air, also the impact parameter of the shower.
84//
85// * The photomultiplier, in turn, has different response properties
86// (quantum efficiencies) for photons of different colour. (Moreover,
87// different pixels have slightly different quantum efficiencies).
88// The resulting number of photo-electrons is then amplified (linearly)
89// with respect to the photo-electron flux.
90//
91// * In the ideal case, one would like to disentagle the effects
92// of the observation conditions from the primary particle energy (which
93// one likes to measure). To do so, one needs:
94//
95// 1) A reliable calibration relating the FADC counts to the photo-electron
96// flux -> This is accomplished with the F-Factor method.
97//
98// 2) A reliable calibration of the wavelength-dependent quantum efficiency
99// -> This is accomplished with the combination of the three methods,
100// together with QE-measurements performed by David in order to do
101// the interpolation.
102//
103// 3) A reliable calibration of the observation conditions. This means:
104// - Tracing the atmospheric conditions -> LIDAR
105// - Tracing the observation zenith angle -> Drive System
106// 4) Some knowlegde about the impact parameter:
107// - This is the only part which cannot be accomplished well with a
108// single telescope. We would thus need to convolute the spectrum
109// over the distribution of impact parameters.
110//
111//
112// How an ideal calibration would look like:
113// =========================================
114//
115// We know from the combined PIN-Diode and Blind-Pixel Method the response of
116// each pixel to well-measured light fluxes in three representative
117// wavelengths (green, blue, UV). We also know the response to these light
118// fluxes in photo-electrons. Thus, we can derive:
119//
120// - conversion factors to photo-electrons
121// - conversion factors to photons in three wavelengths.
122//
123// Together with David's measurements and some MC-simulation, we should be
124// able to derive tables for typical Cherenkov-photon spectra - convoluted
125// with the impact parameters and depending on the athmospheric conditions
126// and the zenith angle (the "outer parameters").
127//
128// From these tables we can create "calibration tables" containing some
129// effective quantum efficiency depending on these outer parameters and which
130// are different for each pixel.
131//
132// In an ideal MCalibrate, one would thus have to convert first the FADC
133// slices to Photo-electrons and then, depending on the outer parameters,
134// look up the effective quantum efficiency and get the mean number of
135// photons which is then used for the further analysis.
136//
137// How the (first) MAGIC calibration should look like:
138// ===================================================
139//
140// For the moment, we have only one reliable calibration method, although
141// with very large systematic errors. This is the F-Factor method. Knowing
142// that the light is uniform over the whole camera (which I would not at all
143// guarantee in the case of the CT1 pulser), one could in principle already
144// perform a relative calibration of the quantum efficiencies in the UV.
145// However, the spread in QE at UV is about 10-15% (according to the plot
146// that Abelardo sent around last time. The spread in photo-electrons is 15%
147// for the inner pixels, but much larger (40%) for the outer ones.
148//
149// I'm not sure if we can already say that we have measured the relative
150// difference in quantum efficiency for the inner pixels and produce a first
151// QE-table for each pixel. To so, I would rather check in other wavelengths
152// (which we can do in about one-two weeks when the optical transmission of
153// the calibration trigger is installed).
154//
155// Thus, for the moment being, I would join Thomas proposal to calibrate in
156// photo-electrons and apply one stupid average quantum efficiency for all
157// pixels. This keeping in mind that we will have much preciser information
158// in about one to two weeks.
159//
160//
161// What MCalibrate should calculate and what should be stored:
162// ===========================================================
163//
164// It is clear that in the end, MCerPhotEvt will store photons.
165// MCalibrationCam stores the conversionfactors to photo-electrons and also
166// some tables of how to apply the conversion to photons, given the outer
167// parameters. This is not yet implemented and not even discussed.
168//
169// To start, I would suggest that we define the "average quantum efficiency"
170// (maybe something like 25+-3%) and apply them equally to all
171// photo-electrons. Later, this average factor can be easily replaced by a
172// pixel-dependent factor and later by a (pixel-dependent) table.
173//
174//
175//
176//////////////////////////////////////////////////////////////////////////////
177#include "MCalibrationChargeCalc.h"
178
179#include <TSystem.h>
180#include <TH1.h>
181
182#include "MLog.h"
183#include "MLogManip.h"
184
185#include "MParList.h"
186
187#include "MGeomCam.h"
188#include "MRawRunHeader.h"
189#include "MRawEvtPixelIter.h"
190
191#include "MPedestalCam.h"
192#include "MPedestalPix.h"
193
194#include "MCalibrationChargeCam.h"
195#include "MCalibrationChargePix.h"
196#include "MCalibrationChargePINDiode.h"
197#include "MCalibrationChargeBlindPix.h"
198
199#include "MExtractedSignalCam.h"
200#include "MExtractedSignalPix.h"
201
202#include "MBadPixelsCam.h"
203#include "MBadPixelsPix.h"
204
205#include "MCalibrationQECam.h"
206#include "MCalibrationQEPix.h"
207
208
209ClassImp(MCalibrationChargeCalc);
210
211using namespace std;
212
213// --------------------------------------------------------------------------
214//
215// Default constructor.
216//
217MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
218 : fPedestals(NULL), fCam(NULL), fQECam(NULL),
219 fRawEvt(NULL), fRunHeader(NULL), fGeom(NULL),
220 fBadPixels(NULL), fEvtTime(NULL),
221 fSignals(NULL), fPINDiode(NULL), fBlindPixel(NULL)
222{
223
224 fName = name ? name : "MCalibrationChargeCalc";
225 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
226
227 AddToBranchList("MRawEvtData.fHiGainPixId");
228 AddToBranchList("MRawEvtData.fLoGainPixId");
229 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
230 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
231
232 Clear();
233}
234
235void MCalibrationChargeCalc::Clear(const Option_t *o)
236{
237
238 SETBIT(fFlags, kUseQualityChecks);
239 SETBIT(fFlags, kHiLoGainCalibration);
240
241 fNumHiGainSamples = 0.;
242 fNumLoGainSamples = 0.;
243 fSqrtHiGainSamples = 0.;
244 fSqrtLoGainSamples = 0.;
245 fConversionHiLo = 0;
246 SkipQualityChecks ( kFALSE );
247 SkipHiLoGainCalibration( kFALSE );
248
249}
250
251
252// --------------------------------------------------------------------------
253//
254// The PreProcess searches for the following input containers:
255// - MRawEvtData
256// - MPedestalCam
257//
258// The following output containers are also searched and created if
259// they were not found:
260//
261// - MCalibrationCam
262// - MCalibrationQECam
263//
264// The following output containers are only searched, but not created
265//
266// - MTime
267//
268Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
269{
270
271 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
272 if (!fRawEvt)
273 {
274 *fLog << err << "MRawEvtData not found... aborting." << endl;
275 return kFALSE;
276 }
277
278 fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
279 if (!fCam)
280 return kFALSE;
281
282 fQECam = (MCalibrationQECam*)pList->FindCreateObj("MCalibrationQECam");
283 if (!fQECam)
284 return kFALSE;
285
286 fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
287 if (!fPINDiode)
288 return kFALSE;
289
290 fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
291 if (!fBlindPixel)
292 return kFALSE;
293
294 fEvtTime = (MTime*)pList->FindObject("MTime");
295
296 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
297 if (!fPedestals)
298 {
299 *fLog << err << "MPedestalCam not found... aborting" << endl;
300 return kFALSE;
301 }
302
303 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
304 if (!fSignals)
305 {
306 *fLog << err << "MExtractedSignalCam not found... aborting" << endl;
307 return kFALSE;
308 }
309
310 return kTRUE;
311}
312
313
314// --------------------------------------------------------------------------
315//
316// The ReInit searches for the following input containers:
317// - MRawRunHeader
318// - MGeomCam
319// - MBadPixelsCam
320//
321// It retrieves the following variables from MExtractedSignalCam:
322//
323// fNumHiGainSamples
324// fNumLoGainSamples
325//
326// fFirstUsedSliceHiGain
327// fLastUsedSliceHiGain
328// fFirstUsedSliceLoGain
329// fLastUsedSliceLoGain
330//
331// It defines the PixId of every pixel in MCalibrationChargeCam and MCalibrationQECam
332// It sets all pixels excluded which have the flag fBadBixelsPix::IsBad() set.
333//
334Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
335{
336
337 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
338 if (!fRunHeader)
339 {
340 *fLog << err << "MRawRunHeader not found... aborting." << endl;
341 return kFALSE;
342 }
343
344 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
345 if (!fGeom)
346 {
347 *fLog << err << "No MGeomCam found... aborting." << endl;
348 return kFALSE;
349 }
350
351 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
352 if (!fBadPixels)
353 {
354 *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
355 return kFALSE;
356 }
357
358 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
359 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
360 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
361 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
362
363 UInt_t npixels = fGeom->GetNumPixels();
364
365 for (UInt_t i=0; i<npixels; i++)
366 {
367 MCalibrationChargePix &pix = (*fCam) [i];
368 MCalibrationQEPix &pqe = (*fQECam) [i];
369 MBadPixelsPix &bad = (*fBadPixels)[i];
370
371 pix.SetPixId(i);
372 pqe.SetPixId(i);
373
374 if (bad.IsBad())
375 {
376 pix.SetExcluded();
377 pqe.SetExcluded();
378 continue;
379 }
380
381 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
382 fSignals->GetLastUsedSliceHiGain());
383 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
384 fSignals->GetLastUsedSliceLoGain());
385 }
386
387 return kTRUE;
388}
389
390
391Int_t MCalibrationChargeCalc::Process()
392{
393 return kTRUE;
394}
395
396// --------------------------------------------------------------------------
397//
398// Finalize pedestals:
399//
400// * Retrieve pedestal and pedestal RMS from MPedestalPix
401// * Retrieve total entries from MPedestalCam
402// * sum up pedestal and pedestalRMS for the average pixel
403// * set pedestal*number of used samples in MCalibrationChargePix
404// * set pedestal RMS * sqrt of number of used samples in MCalibrationChargePix
405//
406//
407void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal,
408 Float_t &avped, Float_t &avrms)
409{
410
411 //
412 // get the pedestals
413 //
414 const Float_t pedes = ped.GetPedestal();
415 const Float_t prms = ped.GetPedestalRms();
416 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
417
418 //
419 // Calculate the average pedestal
420 //
421 avped += pedes;
422 avrms += prms;
423
424 //
425 // set them in the calibration camera
426 //
427 if (cal.IsHiGainSaturation())
428 {
429 cal.SetPedestal(pedes* fNumLoGainSamples,
430 prms * fSqrtLoGainSamples,
431 prms * fNumLoGainSamples / num);
432 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples);
433 }
434 else
435 {
436 cal.SetPedestal(pedes* fNumHiGainSamples,
437 prms * fSqrtHiGainSamples,
438 prms * fNumHiGainSamples / num);
439 }
440
441}
442
443void MCalibrationChargeCalc::FinalizeAvPedestals(MCalibrationChargePix &cal,
444 Float_t avped, Float_t avrms, Int_t avnum)
445{
446
447 //
448 // set the pedestans in the calibration camera
449 //
450 if (cal.IsHiGainSaturation())
451 {
452 cal.SetPedestal(avped/avnum * fNumLoGainSamples,
453 avrms/avnum * fSqrtLoGainSamples,
454 avrms/avnum * fSqrtLoGainSamples/avnum);
455 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples);
456 }
457 else
458 {
459 cal.SetPedestal(avped/avnum * fNumHiGainSamples,
460 avrms/avnum * fSqrtHiGainSamples,
461 avrms/avnum * fSqrtHiGainSamples/avnum);
462 }
463}
464
465//
466// Finalize charges per pixel:
467// * Check chage validity
468// * check absolute time validity
469// * calculate the reduced sigma
470// * caluclate the number of photo-electrons
471//
472//
473Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad)
474{
475
476 cal.CheckChargeValidity (&bad);
477 cal.CheckTimeValidity (&bad);
478
479 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
480 return kFALSE;
481
482 if (!cal.CalcReducedSigma())
483 {
484 *fLog << warn << GetDescriptor()
485 << ": Could not calculate reduced sigmas of pixel: " << cal.GetPixId() << endl;
486 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
487 return kFALSE;
488 }
489
490 if (!cal.CalcFFactorMethod())
491 {
492 *fLog << warn << GetDescriptor()
493 << ": Could not calculate F-Factor of pixel: " << cal.GetPixId() << endl;
494 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
495 return kFALSE;
496 }
497 return kTRUE;
498}
499
500//
501// * Finalize the pedestals
502// * Do the quality checks
503// * Calculate the reduced sigma
504// * Calculate the F-Factor Method
505//
506Int_t MCalibrationChargeCalc::PostProcess()
507{
508
509 //
510 // loop over the pedestal events and check if we have calibration
511 //
512 Int_t nvalid = 0;
513 Float_t avinnerped = 0.;
514 Float_t avinnerprms = 0.;
515 Int_t avinnernum = 0;
516 Float_t avouterped = 0.;
517 Float_t avouterprms = 0.;
518 Int_t avouternum = 0;
519
520 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
521 {
522
523 MCalibrationChargePix &pix = (*fCam)[pixid];
524 //
525 // Check if the pixel has been excluded from the fits
526 //
527 if (pix.IsExcluded())
528 continue;
529
530 MPedestalPix &ped = (*fPedestals)[pixid];
531 MBadPixelsPix &bad = (*fBadPixels)[pixid];
532
533 if (fGeom->GetPixRatio(pixid) == 1.)
534 {
535 FinalizePedestals(ped,pix,avinnerped,avinnerprms);
536 avinnernum++;
537 }
538 else
539 {
540 FinalizePedestals(ped,pix,avouterped,avouterprms);
541 avouternum++;
542 }
543
544 if (FinalizeCharges(pix,bad))
545 nvalid++;
546 }
547
548 //
549 // The Michele check ...
550 //
551 if (nvalid == 0)
552 {
553 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
554 << "Did you forget to fill the histograms "
555 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
556 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
557 << "instead of a calibration run " << endl;
558 return kFALSE;
559 }
560
561 FinalizeAvPedestals(*fCam->GetAverageInnerPix(), avinnerped, avinnerprms,avinnernum);
562 FinalizeAvPedestals(*fCam->GetAverageOuterPix(), avouterped, avouterprms,avouternum);
563
564 FinalizeCharges(*fCam->GetAverageInnerPix(),*fCam->GetAverageInnerBadPix());
565 FinalizeCharges(*fCam->GetAverageOuterPix(),*fCam->GetAverageOuterBadPix());
566
567 //
568 // F-Factor calibration
569 //
570 if (fCam->CalcMeanFluxPhotonsFFactorMethod(*fGeom, *fBadPixels))
571 {
572 fCam->ApplyFFactorCalibration(*fGeom,*fBadPixels);
573 fCam->SetFFactorMethodValid(kTRUE);
574 }
575 else
576 {
577 *fLog << warn << "Could not calculate the flux of photo-electrons from the F-Factor method, " << endl;
578 fCam->SetFFactorMethodValid(kFALSE);
579 }
580
581 //
582 // Blind Pixel calibration
583 //
584 if (!fBlindPixel->CheckChargeFitValidity())
585 {
586 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
587 << "charge fit not valid " << endl;
588 fCam->SetBlindPixelMethodValid(kFALSE);
589 }
590 else
591 {
592 if (!fBlindPixel->CalcFluxInsidePlexiglass())
593 {
594 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
595 << "will skip PIN Diode Calibration " << endl;
596 fCam->SetBlindPixelMethodValid(kFALSE);
597 }
598 else
599 {
600 fCam->SetBlindPixelMethodValid(kTRUE);
601 fCam->ApplyBlindPixelCalibration(*fGeom,*fBadPixels, *fBlindPixel);
602 }
603 }
604
605 //
606 // PIN Diode calibration
607 //
608 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
609 {
610 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
611 << "charge fit not valid " << endl;
612 fCam->SetPINDiodeMethodValid(kFALSE);
613 }
614 else
615 {
616 if (!fPINDiode->CalcFluxOutsidePlexiglass())
617 {
618 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
619 << "will skip PIN Diode Calibration " << endl;
620 fCam->SetPINDiodeMethodValid(kFALSE);
621 }
622 else
623 {
624 fCam->SetPINDiodeMethodValid(kTRUE);
625 fCam->ApplyPINDiodeCalibration(*fGeom,*fBadPixels, *fPINDiode);
626 }
627 }
628 fCam->SetReadyToSave();
629
630return kTRUE;
631}
Note: See TracBrowser for help on using the repository browser.