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

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