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

Last change on this file since 3451 was 3451, checked in by gaug, 21 years ago
*** empty log message ***
File size: 18.3 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(fNumHiGainSamples);
339 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
340
341 UInt_t npixels = fGeom->GetNumPixels();
342
343 for (UInt_t i=0; i<npixels; i++)
344 {
345 MCalibrationChargePix &pix = (*fCam)[i];
346 MBadPixelsPix &bad = (*fBadPixels)[i];
347
348 pix.DefinePixId(i);
349
350 if (bad.IsBad())
351 {
352 pix.SetExcluded();
353 continue;
354 }
355
356 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
357 fSignals->GetLastUsedSliceHiGain());
358 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
359 fSignals->GetLastUsedSliceLoGain());
360 }
361
362 return kTRUE;
363}
364
365
366// --------------------------------------------------------------------------
367//
368// Calculate the integral of the FADC time slices and store them as a new
369// pixel in the MCerPhotEvt container.
370//
371Int_t MCalibrationChargeCalc::Process()
372{
373 return kTRUE;
374}
375
376Int_t MCalibrationChargeCalc::PostProcess()
377{
378
379 //
380 // loop over the pedestal events and check if we have calibration
381 //
382 Int_t nvalid = 0;
383 Float_t avinnerped = 0;
384 Float_t avinnerprms = 0;
385 Float_t avinnernum = 0;
386 Float_t avouterped = 0;
387 Float_t avouterprms = 0;
388 Float_t avouternum = 0;
389
390 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
391 {
392
393 MCalibrationChargePix &pix = (*fCam)[pixid];
394 MBadPixelsPix &bad = (*fBadPixels)[pixid];
395
396 //
397 // Check if the pixel has been excluded from the fits
398 //
399 if (pix.IsExcluded())
400 continue;
401
402 //
403 // get the pedestals
404 //
405 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
406 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
407 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
408
409 if (fGeom->GetPixRatio(pixid) == 1.)
410 {
411 avinnerped += ped;
412 avinnerprms += prms;
413 avinnernum++;
414 }
415 else
416 {
417 avouterped += ped;
418 avouterprms += prms;
419 avouternum++;
420 }
421 //
422 // set them in the calibration camera
423 //
424 if (pix.IsHiGainSaturation())
425 {
426 pix.SetPedestal(ped * fNumLoGainSamples,
427 prms * fSqrtLoGainSamples,
428 prms * fNumLoGainSamples / num);
429 pix.SetNumLoGainSamples(fNumLoGainSamples);
430 pix.ApplyLoGainConversion();
431 }
432 else
433 {
434 pix.SetPedestal(ped * fNumHiGainSamples,
435 prms * fSqrtHiGainSamples,
436 prms * fNumHiGainSamples / num);
437 }
438
439 pix.CheckChargeValidity(&bad);
440 pix.CheckTimeValidity(&bad);
441
442 if (bad.IsUnsuitableRun())
443 continue;
444
445 nvalid++;
446
447 if (!pix.CalcReducedSigma())
448 {
449 bad.SetUnsuitableRun();
450 continue;
451 }
452
453 pix.CalcFFactorMethod();
454 }
455
456
457
458 //
459 // The Michele check ...
460 //
461 if (nvalid == 0)
462 {
463 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
464 << "Did you forget to fill the histograms "
465 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
466 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
467 << "instead of a calibration run " << endl;
468 return kFALSE;
469 }
470
471 MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
472 MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
473 //
474 // set the pedestans in the calibration camera
475 //
476 if (avinnerpix->IsHiGainSaturation())
477 {
478 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumLoGainSamples,
479 avinnerprms/avinnernum * fSqrtLoGainSamples,
480 avinnerprms/avinnernum * fSqrtLoGainSamples/avinnernum);
481 avinnerpix->SetNumLoGainSamples(fNumLoGainSamples);
482 avinnerpix->ApplyLoGainConversion();
483 }
484 else
485 {
486 avinnerpix->SetPedestal(avinnerped/avinnernum * fNumHiGainSamples,
487 avinnerprms/avinnernum * fSqrtHiGainSamples,
488 avinnerprms/avinnernum * fSqrtHiGainSamples/avinnernum);
489 }
490
491 if (avouterpix->IsHiGainSaturation())
492 {
493 avouterpix->SetPedestal(avouterped/avouternum * fNumLoGainSamples,
494 avouterprms/avouternum * fSqrtLoGainSamples,
495 avouterprms/avouternum * fSqrtLoGainSamples/avouternum);
496 avouterpix->SetNumLoGainSamples(fNumLoGainSamples);
497 avouterpix->ApplyLoGainConversion();
498 }
499 else
500 {
501 avouterpix->SetPedestal(avouterped/avouternum * fNumHiGainSamples,
502 avouterprms/avouternum * fSqrtHiGainSamples,
503 avouterprms/avouternum * fSqrtHiGainSamples/avouternum);
504 }
505
506 MBadPixelsPix *bad = fCam->GetAverageInnerBadPix();
507
508 avinnerpix->CheckChargeValidity(bad);
509 avinnerpix->CheckTimeValidity(bad);
510
511 if (bad->IsCalibrationSignalOK())
512 if (!avinnerpix->CalcReducedSigma())
513 avinnerpix->CalcFFactorMethod();
514
515 bad = fCam->GetAverageInnerBadPix();
516
517 avouterpix->CheckChargeValidity(bad);
518 avouterpix->CheckTimeValidity(bad);
519
520 if (bad->IsCalibrationSignalOK())
521 if (!avouterpix->CalcReducedSigma())
522 avouterpix->CalcFFactorMethod();
523
524 if (!fBlindPixel->CheckChargeFitValidity())
525 {
526 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
527 << "charge fit not valid " << endl;
528 fCam->SetBlindPixelMethodValid(kFALSE);
529 }
530 else
531 {
532 if (!fBlindPixel->CalcFluxInsidePlexiglass())
533 {
534 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
535 << "will skip PIN Diode Calibration " << endl;
536 fCam->SetBlindPixelMethodValid(kFALSE);
537 }
538 else
539 {
540 fCam->SetBlindPixelMethodValid(kTRUE);
541 fCam->ApplyBlindPixelCalibration();
542 }
543 }
544
545 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
546 {
547 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
548 << "charge fit not valid " << endl;
549 fCam->SetPINDiodeMethodValid(kFALSE);
550 }
551 else
552 {
553 if (!fPINDiode->CalcFluxOutsidePlexiglass())
554 {
555 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
556 << "will skip PIN Diode Calibration " << endl;
557 fCam->SetPINDiodeMethodValid(kFALSE);
558 }
559 else
560 {
561 fCam->SetPINDiodeMethodValid(kTRUE);
562 fCam->ApplyPINDiodeCalibration();
563 }
564 }
565 fCam->SetReadyToSave();
566
567return kTRUE;
568}
Note: See TracBrowser for help on using the repository browser.