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

Last change on this file since 3551 was 3511, checked in by gaug, 21 years ago
*** empty log message ***
File size: 19.1 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 fEvtTime = (MTime*)pList->FindObject("MTime");
282
283 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
284 if (!fPedestals)
285 {
286 *fLog << err << "MPedestalCam not found... aborting" << endl;
287 return kFALSE;
288 }
289
290 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
291 if (!fSignals)
292 {
293 *fLog << err << "MExtractedSignalCam not found... aborting" << endl;
294 return kFALSE;
295 }
296
297 return kTRUE;
298}
299
300
301// --------------------------------------------------------------------------
302//
303// The ReInit searches for the following input containers:
304// - MRawRunHeader
305//
306Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
307{
308
309 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
310 if (!fRunHeader)
311 {
312 *fLog << err << "MRawRunHeader not found... aborting." << endl;
313 return kFALSE;
314 }
315
316 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
317 if (!fGeom)
318 {
319 *fLog << err << "No MGeomCam found... aborting." << endl;
320 return kFALSE;
321 }
322
323 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
324 if (!fBadPixels)
325 {
326 *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
327 return kFALSE;
328 }
329
330 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
331 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
332 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
333 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
334
335 UInt_t npixels = fGeom->GetNumPixels();
336
337 for (UInt_t i=0; i<npixels; i++)
338 {
339 MCalibrationChargePix &pix = (*fCam)[i];
340 MBadPixelsPix &bad = (*fBadPixels)[i];
341
342 pix.DefinePixId(i);
343
344 if (bad.IsBad())
345 {
346 pix.SetExcluded();
347 continue;
348 }
349
350 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
351 fSignals->GetLastUsedSliceHiGain());
352 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
353 fSignals->GetLastUsedSliceLoGain());
354 }
355
356 return kTRUE;
357}
358
359
360// --------------------------------------------------------------------------
361//
362// Calculate the integral of the FADC time slices and store them as a new
363// pixel in the MCerPhotEvt container.
364//
365Int_t MCalibrationChargeCalc::Process()
366{
367 return kTRUE;
368}
369
370Int_t MCalibrationChargeCalc::PostProcess()
371{
372
373 //
374 // loop over the pedestal events and check if we have calibration
375 //
376 Int_t nvalid = 0;
377 Float_t avinnerped = 0;
378 Float_t avinnerprms = 0;
379 Float_t avinnernum = 0;
380 Float_t avouterped = 0;
381 Float_t avouterprms = 0;
382 Float_t avouternum = 0;
383
384 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
385 {
386
387 MCalibrationChargePix &pix = (*fCam)[pixid];
388 MBadPixelsPix &bad = (*fBadPixels)[pixid];
389
390 //
391 // Check if the pixel has been excluded from the fits
392 //
393 if (pix.IsExcluded())
394 continue;
395
396 //
397 // get the pedestals
398 //
399 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
400 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
401 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
402
403 if (fGeom->GetPixRatio(pixid) == 1.)
404 {
405 avinnerped += ped;
406 avinnerprms += prms;
407 avinnernum++;
408 }
409 else
410 {
411 avouterped += ped;
412 avouterprms += prms;
413 avouternum++;
414 }
415 //
416 // set them in the calibration camera
417 //
418 if (pix.IsHiGainSaturation())
419 {
420 pix.SetPedestal(ped * fNumLoGainSamples,
421 prms * fSqrtLoGainSamples,
422 prms * fNumLoGainSamples / num);
423 pix.SetNumLoGainSamples(fNumLoGainSamples);
424 pix.ApplyLoGainConversion();
425 }
426 else
427 {
428 pix.SetPedestal(ped * fNumHiGainSamples,
429 prms * fSqrtHiGainSamples,
430 prms * fNumHiGainSamples / num);
431 }
432
433 pix.CheckChargeValidity (&bad);
434 pix.CheckTimeValidity (&bad);
435
436 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
437 continue;
438
439 nvalid++;
440
441 if (!pix.CalcReducedSigma())
442 {
443 *fLog << warn << GetDescriptor()
444 << ": Could not calculate reduced sigmas of pixel: " << pix.GetPixId() << endl;
445 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
446 continue;
447 }
448
449 if (!pix.CalcFFactorMethod())
450 {
451 *fLog << warn << GetDescriptor()
452 << ": Could not calculate F-Factor of pixel: " << pix.GetPixId() << endl;
453 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
454 }
455 }
456
457 //
458 // The Michele check ...
459 //
460 if (nvalid == 0)
461 {
462 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
463 << "Did you forget to fill the histograms "
464 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
465 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
466 << "instead of a calibration run " << endl;
467 return kFALSE;
468 }
469
470 MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
471 MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
472 //
473 // set the pedestans in the calibration camera
474 //
475 if (avinnerpix->IsHiGainSaturation())
476 {
477 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumLoGainSamples,
478 avinnerprms/avinnernum * fSqrtLoGainSamples,
479 avinnerprms/avinnernum * fSqrtLoGainSamples/avinnernum);
480 avinnerpix->SetNumLoGainSamples(fNumLoGainSamples);
481 avinnerpix->ApplyLoGainConversion();
482 }
483 else
484 {
485 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumHiGainSamples,
486 avinnerprms/avinnernum * fSqrtHiGainSamples,
487 avinnerprms/avinnernum * fSqrtHiGainSamples/avinnernum);
488 }
489
490 if (avouterpix->IsHiGainSaturation())
491 {
492 avouterpix->SetPedestal(avouterped/avouternum * fNumLoGainSamples,
493 avouterprms/avouternum * fSqrtLoGainSamples,
494 avouterprms/avouternum * fSqrtLoGainSamples/avouternum);
495 avouterpix->SetNumLoGainSamples(fNumLoGainSamples);
496 avouterpix->ApplyLoGainConversion();
497 }
498 else
499 {
500 avouterpix->SetPedestal(avouterped/avouternum * fNumHiGainSamples,
501 avouterprms/avouternum * fSqrtHiGainSamples,
502 avouterprms/avouternum * fSqrtHiGainSamples/avouternum);
503 }
504
505 MBadPixelsPix *bad = fCam->GetAverageInnerBadPix();
506
507 avinnerpix->CheckChargeValidity(bad);
508 avinnerpix->CheckTimeValidity(bad);
509
510 if (bad->IsCalibrationSignalOK())
511 if (!avinnerpix->CalcReducedSigma())
512 avinnerpix->CalcFFactorMethod();
513
514 bad = fCam->GetAverageInnerBadPix();
515
516 avouterpix->CheckChargeValidity(bad);
517 avouterpix->CheckTimeValidity(bad);
518
519 if (bad->IsCalibrationSignalOK())
520 if (!avouterpix->CalcReducedSigma())
521 avouterpix->CalcFFactorMethod();
522
523 //
524 // F-Factor calibration
525 //
526 if (fCam->CalcMeanFluxPhotonsFFactorMethod(*fGeom, *fBadPixels))
527 {
528 fCam->ApplyFFactorCalibration(*fGeom,*fBadPixels);
529 fCam->SetFFactorMethodValid(kTRUE);
530 }
531 else
532 {
533 *fLog << warn << "Could not calculate the flux of photo-electrons from the F-Factor method, " << endl;
534 fCam->SetFFactorMethodValid(kFALSE);
535 }
536
537 //
538 // Blind Pixel calibration
539 //
540 if (!fBlindPixel->CheckChargeFitValidity())
541 {
542 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
543 << "charge fit not valid " << endl;
544 fCam->SetBlindPixelMethodValid(kFALSE);
545 }
546 else
547 {
548 if (!fBlindPixel->CalcFluxInsidePlexiglass())
549 {
550 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
551 << "will skip PIN Diode Calibration " << endl;
552 fCam->SetBlindPixelMethodValid(kFALSE);
553 }
554 else
555 {
556 fCam->SetBlindPixelMethodValid(kTRUE);
557 fCam->ApplyBlindPixelCalibration(*fGeom,*fBadPixels, *fBlindPixel);
558 }
559 }
560
561 //
562 // PIN Diode calibration
563 //
564 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
565 {
566 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
567 << "charge fit not valid " << endl;
568 fCam->SetPINDiodeMethodValid(kFALSE);
569 }
570 else
571 {
572 if (!fPINDiode->CalcFluxOutsidePlexiglass())
573 {
574 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
575 << "will skip PIN Diode Calibration " << endl;
576 fCam->SetPINDiodeMethodValid(kFALSE);
577 }
578 else
579 {
580 fCam->SetPINDiodeMethodValid(kTRUE);
581 fCam->ApplyPINDiodeCalibration(*fGeom,*fBadPixels, *fPINDiode);
582 }
583 }
584 fCam->SetReadyToSave();
585
586return kTRUE;
587}
Note: See TracBrowser for help on using the repository browser.