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

Last change on this file since 3492 was 3479, checked in by gaug, 21 years ago
*** empty log message ***
File size: 18.9 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 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
334 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
335 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
336 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
337
338 UInt_t npixels = fGeom->GetNumPixels();
339
340 for (UInt_t i=0; i<npixels; i++)
341 {
342 MCalibrationChargePix &pix = (*fCam)[i];
343 MBadPixelsPix &bad = (*fBadPixels)[i];
344
345 pix.DefinePixId(i);
346
347 if (bad.IsBad())
348 {
349 pix.SetExcluded();
350 continue;
351 }
352
353 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
354 fSignals->GetLastUsedSliceHiGain());
355 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
356 fSignals->GetLastUsedSliceLoGain());
357 }
358
359 return kTRUE;
360}
361
362
363// --------------------------------------------------------------------------
364//
365// Calculate the integral of the FADC time slices and store them as a new
366// pixel in the MCerPhotEvt container.
367//
368Int_t MCalibrationChargeCalc::Process()
369{
370 return kTRUE;
371}
372
373Int_t MCalibrationChargeCalc::PostProcess()
374{
375
376 //
377 // loop over the pedestal events and check if we have calibration
378 //
379 Int_t nvalid = 0;
380 Float_t avinnerped = 0;
381 Float_t avinnerprms = 0;
382 Float_t avinnernum = 0;
383 Float_t avouterped = 0;
384 Float_t avouterprms = 0;
385 Float_t avouternum = 0;
386
387 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
388 {
389
390 MCalibrationChargePix &pix = (*fCam)[pixid];
391 MBadPixelsPix &bad = (*fBadPixels)[pixid];
392
393 //
394 // Check if the pixel has been excluded from the fits
395 //
396 if (pix.IsExcluded())
397 continue;
398
399 //
400 // get the pedestals
401 //
402 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
403 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
404 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
405
406 if (fGeom->GetPixRatio(pixid) == 1.)
407 {
408 avinnerped += ped;
409 avinnerprms += prms;
410 avinnernum++;
411 }
412 else
413 {
414 avouterped += ped;
415 avouterprms += prms;
416 avouternum++;
417 }
418 //
419 // set them in the calibration camera
420 //
421 if (pix.IsHiGainSaturation())
422 {
423 pix.SetPedestal(ped * fNumLoGainSamples,
424 prms * fSqrtLoGainSamples,
425 prms * fNumLoGainSamples / num);
426 pix.SetNumLoGainSamples(fNumLoGainSamples);
427 pix.ApplyLoGainConversion();
428 }
429 else
430 {
431 pix.SetPedestal(ped * fNumHiGainSamples,
432 prms * fSqrtHiGainSamples,
433 prms * fNumHiGainSamples / num);
434 }
435
436 pix.CheckChargeValidity (&bad);
437 pix.CheckTimeValidity (&bad);
438
439 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
440 continue;
441
442 nvalid++;
443
444 if (!pix.CalcReducedSigma())
445 {
446 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
447 continue;
448 }
449
450 if (!pix.CalcFFactorMethod())
451 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
452
453 }
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);
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);
582 }
583 }
584 fCam->SetReadyToSave();
585
586return kTRUE;
587}
Note: See TracBrowser for help on using the repository browser.