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

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