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

Last change on this file since 5682 was 5682, checked in by gaug, 20 years ago
*** empty log message ***
File size: 74.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// MCalibrationChargeCalc
27//
28// Task to calculate the calibration conversion factors and quantum efficiencies
29// from the fit results to the summed FADC slice distributions delivered by
30// MCalibrationChargeCam, MCalibrationChargePix, MCalibrationChargeBlindPix and
31// MCalibrationChargePINDiode, calculated and filled by MHCalibrationChargeCam,
32// MHCalibrationChargePix, MHCalibrationChargeBlindPix and MHCalibrationChargePINDiode.
33//
34// PreProcess(): Initialize pointers to MCalibrationChargeCam, MCalibrationChargeBlindPix
35// MCalibrationChargePINDiode and MCalibrationQECam
36//
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 to be done, histograms getting filled by MHCalibrationChargeCam
44//
45// PostProcess(): - FinalizePedestals()
46// - FinalizeCharges()
47// - FinalizeFFactorMethod()
48// - FinalizeBadPixels()
49// - FinalizeBlindCam()
50// - FinalizePINDiode()
51// - FinalizeFFactorQECam()
52// - FinalizeBlindPixelQECam()
53// - FinalizePINDiodeQECam()
54//
55// Input Containers:
56// MCalibrationChargeCam
57// MCalibrationChargeBlindPix
58// MCalibrationChargePINDiode
59// MCalibrationQECam
60// MPedestalCam
61// MBadPixelsCam
62// MGeomCam
63// MTime
64//
65// Output Containers:
66// MCalibrationChargeCam
67// MCalibrationChargeBlindPix
68// MCalibrationChargePINDiode
69// MCalibrationQECam
70// MBadPixelsCam
71//
72//
73// Preliminary description of the calibration in photons (email from 12/02/04)
74//
75// Why calibrating in photons:
76// ===========================
77//
78// At the Barcelona meeting in 2002, we decided to calibrate the camera in
79// photons. This for the following reasons:
80//
81// * The physical quantity arriving at the camera are photons. This is
82// the direct physical information from the air shower. The photons
83// have a flux and a spectrum.
84//
85// * The photon fluxes depend mostly on the shower energy (with
86// corrections deriving from the observation conditions), while the photon
87// spectra depend mostly on the observation conditions: zenith angle,
88// quality of the air, also the impact parameter of the shower.
89//
90// * The photomultiplier, in turn, has different response properties
91// (quantum efficiencies) for photons of different colour. (Moreover,
92// different pixels have slightly different quantum efficiencies).
93// The resulting number of photo-electrons is then amplified (linearly)
94// with respect to the photo-electron flux.
95//
96// * In the ideal case, one would like to disentagle the effects
97// of the observation conditions from the primary particle energy (which
98// one likes to measure). To do so, one needs:
99//
100// 1) A reliable calibration relating the FADC counts to the photo-electron
101// flux -> This is accomplished with the F-Factor method.
102//
103// 2) A reliable calibration of the wavelength-dependent quantum efficiency
104// -> This is accomplished with the combination of the three methods,
105// together with QE-measurements performed by David in order to do
106// the interpolation.
107//
108// 3) A reliable calibration of the observation conditions. This means:
109// - Tracing the atmospheric conditions -> LIDAR
110// - Tracing the observation zenith angle -> Drive System
111//
112// 4) Some knowlegde about the impact parameter:
113// - This is the only part which cannot be accomplished well with a
114// single telescope. We would thus need to convolute the spectrum
115// over the distribution of impact parameters.
116//
117//
118// How an ideal calibration would look like:
119// =========================================
120//
121// We know from the combined PIN-Diode and Blind-Pixel Method the response of
122// each pixel to well-measured light fluxes in three representative
123// wavelengths (green, blue, UV). We also know the response to these light
124// fluxes in photo-electrons. Thus, we can derive:
125//
126// - conversion factors to photo-electrons
127// - conversion factors to photons in three wavelengths.
128//
129// Together with David's measurements and some MC-simulation, we should be
130// able to derive tables for typical Cherenkov-photon spectra - convoluted
131// with the impact parameters and depending on the athmospheric conditions
132// and the zenith angle (the "outer parameters").
133//
134// From these tables we can create "calibration tables" containing some
135// effective quantum efficiency depending on these outer parameters and which
136// are different for each pixel.
137//
138// In an ideal MCalibrate, one would thus have to convert first the FADC
139// slices to Photo-electrons and then, depending on the outer parameters,
140// look up the effective quantum efficiency and get the mean number of
141// photons which is then used for the further analysis.
142//
143// How the (first) MAGIC calibration should look like:
144// ===================================================
145//
146// For the moment, we have only one reliable calibration method, although
147// with very large systematic errors. This is the F-Factor method. Knowing
148// that the light is uniform over the whole camera (which I would not at all
149// guarantee in the case of the CT1 pulser), one could in principle already
150// perform a relative calibration of the quantum efficiencies in the UV.
151// However, the spread in QE at UV is about 10-15% (according to the plot
152// that Abelardo sent around last time. The spread in photo-electrons is 15%
153// for the inner pixels, but much larger (40%) for the outer ones.
154//
155// I'm not sure if we can already say that we have measured the relative
156// difference in quantum efficiency for the inner pixels and produce a first
157// QE-table for each pixel. To so, I would rather check in other wavelengths
158// (which we can do in about one-two weeks when the optical transmission of
159// the calibration trigger is installed).
160//
161// Thus, for the moment being, I would join Thomas proposal to calibrate in
162// photo-electrons and apply one stupid average quantum efficiency for all
163// pixels. This keeping in mind that we will have much preciser information
164// in about one to two weeks.
165//
166//
167// What MCalibrate should calculate and what should be stored:
168// ===========================================================
169//
170// It is clear that in the end, MCerPhotEvt will store photons.
171// MCalibrationCam stores the conversionfactors to photo-electrons and also
172// some tables of how to apply the conversion to photons, given the outer
173// parameters. This is not yet implemented and not even discussed.
174//
175// To start, I would suggest that we define the "average quantum efficiency"
176// (maybe something like 25+-3%) and apply them equally to all
177// photo-electrons. Later, this average factor can be easily replaced by a
178// pixel-dependent factor and later by a (pixel-dependent) table.
179//
180//
181//
182//////////////////////////////////////////////////////////////////////////////
183#include "MCalibrationChargeCalc.h"
184
185#include <TSystem.h>
186#include <TH1.h>
187#include <TF1.h>
188
189#include "MLog.h"
190#include "MLogManip.h"
191
192#include "MParList.h"
193
194#include "MStatusDisplay.h"
195
196#include "MRawEvtHeader.h"
197
198#include "MGeomCam.h"
199#include "MGeomPix.h"
200#include "MHCamera.h"
201
202#include "MPedestalCam.h"
203#include "MPedestalPix.h"
204
205#include "MCalibrationIntensityChargeCam.h"
206#include "MCalibrationIntensityQECam.h"
207#include "MCalibrationIntensityBlindCam.h"
208
209#include "MHCalibrationChargeCam.h"
210#include "MHCalibrationChargeBlindCam.h"
211
212#include "MCalibrationChargeCam.h"
213#include "MCalibrationChargePix.h"
214#include "MCalibrationChargePINDiode.h"
215#include "MCalibrationBlindPix.h"
216#include "MCalibrationBlindCam.h"
217
218#include "MExtractedSignalCam.h"
219#include "MExtractedSignalPix.h"
220#include "MExtractedSignalBlindPixel.h"
221#include "MExtractedSignalPINDiode.h"
222
223#include "MBadPixelsIntensityCam.h"
224#include "MBadPixelsCam.h"
225
226#include "MCalibrationQECam.h"
227#include "MCalibrationQEPix.h"
228
229ClassImp(MCalibrationChargeCalc);
230
231using namespace std;
232
233const Float_t MCalibrationChargeCalc::fgChargeLimit = 2.5;
234const Float_t MCalibrationChargeCalc::fgChargeErrLimit = 0.;
235const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit = 1.;
236const Float_t MCalibrationChargeCalc::fgLambdaErrLimit = 0.2;
237const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit = 0.5;
238const Float_t MCalibrationChargeCalc::fgPheErrLimit = 4.5;
239const Float_t MCalibrationChargeCalc::fgFFactorErrLimit = 4.5;
240const TString MCalibrationChargeCalc::fgNamePedestalCam = "MPedestalCam";
241// --------------------------------------------------------------------------
242//
243// Default constructor.
244//
245// Sets the pointer to fQECam and fGeom to NULL
246//
247// Calls AddToBranchList for:
248// - MRawEvtData.fHiGainPixId
249// - MRawEvtData.fLoGainPixId
250// - MRawEvtData.fHiGainFadcSamples
251// - MRawEvtData.fLoGainFadcSamples
252//
253// Initializes:
254// - fChargeLimit to fgChargeLimit
255// - fChargeErrLimit to fgChargeErrLimit
256// - fChargeRelErrLimit to fgChargeRelErrLimit
257// - fFFactorErrLimit to fgFFactorErrLimit
258// - fLambdaCheckLimit to fgLambdaCheckLimit
259// - fLambdaErrLimit to fgLambdaErrLimit
260// - fNamePedestalCam to fgNamePedestalCam
261// - fPheErrLimit to fgPheErrLimit
262// - fPulserColor to MCalibrationCam::kCT1
263// - fOutputPath to "."
264// - fOutputFile to "ChargeCalibStat.txt"
265// - flag debug to kFALSE
266//
267// Sets all checks
268//
269// Calls:
270// - Clear()
271//
272MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
273 : fGeom(NULL), fSignal(NULL), fHeader(NULL)
274{
275
276 fName = name ? name : "MCalibrationChargeCalc";
277 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
278
279 AddToBranchList("MRawEvtData.fHiGainPixId");
280 AddToBranchList("MRawEvtData.fLoGainPixId");
281 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
282 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
283
284 SetChargeLimit ();
285 SetChargeErrLimit ();
286 SetChargeRelErrLimit ();
287 SetFFactorErrLimit ();
288 SetLambdaCheckLimit ();
289 SetLambdaErrLimit ();
290 SetNamePedestalCam ();
291 SetPheErrLimit ();
292 SetOutputPath ();
293 SetOutputFile ();
294 SetDebug ( kFALSE );
295
296 SetCheckDeadPixels ();
297 SetCheckDeviatingBehavior();
298 SetCheckExtractionWindow ();
299 SetCheckHistOverflow ();
300 SetCheckOscillations ();
301
302 Clear();
303
304}
305
306// --------------------------------------------------------------------------
307//
308// Sets:
309// - all variables to 0.,
310// - all flags to kFALSE
311// - all pointers to NULL
312// - the pulser colour to kNONE
313// - fBlindPixelFlags to 0
314// - fPINDiodeFlags to 0
315//
316void MCalibrationChargeCalc::Clear(const Option_t *o)
317{
318
319 fNumHiGainSamples = 0.;
320 fNumLoGainSamples = 0.;
321 fSqrtHiGainSamples = 0.;
322 fSqrtLoGainSamples = 0.;
323 fNumInnerFFactorMethodUsed = 0;
324
325 fNumProcessed = 0;
326
327 fIntensBad = NULL;
328 fBadPixels = NULL;
329 fIntensCam = NULL;
330 fCam = NULL;
331 fHCam = NULL;
332 fIntensQE = NULL;
333 fQECam = NULL;
334 fIntensBlind = NULL;
335 fBlindCam = NULL;
336 fHBlindCam = NULL;
337 fPINDiode = NULL;
338 fPedestals = NULL;
339
340 SetPulserColor ( MCalibrationCam::kNONE );
341
342 fBlindPixelFlags.Set(0);
343 fPINDiodeFlags .Set(0);
344 fResultFlags .Set(0);
345}
346
347
348// -----------------------------------------------------------------------------------
349//
350// The following container are searched for and execution aborted if not in MParList:
351// - MRawEvtHeader
352// - MPedestalCam
353// - MExtractedSignalCam
354//
355Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
356{
357
358 fHeader = (MRawEvtHeader*)pList->FindObject("MRawEvtHeader");
359 if (!fHeader)
360 {
361 *fLog << err << "MRawEvtHeader not found... abort." << endl;
362 return kFALSE;
363 }
364
365 //
366 // Containers that have to be there.
367 //
368 fSignal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
369 if (!fSignal)
370 {
371 *fLog << err << "MExtractedSignalCam not found... aborting" << endl;
372 return kFALSE;
373 }
374
375 if (fPedestals)
376 return kTRUE;
377
378 fPedestals = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCam), "MPedestalCam");
379 if (!fPedestals)
380 {
381 *fLog << err << GetDescriptor() << fNamePedestalCam << " not found... aborting" << endl;
382 return kFALSE;
383 }
384
385 return kTRUE;
386}
387
388
389// --------------------------------------------------------------------------
390//
391// Search for the following input containers and abort if not existing:
392// - MGeomCam
393// - MCalibrationIntensityChargeCam or MCalibrationChargeCam
394// - MCalibrationIntensityQECam or MCalibrationQECam
395// - MBadPixelsIntensityCam or MBadPixelsCam
396//
397// Search for the following input containers and give a warning if not existing:
398// - MCalibrationBlindPix
399// - MCalibrationChargePINDiode
400//
401// It retrieves the following variables from MCalibrationChargeCam:
402//
403// - fNumHiGainSamples
404// - fNumLoGainSamples
405//
406// It defines the PixId of every pixel in:
407//
408// - MCalibrationIntensityChargeCam
409// - MCalibrationChargeCam
410// - MCalibrationQECam
411//
412// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
413//
414// - MCalibrationChargePix
415// - MCalibrationQEPix
416//
417// Sets the pulser colour and tests if it has not changed w.r.t. fPulserColor in:
418//
419// - MCalibrationChargeCam
420// - MCalibrationBlindPix (if existing)
421// - MCalibrationChargePINDiode (if existing)
422//
423Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
424{
425
426 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
427 if (!fGeom)
428 {
429 *fLog << err << "No MGeomCam found... aborting." << endl;
430 return kFALSE;
431 }
432
433 fIntensCam = (MCalibrationIntensityChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityChargeCam"));
434 if (fIntensCam)
435 *fLog << inf << "Found MCalibrationIntensityChargeCam ... " << endl;
436 else
437 {
438 fCam = (MCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
439 if (!fCam)
440 {
441 *fLog << err << "Cannot find MCalibrationChargeCam ... abort." << endl;
442 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationChargeCam before..." << endl;
443 return kFALSE;
444 }
445 }
446
447 fHCam = (MHCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MHCalibrationChargeCam"));
448 if (!fHCam)
449 {
450 *fLog << err << "Cannot find MHCalibrationChargeCam ... abort." << endl;
451 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationChargeCam before..." << endl;
452 return kFALSE;
453 }
454
455 fIntensQE = (MCalibrationIntensityQECam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityQECam"));
456 if (fIntensQE)
457 *fLog << inf << "Found MCalibrationIntensityQECam ... " << endl;
458 else
459 {
460 fQECam = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
461 if (!fQECam)
462 {
463 *fLog << err << "Cannot find MCalibrationQECam ... abort." << endl;
464 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationQECam before..." << endl;
465 return kFALSE;
466 }
467 }
468
469 fIntensBlind = (MCalibrationIntensityBlindCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityBlindCam"));
470 if (fIntensBlind)
471 *fLog << inf << "Found MCalibrationIntensityBlindCam ... " << endl;
472 else
473 {
474 fBlindCam = (MCalibrationBlindCam*)pList->FindObject(AddSerialNumber("MCalibrationBlindCam"));
475 if (!fBlindCam)
476 {
477 *fLog << endl;
478 *fLog << warn << GetDescriptor()
479 << ": No MCalibrationBlindCam found... no Blind Pixel method! " << endl;
480 }
481 }
482
483 fHBlindCam = (MHCalibrationChargeBlindCam*)pList->FindObject(AddSerialNumber("MHCalibrationChargeBlindCam"));
484 if (!fHBlindCam)
485 {
486 *fLog << endl;
487 *fLog << warn << GetDescriptor()
488 << ": No MHCalibrationChargeBlindCam found... no Blind Pixel method! " << endl;
489 }
490
491 fIntensBad = (MBadPixelsIntensityCam*)pList->FindObject(AddSerialNumber("MBadPixelsIntensityCam"));
492 if (fIntensBad)
493 *fLog << inf << "Found MBadPixelsIntensityCam ... " << endl;
494 else
495 {
496 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
497 if (!fBadPixels)
498 {
499 *fLog << err << "Cannot find MBadPixelsCam ... abort." << endl;
500 return kFALSE;
501 }
502 }
503
504 //
505 // Optional Containers
506 //
507 fPINDiode = (MCalibrationChargePINDiode*)pList->FindObject("MCalibrationChargePINDiode");
508 if (!fPINDiode)
509 {
510 *fLog << endl;
511 *fLog << warn << GetDescriptor()
512 << ": MCalibrationChargePINDiode not found... no PIN Diode method! " << endl;
513 }
514
515 MCalibrationQECam *qecam = fIntensQE
516 ? (MCalibrationQECam*) fIntensQE->GetCam() : fQECam;
517 MCalibrationChargeCam *chargecam = fIntensCam
518 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
519 MBadPixelsCam *badcam = fIntensBad
520 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
521
522 UInt_t npixels = fGeom->GetNumPixels();
523
524 for (UInt_t i=0; i<npixels; i++)
525 {
526
527 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
528 MCalibrationQEPix &pqe = (MCalibrationQEPix&) (*qecam) [i];
529 MBadPixelsPix &bad = (*badcam) [i];
530
531 if (bad.IsBad())
532 {
533 pix.SetExcluded();
534 pqe.SetExcluded();
535 continue;
536 }
537
538 if (IsDebug())
539 pix.SetDebug();
540 }
541
542 return kTRUE;
543}
544
545// ----------------------------------------------------------------------------------
546//
547// Set the correct colour to the charge containers
548//
549Int_t MCalibrationChargeCalc::Process()
550{
551
552 const MCalibrationCam::PulserColor_t col = fHeader->GetPulserColor();
553
554 if (col == fPulserColor)
555 {
556 fNumProcessed++;
557 return kTRUE;
558 }
559
560 if (col == MCalibrationCam::kNONE)
561 return kTRUE;
562
563 //
564 // Now retrieve the colour and check if not various colours have been used
565 //
566 if (!fIntensCam)
567 {
568 if (fPulserColor != MCalibrationCam::kNONE)
569 {
570 *fLog << err << "Multiple colours used simultaneously in calibration file. Will reset histograms "
571 << " and merge the MCalibrationQECam!" << endl;
572 fHCam->Finalize();
573 fHBlindCam->Finalize();
574 fHCam->ResetHists();
575 fHBlindCam->ResetHists();
576 Finalize();
577 }
578 }
579
580 fPulserColor = col;
581
582 *fLog << inf << "Found new colour ... " << flush;
583
584 switch (col)
585 {
586 case MCalibrationCam::kGREEN: *fLog << "Green"; break;
587 case MCalibrationCam::kBLUE: *fLog << "Blue"; break;
588 case MCalibrationCam::kUV: *fLog << "UV"; break;
589 case MCalibrationCam::kCT1: *fLog << "CT1"; break;
590 default: break;
591 }
592
593 *fLog << inf << " with strength: " << fHeader->GetPulserStrength() << endl;
594
595 if (fPINDiode)
596 fPINDiode->SetColor( fPulserColor );
597
598 fNumProcessed = 0;
599
600 return kTRUE;
601}
602
603// -----------------------------------------------------------------------
604//
605// Return if number of executions is null.
606//
607Int_t MCalibrationChargeCalc::PostProcess()
608{
609
610 if (GetNumExecutions()==0)
611 return kFALSE;
612
613 if (fPulserColor == MCalibrationCam::kNONE)
614 return kTRUE;
615
616 if (fNumProcessed == 0)
617 return kTRUE;
618
619 *fLog << endl;
620
621 return Finalize();
622}
623
624// -----------------------------------------------------------------------
625//
626// Return kTRUE if fPulserColor is kNONE
627//
628// First loop over pixels, average areas and sectors, call:
629// - FinalizePedestals()
630// - FinalizeCharges()
631// for every entry. Count number of valid pixels in loop and return kFALSE
632// if there are none (the "Michele check").
633//
634// Call FinalizeBadPixels()
635//
636// Call FinalizeFFactorMethod() (second and third loop over pixels and areas)
637//
638// Call FinalizeBlindCam()
639// Call FinalizePINDiode()
640//
641// Call FinalizeFFactorQECam() (fourth loop over pixels and areas)
642// Call FinalizeBlindPixelQECam() (fifth loop over pixels and areas)
643// Call FinalizePINDiodeQECam() (sixth loop over pixels and areas)
644//
645// Call FinalizeUnsuitablePixels()
646//
647// Call MParContainer::SetReadyToSave() for fIntensCam, fCam, fQECam, fBadPixels and
648// fBlindCam and fPINDiode if they exist
649//
650// Print out some statistics
651//
652Int_t MCalibrationChargeCalc::Finalize()
653{
654
655 fNumHiGainSamples = fSignal->GetNumUsedHiGainFADCSlices();
656 fNumLoGainSamples = fSignal->GetNumUsedLoGainFADCSlices();
657
658 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
659 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
660
661 if (fPINDiode)
662 if (!fPINDiode->IsValid())
663 {
664 *fLog << warn << GetDescriptor()
665 << ": MCalibrationChargePINDiode is declared not valid... no PIN Diode method! " << endl;
666 fPINDiode = NULL;
667 }
668
669 MCalibrationBlindCam *blindcam = fIntensBlind
670 ? (MCalibrationBlindCam*)fIntensBlind->GetCam() : fBlindCam;
671 MCalibrationQECam *qecam = fIntensQE
672 ? (MCalibrationQECam*) fIntensQE->GetCam() : fQECam;
673 MCalibrationChargeCam *chargecam = fIntensCam
674 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
675 MBadPixelsCam *badcam = fIntensBad
676 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
677
678 //
679 // First loop over pixels, call FinalizePedestals and FinalizeCharges
680 //
681 Int_t nvalid = 0;
682
683 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
684 {
685
686 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[pixid];
687 //
688 // Check if the pixel has been excluded from the fits
689 //
690 if (pix.IsExcluded())
691 continue;
692
693 MPedestalPix &ped = (*fPedestals)[pixid];
694 MBadPixelsPix &bad = (*badcam) [pixid];
695
696 const Int_t aidx = (*fGeom)[pixid].GetAidx();
697
698 FinalizePedestals(ped,pix,aidx);
699
700 if (FinalizeCharges(pix,bad,"pixel "))
701 nvalid++;
702 }
703
704 *fLog << endl;
705 //
706 // The Michele check ...
707 //
708 if (nvalid == 0)
709 {
710 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
711 << "Did you forget to fill the histograms "
712 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
713 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
714 << "instead of a calibration run " << endl;
715 return kFALSE;
716 }
717
718 for (UInt_t aidx=0; aidx<fGeom->GetNumAreas(); aidx++)
719 {
720
721 const MPedestalPix &ped = fPedestals->GetAverageArea(aidx);
722 MCalibrationChargePix &pix = (MCalibrationChargePix&)chargecam->GetAverageArea(aidx);
723
724 FinalizePedestals(ped,pix,aidx);
725 FinalizeCharges(pix,
726 fIntensCam ? fIntensCam->GetAverageBadArea(aidx) : fCam->GetAverageBadArea(aidx),
727 "area id");
728 }
729
730 *fLog << endl;
731
732 for (UInt_t sector=0; sector<fGeom->GetNumSectors(); sector++)
733 {
734
735 const MPedestalPix &ped = fPedestals->GetAverageSector(sector);
736
737 MCalibrationChargePix &pix = (MCalibrationChargePix&)chargecam->GetAverageSector(sector);
738 FinalizePedestals(ped,pix, 0);
739 }
740
741 *fLog << endl;
742
743 //
744 // Finalize Bad Pixels
745 //
746 FinalizeBadPixels();
747
748 //
749 // Finalize F-Factor method
750 //
751 if (FinalizeFFactorMethod())
752 chargecam->SetFFactorMethodValid(kTRUE);
753 else
754 {
755 *fLog << warn << "Could not calculate the photons flux from the F-Factor method " << endl;
756 chargecam->SetFFactorMethodValid(kFALSE);
757 return kFALSE;
758 }
759
760 *fLog << endl;
761 //
762 // Finalize Blind Pixel
763 //
764 qecam->SetBlindPixelMethodValid(FinalizeBlindCam());
765
766 //
767 // Finalize PIN Diode
768 //
769 qecam->SetBlindPixelMethodValid(FinalizePINDiode());
770
771 //
772 // Finalize QE Cam
773 //
774 FinalizeFFactorQECam();
775 FinalizeBlindPixelQECam();
776 FinalizePINDiodeQECam();
777 FinalizeCombinedQECam();
778
779 //
780 // Re-direct the output to an ascii-file from now on:
781 //
782 MLog *oldlog = fLog;
783 MLog asciilog;
784 if (!fOutputFile.IsNull())
785 {
786 asciilog.SetOutputFile(GetOutputFile(),kTRUE);
787 SetLogStream(&asciilog);
788 }
789
790 //
791 // Finalize calibration statistics
792 //
793 FinalizeUnsuitablePixels();
794
795 chargecam->SetReadyToSave();
796 qecam ->SetReadyToSave();
797 badcam ->SetReadyToSave();
798
799 if (blindcam)
800 blindcam->SetReadyToSave();
801 if (fPINDiode)
802 fPINDiode->SetReadyToSave();
803
804 *fLog << inf << endl;
805 *fLog << GetDescriptor() << ": Fatal errors statistics:" << endl;
806
807 PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal,
808 Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS: "));
809 PrintUncalibrated(MBadPixelsPix::kChargeRelErrNotValid,
810 Form("%s%2.1f%s","Signal Error bigger than ",fChargeRelErrLimit," times Mean Signal: "));
811 PrintUncalibrated(MBadPixelsPix::kLoGainSaturation,
812 "Pixels with Low Gain Saturation: ");
813 PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin,
814 Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s): "));
815 PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins,
816 Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s): "));
817 PrintUncalibrated(MBadPixelsPix::kHiGainOverFlow,
818 "Pixels with High Gain Overflow: ");
819 PrintUncalibrated(MBadPixelsPix::kLoGainOverFlow,
820 "Pixels with Low Gain Overflow : ");
821
822 *fLog << inf << endl;
823 *fLog << GetDescriptor() << ": Unreliable errors statistics:" << endl;
824
825 PrintUncalibrated(MBadPixelsPix::kChargeSigmaNotValid,
826 "Signal Sigma smaller than Pedestal RMS: ");
827 PrintUncalibrated(MBadPixelsPix::kHiGainOscillating,
828 "Pixels with changing Hi Gain signal over time: ");
829 PrintUncalibrated(MBadPixelsPix::kLoGainOscillating,
830 "Pixels with changing Lo Gain signal over time: ");
831 PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted,
832 "Pixels with unsuccesful Gauss fit to the Hi Gain: ");
833 PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted,
834 "Pixels with unsuccesful Gauss fit to the Lo Gain: ");
835 PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes,
836 "Pixels with deviating number of phes: ");
837 PrintUncalibrated(MBadPixelsPix::kDeviatingFFactor,
838 "Pixels with deviating F-Factor: ");
839
840 if (!fOutputFile.IsNull())
841 SetLogStream(oldlog);
842
843 return kTRUE;
844}
845
846// ----------------------------------------------------------------------------------
847//
848// Retrieves pedestal and pedestal RMS from MPedestalPix
849// Retrieves total entries from MPedestalCam
850// Sets pedestal*fNumHiGainSamples and pedestal*fNumLoGainSamples in MCalibrationChargePix
851// Sets pedRMS *fSqrtHiGainSamples and pedRMS *fSqrtLoGainSamples in MCalibrationChargePix
852//
853// If the flag MCalibrationPix::IsHiGainSaturation() is set, call also:
854// - MCalibrationChargePix::CalcLoGainPedestal()
855//
856void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal, const Int_t aidx)
857{
858
859 //
860 // get the pedestals
861 //
862 const Float_t pedes = ped.GetPedestal();
863 const Float_t prms = ped.GetPedestalRms();
864 const Int_t num = fPedestals->GetTotalEntries();
865
866 //
867 // RMS error set by PedCalcFromLoGain, 0 in case MPedCalcPedRun was used.
868 //
869 const Float_t prmserr = num>0 ? prms/TMath::Sqrt(2.*num) : ped.GetPedestalRmsError();
870
871 //
872 // set them in the calibration camera
873 //
874 if (cal.IsHiGainSaturation())
875 {
876 cal.SetPedestal(pedes * fNumLoGainSamples,
877 prms * fSqrtLoGainSamples,
878 prmserr * fSqrtLoGainSamples);
879 cal.CalcLoGainPedestal(fNumLoGainSamples);
880 }
881 else
882 {
883
884 cal.SetPedestal(pedes * fNumHiGainSamples,
885 prms * fSqrtHiGainSamples,
886 prmserr * fSqrtHiGainSamples);
887 }
888
889}
890
891// ----------------------------------------------------------------------------------------------------
892//
893// Check fit results validity. Bad Pixels flags are set if:
894//
895// 1) Pixel has a mean smaller than fChargeLimit*PedRMS ( Flag: MBadPixelsPix::kChargeIsPedestal)
896// 2) Pixel has a mean error smaller than fChargeErrLimit ( Flag: MBadPixelsPix::kChargeErrNotValid)
897// 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error
898// ( Flag: MBadPixelsPix::kChargeRelErrNotValid)
899// 4) Pixel has a sigma bigger than its Pedestal RMS ( Flag: MBadPixelsPix::kChargeSigmaNotValid )
900//
901// Further returns if flags: MBadPixelsPix::kUnsuitableRun is set
902//
903// Calls MCalibrationChargePix::CalcReducedSigma() and sets flag: MBadPixelsPix::kChargeIsPedestal
904// and returns kFALSE if not succesful.
905//
906// Calls MCalibrationChargePix::CalcFFactor() and sets flag: MBadPixelsPix::kDeviatingNumPhes)
907// and returns kFALSE if not succesful.
908//
909// Calls MCalibrationChargePix::CalcConvFFactor()and sets flag: MBadPixelsPix::kDeviatingNumPhes)
910// and returns kFALSE if not succesful.
911//
912Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what)
913{
914
915 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
916 return kFALSE;
917
918 if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
919 {
920 *fLog << warn
921 << Form("Fitted Charge: %5.2f < %2.1f",cal.GetMean(),fChargeLimit)
922 << Form(" * Pedestal RMS %5.2f in %s%3i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
923 bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
924 }
925
926 if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr())
927 {
928 *fLog << warn
929 << Form("Fitted Charge: %4.2f < %2.1f",cal.GetMean(),fChargeRelErrLimit)
930 << Form(" * its error %4.2f in %s%3i",cal.GetMeanErr(),what,cal.GetPixId()) << endl;
931 bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
932 }
933
934 if (cal.GetSigma() < cal.GetPedRms())
935 {
936 *fLog << warn
937 << Form("Sigma of Fitted Charge: %6.2f <",cal.GetSigma())
938 << Form(" Ped. RMS=%5.2f in %s%3i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
939 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
940 return kFALSE;
941 }
942
943 if (!cal.CalcReducedSigma())
944 {
945 *fLog << warn
946 << Form("Could not calculate the reduced sigma in %s: ",what)
947 << Form(" %4i",cal.GetPixId())
948 << endl;
949 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
950 return kFALSE;
951 }
952
953 if (!cal.CalcFFactor())
954 {
955 *fLog << warn
956 << Form("Could not calculate the F-Factor in %s: ",what)
957 << Form(" %4i",cal.GetPixId())
958 << endl;
959 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
960 return kFALSE;
961 }
962
963 if (cal.GetPheFFactorMethod() < 0.)
964 {
965 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
966 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
967 cal.SetFFactorMethodValid(kFALSE);
968 return kFALSE;
969 }
970
971 if (!cal.CalcConvFFactor())
972 {
973 *fLog << warn
974 << Form("Could not calculate the Conv. FADC counts to Phes in %s: ",what)
975 << Form(" %4i",cal.GetPixId())
976 << endl;
977 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
978 return kFALSE;
979 }
980
981 return kTRUE;
982}
983
984
985
986// -----------------------------------------------------------------------------------
987//
988// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
989// - MBadPixelsPix::kChargeIsPedestal
990// - MBadPixelsPix::kChargeRelErrNotValid
991// - MBadPixelsPix::kMeanTimeInFirstBin
992// - MBadPixelsPix::kMeanTimeInLast2Bins
993// - MBadPixelsPix::kDeviatingNumPhes
994// - MBadPixelsPix::kHiGainOverFlow
995// - MBadPixelsPix::kLoGainOverFlow
996//
997// - Call MCalibrationPix::SetExcluded() for the bad pixels
998//
999// Sets pixel to MBadPixelsPix::kUnreliableRun, if one of the following flags is set:
1000// - MBadPixelsPix::kChargeSigmaNotValid
1001//
1002void MCalibrationChargeCalc::FinalizeBadPixels()
1003{
1004
1005 MBadPixelsCam *badcam = fIntensBad
1006 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
1007 MCalibrationChargeCam *chargecam = fIntensCam
1008 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
1009
1010 for (Int_t i=0; i<badcam->GetSize(); i++)
1011 {
1012
1013 MBadPixelsPix &bad = (*badcam) [i];
1014 MCalibrationPix &pix = (*chargecam)[i];
1015
1016 if (IsCheckDeadPixels())
1017 {
1018 if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal))
1019 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1020
1021 if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid ))
1022 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1023
1024 if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ))
1025 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1026 }
1027
1028 if (IsCheckExtractionWindow())
1029 {
1030 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ))
1031 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1032
1033 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ))
1034 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1035 }
1036
1037 if (IsCheckDeviatingBehavior())
1038 {
1039 if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes ))
1040 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1041 }
1042
1043 if (IsCheckHistOverflow())
1044 {
1045 if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOverFlow ))
1046 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1047
1048 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOverFlow ))
1049 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
1050 }
1051
1052 if (bad.IsUnsuitable( MBadPixelsPix::kUnsuitableRun ))
1053 pix.SetExcluded();
1054
1055 if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ))
1056 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
1057 }
1058}
1059
1060// ------------------------------------------------------------------------
1061//
1062//
1063// First loop: Calculate a mean and mean RMS of photo-electrons per area index
1064// Include only pixels which are not MBadPixelsPix::kUnsuitableRun nor
1065// MBadPixelsPix::kChargeSigmaNotValid (see FinalizeBadPixels()) and set
1066// MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case.
1067//
1068// Second loop: Get mean number of photo-electrons and its RMS including
1069// only pixels with flag MCalibrationChargePix::IsFFactorMethodValid()
1070// and further exclude those deviating by more than fPheErrLimit mean
1071// sigmas from the mean (obtained in first loop). Set
1072// MBadPixelsPix::kDeviatingNumPhes if excluded.
1073//
1074// For the suitable pixels with flag MBadPixelsPix::kChargeSigmaNotValid
1075// set the number of photo-electrons as the mean number of photo-electrons
1076// calculated in that area index.
1077//
1078// Set weighted mean and variance of photo-electrons per area index in:
1079// average area pixels of MCalibrationChargeCam (obtained from:
1080// MCalibrationChargeCam::GetAverageArea() )
1081//
1082// Set weighted mean and variance of photo-electrons per sector in:
1083// average sector pixels of MCalibrationChargeCam (obtained from:
1084// MCalibrationChargeCam::GetAverageSector() )
1085//
1086//
1087// Third loop: Set mean number of photo-electrons and its RMS in the pixels
1088// only excluded as: MBadPixelsPix::kChargeSigmaNotValid
1089//
1090Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod()
1091{
1092
1093 MBadPixelsCam *badcam = fIntensBad
1094 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
1095 MCalibrationChargeCam *chargecam = fIntensCam
1096 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
1097
1098 const Int_t npixels = fGeom->GetNumPixels();
1099 const Int_t nareas = fGeom->GetNumAreas();
1100 const Int_t nsectors = fGeom->GetNumSectors();
1101
1102 TArrayF lowlim (nareas);
1103 TArrayF upplim (nareas);
1104 TArrayD areavars (nareas);
1105 TArrayD areaweights (nareas);
1106 TArrayD sectorweights (nsectors);
1107 TArrayD areaphes (nareas);
1108 TArrayD sectorphes (nsectors);
1109 TArrayI numareavalid (nareas);
1110 TArrayI numsectorvalid(nsectors);
1111
1112 //
1113 // First loop: Get mean number of photo-electrons and the RMS
1114 // The loop is only to recognize later pixels with very deviating numbers
1115 //
1116 MHCamera camphes(*fGeom,"Camphes","Phes in Camera");
1117
1118 for (Int_t i=0; i<npixels; i++)
1119 {
1120
1121 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1122 MBadPixelsPix &bad = (*badcam)[i];
1123
1124 if (!pix.IsFFactorMethodValid())
1125 continue;
1126
1127 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1128 {
1129 pix.SetFFactorMethodValid(kFALSE);
1130 continue;
1131 }
1132
1133 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1134 continue;
1135
1136 const Float_t nphe = pix.GetPheFFactorMethod();
1137 const Int_t aidx = (*fGeom)[i].GetAidx();
1138
1139 camphes.Fill(i,nphe);
1140 camphes.SetUsed(i);
1141
1142 areaphes [aidx] += nphe;
1143 areavars [aidx] += nphe*nphe;
1144 numareavalid[aidx] ++;
1145 }
1146
1147 for (Int_t i=0; i<nareas; i++)
1148 {
1149 if (numareavalid[i] == 0)
1150 {
1151 *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
1152 << "in area index: " << i << endl;
1153 continue;
1154 }
1155
1156 if (numareavalid[i] == 1)
1157 areavars[i] = 0.;
1158 else if (numareavalid[i] == 0)
1159 {
1160 areaphes[i] = -1.;
1161 areaweights[i] = -1.;
1162 }
1163 else
1164 {
1165 areavars[i] = (areavars[i] - areaphes[i]*areaphes[i]/numareavalid[i]) / (numareavalid[i]-1);
1166 areaphes[i] = areaphes[i] / numareavalid[i];
1167 }
1168
1169 if (areavars[i] < 0.)
1170 {
1171 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of photo-electrons found "
1172 << "in area index: " << i << endl;
1173 continue;
1174 }
1175
1176 lowlim [i] = areaphes[i] - fPheErrLimit*TMath::Sqrt(areavars[i]);
1177 upplim [i] = areaphes[i] + fPheErrLimit*TMath::Sqrt(areavars[i]);
1178
1179 TH1D *hist = camphes.ProjectionS(TArrayI(),TArrayI(1,&i),"_py",100);
1180 hist->Fit("gaus","Q");
1181 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1182 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1183 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1184
1185 if (IsDebug())
1186 hist->DrawClone();
1187
1188 if (ndf < 2)
1189 {
1190 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1191 << "in the camera with area index: " << i << endl;
1192 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1193 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1194 delete hist;
1195 continue;
1196 }
1197
1198 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1199
1200 if (prob < 0.001)
1201 {
1202 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1203 << "in the camera with area index: " << i << endl;
1204 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1205 << " is smaller than 0.001 " << endl;
1206 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1207 delete hist;
1208 continue;
1209 }
1210
1211 if (mean < 0.)
1212 {
1213 *fLog << inf << GetDescriptor() << ": Fitted mean number of photo-electrons "
1214 << "with area idx " << i << ": " << mean << " is smaller than 0. " << endl;
1215 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1216 delete hist;
1217 continue;
1218 }
1219
1220 *fLog << inf << GetDescriptor() << ": Mean number of photo-electrons "
1221 << "with area idx " << i << ": "
1222 << Form("%7.2f+-%6.2f",mean,sigma) << endl;
1223
1224 lowlim [i] = mean - fPheErrLimit*sigma;
1225 upplim [i] = mean + fPheErrLimit*sigma;
1226
1227 delete hist;
1228 }
1229
1230 *fLog << endl;
1231
1232 numareavalid.Reset();
1233 areaphes .Reset();
1234 areavars .Reset();
1235 //
1236 // Second loop: Get mean number of photo-electrons and its RMS excluding
1237 // pixels deviating by more than fPheErrLimit sigma.
1238 // Set the conversion factor FADC counts to photo-electrons
1239 //
1240 for (Int_t i=0; i<npixels; i++)
1241 {
1242
1243 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1244
1245 if (!pix.IsFFactorMethodValid())
1246 continue;
1247
1248 MBadPixelsPix &bad = (*badcam)[i];
1249
1250 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1251 continue;
1252
1253 const Float_t nvar = pix.GetPheFFactorMethodVar();
1254 if (nvar <= 0.)
1255 {
1256 pix.SetFFactorMethodValid(kFALSE);
1257 continue;
1258 }
1259
1260 const Int_t aidx = (*fGeom)[i].GetAidx();
1261 const Int_t sector = (*fGeom)[i].GetSector();
1262 const Float_t area = (*fGeom)[i].GetA();
1263 const Float_t nphe = pix.GetPheFFactorMethod();
1264
1265 if ( nphe < lowlim[aidx] || nphe > upplim[aidx] )
1266 {
1267 *fLog << warn << GetDescriptor() << ": Number of phes: "
1268 << Form("%7.2f out of %3.1f sigma limit: ",nphe,fPheErrLimit)
1269 << Form("[%7.2f,%7.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
1270 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1271 if (IsCheckDeviatingBehavior())
1272 {
1273 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1274 pix.SetFFactorMethodValid(kFALSE);
1275 }
1276 continue;
1277 }
1278
1279 areaweights [aidx] += nphe*nphe;
1280 areaphes [aidx] += nphe;
1281 numareavalid [aidx] ++;
1282
1283 if (aidx == 0)
1284 fNumInnerFFactorMethodUsed++;
1285
1286 sectorweights [sector] += nphe*nphe/area/area;
1287 sectorphes [sector] += nphe/area;
1288 numsectorvalid[sector] ++;
1289 }
1290
1291 *fLog << endl;
1292
1293 for (Int_t aidx=0; aidx<nareas; aidx++)
1294 {
1295
1296 MCalibrationChargePix &apix = (MCalibrationChargePix&)chargecam->GetAverageArea(aidx);
1297
1298 if (numareavalid[aidx] == 1)
1299 areaweights[aidx] = 0.;
1300 else if (numareavalid[aidx] == 0)
1301 {
1302 areaphes[aidx] = -1.;
1303 areaweights[aidx] = -1.;
1304 }
1305 else
1306 {
1307 areaweights[aidx] = (areaweights[aidx] - areaphes[aidx]*areaphes[aidx]/numareavalid[aidx])
1308 / (numareavalid[aidx]-1);
1309 areaphes[aidx] /= numareavalid[aidx];
1310 }
1311
1312 if (areaweights[aidx] < 0. || areaphes[aidx] <= 0.)
1313 {
1314 *fLog << warn << GetDescriptor()
1315 << ": Mean number phes from area index " << aidx << " could not be calculated: "
1316 << " Mean: " << areaphes[aidx]
1317 << " Variance: " << areaweights[aidx] << endl;
1318 apix.SetFFactorMethodValid(kFALSE);
1319 continue;
1320 }
1321
1322 *fLog << inf << GetDescriptor()
1323 << ": Average total number phes in area idx " << aidx << ": "
1324 << Form("%7.2f +- 6.2f",areaphes[aidx],TMath::Sqrt(areaweights[aidx])) << endl;
1325
1326 apix.SetPheFFactorMethod ( areaphes[aidx] );
1327 apix.SetPheFFactorMethodVar( areaweights[aidx] / numareavalid[aidx] );
1328 apix.SetFFactorMethodValid ( kTRUE );
1329
1330 }
1331
1332 *fLog << endl;
1333
1334 for (Int_t sector=0; sector<nsectors; sector++)
1335 {
1336
1337 if (numsectorvalid[sector] == 1)
1338 sectorweights[sector] = 0.;
1339 else if (numsectorvalid[sector] == 0)
1340 {
1341 sectorphes[sector] = -1.;
1342 sectorweights[sector] = -1.;
1343 }
1344 else
1345 {
1346 sectorweights[sector] = (sectorweights[sector]
1347 - sectorphes[sector]*sectorphes[sector]/numsectorvalid[sector]
1348 )
1349 / (numsectorvalid[sector]-1.);
1350 sectorphes[sector] /= numsectorvalid[sector];
1351 }
1352
1353 MCalibrationChargePix &spix = (MCalibrationChargePix&)chargecam->GetAverageSector(sector);
1354
1355 if (sectorweights[sector] < 0. || sectorphes[sector] <= 0.)
1356 {
1357 *fLog << warn << GetDescriptor()
1358 <<": Mean number phes per area for sector " << sector << " could not be calculated: "
1359 << " Mean: " << sectorphes[sector]
1360 << " Variance: " << sectorweights[sector] << endl;
1361 spix.SetFFactorMethodValid(kFALSE);
1362 continue;
1363 }
1364
1365 *fLog << inf << GetDescriptor()
1366 << ": Average number phes per area in sector " << sector << ": "
1367 << Form("%5.3f+-%4.3f [phe/mm^2]",sectorphes[sector],TMath::Sqrt(sectorweights[sector]))
1368 << endl;
1369
1370 spix.SetPheFFactorMethod ( sectorphes[sector] );
1371 spix.SetPheFFactorMethodVar( sectorweights[sector] / numsectorvalid[sector]);
1372 spix.SetFFactorMethodValid ( kTRUE );
1373
1374 }
1375
1376 //
1377 // Third loop: Set mean number of photo-electrons and its RMS in the pixels
1378 // only excluded as: MBadPixelsPix::kChargeSigmaNotValid
1379 //
1380 for (Int_t i=0; i<npixels; i++)
1381 {
1382
1383 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1384 MBadPixelsPix &bad = (*badcam)[i];
1385
1386 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1387 continue;
1388
1389 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1390 {
1391 const Int_t aidx = (*fGeom)[i].GetAidx();
1392 MCalibrationChargePix &apix = (MCalibrationChargePix&)chargecam->GetAverageArea(aidx);
1393
1394 pix.SetPheFFactorMethod ( apix.GetPheFFactorMethod() );
1395 pix.SetPheFFactorMethodVar( apix.GetPheFFactorMethodVar() );
1396 if (!pix.CalcConvFFactor())
1397 {
1398 *fLog << warn << GetDescriptor()
1399 << ": Could not calculate the Conv. FADC counts to Phes in pixel: "
1400 << Form(" %4i",pix.GetPixId())
1401 << endl;
1402 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1403 if (IsCheckDeviatingBehavior())
1404 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1405 }
1406 }
1407 }
1408
1409 return kTRUE;
1410}
1411
1412
1413
1414// ------------------------------------------------------------------------
1415//
1416// Returns kFALSE if pointer to MCalibrationBlindCam is NULL
1417//
1418// The check returns kFALSE if:
1419//
1420// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
1421// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
1422//
1423// Calls:
1424// - MCalibrationBlindPix::CalcFluxInsidePlexiglass()
1425//
1426Bool_t MCalibrationChargeCalc::FinalizeBlindCam()
1427{
1428
1429 MCalibrationBlindCam *blindcam = fIntensBlind
1430 ? (MCalibrationBlindCam*)fIntensBlind->GetCam() : fBlindCam;
1431
1432 if (!blindcam)
1433 return kFALSE;
1434
1435 Int_t nvalid = 0;
1436
1437 for (Int_t i=0; i<blindcam->GetSize(); i++)
1438 {
1439
1440 MCalibrationBlindPix &blindpix = (MCalibrationBlindPix&)(*blindcam)[i];
1441
1442 if (!blindpix.IsValid())
1443 continue;
1444
1445 const Float_t lambda = blindpix.GetLambda();
1446 const Float_t lambdaerr = blindpix.GetLambdaErr();
1447 const Float_t lambdacheck = blindpix.GetLambdaCheck();
1448
1449 if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
1450 {
1451 *fLog << warn << GetDescriptor()
1452 << Form("%s%4.2f%s%4.2f%s%4.2f%s%2i",": Lambda: ",lambda," and Lambda-Check: ",
1453 lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel Nr.",i)
1454 << endl;
1455 blindpix.SetValid(kFALSE);
1456 continue;
1457 }
1458
1459 if (lambdaerr > fLambdaErrLimit)
1460 {
1461 *fLog << warn << GetDescriptor()
1462 << Form("%s%4.2f%s%4.2f%s%2i",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
1463 fLambdaErrLimit," in Blind Pixel Nr.",i) << endl;
1464 blindpix.SetValid(kFALSE);
1465 continue;
1466 }
1467
1468 if (!blindpix.CalcFluxInsidePlexiglass())
1469 {
1470 *fLog << warn << "Could not calculate the flux of photons from Blind Pixel Nr." << i << endl;
1471 blindpix.SetValid(kFALSE);
1472 continue;
1473 }
1474
1475 nvalid++;
1476 }
1477
1478 if (!nvalid)
1479 return kFALSE;
1480
1481 return kTRUE;
1482}
1483
1484// ------------------------------------------------------------------------
1485//
1486// Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL
1487//
1488// The check returns kFALSE if:
1489//
1490// 1) PINDiode has a fitted charge smaller than fChargeLimit*PedRMS
1491// 2) PINDiode has a fit error smaller than fChargeErrLimit
1492// 3) PINDiode has a fitted charge smaller its fChargeRelErrLimit times its charge error
1493// 4) PINDiode has a charge sigma smaller than its Pedestal RMS
1494//
1495// Calls:
1496// - MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass()
1497//
1498Bool_t MCalibrationChargeCalc::FinalizePINDiode()
1499{
1500
1501 if (!fPINDiode)
1502 return kFALSE;
1503
1504 if (fPINDiode->GetMean() < fChargeLimit*fPINDiode->GetPedRms())
1505 {
1506 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1507 << fChargeLimit << " Pedestal RMS in PINDiode " << endl;
1508 return kFALSE;
1509 }
1510
1511 if (fPINDiode->GetMeanErr() < fChargeErrLimit)
1512 {
1513 *fLog << warn << GetDescriptor() << ": Error of Fitted Charge is smaller than "
1514 << fChargeErrLimit << " in PINDiode " << endl;
1515 return kFALSE;
1516 }
1517
1518 if (fPINDiode->GetMean() < fChargeRelErrLimit*fPINDiode->GetMeanErr())
1519 {
1520 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1521 << fChargeRelErrLimit << "* its error in PINDiode " << endl;
1522 return kFALSE;
1523 }
1524
1525 if (fPINDiode->GetSigma() < fPINDiode->GetPedRms())
1526 {
1527 *fLog << warn << GetDescriptor()
1528 << ": Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl;
1529 return kFALSE;
1530 }
1531
1532
1533 if (!fPINDiode->CalcFluxOutsidePlexiglass())
1534 {
1535 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
1536 << "will skip PIN Diode Calibration " << endl;
1537 return kFALSE;
1538 }
1539
1540 return kTRUE;
1541}
1542
1543// ------------------------------------------------------------------------
1544//
1545// Calculate the average number of photons outside the plexiglass with the
1546// formula:
1547//
1548// av.Num.photons(area index) = av.Num.Phes(area index)
1549// / MCalibrationQEPix::GetDefaultQE(fPulserColor)
1550// / MCalibrationQEPix::GetPMTCollectionEff()
1551// / MCalibrationQEPix::GetLightGuidesEff(fPulserColor)
1552// / MCalibrationQECam::GetPlexiglassQE()
1553//
1554// Calculate the variance on the average number of photons assuming that the error on the
1555// Quantum efficiency is reduced by the number of used inner pixels, but the rest of the
1556// values keeps it ordinary error since it is systematic.
1557//
1558// Loop over pixels:
1559//
1560// - Continue, if not MCalibrationChargePix::IsFFactorMethodValid() and set:
1561// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor)
1562//
1563// - Call MCalibrationChargePix::CalcMeanFFactor(av.Num.photons) and set:
1564// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) if not succesful
1565//
1566// - Calculate the quantum efficiency with the formula:
1567//
1568// QE = ( Num.Phes / av.Num.photons ) * MGeomCam::GetPixRatio()
1569//
1570// - Set QE in MCalibrationQEPix::SetQEFFactor ( QE, fPulserColor );
1571//
1572// - Set Variance of QE in MCalibrationQEPix::SetQEFFactorVar ( Variance, fPulserColor );
1573// - Set bit MCalibrationQEPix::SetFFactorMethodValid(kTRUE,fPulserColor)
1574//
1575// - Call MCalibrationQEPix::UpdateFFactorMethod()
1576//
1577void MCalibrationChargeCalc::FinalizeFFactorQECam()
1578{
1579
1580 if (fNumInnerFFactorMethodUsed < 2)
1581 {
1582 *fLog << warn << GetDescriptor()
1583 << ": Could not calculate F-Factor Method: Less than 2 inner pixels valid! " << endl;
1584 return;
1585 }
1586
1587 MCalibrationQECam *qecam = fIntensQE
1588 ? (MCalibrationQECam*) fIntensQE->GetCam() : fQECam;
1589 MCalibrationChargeCam *chargecam = fIntensCam
1590 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
1591 MBadPixelsCam *badcam = fIntensBad
1592 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
1593
1594 MCalibrationChargePix &avpix = (MCalibrationChargePix&)chargecam->GetAverageArea(0);
1595 MCalibrationQEPix &qepix = (MCalibrationQEPix&) qecam->GetAverageArea(0);
1596
1597 const Float_t avphotons = avpix.GetPheFFactorMethod()
1598 / qepix.GetDefaultQE(fPulserColor)
1599 / qepix.GetPMTCollectionEff()
1600 / qepix.GetLightGuidesEff(fPulserColor)
1601 / qecam->GetPlexiglassQE();
1602
1603 const Float_t avphotrelvar = avpix.GetPheFFactorMethodRelVar()
1604 + qepix.GetDefaultQERelVar(fPulserColor) / fNumInnerFFactorMethodUsed
1605 + qepix.GetPMTCollectionEffRelVar()
1606 + qepix.GetLightGuidesEffRelVar(fPulserColor)
1607 + qecam->GetPlexiglassQERelVar();
1608
1609 const UInt_t nareas = fGeom->GetNumAreas();
1610
1611 //
1612 // Set the results in the MCalibrationChargeCam
1613 //
1614 chargecam->SetNumPhotonsFFactorMethod (avphotons);
1615
1616 if (avphotrelvar > 0.)
1617 chargecam->SetNumPhotonsFFactorMethodErr(TMath::Sqrt( avphotrelvar * avphotons * avphotons));
1618
1619 TArrayF lowlim (nareas);
1620 TArrayF upplim (nareas);
1621 TArrayD avffactorphotons (nareas);
1622 TArrayD avffactorphotvar (nareas);
1623 TArrayI numffactor (nareas);
1624
1625 const UInt_t npixels = fGeom->GetNumPixels();
1626
1627 MHCamera camffactor(*fGeom,"Camffactor","F-Factor in Camera");
1628
1629 for (UInt_t i=0; i<npixels; i++)
1630 {
1631
1632 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1633 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*qecam) [i];
1634 MBadPixelsPix &bad = (*badcam) [i];
1635
1636 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1637 continue;
1638
1639 const Float_t photons = avphotons / fGeom->GetPixRatio(i);
1640 const Float_t qe = pix.GetPheFFactorMethod() / photons ;
1641
1642 const Float_t qerelvar = avphotrelvar + pix.GetPheFFactorMethodRelVar();
1643
1644 qepix.SetQEFFactor ( qe , fPulserColor );
1645 qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor );
1646 qepix.SetFFactorMethodValid( kTRUE , fPulserColor );
1647
1648 if (!qepix.UpdateFFactorMethod( qecam->GetPlexiglassQE() ))
1649 *fLog << warn << GetDescriptor()
1650 << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl;
1651
1652 //
1653 // The following pixels are those with deviating sigma, but otherwise OK,
1654 // probably those with stars during the pedestal run, but not the cal. run.
1655 //
1656 if (!pix.CalcMeanFFactor( photons , avphotrelvar ))
1657 {
1658 bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1659 if (IsCheckDeviatingBehavior())
1660 bad.SetUnsuitable ( MBadPixelsPix::kUnreliableRun );
1661 continue;
1662 }
1663
1664 const Int_t aidx = (*fGeom)[i].GetAidx();
1665 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1666
1667 camffactor.Fill(i,ffactor);
1668 camffactor.SetUsed(i);
1669
1670 avffactorphotons[aidx] += ffactor;
1671 avffactorphotvar[aidx] += ffactor*ffactor;
1672 numffactor[aidx]++;
1673 }
1674
1675 for (UInt_t i=0; i<nareas; i++)
1676 {
1677
1678 if (numffactor[i] == 0)
1679 {
1680 *fLog << warn << GetDescriptor() << ": No pixels with valid total F-Factor found "
1681 << "in area index: " << i << endl;
1682 continue;
1683 }
1684
1685 avffactorphotvar[i] = (avffactorphotvar[i] - avffactorphotons[i]*avffactorphotons[i]/numffactor[i])
1686 / (numffactor[i]-1.);
1687 avffactorphotons[i] = avffactorphotons[i] / numffactor[i];
1688
1689 if (avffactorphotvar[i] < 0.)
1690 {
1691 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of total F-Factor found "
1692 << "in area index: " << i << endl;
1693 continue;
1694 }
1695
1696 lowlim [i] = 1.; // Lowest known F-Factor of a PMT
1697 upplim [i] = avffactorphotons[i] + fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);
1698
1699 TArrayI area(1);
1700 area[0] = i;
1701
1702 TH1D *hist = camffactor.ProjectionS(TArrayI(),area,"_py",100);
1703 hist->Fit("gaus","Q");
1704 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1705 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1706 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1707
1708 if (IsDebug())
1709 camffactor.DrawClone();
1710
1711 if (ndf < 2)
1712 {
1713 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1714 << "in the camera with area index: " << i << endl;
1715 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1716 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1717 delete hist;
1718 continue;
1719 }
1720
1721 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1722
1723 if (prob < 0.001)
1724 {
1725 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1726 << "in the camera with area index: " << i << endl;
1727 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1728 << " is smaller than 0.001 " << endl;
1729 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1730 delete hist;
1731 continue;
1732 }
1733
1734 *fLog << inf << GetDescriptor() << ": Mean F-Factor "
1735 << "with area index #" << i << ": "
1736 << Form("%4.2f+-%4.2f",mean,sigma) << endl;
1737
1738 lowlim [i] = 1.;
1739 upplim [i] = mean + fFFactorErrLimit*sigma;
1740
1741 delete hist;
1742 }
1743
1744 *fLog << endl;
1745
1746 for (UInt_t i=0; i<npixels; i++)
1747 {
1748
1749 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1750 MBadPixelsPix &bad = (*badcam) [i];
1751
1752 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1753 continue;
1754
1755 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1756 const Int_t aidx = (*fGeom)[i].GetAidx();
1757
1758 if ( ffactor < lowlim[aidx] || ffactor > upplim[aidx] )
1759 {
1760 *fLog << warn << GetDescriptor() << ": Overall F-Factor "
1761 << Form("%5.2f",ffactor) << " out of range ["
1762 << Form("%5.2f,%5.2f",lowlim[aidx],upplim[aidx]) << "] pixel " << i << endl;
1763
1764 bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1765 if (IsCheckDeviatingBehavior())
1766 bad.SetUnsuitable ( MBadPixelsPix::kUnreliableRun );
1767 }
1768 }
1769
1770 for (UInt_t i=0; i<npixels; i++)
1771 {
1772
1773 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1774 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*qecam) [i];
1775 MBadPixelsPix &bad = (*badcam) [i];
1776
1777 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1778 {
1779 qepix.SetFFactorMethodValid(kFALSE,fPulserColor);
1780 pix.SetFFactorMethodValid(kFALSE);
1781 pix.SetExcluded();
1782 continue;
1783 }
1784 }
1785}
1786
1787
1788// ------------------------------------------------------------------------
1789//
1790// Loop over pixels:
1791//
1792// - Continue, if not MCalibrationBlindPix::IsFluxInsidePlexiglassAvailable() and set:
1793// MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor)
1794//
1795// - Calculate the quantum efficiency with the formula:
1796//
1797// QE = Num.Phes / MCalibrationBlindPix::GetFluxInsidePlexiglass()
1798// / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE()
1799//
1800// - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor );
1801// - Set Variance of QE in MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor );
1802// - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor)
1803//
1804// - Call MCalibrationQEPix::UpdateBlindPixelMethod()
1805//
1806void MCalibrationChargeCalc::FinalizeBlindPixelQECam()
1807{
1808
1809
1810 MCalibrationBlindCam *blindcam = fIntensBlind
1811 ? (MCalibrationBlindCam*) fIntensBlind->GetCam(): fBlindCam;
1812 MBadPixelsCam *badcam = fIntensBad
1813 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
1814 MCalibrationQECam *qecam = fIntensQE
1815 ? (MCalibrationQECam*) fIntensQE->GetCam() : fQECam;
1816 MCalibrationChargeCam *chargecam = fIntensCam
1817 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
1818
1819 //
1820 // Set the results in the MCalibrationChargeCam
1821 //
1822 if (!blindcam || !(blindcam->IsFluxInsidePlexiglassAvailable()))
1823 {
1824
1825 const Float_t photons = blindcam->GetFluxInsidePlexiglass() * (*fGeom)[0].GetA()
1826 / qecam->GetPlexiglassQE();
1827 chargecam->SetNumPhotonsBlindPixelMethod(photons);
1828
1829 const Float_t photrelvar = blindcam->GetFluxInsidePlexiglassRelVar()
1830 + qecam->GetPlexiglassQERelVar();
1831
1832 if (photrelvar > 0.)
1833 chargecam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1834 }
1835
1836 //
1837 // With the knowledge of the overall photon flux, calculate the
1838 // quantum efficiencies after the Blind Pixel and PIN Diode method
1839 //
1840 const UInt_t npixels = fGeom->GetNumPixels();
1841 for (UInt_t i=0; i<npixels; i++)
1842 {
1843
1844 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*qecam) [i];
1845
1846 if (!blindcam || !(blindcam->IsFluxInsidePlexiglassAvailable()))
1847 {
1848 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1849 continue;
1850 }
1851
1852 MBadPixelsPix &bad = (*badcam) [i];
1853
1854 if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1855 {
1856 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1857 continue;
1858 }
1859
1860 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1861 MGeomPix &geo = (*fGeom) [i];
1862
1863 const Float_t qe = pix.GetPheFFactorMethod()
1864 / blindcam->GetFluxInsidePlexiglass()
1865 / geo.GetA()
1866 * qecam->GetPlexiglassQE();
1867
1868 const Float_t qerelvar = blindcam->GetFluxInsidePlexiglassRelVar()
1869 + qecam->GetPlexiglassQERelVar()
1870 + pix.GetPheFFactorMethodRelVar();
1871
1872 qepix.SetQEBlindPixel ( qe , fPulserColor );
1873 qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor );
1874 qepix.SetBlindPixelMethodValid( kTRUE , fPulserColor );
1875
1876 if (!qepix.UpdateBlindPixelMethod( qecam->GetPlexiglassQE()))
1877 *fLog << warn << GetDescriptor()
1878 << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl;
1879 }
1880}
1881
1882// ------------------------------------------------------------------------
1883//
1884// Loop over pixels:
1885//
1886// - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set:
1887// MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor)
1888//
1889// - Calculate the quantum efficiency with the formula:
1890//
1891// QE = Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA()
1892//
1893// - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor );
1894// - Set Variance of QE in MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor );
1895// - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor)
1896//
1897// - Call MCalibrationQEPix::UpdatePINDiodeMethod()
1898//
1899void MCalibrationChargeCalc::FinalizePINDiodeQECam()
1900{
1901
1902 const UInt_t npixels = fGeom->GetNumPixels();
1903
1904 MCalibrationQECam *qecam = fIntensQE
1905 ? (MCalibrationQECam*) fIntensQE->GetCam() : fQECam;
1906 MCalibrationChargeCam *chargecam = fIntensCam
1907 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
1908 MBadPixelsCam *badcam = fIntensBad
1909 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
1910 //
1911 // With the knowledge of the overall photon flux, calculate the
1912 // quantum efficiencies after the PIN Diode method
1913 //
1914 for (UInt_t i=0; i<npixels; i++)
1915 {
1916
1917 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*qecam) [i];
1918
1919 if (!fPINDiode)
1920 {
1921 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1922 continue;
1923 }
1924
1925 if (!fPINDiode->IsFluxOutsidePlexiglassAvailable())
1926 {
1927 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1928 continue;
1929 }
1930
1931 MBadPixelsPix &bad = (*badcam) [i];
1932
1933 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1934 {
1935 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1936 continue;
1937 }
1938
1939 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i];
1940 MGeomPix &geo = (*fGeom) [i];
1941
1942 const Float_t qe = pix.GetPheFFactorMethod()
1943 / fPINDiode->GetFluxOutsidePlexiglass()
1944 / geo.GetA();
1945
1946 const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar();
1947
1948 qepix.SetQEPINDiode ( qe , fPulserColor );
1949 qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor );
1950 qepix.SetPINDiodeMethodValid( kTRUE , fPulserColor );
1951
1952 if (!qepix.UpdatePINDiodeMethod())
1953 *fLog << warn << GetDescriptor()
1954 << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl;
1955 }
1956}
1957
1958// ------------------------------------------------------------------------
1959//
1960// Loop over pixels:
1961//
1962// - Call MCalibrationQEPix::UpdateCombinedMethod()
1963//
1964void MCalibrationChargeCalc::FinalizeCombinedQECam()
1965{
1966
1967 const UInt_t npixels = fGeom->GetNumPixels();
1968
1969 MCalibrationQECam *qecam = fIntensQE
1970 ? (MCalibrationQECam*) fIntensQE->GetCam() : fQECam;
1971 MBadPixelsCam *badcam = fIntensBad
1972 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
1973
1974 for (UInt_t i=0; i<npixels; i++)
1975 {
1976
1977 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*qecam) [i];
1978 MBadPixelsPix &bad = (*badcam) [i];
1979
1980 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1981 {
1982 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1983 continue;
1984 }
1985
1986 qepix.UpdateCombinedMethod();
1987 }
1988}
1989
1990// -----------------------------------------------------------------------------------------------
1991//
1992// - Print out statistics about BadPixels of type UnsuitableType_t
1993// - store numbers of bad pixels of each type in fCam or fIntensCam
1994//
1995void MCalibrationChargeCalc::FinalizeUnsuitablePixels()
1996{
1997
1998 *fLog << inf << endl;
1999 *fLog << GetDescriptor() << ": Charge Calibration status:" << endl;
2000 *fLog << dec << setfill(' ');
2001
2002 const Int_t nareas = fGeom->GetNumAreas();
2003
2004 TArrayI counts(nareas);
2005
2006 MBadPixelsCam *badcam = fIntensBad
2007 ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels;
2008 MCalibrationChargeCam *chargecam = fIntensCam
2009 ? (MCalibrationChargeCam*)fIntensCam->GetCam() : fCam;
2010
2011 for (Int_t i=0; i<badcam->GetSize(); i++)
2012 {
2013 MBadPixelsPix &bad = (*badcam)[i];
2014 if (!bad.IsBad())
2015 {
2016 const Int_t aidx = (*fGeom)[i].GetAidx();
2017 counts[aidx]++;
2018 }
2019 }
2020
2021 if (fGeom->InheritsFrom("MGeomCamMagic"))
2022 *fLog << " " << setw(7) << "Successfully calibrated Pixels: "
2023 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
2024
2025 counts.Reset();
2026
2027 for (Int_t i=0; i<badcam->GetSize(); i++)
2028 {
2029 MBadPixelsPix &bad = (*badcam)[i];
2030
2031 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
2032 {
2033 const Int_t aidx = (*fGeom)[i].GetAidx();
2034 counts[aidx]++;
2035 }
2036 }
2037
2038 for (Int_t aidx=0; aidx<nareas; aidx++)
2039 chargecam->SetNumUnsuitable(counts[aidx], aidx);
2040
2041 if (fGeom->InheritsFrom("MGeomCamMagic"))
2042 *fLog << " " << setw(7) << "Uncalibrated Pixels: "
2043 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
2044
2045 counts.Reset();
2046
2047 for (Int_t i=0; i<badcam->GetSize(); i++)
2048 {
2049
2050 MBadPixelsPix &bad = (*badcam)[i];
2051
2052 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
2053 {
2054 const Int_t aidx = (*fGeom)[i].GetAidx();
2055 counts[aidx]++;
2056 }
2057 }
2058
2059 for (Int_t aidx=0; aidx<nareas; aidx++)
2060 chargecam->SetNumUnreliable(counts[aidx], aidx);
2061
2062 *fLog << " " << setw(7) << "Unreliable Pixels: "
2063 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
2064
2065}
2066
2067// -----------------------------------------------------------------------------------------------
2068//
2069// Print out statistics about BadPixels of type UncalibratedType_t
2070//
2071void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
2072{
2073
2074 UInt_t countinner = 0;
2075 UInt_t countouter = 0;
2076
2077 MBadPixelsCam *badcam = fIntensBad
2078 ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels;
2079
2080 for (Int_t i=0; i<badcam->GetSize(); i++)
2081 {
2082 MBadPixelsPix &bad = (*badcam)[i];
2083
2084 if (bad.IsUncalibrated(typ))
2085 {
2086 if (fGeom->GetPixRatio(i) == 1.)
2087 countinner++;
2088 else
2089 countouter++;
2090 }
2091 }
2092
2093 *fLog << " " << setw(7) << text
2094 << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
2095}
2096
2097// --------------------------------------------------------------------------
2098//
2099// Set the path for output file
2100//
2101void MCalibrationChargeCalc::SetOutputPath(TString path)
2102{
2103 fOutputPath = path;
2104 if (fOutputPath.EndsWith("/"))
2105 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
2106}
2107
2108// --------------------------------------------------------------------------
2109//
2110// Set the output file
2111//
2112void MCalibrationChargeCalc::SetOutputFile(TString file)
2113{
2114 fOutputFile = file;
2115}
2116
2117// --------------------------------------------------------------------------
2118//
2119// Get the output file
2120//
2121const char* MCalibrationChargeCalc::GetOutputFile()
2122{
2123 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
2124}
2125
2126// --------------------------------------------------------------------------
2127//
2128// Read the environment for the following data members:
2129// - fChargeLimit
2130// - fChargeErrLimit
2131// - fChargeRelErrLimit
2132// - fDebug
2133// - fFFactorErrLimit
2134// - fLambdaErrLimit
2135// - fLambdaCheckErrLimit
2136// - fPheErrLimit
2137//
2138Int_t MCalibrationChargeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
2139{
2140
2141 Bool_t rc = kFALSE;
2142 if (IsEnvDefined(env, prefix, "ChargeLimit", print))
2143 {
2144 SetChargeLimit(GetEnvValue(env, prefix, "ChargeLimit", fChargeLimit));
2145 rc = kTRUE;
2146 }
2147 if (IsEnvDefined(env, prefix, "ChargeErrLimit", print))
2148 {
2149 SetChargeErrLimit(GetEnvValue(env, prefix, "ChargeErrLimit", fChargeErrLimit));
2150 rc = kTRUE;
2151 }
2152 if (IsEnvDefined(env, prefix, "ChargeRelErrLimit", print))
2153 {
2154 SetChargeRelErrLimit(GetEnvValue(env, prefix, "ChargeRelErrLimit", fChargeRelErrLimit));
2155 rc = kTRUE;
2156 }
2157 if (IsEnvDefined(env, prefix, "Debug", print))
2158 {
2159 SetDebug(GetEnvValue(env, prefix, "Debug", IsDebug()));
2160 rc = kTRUE;
2161 }
2162 if (IsEnvDefined(env, prefix, "FFactorErrLimit", print))
2163 {
2164 SetFFactorErrLimit(GetEnvValue(env, prefix, "FFactorErrLimit", fFFactorErrLimit));
2165 rc = kTRUE;
2166 }
2167 if (IsEnvDefined(env, prefix, "LambdaErrLimit", print))
2168 {
2169 SetLambdaErrLimit(GetEnvValue(env, prefix, "LambdaErrLimit", fLambdaErrLimit));
2170 rc = kTRUE;
2171 }
2172 if (IsEnvDefined(env, prefix, "LambdaCheckLimit", print))
2173 {
2174 SetLambdaCheckLimit(GetEnvValue(env, prefix, "LambdaCheckLimit", fLambdaCheckLimit));
2175 rc = kTRUE;
2176 }
2177 if (IsEnvDefined(env, prefix, "PheErrLimit", print))
2178 {
2179 SetPheErrLimit(GetEnvValue(env, prefix, "PheErrLimit", fPheErrLimit));
2180 rc = kTRUE;
2181 }
2182 if (IsEnvDefined(env, prefix, "CheckDeadPixels", print))
2183 {
2184 SetCheckDeadPixels(GetEnvValue(env, prefix, "CheckDeadPixels", IsCheckDeadPixels()));
2185 rc = kTRUE;
2186 }
2187 if (IsEnvDefined(env, prefix, "CheckDeviatingBehavior", print))
2188 {
2189 SetCheckDeviatingBehavior(GetEnvValue(env, prefix, "CheckDeviatingBehavior", IsCheckDeviatingBehavior()));
2190 rc = kTRUE;
2191 }
2192 if (IsEnvDefined(env, prefix, "CheckExtractionWindow", print))
2193 {
2194 SetCheckExtractionWindow(GetEnvValue(env, prefix, "CheckExtractionWindow", IsCheckExtractionWindow()));
2195 rc = kTRUE;
2196 }
2197 if (IsEnvDefined(env, prefix, "CheckHistOverflow", print))
2198 {
2199 SetCheckHistOverflow(GetEnvValue(env, prefix, "CheckHistOverflow", IsCheckHistOverflow()));
2200 rc = kTRUE;
2201 }
2202 if (IsEnvDefined(env, prefix, "CheckOscillations", print))
2203 {
2204 SetCheckOscillations(GetEnvValue(env, prefix, "CheckOscillations", IsCheckOscillations()));
2205 rc = kTRUE;
2206 }
2207
2208 return rc;
2209}
Note: See TracBrowser for help on using the repository browser.