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

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