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

Last change on this file since 3555 was 3555, checked in by gaug, 21 years ago
*** empty log message ***
File size: 20.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: 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(MPedestalPix &ped, MCalibrationChargePix &cal,
408 Float_t &avped, Float_t &avrms, Float_t &avnum)
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 avnum++;
424
425 //
426 // set them in the calibration camera
427 //
428 if (cal.IsHiGainSaturation())
429 {
430 cal.SetPedestal(pedes* fNumLoGainSamples,
431 prms * fSqrtLoGainSamples,
432 prms * fNumLoGainSamples / num);
433 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples);
434 }
435 else
436 {
437 cal.SetPedestal(pedes* fNumHiGainSamples,
438 prms * fSqrtHiGainSamples,
439 prms * fNumHiGainSamples / num);
440 }
441
442}
443
444
445
446Int_t MCalibrationChargeCalc::PostProcess()
447{
448
449 //
450 // loop over the pedestal events and check if we have calibration
451 //
452 Int_t nvalid = 0;
453 Float_t avinnerped = 0;
454 Float_t avinnerprms = 0;
455 Float_t avinnernum = 0;
456 Float_t avouterped = 0;
457 Float_t avouterprms = 0;
458 Float_t avouternum = 0;
459
460 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
461 {
462
463 MCalibrationChargePix &pix = (*fCam)[pixid];
464 //
465 // Check if the pixel has been excluded from the fits
466 //
467 if (pix.IsExcluded())
468 continue;
469
470 MPedestalPix &ped = (*fPedestals)[pixid];
471
472
473 if (fGeom->GetPixRatio(pixid) == 1.)
474 FinalizePedestals(ped,pix,avinnerped,avinnerprms,avinnernum);
475 else
476 FinalizePedestals(ped,pix,avouterped,avouterprms,avouternum);
477
478
479 MBadPixelsPix &bad = (*fBadPixels)[pixid];
480
481 pix.CheckChargeValidity (&bad);
482 pix.CheckTimeValidity (&bad);
483
484 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
485 continue;
486
487 nvalid++;
488
489 if (!pix.CalcReducedSigma())
490 {
491 *fLog << warn << GetDescriptor()
492 << ": Could not calculate reduced sigmas of pixel: " << pix.GetPixId() << endl;
493 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
494 continue;
495 }
496
497 if (!pix.CalcFFactorMethod())
498 {
499 *fLog << warn << GetDescriptor()
500 << ": Could not calculate F-Factor of pixel: " << pix.GetPixId() << endl;
501 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
502 }
503 }
504
505 //
506 // The Michele check ...
507 //
508 if (nvalid == 0)
509 {
510 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
511 << "Did you forget to fill the histograms "
512 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
513 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
514 << "instead of a calibration run " << endl;
515 return kFALSE;
516 }
517
518 MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
519 MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
520 //
521 // set the pedestans in the calibration camera
522 //
523 if (avinnerpix->IsHiGainSaturation())
524 {
525 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumLoGainSamples,
526 avinnerprms/avinnernum * fSqrtLoGainSamples,
527 avinnerprms/avinnernum * fSqrtLoGainSamples/avinnernum);
528 avinnerpix->CalcLoGainPedestal((Float_t)fNumLoGainSamples);
529 }
530 else
531 {
532 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumHiGainSamples,
533 avinnerprms/avinnernum * fSqrtHiGainSamples,
534 avinnerprms/avinnernum * fSqrtHiGainSamples/avinnernum);
535 }
536
537 if (avouterpix->IsHiGainSaturation())
538 {
539 avouterpix->SetPedestal(avouterped/avouternum * fNumLoGainSamples,
540 avouterprms/avouternum * fSqrtLoGainSamples,
541 avouterprms/avouternum * fSqrtLoGainSamples/avouternum);
542 avouterpix->CalcLoGainPedestal((Float_t)fNumLoGainSamples);
543 }
544 else
545 {
546 avouterpix->SetPedestal(avouterped/avouternum * fNumHiGainSamples,
547 avouterprms/avouternum * fSqrtHiGainSamples,
548 avouterprms/avouternum * fSqrtHiGainSamples/avouternum);
549 }
550
551 MBadPixelsPix *bad = fCam->GetAverageInnerBadPix();
552
553 avinnerpix->CheckChargeValidity(bad);
554 avinnerpix->CheckTimeValidity(bad);
555
556 if (bad->IsCalibrationSignalOK())
557 if (!avinnerpix->CalcReducedSigma())
558 avinnerpix->CalcFFactorMethod();
559
560 bad = fCam->GetAverageInnerBadPix();
561
562 avouterpix->CheckChargeValidity(bad);
563 avouterpix->CheckTimeValidity(bad);
564
565 if (bad->IsCalibrationSignalOK())
566 if (!avouterpix->CalcReducedSigma())
567 avouterpix->CalcFFactorMethod();
568
569 //
570 // F-Factor calibration
571 //
572 if (fCam->CalcMeanFluxPhotonsFFactorMethod(*fGeom, *fBadPixels))
573 {
574 fCam->ApplyFFactorCalibration(*fGeom,*fBadPixels);
575 fCam->SetFFactorMethodValid(kTRUE);
576 }
577 else
578 {
579 *fLog << warn << "Could not calculate the flux of photo-electrons from the F-Factor method, " << endl;
580 fCam->SetFFactorMethodValid(kFALSE);
581 }
582
583 //
584 // Blind Pixel calibration
585 //
586 if (!fBlindPixel->CheckChargeFitValidity())
587 {
588 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
589 << "charge fit not valid " << endl;
590 fCam->SetBlindPixelMethodValid(kFALSE);
591 }
592 else
593 {
594 if (!fBlindPixel->CalcFluxInsidePlexiglass())
595 {
596 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
597 << "will skip PIN Diode Calibration " << endl;
598 fCam->SetBlindPixelMethodValid(kFALSE);
599 }
600 else
601 {
602 fCam->SetBlindPixelMethodValid(kTRUE);
603 fCam->ApplyBlindPixelCalibration(*fGeom,*fBadPixels, *fBlindPixel);
604 }
605 }
606
607 //
608 // PIN Diode calibration
609 //
610 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
611 {
612 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
613 << "charge fit not valid " << endl;
614 fCam->SetPINDiodeMethodValid(kFALSE);
615 }
616 else
617 {
618 if (!fPINDiode->CalcFluxOutsidePlexiglass())
619 {
620 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
621 << "will skip PIN Diode Calibration " << endl;
622 fCam->SetPINDiodeMethodValid(kFALSE);
623 }
624 else
625 {
626 fCam->SetPINDiodeMethodValid(kTRUE);
627 fCam->ApplyPINDiodeCalibration(*fGeom,*fBadPixels, *fPINDiode);
628 }
629 }
630 fCam->SetReadyToSave();
631
632return kTRUE;
633}
Note: See TracBrowser for help on using the repository browser.