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

Last change on this file since 4778 was 4774, checked in by gaug, 20 years ago
*** empty log message ***
File size: 69.4 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 "MRawRunHeader.h"
197#include "MRawEvtPixelIter.h"
198
199#include "MGeomCam.h"
200#include "MGeomPix.h"
201#include "MHCamera.h"
202
203#include "MPedestalCam.h"
204#include "MPedestalPix.h"
205
206#include "MCalibrationChargeCam.h"
207#include "MCalibrationChargePix.h"
208#include "MCalibrationChargePINDiode.h"
209#include "MCalibrationChargeBlindPix.h"
210#include "MCalibrationChargeBlindCam.h"
211
212#include "MExtractedSignalCam.h"
213#include "MExtractedSignalPix.h"
214#include "MExtractedSignalBlindPixel.h"
215#include "MExtractedSignalPINDiode.h"
216
217#include "MBadPixelsCam.h"
218#include "MBadPixelsPix.h"
219
220#include "MCalibrationQECam.h"
221#include "MCalibrationQEPix.h"
222
223#include "MCalibrationCam.h"
224
225ClassImp(MCalibrationChargeCalc);
226
227using namespace std;
228
229const Float_t MCalibrationChargeCalc::fgChargeLimit = 2.5;
230const Float_t MCalibrationChargeCalc::fgChargeErrLimit = 0.;
231const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit = 1.;
232const Float_t MCalibrationChargeCalc::fgLambdaErrLimit = 0.2;
233const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit = 0.5;
234const Float_t MCalibrationChargeCalc::fgPheErrLimit = 4.5;
235const Float_t MCalibrationChargeCalc::fgFFactorErrLimit = 4.5;
236// --------------------------------------------------------------------------
237//
238// Default constructor.
239//
240// Sets the pointer to fQECam and fGeom to NULL
241//
242// Calls AddToBranchList for:
243// - MRawEvtData.fHiGainPixId
244// - MRawEvtData.fLoGainPixId
245// - MRawEvtData.fHiGainFadcSamples
246// - MRawEvtData.fLoGainFadcSamples
247//
248// Initializes:
249// - fChargeLimit to fgChargeLimit
250// - fChargeErrLimit to fgChargeErrLimit
251// - fChargeRelErrLimit to fgChargeRelErrLimit
252// - fFFactorErrLimit to fgFFactorErrLimit
253// - fLambdaCheckLimit to fgLambdaCheckLimit
254// - fLambdaErrLimit to fgLambdaErrLimit
255// - fPheErrLimit to fgPheErrLimit
256// - fPulserColor to MCalibrationCam::kCT1
257// - fOutputPath to "."
258// - fOutputFile to "ChargeCalibStat.txt"
259// - flag debug to kFALSE
260//
261// Calls:
262// - Clear()
263//
264MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
265 : fQECam(NULL), fGeom(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 Clear();
288
289}
290
291// --------------------------------------------------------------------------
292//
293// Sets:
294// - all variables to 0.,
295// - all flags to kFALSE
296// - all pointers to NULL
297// - the pulser colour to kNONE
298// - fBlindPixelFlags to 0
299// - fPINDiodeFlags to 0
300//
301void MCalibrationChargeCalc::Clear(const Option_t *o)
302{
303
304 fNumHiGainSamples = 0.;
305 fNumLoGainSamples = 0.;
306 fSqrtHiGainSamples = 0.;
307 fSqrtLoGainSamples = 0.;
308 fNumInnerFFactorMethodUsed = 0;
309
310 fBadPixels = NULL;
311 fCam = NULL;
312 fBlindPixel = NULL;
313 fBlindCam = NULL;
314 fPINDiode = NULL;
315 fPedestals = NULL;
316
317 SetPulserColor ( MCalibrationCam::kNONE );
318
319 fBlindPixelFlags.Set(0);
320 fPINDiodeFlags .Set(0);
321 fResultFlags .Set(0);
322}
323
324
325// -----------------------------------------------------------------------------------
326//
327// The following container are searched for and execution aborted if not in MParList:
328// - MPedestalCam
329//
330// The following containers are searched and created if they were not found:
331//
332// - MCalibrationQECam
333// - MBadPixelsCam
334//
335Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
336{
337
338 //
339 // Containers that have to be there.
340 //
341 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
342 if (!fPedestals)
343 {
344 *fLog << err << "MPedestalCam not found... aborting" << endl;
345 return kFALSE;
346 }
347
348 //
349 // Containers that are created in case that they are not there.
350 //
351 fQECam = (MCalibrationQECam*)pList->FindCreateObj("MCalibrationQECam");
352 if (!fQECam)
353 {
354 *fLog << err << "Cannot find nor create MCalibrationQECam... aborting" << endl;
355 return kFALSE;
356 }
357
358 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
359 if (!fBadPixels)
360 {
361 *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
362 return kFALSE;
363 }
364
365
366 //
367 // Check the pulser colour --> FIXME: this solution is only valid until the arrival of the DM's
368 //
369 if (fPulserColor == MCalibrationCam::kNONE)
370 {
371 *fLog << endl;
372 *fLog << err << GetDescriptor()
373 << ": No Pulser colour has been chosen. Since the installation of the IFAE pulser box,"
374 << " you HAVE to provide the LEDs colour, otherwise there is no calibration. " << endl;
375 *fLog << "See e.g. the macro calibration.C " << endl;
376 return kFALSE;
377 }
378
379 return kTRUE;
380}
381
382
383// --------------------------------------------------------------------------
384//
385// Search for the following input containers and abort if not existing:
386// - MGeomCam
387// - MCalibrationChargeCam
388//
389// Search for the following input containers and give a warning if not existing:
390// - MCalibrationChargeBlindPix
391// - MCalibrationChargePINDiode
392//
393// It retrieves the following variables from MCalibrationChargeCam:
394//
395// - fNumHiGainSamples
396// - fNumLoGainSamples
397//
398// It defines the PixId of every pixel in:
399//
400// - MCalibrationChargeCam
401// - MCalibrationQECam
402//
403// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
404//
405// - MCalibrationChargePix
406// - MCalibrationQEPix
407//
408// Sets the pulser colour and tests if it has not changed w.r.t. fPulserColor in:
409//
410// - MCalibrationChargeCam
411// - MCalibrationChargeBlindPix (if existing)
412// - MCalibrationChargePINDiode (if existing)
413//
414Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
415{
416
417 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
418 if (!fGeom)
419 {
420 *fLog << err << "No MGeomCam found... aborting." << endl;
421 return kFALSE;
422 }
423
424 fCam = (MCalibrationChargeCam*)pList->FindObject("MCalibrationChargeCam");
425 if (!fCam)
426 {
427 *fLog << err << "Cannot find MCalibrationChargeCam... aborting" << endl;
428 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationChargeCam before..." << endl;
429 return kFALSE;
430 }
431
432 //
433 // Optional Containers
434 //
435 fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindObject("MCalibrationChargeBlindPix");
436 if (!fBlindPixel)
437 {
438 fBlindCam = (MCalibrationChargeBlindCam*)pList->FindObject("MCalibrationChargeBlindCam");
439 if (!fBlindCam)
440 {
441 *fLog << endl;
442 *fLog << warn << GetDescriptor()
443 << ": MCalibrationChargeBlindPix nor MCalibrationChargeBlindCam "
444 << " found... no Blind Pixel method! " << endl;
445 }
446 }
447
448 fPINDiode = (MCalibrationChargePINDiode*)pList->FindObject("MCalibrationChargePINDiode");
449 if (!fPINDiode)
450 {
451 *fLog << endl;
452 *fLog << warn << GetDescriptor()
453 << ": MCalibrationChargePINDiode not found... no PIN Diode method! " << endl;
454 }
455
456
457 //
458 // Initialize the pulser colours
459 //
460 if (fCam->GetPulserColor() == MCalibrationCam::kNONE)
461 {
462 fCam->SetPulserColor( fPulserColor );
463
464 if (fBlindPixel)
465 fBlindPixel->SetColor( fPulserColor );
466
467 if (fBlindCam)
468 fBlindCam->SetColor( fPulserColor );
469
470 if (fPINDiode)
471 fPINDiode->SetColor( fPulserColor );
472 }
473
474 if (fPulserColor != fCam->GetPulserColor())
475 {
476 *fLog << err << GetDescriptor()
477 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeCam" << endl;
478 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
479 return kFALSE;
480 }
481
482 if (fBlindPixel)
483 if (fPulserColor != fBlindPixel->GetColor())
484 {
485 *fLog << err << GetDescriptor()
486 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeBlindPix." << endl;
487 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
488 return kFALSE;
489 }
490
491 if (fBlindCam)
492 if (fPulserColor != fBlindCam->GetColor())
493 {
494 *fLog << err << GetDescriptor()
495 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeBlindCam." << endl;
496 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
497 return kFALSE;
498 }
499
500 if (fPINDiode)
501 if (fPulserColor != fPINDiode->GetColor())
502 {
503 *fLog << err << GetDescriptor()
504 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargePINDiode." << endl;
505 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
506 return kFALSE;
507 }
508
509
510 fNumHiGainSamples = fCam->GetNumHiGainFADCSlices();
511 fNumLoGainSamples = fCam->GetNumLoGainFADCSlices();
512
513 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
514 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
515
516 UInt_t npixels = fGeom->GetNumPixels();
517
518 for (UInt_t i=0; i<npixels; i++)
519 {
520
521 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam) [i];
522 MCalibrationQEPix &pqe = (MCalibrationQEPix&) (*fQECam)[i];
523 MBadPixelsPix &bad = (*fBadPixels)[i];
524
525 pix.SetPixId(i);
526 pqe.SetPixId(i);
527
528 if (bad.IsBad())
529 {
530 pix.SetExcluded();
531 pqe.SetExcluded();
532 continue;
533 }
534
535 if (IsDebug())
536 pix.SetDebug();
537 }
538
539 return kTRUE;
540}
541
542// ----------------------------------------------------------------------------------
543//
544// Nothing to be done in Process, but have a look at MHCalibrationChargeCam, instead
545//
546Int_t MCalibrationChargeCalc::Process()
547{
548 return kTRUE;
549}
550
551// -----------------------------------------------------------------------
552//
553// Return if number of executions is null.
554//
555// First loop over pixels, average areas and sectors, call:
556// - FinalizePedestals()
557// - FinalizeCharges()
558// for every entry. Count number of valid pixels in loop and return kFALSE
559// if there are none (the "Michele check").
560//
561// Call FinalizeBadPixels()
562//
563// Call FinalizeFFactorMethod() (second and third loop over pixels and areas)
564//
565// Call FinalizeBlindPixel()
566// Call FinalizeBlindCam()
567// Call FinalizePINDiode()
568//
569// Call FinalizeFFactorQECam() (fourth loop over pixels and areas)
570// Call FinalizeBlindPixelQECam() (fifth loop over pixels and areas)
571// Call FinalizePINDiodeQECam() (sixth loop over pixels and areas)
572//
573// Call FinalizeUnsuitablePixels()
574//
575// Call MParContainer::SetReadyToSave() for fCam, fQECam, fBadPixels and
576// fBlindPixel and fPINDiode if they exist
577//
578// Print out some statistics
579//
580Int_t MCalibrationChargeCalc::PostProcess()
581{
582
583 if (GetNumExecutions()==0)
584 return kFALSE;
585
586 *fLog << endl;
587
588 if (fPINDiode)
589 if (!fPINDiode->IsValid())
590 {
591 *fLog << warn << GetDescriptor()
592 << ": MCalibrationChargePINDiode is declared not valid... no PIN Diode method! " << endl;
593 fPINDiode = NULL;
594 }
595
596 if (fBlindPixel)
597 if (!fBlindPixel->IsValid())
598 {
599 *fLog << warn << GetDescriptor()
600 << ": MCalibrationChargeBlindPix is declared not valid... no Blind Pixel method! " << endl;
601 fBlindPixel = NULL;
602 }
603
604 *fLog << endl;
605 //
606 // First loop over pixels, call FinalizePedestals and FinalizeCharges
607 //
608 Int_t nvalid = 0;
609
610 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
611 {
612
613 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[pixid];
614 //
615 // Check if the pixel has been excluded from the fits
616 //
617 if (pix.IsExcluded())
618 continue;
619
620 MPedestalPix &ped = (*fPedestals)[pixid];
621 MBadPixelsPix &bad = (*fBadPixels)[pixid];
622 const Int_t aidx = (*fGeom)[pixid].GetAidx();
623
624 FinalizePedestals(ped,pix,aidx);
625
626 if (FinalizeCharges(pix,bad,"pixel "))
627 nvalid++;
628 }
629
630 *fLog << endl;
631 //
632 // The Michele check ...
633 //
634 if (nvalid == 0)
635 {
636 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
637 << "Did you forget to fill the histograms "
638 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
639 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
640 << "instead of a calibration run " << endl;
641 return kFALSE;
642 }
643
644 for (UInt_t aidx=0; aidx<fGeom->GetNumAreas(); aidx++)
645 {
646
647 const MPedestalPix &ped = fPedestals->GetAverageArea(aidx);
648 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
649
650 FinalizePedestals(ped,pix,aidx);
651 FinalizeCharges(pix, fCam->GetAverageBadArea(aidx),"area id");
652 }
653
654 *fLog << endl;
655
656 for (UInt_t sector=0; sector<fGeom->GetNumSectors(); sector++)
657 {
658
659 const MPedestalPix &ped = fPedestals->GetAverageSector(sector);
660
661 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
662 FinalizePedestals(ped,pix, 0);
663 }
664
665 *fLog << endl;
666
667 //
668 // Finalize Bad Pixels
669 //
670 FinalizeBadPixels();
671
672 //
673 // Finalize F-Factor method
674 //
675 if (!FinalizeFFactorMethod())
676 {
677 *fLog << warn << "Could not calculate the photons flux from the F-Factor method " << endl;
678 fCam->SetFFactorMethodValid(kFALSE);
679 return kFALSE;
680 }
681 else
682 fCam->SetFFactorMethodValid(kTRUE);
683
684 *fLog << endl;
685 //
686 // Finalize Blind Pixel
687 //
688 if (fBlindPixel)
689 if (FinalizeBlindPixel())
690 fQECam->SetBlindPixelMethodValid(kTRUE);
691 else
692 fQECam->SetBlindPixelMethodValid(kFALSE);
693 else
694 if (FinalizeBlindCam())
695 fQECam->SetBlindPixelMethodValid(kTRUE);
696 else
697 fQECam->SetBlindPixelMethodValid(kFALSE);
698
699 //
700 // Finalize PIN Diode
701 //
702 if (FinalizePINDiode())
703 fQECam->SetPINDiodeMethodValid(kTRUE);
704 else
705 fQECam->SetPINDiodeMethodValid(kFALSE);
706
707 //
708 // Finalize QE Cam
709 //
710 FinalizeFFactorQECam();
711 FinalizeBlindPixelQECam();
712 FinalizePINDiodeQECam();
713
714 //
715 // Re-direct the output to an ascii-file from now on:
716 //
717 MLog asciilog;
718 asciilog.SetOutputFile(GetOutputFile(),kTRUE);
719 SetLogStream(&asciilog);
720 //
721 // Finalize calibration statistics
722 //
723 FinalizeUnsuitablePixels();
724
725 fCam ->SetReadyToSave();
726 fQECam ->SetReadyToSave();
727 fBadPixels->SetReadyToSave();
728
729 if (fBlindPixel)
730 fBlindPixel->SetReadyToSave();
731 if (fBlindCam)
732 fBlindCam->SetReadyToSave();
733 if (fPINDiode)
734 fPINDiode->SetReadyToSave();
735
736 *fLog << inf << endl;
737 *fLog << GetDescriptor() << ": Fatal errors statistics:" << endl;
738
739 PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal,
740 Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS: "));
741 PrintUncalibrated(MBadPixelsPix::kChargeRelErrNotValid,
742 Form("%s%2.1f%s","Signal Error bigger than ",fChargeRelErrLimit," times Mean Signal: "));
743 PrintUncalibrated(MBadPixelsPix::kLoGainSaturation,
744 "Pixels with Low Gain Saturation: ");
745 PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin,
746 Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s): "));
747 PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins,
748 Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s): "));
749 PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes,
750 "Pixels with deviating number of phes: ");
751 PrintUncalibrated(MBadPixelsPix::kDeviatingFFactor,
752 "Pixels with deviating F-Factor: ");
753
754 *fLog << inf << endl;
755 *fLog << GetDescriptor() << ": Unreliable errors statistics:" << endl;
756
757 PrintUncalibrated(MBadPixelsPix::kChargeSigmaNotValid,
758 "Signal Sigma smaller than Pedestal RMS: ");
759 PrintUncalibrated(MBadPixelsPix::kHiGainOscillating,
760 "Pixels with changing Hi Gain signal over time: ");
761 PrintUncalibrated(MBadPixelsPix::kLoGainOscillating,
762 "Pixels with changing Lo Gain signal over time: ");
763 PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted,
764 "Pixels with unsuccesful Gauss fit to the Hi Gain: ");
765 PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted,
766 "Pixels with unsuccesful Gauss fit to the Lo Gain: ");
767
768 SetLogStream(&gLog);
769
770 return kTRUE;
771}
772
773// ----------------------------------------------------------------------------------
774//
775// Retrieves pedestal and pedestal RMS from MPedestalPix
776// Retrieves total entries from MPedestalCam
777// Sets pedestal*fNumHiGainSamples and pedestal*fNumLoGainSamples in MCalibrationChargePix
778// Sets pedRMS *fSqrtHiGainSamples and pedRMS *fSqrtLoGainSamples in MCalibrationChargePix
779//
780// If the flag MCalibrationPix::IsHiGainSaturation() is set, call also:
781// - MCalibrationChargePix::CalcLoGainPedestal()
782//
783void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal, const Int_t aidx)
784{
785
786 //
787 // get the pedestals
788 //
789 const Float_t pedes = ped.GetPedestal();
790 const Float_t prms = ped.GetPedestalRms();
791 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
792
793
794 //
795 // set them in the calibration camera
796 //
797 if (cal.IsHiGainSaturation())
798 {
799 cal.SetPedestal(pedes* fNumLoGainSamples,
800 prms * fSqrtLoGainSamples,
801 prms * fNumLoGainSamples / num);
802 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples, aidx);
803 }
804 else
805 {
806 cal.SetPedestal(pedes* fNumHiGainSamples,
807 prms * fSqrtHiGainSamples,
808 prms * fNumHiGainSamples / num);
809 }
810
811}
812
813// ----------------------------------------------------------------------------------------------------
814//
815// Check fit results validity. Bad Pixels flags are set if:
816//
817// 1) Pixel has a mean smaller than fChargeLimit*PedRMS ( Flag: MBadPixelsPix::kChargeIsPedestal)
818// 2) Pixel has a mean error smaller than fChargeErrLimit ( Flag: MBadPixelsPix::kChargeErrNotValid)
819// 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error
820// ( Flag: MBadPixelsPix::kChargeRelErrNotValid)
821// 4) Pixel has a sigma bigger than its Pedestal RMS ( Flag: MBadPixelsPix::kChargeSigmaNotValid )
822//
823// Further returns if flags: MBadPixelsPix::kUnsuitableRun is set
824//
825// Calls MCalibrationChargePix::CalcReducedSigma() and sets flag: MBadPixelsPix::kChargeIsPedestal
826// and returns kFALSE if not succesful.
827//
828// Calls MCalibrationChargePix::CalcFFactor() and sets flag: MBadPixelsPix::kDeviatingNumPhes)
829// and returns kFALSE if not succesful.
830//
831// Calls MCalibrationChargePix::CalcConvFFactor()and sets flag: MBadPixelsPix::kDeviatingNumPhes)
832// and returns kFALSE if not succesful.
833//
834Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what)
835{
836
837 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
838 return kFALSE;
839
840 if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
841 {
842 *fLog << warn << GetDescriptor()
843 << Form(": Fitted Charge: %5.2f is smaller than %2.1f",cal.GetMean(),fChargeLimit)
844 << Form(" Pedestal RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
845 bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
846 }
847
848 if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr())
849 {
850 *fLog << warn << GetDescriptor()
851 << Form(": Fitted Charge: %4.2f is smaller than %2.1f",cal.GetMean(),fChargeRelErrLimit)
852 << Form(" times its error: %4.2f in %s%4i",cal.GetMeanErr(),what,cal.GetPixId()) << endl;
853 bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
854 }
855
856 if (cal.GetSigma() < cal.GetPedRms())
857 {
858 *fLog << warn << GetDescriptor()
859 << Form(": Sigma of Fitted Charge: %6.2f is smaller than",cal.GetSigma())
860 << Form(" Ped. RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
861 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
862 return kFALSE;
863 }
864
865 if (!cal.CalcReducedSigma())
866 {
867 *fLog << warn << GetDescriptor()
868 << Form(": Could not calculate the reduced sigma in %s: ",what)
869 << Form(" %4i",cal.GetPixId())
870 << endl;
871 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
872 return kFALSE;
873 }
874
875 if (!cal.CalcFFactor())
876 {
877 *fLog << warn << GetDescriptor()
878 << Form(": Could not calculate the F-Factor in %s: ",what)
879 << Form(" %4i",cal.GetPixId())
880 << endl;
881 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
882 return kFALSE;
883 }
884
885 if (!cal.CalcConvFFactor())
886 {
887 *fLog << warn << GetDescriptor()
888 << Form(": Could not calculate the Conv. FADC counts to Phes in %s: ",what)
889 << Form(" %4i",cal.GetPixId())
890 << endl;
891 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
892 return kFALSE;
893 }
894
895 return kTRUE;
896}
897
898
899
900// -----------------------------------------------------------------------------------
901//
902// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
903// - MBadPixelsPix::kChargeIsPedestal
904// - MBadPixelsPix::kChargeRelErrNotValid
905// - MBadPixelsPix::kMeanTimeInFirstBin
906// - MBadPixelsPix::kMeanTimeInLast2Bins
907// - MBadPixelsPix::kDeviatingNumPhes
908//
909// - Call MCalibrationPix::SetExcluded() for the bad pixels
910//
911// Sets pixel to MBadPixelsPix::kUnreliableRun, if one of the following flags is set:
912// - MBadPixelsPix::kChargeSigmaNotValid
913//
914void MCalibrationChargeCalc::FinalizeBadPixels()
915{
916
917 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
918 {
919
920 MBadPixelsPix &bad = (*fBadPixels)[i];
921 MCalibrationPix &pix = (*fCam)[i];
922
923 if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal))
924 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
925
926 if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid ))
927 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
928
929 if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ))
930 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
931
932 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ))
933 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
934
935 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ))
936 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
937
938 if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes ))
939 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
940
941 if (bad.IsUnsuitable( MBadPixelsPix::kUnsuitableRun ))
942 pix.SetExcluded();
943
944 if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ))
945 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
946 }
947}
948
949// ------------------------------------------------------------------------
950//
951//
952// First loop: Calculate a mean and mean RMS of photo-electrons per area index
953// Include only pixels which are not MBadPixelsPix::kUnsuitableRun nor
954// MBadPixelsPix::kChargeSigmaNotValid (see FinalizeBadPixels()) and set
955// MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case.
956//
957// Second loop: Get mean number of photo-electrons and its RMS including
958// only pixels with flag MCalibrationChargePix::IsFFactorMethodValid()
959// and further exclude those deviating by more than fPheErrLimit mean
960// sigmas from the mean (obtained in first loop). Set
961// MBadPixelsPix::kDeviatingNumPhes if excluded.
962//
963// For the suitable pixels with flag MBadPixelsPix::kChargeSigmaNotValid
964// set the number of photo-electrons as the mean number of photo-electrons
965// calculated in that area index.
966//
967// Set weighted mean and variance of photo-electrons per area index in:
968// average area pixels of MCalibrationChargeCam (obtained from:
969// MCalibrationChargeCam::GetAverageArea() )
970//
971// Set weighted mean and variance of photo-electrons per sector in:
972// average sector pixels of MCalibrationChargeCam (obtained from:
973// MCalibrationChargeCam::GetAverageSector() )
974//
975//
976// Third loop: Set mean number of photo-electrons and its RMS in the pixels
977// only excluded as: MBadPixelsPix::kChargeSigmaNotValid
978//
979Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod()
980{
981
982 const UInt_t npixels = fGeom->GetNumPixels();
983 const UInt_t nareas = fGeom->GetNumAreas();
984 const UInt_t nsectors = fGeom->GetNumSectors();
985
986 Float_t lowlim [nareas];
987 Float_t upplim [nareas];
988 Double_t areavars [nareas];
989 Double_t areaweights[nareas], sectorweights [nsectors];
990 Double_t areaphes [nareas], sectorphes [nsectors];
991 Int_t numareavalid[nareas], numsectorvalid[nsectors];
992
993 memset(lowlim ,0, nareas * sizeof(Float_t));
994 memset(upplim ,0, nareas * sizeof(Float_t));
995 memset(areaphes ,0, nareas * sizeof(Double_t));
996 memset(areavars ,0, nareas * sizeof(Double_t));
997 memset(areaweights ,0, nareas * sizeof(Double_t));
998 memset(numareavalid ,0, nareas * sizeof(Int_t ));
999 memset(sectorweights ,0, nsectors * sizeof(Double_t));
1000 memset(sectorphes ,0, nsectors * sizeof(Double_t));
1001 memset(numsectorvalid,0, nsectors * sizeof(Int_t ));
1002
1003 //
1004 // First loop: Get mean number of photo-electrons and the RMS
1005 // The loop is only to recognize later pixels with very deviating numbers
1006 //
1007 MHCamera camphes(*fGeom,"Camphes","Phes in Camera");
1008
1009 for (UInt_t i=0; i<npixels; i++)
1010 {
1011
1012 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam) [i];
1013 MBadPixelsPix &bad = (*fBadPixels)[i];
1014
1015 if (!pix.IsFFactorMethodValid())
1016 continue;
1017
1018 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1019 {
1020 pix.SetFFactorMethodValid(kFALSE);
1021 continue;
1022 }
1023
1024 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1025 continue;
1026
1027 const Float_t nphe = pix.GetPheFFactorMethod();
1028 const Int_t aidx = (*fGeom)[i].GetAidx();
1029
1030 camphes.Fill(i,nphe);
1031 camphes.SetUsed(i);
1032
1033 areaphes [aidx] += nphe;
1034 areavars [aidx] += nphe*nphe;
1035 numareavalid[aidx] ++;
1036 }
1037
1038 for (UInt_t i=0; i<nareas; i++)
1039 {
1040 if (numareavalid[i] == 0)
1041 {
1042 *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
1043 << "in area index: " << i << endl;
1044 continue;
1045 }
1046
1047 if (numareavalid[i] == 1)
1048 areavars[i] = 0.;
1049 else if (numareavalid[i] == 0)
1050 {
1051 areaphes[i] = -1.;
1052 areaweights[i] = -1.;
1053 }
1054 else
1055 {
1056 areavars[i] = (areavars[i] - areaphes[i]*areaphes[i]/numareavalid[i]) / (numareavalid[i]-1);
1057 areaphes[i] = areaphes[i] / numareavalid[i];
1058 }
1059
1060 if (areavars[i] < 0.)
1061 {
1062 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of photo-electrons found "
1063 << "in area index: " << i << endl;
1064 continue;
1065 }
1066
1067 lowlim [i] = areaphes[i] - fPheErrLimit*TMath::Sqrt(areavars[i]);
1068 upplim [i] = areaphes[i] + fPheErrLimit*TMath::Sqrt(areavars[i]);
1069
1070 TArrayI area(1);
1071 area[0] = i;
1072
1073 TH1D *hist = camphes.ProjectionS(TArrayI(),area,"_py",100);
1074 hist->Fit("gaus","Q");
1075 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1076 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1077 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1078
1079 if (IsDebug())
1080 camphes.DrawClone();
1081
1082 if (ndf < 2)
1083 {
1084 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1085 << "in the camera with area index: " << i << endl;
1086 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1087 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1088 delete hist;
1089 continue;
1090 }
1091
1092 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1093
1094 if (prob < 0.001)
1095 {
1096 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1097 << "in the camera with area index: " << i << endl;
1098 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1099 << " is smaller than 0.001 " << endl;
1100 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1101 delete hist;
1102 continue;
1103 }
1104
1105 *fLog << inf << GetDescriptor() << ": Mean number of photo-electrons "
1106 << "with area idx " << i << ": "
1107 << Form("%7.2f+-%6.2f",mean,sigma) << endl;
1108
1109 lowlim [i] = mean - fPheErrLimit*sigma;
1110 upplim [i] = mean + fPheErrLimit*sigma;
1111
1112 delete hist;
1113 }
1114
1115 *fLog << endl;
1116
1117 memset(numareavalid,0,nareas*sizeof(Int_t));
1118 memset(areaphes ,0,nareas*sizeof(Double_t));
1119 memset(areavars ,0,nareas*sizeof(Double_t));
1120
1121 //
1122 // Second loop: Get mean number of photo-electrons and its RMS excluding
1123 // pixels deviating by more than fPheErrLimit sigma.
1124 // Set the conversion factor FADC counts to photo-electrons
1125 //
1126 for (UInt_t i=0; i<npixels; i++)
1127 {
1128
1129 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1130
1131 if (!pix.IsFFactorMethodValid())
1132 continue;
1133
1134 MBadPixelsPix &bad = (*fBadPixels)[i];
1135
1136 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1137 continue;
1138
1139 const Float_t nvar = pix.GetPheFFactorMethodVar();
1140
1141 if (nvar <= 0.)
1142 {
1143 pix.SetFFactorMethodValid(kFALSE);
1144 continue;
1145 }
1146
1147 const Int_t aidx = (*fGeom)[i].GetAidx();
1148 const Int_t sector = (*fGeom)[i].GetSector();
1149 const Float_t area = (*fGeom)[i].GetA();
1150 const Float_t nphe = pix.GetPheFFactorMethod();
1151
1152 if ( nphe < lowlim[aidx] || nphe > upplim[aidx] )
1153 {
1154 *fLog << warn << GetDescriptor() << ": Number of phes: "
1155 << Form("%7.2f out of %3.1f sigma limit: ",nphe,fPheErrLimit)
1156 << Form("[%7.2f,%7.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
1157 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1158 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1159 pix.SetFFactorMethodValid(kFALSE);
1160 continue;
1161 }
1162
1163 areaweights [aidx] += nphe*nphe;
1164 areaphes [aidx] += nphe;
1165 numareavalid [aidx] ++;
1166
1167 if (aidx == 0)
1168 fNumInnerFFactorMethodUsed++;
1169
1170 sectorweights [sector] += nphe*nphe/area/area;
1171 sectorphes [sector] += nphe/area;
1172 numsectorvalid[sector] ++;
1173 }
1174
1175 *fLog << endl;
1176
1177 for (UInt_t aidx=0; aidx<nareas; aidx++)
1178 {
1179
1180 MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
1181
1182 if (numareavalid[aidx] == 1)
1183 areaweights[aidx] = 0.;
1184 else if (numareavalid[aidx] == 0)
1185 {
1186 areaphes[aidx] = -1.;
1187 areaweights[aidx] = -1.;
1188 }
1189 else
1190 {
1191 areaweights[aidx] = (areaweights[aidx] - areaphes[aidx]*areaphes[aidx]/numareavalid[aidx])
1192 / (numareavalid[aidx]-1);
1193 areaphes[aidx] /= numareavalid[aidx];
1194 }
1195
1196 if (areaweights[aidx] < 0. || areaphes[aidx] <= 0.)
1197 {
1198 *fLog << warn << GetDescriptor()
1199 << ": Mean number phes from area index " << aidx << " could not be calculated: "
1200 << " Mean: " << areaphes[aidx]
1201 << " Variance: " << areaweights[aidx] << endl;
1202 apix.SetFFactorMethodValid(kFALSE);
1203 continue;
1204 }
1205
1206 *fLog << inf << GetDescriptor()
1207 << ": Average total number phes in area idx " << aidx << ": "
1208 << Form("%7.2f%s%6.2f",areaphes[aidx]," +- ",TMath::Sqrt(areaweights[aidx])) << endl;
1209
1210 apix.SetPheFFactorMethod ( areaphes[aidx] );
1211 apix.SetPheFFactorMethodVar( areaweights[aidx] / numareavalid[aidx] );
1212 apix.SetFFactorMethodValid ( kTRUE );
1213
1214 }
1215
1216 *fLog << endl;
1217
1218 for (UInt_t sector=0; sector<nsectors; sector++)
1219 {
1220
1221 if (numsectorvalid[sector] == 1)
1222 sectorweights[sector] = 0.;
1223 else if (numsectorvalid[sector] == 0)
1224 {
1225 sectorphes[sector] = -1.;
1226 sectorweights[sector] = -1.;
1227 }
1228 else
1229 {
1230 sectorweights[sector] = (sectorweights[sector]
1231 - sectorphes[sector]*sectorphes[sector]/numsectorvalid[sector]
1232 )
1233 / (numsectorvalid[sector]-1.);
1234 sectorphes[sector] /= numsectorvalid[sector];
1235 }
1236
1237 MCalibrationChargePix &spix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
1238
1239 if (sectorweights[sector] < 0. || sectorphes[sector] <= 0.)
1240 {
1241 *fLog << warn << GetDescriptor()
1242 <<": Mean number phes per area for sector " << sector << " could not be calculated: "
1243 << " Mean: " << sectorphes[sector]
1244 << " Variance: " << sectorweights[sector] << endl;
1245 spix.SetFFactorMethodValid(kFALSE);
1246 continue;
1247 }
1248
1249 *fLog << inf << GetDescriptor()
1250 << ": Average number phes per area in sector " << sector << ": "
1251 << Form("%5.2f+-%4.2f [phe/mm^2]",sectorphes[sector],TMath::Sqrt(sectorweights[sector]))
1252 << endl;
1253
1254 spix.SetPheFFactorMethod ( sectorphes[sector] );
1255 spix.SetPheFFactorMethodVar( sectorweights[sector] / numsectorvalid[sector]);
1256 spix.SetFFactorMethodValid ( kTRUE );
1257
1258 }
1259
1260 //
1261 // Third loop: Set mean number of photo-electrons and its RMS in the pixels
1262 // only excluded as: MBadPixelsPix::kChargeSigmaNotValid
1263 //
1264 for (UInt_t i=0; i<npixels; i++)
1265 {
1266
1267 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1268 MBadPixelsPix &bad = (*fBadPixels)[i];
1269
1270 if (!pix.IsFFactorMethodValid())
1271 continue;
1272
1273 const Int_t aidx = (*fGeom)[i].GetAidx();
1274
1275 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1276 {
1277 MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
1278 pix.SetPheFFactorMethod ( apix.GetPheFFactorMethod() );
1279 pix.SetPheFFactorMethodVar( apix.GetPheFFactorMethodVar() );
1280 if (!pix.CalcConvFFactor())
1281 {
1282 *fLog << warn << GetDescriptor()
1283 << ": Could not calculate the Conv. FADC counts to Phes in pixel: "
1284 << Form(" %4i",pix.GetPixId())
1285 << endl;
1286 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1287 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1288 }
1289 }
1290 }
1291
1292 return kTRUE;
1293}
1294
1295
1296// ------------------------------------------------------------------------
1297//
1298// Returns kFALSE if pointer to MCalibrationChargeBlindPix is NULL
1299//
1300// The check returns kFALSE if:
1301//
1302// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
1303// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
1304//
1305// Calls:
1306// - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
1307//
1308Bool_t MCalibrationChargeCalc::FinalizeBlindPixel()
1309{
1310
1311 if (!fBlindPixel)
1312 return kFALSE;
1313
1314 const Float_t lambda = fBlindPixel->GetLambda();
1315 const Float_t lambdaerr = fBlindPixel->GetLambdaErr();
1316 const Float_t lambdacheck = fBlindPixel->GetLambdaCheck();
1317
1318 if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
1319 {
1320 *fLog << warn << GetDescriptor()
1321 << Form("%s%4.2f%s%4.2f%s%4.2f%s",": Lambda: ",lambda," and Lambda-Check: ",
1322 lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel ")
1323 << endl;
1324 return kFALSE;
1325 }
1326
1327 if (lambdaerr > fLambdaErrLimit)
1328 {
1329 *fLog << warn << GetDescriptor()
1330 << Form("%s%4.2f%s%4.2f%s",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
1331 fLambdaErrLimit," in Blind Pixel ") << endl;
1332 return kFALSE;
1333 }
1334
1335 if (!fBlindPixel->CalcFluxInsidePlexiglass())
1336 {
1337 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
1338 << "will skip Blind Pixel Calibration " << endl;
1339 return kFALSE;
1340 }
1341
1342 return kTRUE;
1343}
1344
1345// ------------------------------------------------------------------------
1346//
1347// Returns kFALSE if pointer to MCalibrationChargeBlindCam is NULL
1348//
1349// The check returns kFALSE if:
1350//
1351// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
1352// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
1353//
1354// Calls:
1355// - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
1356//
1357Bool_t MCalibrationChargeCalc::FinalizeBlindCam()
1358{
1359
1360 if (!fBlindCam)
1361 return kFALSE;
1362
1363 Float_t flux = 0.;
1364 Float_t fluxvar = 0.;
1365 Int_t nvalid = 0;
1366
1367 for (UInt_t i=0; i<fBlindCam->GetNumBlindPixels(); i++)
1368 {
1369
1370 MCalibrationChargeBlindPix &blindpix = (*fBlindCam)[i];
1371
1372 if (!blindpix.IsValid())
1373 continue;
1374
1375 const Float_t lambda = blindpix.GetLambda();
1376 const Float_t lambdaerr = blindpix.GetLambdaErr();
1377 const Float_t lambdacheck = blindpix.GetLambdaCheck();
1378
1379 if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
1380 {
1381 *fLog << warn << GetDescriptor()
1382 << Form("%s%4.2f%s%4.2f%s%4.2f%s%2i",": Lambda: ",lambda," and Lambda-Check: ",
1383 lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel Nr.",i)
1384 << endl;
1385 blindpix.SetValid(kFALSE);
1386 continue;
1387 }
1388
1389 if (lambdaerr > fLambdaErrLimit)
1390 {
1391 *fLog << warn << GetDescriptor()
1392 << Form("%s%4.2f%s%4.2f%s%2i",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
1393 fLambdaErrLimit," in Blind Pixel Nr.",i) << endl;
1394 blindpix.SetValid(kFALSE);
1395 continue;
1396 }
1397
1398 if (!blindpix.CalcFluxInsidePlexiglass())
1399 {
1400 *fLog << warn << "Could not calculate the flux of photons from Blind Pixel Nr." << i << endl;
1401 blindpix.SetValid(kFALSE);
1402 continue;
1403 }
1404
1405 nvalid++;
1406 const Float_t weight = 1./ blindpix.GetFluxInsidePlexiglassErr() / blindpix.GetFluxInsidePlexiglassErr();
1407 flux += weight * blindpix.GetFluxInsidePlexiglass();
1408 fluxvar += weight;
1409 }
1410
1411 if (!nvalid)
1412 return kFALSE;
1413
1414 flux /= fluxvar;
1415 fluxvar /= 1./fluxvar;
1416
1417 const Float_t photons = flux * (*fGeom)[0].GetA() / fQECam->GetPlexiglassQE();
1418 fCam->SetNumPhotonsBlindPixelMethod(photons);
1419
1420 const Float_t photrelvar = fluxvar / flux / flux + fQECam->GetPlexiglassQERelVar();
1421 if (photrelvar > 0.)
1422 fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1423
1424 return kTRUE;
1425}
1426
1427// ------------------------------------------------------------------------
1428//
1429// Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL
1430//
1431// The check returns kFALSE if:
1432//
1433// 1) PINDiode has a fitted charge smaller than fChargeLimit*PedRMS
1434// 2) PINDiode has a fit error smaller than fChargeErrLimit
1435// 3) PINDiode has a fitted charge smaller its fChargeRelErrLimit times its charge error
1436// 4) PINDiode has a charge sigma smaller than its Pedestal RMS
1437//
1438// Calls:
1439// - MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass()
1440//
1441Bool_t MCalibrationChargeCalc::FinalizePINDiode()
1442{
1443
1444 if (!fPINDiode)
1445 return kFALSE;
1446
1447 if (fPINDiode->GetMean() < fChargeLimit*fPINDiode->GetPedRms())
1448 {
1449 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1450 << fChargeLimit << " Pedestal RMS in PINDiode " << endl;
1451 return kFALSE;
1452 }
1453
1454 if (fPINDiode->GetMeanErr() < fChargeErrLimit)
1455 {
1456 *fLog << warn << GetDescriptor() << ": Error of Fitted Charge is smaller than "
1457 << fChargeErrLimit << " in PINDiode " << endl;
1458 return kFALSE;
1459 }
1460
1461 if (fPINDiode->GetMean() < fChargeRelErrLimit*fPINDiode->GetMeanErr())
1462 {
1463 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1464 << fChargeRelErrLimit << "* its error in PINDiode " << endl;
1465 return kFALSE;
1466 }
1467
1468 if (fPINDiode->GetSigma() < fPINDiode->GetPedRms())
1469 {
1470 *fLog << warn << GetDescriptor()
1471 << ": Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl;
1472 return kFALSE;
1473 }
1474
1475
1476 if (!fPINDiode->CalcFluxOutsidePlexiglass())
1477 {
1478 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
1479 << "will skip PIN Diode Calibration " << endl;
1480 return kFALSE;
1481 }
1482
1483 return kTRUE;
1484}
1485
1486// ------------------------------------------------------------------------
1487//
1488// Calculate the average number of photons outside the plexiglass with the
1489// formula:
1490//
1491// av.Num.photons(area index) = av.Num.Phes(area index)
1492// / MCalibrationQEPix::GetDefaultQE(fPulserColor)
1493// / MCalibrationQEPix::GetPMTCollectionEff()
1494// / MCalibrationQEPix::GetLightGuidesEff(fPulserColor)
1495// / MCalibrationQECam::GetPlexiglassQE()
1496//
1497// Calculate the variance on the average number of photons assuming that the error on the
1498// Quantum efficiency is reduced by the number of used inner pixels, but the rest of the
1499// values keeps it ordinary error since it is systematic.
1500//
1501// Loop over pixels:
1502//
1503// - Continue, if not MCalibrationChargePix::IsFFactorMethodValid() and set:
1504// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor)
1505//
1506// - Call MCalibrationChargePix::CalcMeanFFactor(av.Num.photons) and set:
1507// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) if not succesful
1508//
1509// - Calculate the quantum efficiency with the formula:
1510//
1511// QE = ( Num.Phes / av.Num.photons ) * MGeomCam::GetPixRatio()
1512//
1513// - Set QE in MCalibrationQEPix::SetQEFFactor ( QE, fPulserColor );
1514//
1515// - Set Variance of QE in MCalibrationQEPix::SetQEFFactorVar ( Variance, fPulserColor );
1516// - Set bit MCalibrationQEPix::SetFFactorMethodValid(kTRUE,fPulserColor)
1517//
1518// - Call MCalibrationQEPix::UpdateFFactorMethod()
1519//
1520void MCalibrationChargeCalc::FinalizeFFactorQECam()
1521{
1522
1523 if (fNumInnerFFactorMethodUsed < 2)
1524 {
1525 *fLog << warn << GetDescriptor()
1526 << ": Could not calculate F-Factor Method: Less than 2 inner pixels valid! " << endl;
1527 return;
1528 }
1529
1530 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCam->GetAverageArea(0);
1531 MCalibrationQEPix &qepix = (MCalibrationQEPix&) fQECam->GetAverageArea(0);
1532
1533 const Float_t avphotons = avpix.GetPheFFactorMethod()
1534 / qepix.GetDefaultQE(fPulserColor)
1535 / qepix.GetPMTCollectionEff()
1536 / qepix.GetLightGuidesEff(fPulserColor)
1537 / fQECam->GetPlexiglassQE();
1538
1539 const Float_t avphotrelvar = avpix.GetPheFFactorMethodRelVar()
1540 + qepix.GetDefaultQERelVar(fPulserColor) / fNumInnerFFactorMethodUsed
1541 + qepix.GetPMTCollectionEffRelVar()
1542 + qepix.GetLightGuidesEffRelVar(fPulserColor)
1543 + fQECam->GetPlexiglassQERelVar();
1544
1545 const UInt_t nareas = fGeom->GetNumAreas();
1546
1547 //
1548 // Set the results in the MCalibrationChargeCam
1549 //
1550 fCam->SetNumPhotonsFFactorMethod (avphotons);
1551 if (avphotrelvar > 0.)
1552 fCam->SetNumPhotonsFFactorMethodErr(TMath::Sqrt( avphotrelvar * avphotons * avphotons));
1553
1554 Float_t lowlim [nareas];
1555 Float_t upplim [nareas];
1556 Double_t avffactorphotons [nareas];
1557 Double_t avffactorphotvar [nareas];
1558 Int_t numffactor [nareas];
1559
1560 memset(lowlim ,0, nareas * sizeof(Float_t));
1561 memset(upplim ,0, nareas * sizeof(Float_t));
1562 memset(avffactorphotons,0, nareas * sizeof(Double_t));
1563 memset(avffactorphotvar,0, nareas * sizeof(Double_t));
1564 memset(numffactor ,0, nareas * sizeof(Int_t));
1565
1566 const UInt_t npixels = fGeom->GetNumPixels();
1567
1568 MHCamera camffactor(*fGeom,"Camffactor","F-Factor in Camera");
1569
1570 for (UInt_t i=0; i<npixels; i++)
1571 {
1572
1573 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1574 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1575 MBadPixelsPix &bad = (*fBadPixels)[i];
1576
1577 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1578 continue;
1579
1580 const Float_t photons = avphotons / fGeom->GetPixRatio(i);
1581 const Float_t qe = pix.GetPheFFactorMethod() / photons ;
1582
1583 const Float_t qerelvar = avphotrelvar + pix.GetPheFFactorMethodRelVar();
1584
1585 qepix.SetQEFFactor ( qe , fPulserColor );
1586 qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor );
1587 qepix.SetFFactorMethodValid( kTRUE , fPulserColor );
1588
1589 if (!qepix.UpdateFFactorMethod())
1590 *fLog << warn << GetDescriptor()
1591 << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl;
1592
1593 //
1594 // The following pixels are those with deviating sigma, but otherwise OK,
1595 // probably those with stars during the pedestal run, but not the cal. run.
1596 //
1597 if (!pix.CalcMeanFFactor( photons , avphotrelvar ))
1598 {
1599 (*fBadPixels)[i].SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1600 (*fBadPixels)[i].SetUnsuitable ( MBadPixelsPix::kUnreliableRun );
1601 continue;
1602 }
1603
1604 const Int_t aidx = (*fGeom)[i].GetAidx();
1605 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1606
1607 camffactor.Fill(i,ffactor);
1608 camffactor.SetUsed(i);
1609
1610 avffactorphotons[aidx] += ffactor;
1611 avffactorphotvar[aidx] += ffactor*ffactor;
1612 numffactor[aidx]++;
1613 }
1614
1615 for (UInt_t i=0; i<nareas; i++)
1616 {
1617
1618 if (numffactor[i] == 0)
1619 {
1620 *fLog << warn << GetDescriptor() << ": No pixels with valid total F-Factor found "
1621 << "in area index: " << i << endl;
1622 continue;
1623 }
1624
1625 avffactorphotvar[i] = (avffactorphotvar[i] - avffactorphotons[i]*avffactorphotons[i]/numffactor[i]) / (numffactor[i]-1.);
1626 avffactorphotons[i] = avffactorphotons[i] / numffactor[i];
1627
1628 if (avffactorphotvar[i] < 0.)
1629 {
1630 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of total F-Factor found "
1631 << "in area index: " << i << endl;
1632 continue;
1633 }
1634
1635 lowlim [i] = 1.1; // Lowest known F-Factor of a PMT
1636 upplim [i] = avffactorphotons[i] + fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);
1637
1638 TArrayI area(1);
1639 area[0] = i;
1640
1641 TH1D *hist = camffactor.ProjectionS(TArrayI(),area,"_py",100);
1642 hist->Fit("gaus","Q");
1643 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1644 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1645 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1646
1647 if (IsDebug())
1648 camffactor.DrawClone();
1649
1650 if (ndf < 2)
1651 {
1652 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1653 << "in the camera with area index: " << i << endl;
1654 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1655 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1656 delete hist;
1657 continue;
1658 }
1659
1660 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1661
1662 if (prob < 0.001)
1663 {
1664 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1665 << "in the camera with area index: " << i << endl;
1666 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1667 << " is smaller than 0.001 " << endl;
1668 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1669 delete hist;
1670 continue;
1671 }
1672
1673 *fLog << inf << GetDescriptor() << ": Mean F-Factor "
1674 << "with area index #" << i << ": "
1675 << Form("%4.2f+-%4.2f",mean,sigma) << endl;
1676
1677 lowlim [i] = 1.1;
1678 upplim [i] = mean + fFFactorErrLimit*sigma;
1679
1680 delete hist;
1681 }
1682
1683 *fLog << endl;
1684
1685 for (UInt_t i=0; i<npixels; i++)
1686 {
1687
1688 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1689 MBadPixelsPix &bad = (*fBadPixels)[i];
1690
1691 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1692 continue;
1693
1694 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1695 const Int_t aidx = (*fGeom)[i].GetAidx();
1696
1697 if ( ffactor < lowlim[aidx] || ffactor > upplim[aidx] )
1698 {
1699 *fLog << warn << GetDescriptor() << ": Overall F-Factor "
1700 << Form("%5.2f",ffactor) << " out of range ["
1701 << Form("%5.2f,%5.2f",lowlim[aidx],upplim[aidx]) << "] pixel " << i << endl;
1702
1703 bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1704 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1705 }
1706 }
1707
1708 for (UInt_t i=0; i<npixels; i++)
1709 {
1710
1711 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1712 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1713 MBadPixelsPix &bad = (*fBadPixels)[i];
1714
1715 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1716 {
1717 qepix.SetFFactorMethodValid(kFALSE,fPulserColor);
1718 pix.SetFFactorMethodValid(kFALSE);
1719 pix.SetExcluded();
1720 continue;
1721 }
1722 }
1723}
1724
1725
1726// ------------------------------------------------------------------------
1727//
1728// Loop over pixels:
1729//
1730// - Continue, if not MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() and set:
1731// MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor)
1732//
1733// - Calculate the quantum efficiency with the formula:
1734//
1735// QE = Num.Phes / MCalibrationChargeBlindPix::GetFluxInsidePlexiglass()
1736// / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE()
1737//
1738// - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor );
1739// - Set Variance of QE in MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor );
1740// - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor)
1741//
1742// - Call MCalibrationQEPix::UpdateBlindPixelMethod()
1743//
1744void MCalibrationChargeCalc::FinalizeBlindPixelQECam()
1745{
1746
1747 const UInt_t npixels = fGeom->GetNumPixels();
1748
1749 //
1750 // Set the results in the MCalibrationChargeCam
1751 //
1752 if (fBlindPixel)
1753 {
1754 if (fBlindPixel->IsFluxInsidePlexiglassAvailable())
1755 {
1756
1757 const Float_t photons = fBlindPixel->GetFluxInsidePlexiglass() * (*fGeom)[0].GetA()
1758 / fQECam->GetPlexiglassQE();
1759 fCam->SetNumPhotonsBlindPixelMethod(photons);
1760
1761 const Float_t photrelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1762 + fQECam->GetPlexiglassQERelVar();
1763 if (photrelvar > 0.)
1764 fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1765 }
1766 }
1767 //
1768 // With the knowledge of the overall photon flux, calculate the
1769 // quantum efficiencies after the Blind Pixel and PIN Diode method
1770 //
1771 for (UInt_t i=0; i<npixels; i++)
1772 {
1773
1774 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1775
1776 if (!fBlindPixel)
1777 {
1778 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1779 continue;
1780 }
1781
1782 if (!fBlindPixel->IsFluxInsidePlexiglassAvailable())
1783 {
1784 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1785 continue;
1786 }
1787
1788 MBadPixelsPix &bad = (*fBadPixels)[i];
1789 if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1790 {
1791 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1792 continue;
1793 }
1794
1795 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1796 MGeomPix &geo = (*fGeom)[i];
1797
1798 const Float_t qe = pix.GetPheFFactorMethod()
1799 / fBlindPixel->GetFluxInsidePlexiglass()
1800 / geo.GetA()
1801 * fQECam->GetPlexiglassQE();
1802
1803 const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1804 + fQECam->GetPlexiglassQERelVar()
1805 + pix.GetPheFFactorMethodRelVar();
1806
1807 qepix.SetQEBlindPixel ( qe , fPulserColor );
1808 qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor );
1809 qepix.SetBlindPixelMethodValid( kTRUE , fPulserColor );
1810
1811 if (!qepix.UpdateBlindPixelMethod())
1812 *fLog << warn << GetDescriptor()
1813 << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl;
1814 }
1815}
1816
1817// ------------------------------------------------------------------------
1818//
1819// Loop over pixels:
1820//
1821// - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set:
1822// MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor)
1823//
1824// - Calculate the quantum efficiency with the formula:
1825//
1826// QE = Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA()
1827//
1828// - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor );
1829// - Set Variance of QE in MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor );
1830// - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor)
1831//
1832// - Call MCalibrationQEPix::UpdatePINDiodeMethod()
1833//
1834void MCalibrationChargeCalc::FinalizePINDiodeQECam()
1835{
1836
1837 const UInt_t npixels = fGeom->GetNumPixels();
1838
1839 //
1840 // With the knowledge of the overall photon flux, calculate the
1841 // quantum efficiencies after the PIN Diode method
1842 //
1843 for (UInt_t i=0; i<npixels; i++)
1844 {
1845
1846 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1847
1848 if (!fPINDiode)
1849 {
1850 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1851 continue;
1852 }
1853
1854 if (!fPINDiode->IsFluxOutsidePlexiglassAvailable())
1855 {
1856 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1857 continue;
1858 }
1859
1860 MBadPixelsPix &bad = (*fBadPixels)[i];
1861
1862 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1863 {
1864 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1865 continue;
1866 }
1867
1868 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1869 MGeomPix &geo = (*fGeom)[i];
1870
1871 const Float_t qe = pix.GetPheFFactorMethod()
1872 / fPINDiode->GetFluxOutsidePlexiglass()
1873 / geo.GetA();
1874
1875 const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar();
1876
1877 qepix.SetQEPINDiode ( qe , fPulserColor );
1878 qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor );
1879 qepix.SetPINDiodeMethodValid( kTRUE , fPulserColor );
1880
1881 if (!qepix.UpdatePINDiodeMethod())
1882 *fLog << warn << GetDescriptor()
1883 << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl;
1884 }
1885}
1886
1887// -----------------------------------------------------------------------------------------------
1888//
1889// - Print out statistics about BadPixels of type UnsuitableType_t
1890// - store numbers of bad pixels of each type in fCam
1891//
1892void MCalibrationChargeCalc::FinalizeUnsuitablePixels()
1893{
1894
1895 *fLog << inf << endl;
1896 *fLog << GetDescriptor() << ": Charge Calibration status:" << endl;
1897 *fLog << dec << setfill(' ');
1898
1899 const Int_t nareas = fGeom->GetNumAreas();
1900
1901 Int_t counts[nareas];
1902 memset(counts,0,nareas*sizeof(Int_t));
1903
1904 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1905 {
1906 MBadPixelsPix &bad = (*fBadPixels)[i];
1907 if (!bad.IsBad())
1908 {
1909 const Int_t aidx = (*fGeom)[i].GetAidx();
1910 counts[aidx]++;
1911 }
1912 }
1913
1914 if (fGeom->InheritsFrom("MGeomCamMagic"))
1915 *fLog << " " << setw(7) << "Successfully calibrated Pixels: "
1916 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1917
1918 memset(counts,0,nareas*sizeof(Int_t));
1919
1920 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1921 {
1922 MBadPixelsPix &bad = (*fBadPixels)[i];
1923 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1924 {
1925 const Int_t aidx = (*fGeom)[i].GetAidx();
1926 counts[aidx]++;
1927 }
1928 }
1929
1930 for (Int_t aidx=0; aidx<nareas; aidx++)
1931 fCam->SetNumUnsuitable(counts[aidx], aidx);
1932
1933 if (fGeom->InheritsFrom("MGeomCamMagic"))
1934 *fLog << " " << setw(7) << "Uncalibrated Pixels: "
1935 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1936
1937 memset(counts,0,nareas*sizeof(Int_t));
1938
1939 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1940 {
1941 MBadPixelsPix &bad = (*fBadPixels)[i];
1942 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
1943 {
1944 const Int_t aidx = (*fGeom)[i].GetAidx();
1945 counts[aidx]++;
1946 }
1947 }
1948
1949 for (Int_t aidx=0; aidx<nareas; aidx++)
1950 fCam->SetNumUnreliable(counts[aidx], aidx);
1951
1952 *fLog << " " << setw(7) << "Unreliable Pixels: "
1953 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1954
1955}
1956
1957// -----------------------------------------------------------------------------------------------
1958//
1959// Print out statistics about BadPixels of type UncalibratedType_t
1960//
1961void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
1962{
1963
1964 UInt_t countinner = 0;
1965 UInt_t countouter = 0;
1966
1967 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1968 {
1969 MBadPixelsPix &bad = (*fBadPixels)[i];
1970 if (bad.IsUncalibrated(typ))
1971 {
1972 if (fGeom->GetPixRatio(i) == 1.)
1973 countinner++;
1974 else
1975 countouter++;
1976 }
1977 }
1978
1979 *fLog << " " << setw(7) << text
1980 << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
1981}
1982
1983// --------------------------------------------------------------------------
1984//
1985// Set the path for output file
1986//
1987void MCalibrationChargeCalc::SetOutputPath(TString path)
1988{
1989 fOutputPath = path;
1990 if (fOutputPath.EndsWith("/"))
1991 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
1992}
1993
1994// --------------------------------------------------------------------------
1995//
1996// Set the output file
1997//
1998void MCalibrationChargeCalc::SetOutputFile(TString file)
1999{
2000 fOutputFile = file;
2001}
2002
2003// --------------------------------------------------------------------------
2004//
2005// Get the output file
2006//
2007const char* MCalibrationChargeCalc::GetOutputFile()
2008{
2009 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
2010}
2011
2012// --------------------------------------------------------------------------
2013//
2014// Read the environment for the following data members:
2015// - fChargeLimit
2016// - fChargeErrLimit
2017// - fChargeRelErrLimit
2018// - fDebug
2019// - fFFactorErrLimit
2020// - fLambdaErrLimit
2021// - fLambdaCheckErrLimit
2022// - fPheErrLimit
2023//
2024Int_t MCalibrationChargeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
2025{
2026
2027 Bool_t rc = kFALSE;
2028 if (IsEnvDefined(env, prefix, "ChargeLimit", print))
2029 {
2030 SetChargeLimit(GetEnvValue(env, prefix, "ChargeLimit", fChargeLimit));
2031 rc = kTRUE;
2032 }
2033 if (IsEnvDefined(env, prefix, "ChargeErrLimit", print))
2034 {
2035 SetChargeErrLimit(GetEnvValue(env, prefix, "ChargeErrLimit", fChargeErrLimit));
2036 rc = kTRUE;
2037 }
2038 if (IsEnvDefined(env, prefix, "ChargeRelErrLimit", print))
2039 {
2040 SetChargeRelErrLimit(GetEnvValue(env, prefix, "ChargeRelErrLimit", fChargeRelErrLimit));
2041 rc = kTRUE;
2042 }
2043 if (IsEnvDefined(env, prefix, "Debug", print))
2044 {
2045 SetDebug(GetEnvValue(env, prefix, "Debug", IsDebug()));
2046 rc = kTRUE;
2047 }
2048 if (IsEnvDefined(env, prefix, "FFactorErrLimit", print))
2049 {
2050 SetFFactorErrLimit(GetEnvValue(env, prefix, "FFactorErrLimit", fFFactorErrLimit));
2051 rc = kTRUE;
2052 }
2053 if (IsEnvDefined(env, prefix, "LambdaErrLimit", print))
2054 {
2055 SetLambdaErrLimit(GetEnvValue(env, prefix, "LambdaErrLimit", fLambdaErrLimit));
2056 rc = kTRUE;
2057 }
2058 if (IsEnvDefined(env, prefix, "LambdaCheckLimit", print))
2059 {
2060 SetLambdaCheckLimit(GetEnvValue(env, prefix, "LambdaCheckLimit", fLambdaCheckLimit));
2061 rc = kTRUE;
2062 }
2063 if (IsEnvDefined(env, prefix, "PheErrLimit", print))
2064 {
2065 SetPheErrLimit(GetEnvValue(env, prefix, "PheErrLimit", fPheErrLimit));
2066 rc = kTRUE;
2067 }
2068 // void SetPulserColor(const MCalibrationCam::PulserColor_t col) { fPulserColor = col; }
2069
2070 return rc;
2071}
Note: See TracBrowser for help on using the repository browser.