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

Last change on this file since 4875 was 4853, checked in by gaug, 20 years ago
*** empty log message ***
File size: 69.7 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 FinalizeCombinedQECam();
714
715 //
716 // Re-direct the output to an ascii-file from now on:
717 //
718 MLog asciilog;
719 asciilog.SetOutputFile(GetOutputFile(),kTRUE);
720 SetLogStream(&asciilog);
721 //
722 // Finalize calibration statistics
723 //
724 FinalizeUnsuitablePixels();
725
726 fCam ->SetReadyToSave();
727 fQECam ->SetReadyToSave();
728 fBadPixels->SetReadyToSave();
729
730 if (fBlindPixel)
731 fBlindPixel->SetReadyToSave();
732 if (fBlindCam)
733 fBlindCam->SetReadyToSave();
734 if (fPINDiode)
735 fPINDiode->SetReadyToSave();
736
737 *fLog << inf << endl;
738 *fLog << GetDescriptor() << ": Fatal errors statistics:" << endl;
739
740 PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal,
741 Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS: "));
742 PrintUncalibrated(MBadPixelsPix::kChargeRelErrNotValid,
743 Form("%s%2.1f%s","Signal Error bigger than ",fChargeRelErrLimit," times Mean Signal: "));
744 PrintUncalibrated(MBadPixelsPix::kLoGainSaturation,
745 "Pixels with Low Gain Saturation: ");
746 PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin,
747 Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s): "));
748 PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins,
749 Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s): "));
750 PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes,
751 "Pixels with deviating number of phes: ");
752 PrintUncalibrated(MBadPixelsPix::kDeviatingFFactor,
753 "Pixels with deviating F-Factor: ");
754
755 *fLog << inf << endl;
756 *fLog << GetDescriptor() << ": Unreliable errors statistics:" << endl;
757
758 PrintUncalibrated(MBadPixelsPix::kChargeSigmaNotValid,
759 "Signal Sigma smaller than Pedestal RMS: ");
760 PrintUncalibrated(MBadPixelsPix::kHiGainOscillating,
761 "Pixels with changing Hi Gain signal over time: ");
762 PrintUncalibrated(MBadPixelsPix::kLoGainOscillating,
763 "Pixels with changing Lo Gain signal over time: ");
764 PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted,
765 "Pixels with unsuccesful Gauss fit to the Hi Gain: ");
766 PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted,
767 "Pixels with unsuccesful Gauss fit to the Lo Gain: ");
768
769 SetLogStream(&gLog);
770
771 return kTRUE;
772}
773
774// ----------------------------------------------------------------------------------
775//
776// Retrieves pedestal and pedestal RMS from MPedestalPix
777// Retrieves total entries from MPedestalCam
778// Sets pedestal*fNumHiGainSamples and pedestal*fNumLoGainSamples in MCalibrationChargePix
779// Sets pedRMS *fSqrtHiGainSamples and pedRMS *fSqrtLoGainSamples in MCalibrationChargePix
780//
781// If the flag MCalibrationPix::IsHiGainSaturation() is set, call also:
782// - MCalibrationChargePix::CalcLoGainPedestal()
783//
784void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal, const Int_t aidx)
785{
786
787 //
788 // get the pedestals
789 //
790 const Float_t pedes = ped.GetPedestal();
791 const Float_t prms = ped.GetPedestalRms();
792 const Int_t num = fPedestals->GetTotalEntries();
793
794 //
795 // RMS error set by PedCalcFromLoGain, 0 in case MPedCalcPedRun was used.
796 //
797 const Float_t prmserr = num>0 ? prms/TMath::Sqrt(2.*num) : ped.GetPedestalRmsError();
798
799 //
800 // set them in the calibration camera
801 //
802 if (cal.IsHiGainSaturation())
803 {
804 cal.SetPedestal(pedes * fNumLoGainSamples,
805 prms * fSqrtLoGainSamples,
806 prmserr * fSqrtLoGainSamples);
807 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples, aidx);
808 }
809 else
810 {
811 cal.SetPedestal(pedes * fNumHiGainSamples,
812 prms * fSqrtHiGainSamples,
813 prmserr * fSqrtHiGainSamples);
814 }
815
816}
817
818// ----------------------------------------------------------------------------------------------------
819//
820// Check fit results validity. Bad Pixels flags are set if:
821//
822// 1) Pixel has a mean smaller than fChargeLimit*PedRMS ( Flag: MBadPixelsPix::kChargeIsPedestal)
823// 2) Pixel has a mean error smaller than fChargeErrLimit ( Flag: MBadPixelsPix::kChargeErrNotValid)
824// 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error
825// ( Flag: MBadPixelsPix::kChargeRelErrNotValid)
826// 4) Pixel has a sigma bigger than its Pedestal RMS ( Flag: MBadPixelsPix::kChargeSigmaNotValid )
827//
828// Further returns if flags: MBadPixelsPix::kUnsuitableRun is set
829//
830// Calls MCalibrationChargePix::CalcReducedSigma() and sets flag: MBadPixelsPix::kChargeIsPedestal
831// and returns kFALSE if not succesful.
832//
833// Calls MCalibrationChargePix::CalcFFactor() and sets flag: MBadPixelsPix::kDeviatingNumPhes)
834// and returns kFALSE if not succesful.
835//
836// Calls MCalibrationChargePix::CalcConvFFactor()and sets flag: MBadPixelsPix::kDeviatingNumPhes)
837// and returns kFALSE if not succesful.
838//
839Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what)
840{
841
842 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
843 return kFALSE;
844
845 if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
846 {
847 *fLog << warn << GetDescriptor()
848 << Form(": Fitted Charge: %5.2f is smaller than %2.1f",cal.GetMean(),fChargeLimit)
849 << Form(" Pedestal RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
850 bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
851 }
852
853 if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr())
854 {
855 *fLog << warn << GetDescriptor()
856 << Form(": Fitted Charge: %4.2f is smaller than %2.1f",cal.GetMean(),fChargeRelErrLimit)
857 << Form(" times its error: %4.2f in %s%4i",cal.GetMeanErr(),what,cal.GetPixId()) << endl;
858 bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
859 }
860
861 if (cal.GetSigma() < cal.GetPedRms())
862 {
863 *fLog << warn << GetDescriptor()
864 << Form(": Sigma of Fitted Charge: %6.2f is smaller than",cal.GetSigma())
865 << Form(" Ped. RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
866 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
867 return kFALSE;
868 }
869
870 if (!cal.CalcReducedSigma())
871 {
872 *fLog << warn << GetDescriptor()
873 << Form(": Could not calculate the reduced sigma in %s: ",what)
874 << Form(" %4i",cal.GetPixId())
875 << endl;
876 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
877 return kFALSE;
878 }
879
880 if (!cal.CalcFFactor())
881 {
882 *fLog << warn << GetDescriptor()
883 << Form(": Could not calculate the F-Factor in %s: ",what)
884 << Form(" %4i",cal.GetPixId())
885 << endl;
886 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
887 return kFALSE;
888 }
889
890 if (!cal.CalcConvFFactor())
891 {
892 *fLog << warn << GetDescriptor()
893 << Form(": Could not calculate the Conv. FADC counts to Phes in %s: ",what)
894 << Form(" %4i",cal.GetPixId())
895 << endl;
896 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
897 return kFALSE;
898 }
899
900 return kTRUE;
901}
902
903
904
905// -----------------------------------------------------------------------------------
906//
907// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
908// - MBadPixelsPix::kChargeIsPedestal
909// - MBadPixelsPix::kChargeRelErrNotValid
910// - MBadPixelsPix::kMeanTimeInFirstBin
911// - MBadPixelsPix::kMeanTimeInLast2Bins
912// - MBadPixelsPix::kDeviatingNumPhes
913// - MBadPixelsPix::kHiGainOverFlow
914// - MBadPixelsPix::kLoGainOverFlow
915//
916// - Call MCalibrationPix::SetExcluded() for the bad pixels
917//
918// Sets pixel to MBadPixelsPix::kUnreliableRun, if one of the following flags is set:
919// - MBadPixelsPix::kChargeSigmaNotValid
920//
921void MCalibrationChargeCalc::FinalizeBadPixels()
922{
923
924 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
925 {
926
927 MBadPixelsPix &bad = (*fBadPixels)[i];
928 MCalibrationPix &pix = (*fCam)[i];
929
930 if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal))
931 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
932
933 if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid ))
934 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
935
936 if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ))
937 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
938
939 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ))
940 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
941
942 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ))
943 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
944
945 if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes ))
946 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
947
948 if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOverFlow ))
949 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
950
951 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOverFlow ))
952 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
953
954 if (bad.IsUnsuitable( MBadPixelsPix::kUnsuitableRun ))
955 pix.SetExcluded();
956
957 if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ))
958 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
959 }
960}
961
962// ------------------------------------------------------------------------
963//
964//
965// First loop: Calculate a mean and mean RMS of photo-electrons per area index
966// Include only pixels which are not MBadPixelsPix::kUnsuitableRun nor
967// MBadPixelsPix::kChargeSigmaNotValid (see FinalizeBadPixels()) and set
968// MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case.
969//
970// Second loop: Get mean number of photo-electrons and its RMS including
971// only pixels with flag MCalibrationChargePix::IsFFactorMethodValid()
972// and further exclude those deviating by more than fPheErrLimit mean
973// sigmas from the mean (obtained in first loop). Set
974// MBadPixelsPix::kDeviatingNumPhes if excluded.
975//
976// For the suitable pixels with flag MBadPixelsPix::kChargeSigmaNotValid
977// set the number of photo-electrons as the mean number of photo-electrons
978// calculated in that area index.
979//
980// Set weighted mean and variance of photo-electrons per area index in:
981// average area pixels of MCalibrationChargeCam (obtained from:
982// MCalibrationChargeCam::GetAverageArea() )
983//
984// Set weighted mean and variance of photo-electrons per sector in:
985// average sector pixels of MCalibrationChargeCam (obtained from:
986// MCalibrationChargeCam::GetAverageSector() )
987//
988//
989// Third loop: Set mean number of photo-electrons and its RMS in the pixels
990// only excluded as: MBadPixelsPix::kChargeSigmaNotValid
991//
992Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod()
993{
994
995 const UInt_t npixels = fGeom->GetNumPixels();
996 const UInt_t nareas = fGeom->GetNumAreas();
997 const UInt_t nsectors = fGeom->GetNumSectors();
998
999 TArrayF lowlim (nareas);
1000 TArrayF upplim (nareas);
1001 TArrayD areavars (nareas);
1002 TArrayD areaweights (nareas);
1003 TArrayD sectorweights (nsectors);
1004 TArrayD areaphes (nareas);
1005 TArrayD sectorphes (nsectors);
1006 TArrayI numareavalid (nareas);
1007 TArrayI numsectorvalid(nsectors);
1008
1009 //
1010 // First loop: Get mean number of photo-electrons and the RMS
1011 // The loop is only to recognize later pixels with very deviating numbers
1012 //
1013 MHCamera camphes(*fGeom,"Camphes","Phes in Camera");
1014
1015 for (UInt_t i=0; i<npixels; i++)
1016 {
1017
1018 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam) [i];
1019 MBadPixelsPix &bad = (*fBadPixels)[i];
1020
1021 if (!pix.IsFFactorMethodValid())
1022 continue;
1023
1024 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1025 {
1026 pix.SetFFactorMethodValid(kFALSE);
1027 continue;
1028 }
1029
1030 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1031 continue;
1032
1033 const Float_t nphe = pix.GetPheFFactorMethod();
1034 const Int_t aidx = (*fGeom)[i].GetAidx();
1035
1036 camphes.Fill(i,nphe);
1037 camphes.SetUsed(i);
1038
1039 areaphes [aidx] += nphe;
1040 areavars [aidx] += nphe*nphe;
1041 numareavalid[aidx] ++;
1042 }
1043
1044 for (UInt_t i=0; i<nareas; i++)
1045 {
1046 if (numareavalid[i] == 0)
1047 {
1048 *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
1049 << "in area index: " << i << endl;
1050 continue;
1051 }
1052
1053 if (numareavalid[i] == 1)
1054 areavars[i] = 0.;
1055 else if (numareavalid[i] == 0)
1056 {
1057 areaphes[i] = -1.;
1058 areaweights[i] = -1.;
1059 }
1060 else
1061 {
1062 areavars[i] = (areavars[i] - areaphes[i]*areaphes[i]/numareavalid[i]) / (numareavalid[i]-1);
1063 areaphes[i] = areaphes[i] / numareavalid[i];
1064 }
1065
1066 if (areavars[i] < 0.)
1067 {
1068 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of photo-electrons found "
1069 << "in area index: " << i << endl;
1070 continue;
1071 }
1072
1073 lowlim [i] = areaphes[i] - fPheErrLimit*TMath::Sqrt(areavars[i]);
1074 upplim [i] = areaphes[i] + fPheErrLimit*TMath::Sqrt(areavars[i]);
1075
1076 TArrayI area(1);
1077 area[0] = i;
1078
1079 TH1D *hist = camphes.ProjectionS(TArrayI(),area,"_py",100);
1080 hist->Fit("gaus","Q");
1081 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1082 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1083 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1084
1085 if (IsDebug())
1086 camphes.DrawClone();
1087
1088 if (ndf < 2)
1089 {
1090 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1091 << "in the camera with area index: " << i << endl;
1092 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1093 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1094 delete hist;
1095 continue;
1096 }
1097
1098 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1099
1100 if (prob < 0.001)
1101 {
1102 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1103 << "in the camera with area index: " << i << endl;
1104 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1105 << " is smaller than 0.001 " << endl;
1106 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1107 delete hist;
1108 continue;
1109 }
1110
1111 *fLog << inf << GetDescriptor() << ": Mean number of photo-electrons "
1112 << "with area idx " << i << ": "
1113 << Form("%7.2f+-%6.2f",mean,sigma) << endl;
1114
1115 lowlim [i] = mean - fPheErrLimit*sigma;
1116 upplim [i] = mean + fPheErrLimit*sigma;
1117
1118 delete hist;
1119 }
1120
1121 *fLog << endl;
1122
1123 numareavalid.Reset();
1124 areaphes .Reset();
1125 areavars .Reset();
1126 //
1127 // Second loop: Get mean number of photo-electrons and its RMS excluding
1128 // pixels deviating by more than fPheErrLimit sigma.
1129 // Set the conversion factor FADC counts to photo-electrons
1130 //
1131 for (UInt_t i=0; i<npixels; i++)
1132 {
1133
1134 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1135
1136 if (!pix.IsFFactorMethodValid())
1137 continue;
1138
1139 MBadPixelsPix &bad = (*fBadPixels)[i];
1140
1141 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1142 continue;
1143
1144 const Float_t nvar = pix.GetPheFFactorMethodVar();
1145
1146 if (nvar <= 0.)
1147 {
1148 pix.SetFFactorMethodValid(kFALSE);
1149 continue;
1150 }
1151
1152 const Int_t aidx = (*fGeom)[i].GetAidx();
1153 const Int_t sector = (*fGeom)[i].GetSector();
1154 const Float_t area = (*fGeom)[i].GetA();
1155 const Float_t nphe = pix.GetPheFFactorMethod();
1156
1157 if ( nphe < lowlim[aidx] || nphe > upplim[aidx] )
1158 {
1159 *fLog << warn << GetDescriptor() << ": Number of phes: "
1160 << Form("%7.2f out of %3.1f sigma limit: ",nphe,fPheErrLimit)
1161 << Form("[%7.2f,%7.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
1162 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1163 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1164 pix.SetFFactorMethodValid(kFALSE);
1165 continue;
1166 }
1167
1168 areaweights [aidx] += nphe*nphe;
1169 areaphes [aidx] += nphe;
1170 numareavalid [aidx] ++;
1171
1172 if (aidx == 0)
1173 fNumInnerFFactorMethodUsed++;
1174
1175 sectorweights [sector] += nphe*nphe/area/area;
1176 sectorphes [sector] += nphe/area;
1177 numsectorvalid[sector] ++;
1178 }
1179
1180 *fLog << endl;
1181
1182 for (UInt_t aidx=0; aidx<nareas; aidx++)
1183 {
1184
1185 MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
1186
1187 if (numareavalid[aidx] == 1)
1188 areaweights[aidx] = 0.;
1189 else if (numareavalid[aidx] == 0)
1190 {
1191 areaphes[aidx] = -1.;
1192 areaweights[aidx] = -1.;
1193 }
1194 else
1195 {
1196 areaweights[aidx] = (areaweights[aidx] - areaphes[aidx]*areaphes[aidx]/numareavalid[aidx])
1197 / (numareavalid[aidx]-1);
1198 areaphes[aidx] /= numareavalid[aidx];
1199 }
1200
1201 if (areaweights[aidx] < 0. || areaphes[aidx] <= 0.)
1202 {
1203 *fLog << warn << GetDescriptor()
1204 << ": Mean number phes from area index " << aidx << " could not be calculated: "
1205 << " Mean: " << areaphes[aidx]
1206 << " Variance: " << areaweights[aidx] << endl;
1207 apix.SetFFactorMethodValid(kFALSE);
1208 continue;
1209 }
1210
1211 *fLog << inf << GetDescriptor()
1212 << ": Average total number phes in area idx " << aidx << ": "
1213 << Form("%7.2f%s%6.2f",areaphes[aidx]," +- ",TMath::Sqrt(areaweights[aidx])) << endl;
1214
1215 apix.SetPheFFactorMethod ( areaphes[aidx] );
1216 apix.SetPheFFactorMethodVar( areaweights[aidx] / numareavalid[aidx] );
1217 apix.SetFFactorMethodValid ( kTRUE );
1218
1219 }
1220
1221 *fLog << endl;
1222
1223 for (UInt_t sector=0; sector<nsectors; sector++)
1224 {
1225
1226 if (numsectorvalid[sector] == 1)
1227 sectorweights[sector] = 0.;
1228 else if (numsectorvalid[sector] == 0)
1229 {
1230 sectorphes[sector] = -1.;
1231 sectorweights[sector] = -1.;
1232 }
1233 else
1234 {
1235 sectorweights[sector] = (sectorweights[sector]
1236 - sectorphes[sector]*sectorphes[sector]/numsectorvalid[sector]
1237 )
1238 / (numsectorvalid[sector]-1.);
1239 sectorphes[sector] /= numsectorvalid[sector];
1240 }
1241
1242 MCalibrationChargePix &spix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
1243
1244 if (sectorweights[sector] < 0. || sectorphes[sector] <= 0.)
1245 {
1246 *fLog << warn << GetDescriptor()
1247 <<": Mean number phes per area for sector " << sector << " could not be calculated: "
1248 << " Mean: " << sectorphes[sector]
1249 << " Variance: " << sectorweights[sector] << endl;
1250 spix.SetFFactorMethodValid(kFALSE);
1251 continue;
1252 }
1253
1254 *fLog << inf << GetDescriptor()
1255 << ": Average number phes per area in sector " << sector << ": "
1256 << Form("%5.2f+-%4.2f [phe/mm^2]",sectorphes[sector],TMath::Sqrt(sectorweights[sector]))
1257 << endl;
1258
1259 spix.SetPheFFactorMethod ( sectorphes[sector] );
1260 spix.SetPheFFactorMethodVar( sectorweights[sector] / numsectorvalid[sector]);
1261 spix.SetFFactorMethodValid ( kTRUE );
1262
1263 }
1264
1265 //
1266 // Third loop: Set mean number of photo-electrons and its RMS in the pixels
1267 // only excluded as: MBadPixelsPix::kChargeSigmaNotValid
1268 //
1269 for (UInt_t i=0; i<npixels; i++)
1270 {
1271
1272 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1273 MBadPixelsPix &bad = (*fBadPixels)[i];
1274
1275 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1276 continue;
1277
1278 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1279 {
1280 const Int_t aidx = (*fGeom)[i].GetAidx();
1281 MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
1282 pix.SetPheFFactorMethod ( apix.GetPheFFactorMethod() );
1283 pix.SetPheFFactorMethodVar( apix.GetPheFFactorMethodVar() );
1284 if (!pix.CalcConvFFactor())
1285 {
1286 *fLog << warn << GetDescriptor()
1287 << ": Could not calculate the Conv. FADC counts to Phes in pixel: "
1288 << Form(" %4i",pix.GetPixId())
1289 << endl;
1290 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1291 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1292 }
1293 }
1294 }
1295
1296 return kTRUE;
1297}
1298
1299
1300// ------------------------------------------------------------------------
1301//
1302// Returns kFALSE if pointer to MCalibrationChargeBlindPix is NULL
1303//
1304// The check returns kFALSE if:
1305//
1306// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
1307// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
1308//
1309// Calls:
1310// - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
1311//
1312Bool_t MCalibrationChargeCalc::FinalizeBlindPixel()
1313{
1314
1315 if (!fBlindPixel)
1316 return kFALSE;
1317
1318 const Float_t lambda = fBlindPixel->GetLambda();
1319 const Float_t lambdaerr = fBlindPixel->GetLambdaErr();
1320 const Float_t lambdacheck = fBlindPixel->GetLambdaCheck();
1321
1322 if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
1323 {
1324 *fLog << warn << GetDescriptor()
1325 << Form("%s%4.2f%s%4.2f%s%4.2f%s",": Lambda: ",lambda," and Lambda-Check: ",
1326 lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel ")
1327 << endl;
1328 return kFALSE;
1329 }
1330
1331 if (lambdaerr > fLambdaErrLimit)
1332 {
1333 *fLog << warn << GetDescriptor()
1334 << Form("%s%4.2f%s%4.2f%s",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
1335 fLambdaErrLimit," in Blind Pixel ") << endl;
1336 return kFALSE;
1337 }
1338
1339 if (!fBlindPixel->CalcFluxInsidePlexiglass())
1340 {
1341 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
1342 << "will skip Blind Pixel Calibration " << endl;
1343 return kFALSE;
1344 }
1345
1346 return kTRUE;
1347}
1348
1349// ------------------------------------------------------------------------
1350//
1351// Returns kFALSE if pointer to MCalibrationChargeBlindCam is NULL
1352//
1353// The check returns kFALSE if:
1354//
1355// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
1356// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
1357//
1358// Calls:
1359// - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
1360//
1361Bool_t MCalibrationChargeCalc::FinalizeBlindCam()
1362{
1363
1364 if (!fBlindCam)
1365 return kFALSE;
1366
1367 Float_t flux = 0.;
1368 Float_t fluxvar = 0.;
1369 Int_t nvalid = 0;
1370
1371 for (UInt_t i=0; i<fBlindCam->GetNumBlindPixels(); i++)
1372 {
1373
1374 MCalibrationChargeBlindPix &blindpix = (*fBlindCam)[i];
1375
1376 if (!blindpix.IsValid())
1377 continue;
1378
1379 const Float_t lambda = blindpix.GetLambda();
1380 const Float_t lambdaerr = blindpix.GetLambdaErr();
1381 const Float_t lambdacheck = blindpix.GetLambdaCheck();
1382
1383 if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
1384 {
1385 *fLog << warn << GetDescriptor()
1386 << Form("%s%4.2f%s%4.2f%s%4.2f%s%2i",": Lambda: ",lambda," and Lambda-Check: ",
1387 lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel Nr.",i)
1388 << endl;
1389 blindpix.SetValid(kFALSE);
1390 continue;
1391 }
1392
1393 if (lambdaerr > fLambdaErrLimit)
1394 {
1395 *fLog << warn << GetDescriptor()
1396 << Form("%s%4.2f%s%4.2f%s%2i",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
1397 fLambdaErrLimit," in Blind Pixel Nr.",i) << endl;
1398 blindpix.SetValid(kFALSE);
1399 continue;
1400 }
1401
1402 if (!blindpix.CalcFluxInsidePlexiglass())
1403 {
1404 *fLog << warn << "Could not calculate the flux of photons from Blind Pixel Nr." << i << endl;
1405 blindpix.SetValid(kFALSE);
1406 continue;
1407 }
1408
1409 nvalid++;
1410 const Float_t weight = 1./ blindpix.GetFluxInsidePlexiglassErr() / blindpix.GetFluxInsidePlexiglassErr();
1411 flux += weight * blindpix.GetFluxInsidePlexiglass();
1412 fluxvar += weight;
1413 }
1414
1415 if (!nvalid)
1416 return kFALSE;
1417
1418 flux /= fluxvar;
1419 fluxvar /= 1./fluxvar;
1420
1421 const Float_t photons = flux * (*fGeom)[0].GetA() / fQECam->GetPlexiglassQE();
1422 fCam->SetNumPhotonsBlindPixelMethod(photons);
1423
1424 const Float_t photrelvar = fluxvar / flux / flux + fQECam->GetPlexiglassQERelVar();
1425 if (photrelvar > 0.)
1426 fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1427
1428 return kTRUE;
1429}
1430
1431// ------------------------------------------------------------------------
1432//
1433// Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL
1434//
1435// The check returns kFALSE if:
1436//
1437// 1) PINDiode has a fitted charge smaller than fChargeLimit*PedRMS
1438// 2) PINDiode has a fit error smaller than fChargeErrLimit
1439// 3) PINDiode has a fitted charge smaller its fChargeRelErrLimit times its charge error
1440// 4) PINDiode has a charge sigma smaller than its Pedestal RMS
1441//
1442// Calls:
1443// - MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass()
1444//
1445Bool_t MCalibrationChargeCalc::FinalizePINDiode()
1446{
1447
1448 if (!fPINDiode)
1449 return kFALSE;
1450
1451 if (fPINDiode->GetMean() < fChargeLimit*fPINDiode->GetPedRms())
1452 {
1453 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1454 << fChargeLimit << " Pedestal RMS in PINDiode " << endl;
1455 return kFALSE;
1456 }
1457
1458 if (fPINDiode->GetMeanErr() < fChargeErrLimit)
1459 {
1460 *fLog << warn << GetDescriptor() << ": Error of Fitted Charge is smaller than "
1461 << fChargeErrLimit << " in PINDiode " << endl;
1462 return kFALSE;
1463 }
1464
1465 if (fPINDiode->GetMean() < fChargeRelErrLimit*fPINDiode->GetMeanErr())
1466 {
1467 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1468 << fChargeRelErrLimit << "* its error in PINDiode " << endl;
1469 return kFALSE;
1470 }
1471
1472 if (fPINDiode->GetSigma() < fPINDiode->GetPedRms())
1473 {
1474 *fLog << warn << GetDescriptor()
1475 << ": Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl;
1476 return kFALSE;
1477 }
1478
1479
1480 if (!fPINDiode->CalcFluxOutsidePlexiglass())
1481 {
1482 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
1483 << "will skip PIN Diode Calibration " << endl;
1484 return kFALSE;
1485 }
1486
1487 return kTRUE;
1488}
1489
1490// ------------------------------------------------------------------------
1491//
1492// Calculate the average number of photons outside the plexiglass with the
1493// formula:
1494//
1495// av.Num.photons(area index) = av.Num.Phes(area index)
1496// / MCalibrationQEPix::GetDefaultQE(fPulserColor)
1497// / MCalibrationQEPix::GetPMTCollectionEff()
1498// / MCalibrationQEPix::GetLightGuidesEff(fPulserColor)
1499// / MCalibrationQECam::GetPlexiglassQE()
1500//
1501// Calculate the variance on the average number of photons assuming that the error on the
1502// Quantum efficiency is reduced by the number of used inner pixels, but the rest of the
1503// values keeps it ordinary error since it is systematic.
1504//
1505// Loop over pixels:
1506//
1507// - Continue, if not MCalibrationChargePix::IsFFactorMethodValid() and set:
1508// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor)
1509//
1510// - Call MCalibrationChargePix::CalcMeanFFactor(av.Num.photons) and set:
1511// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) if not succesful
1512//
1513// - Calculate the quantum efficiency with the formula:
1514//
1515// QE = ( Num.Phes / av.Num.photons ) * MGeomCam::GetPixRatio()
1516//
1517// - Set QE in MCalibrationQEPix::SetQEFFactor ( QE, fPulserColor );
1518//
1519// - Set Variance of QE in MCalibrationQEPix::SetQEFFactorVar ( Variance, fPulserColor );
1520// - Set bit MCalibrationQEPix::SetFFactorMethodValid(kTRUE,fPulserColor)
1521//
1522// - Call MCalibrationQEPix::UpdateFFactorMethod()
1523//
1524void MCalibrationChargeCalc::FinalizeFFactorQECam()
1525{
1526
1527 if (fNumInnerFFactorMethodUsed < 2)
1528 {
1529 *fLog << warn << GetDescriptor()
1530 << ": Could not calculate F-Factor Method: Less than 2 inner pixels valid! " << endl;
1531 return;
1532 }
1533
1534 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCam->GetAverageArea(0);
1535 MCalibrationQEPix &qepix = (MCalibrationQEPix&) fQECam->GetAverageArea(0);
1536
1537 const Float_t avphotons = avpix.GetPheFFactorMethod()
1538 / qepix.GetDefaultQE(fPulserColor)
1539 / qepix.GetPMTCollectionEff()
1540 / qepix.GetLightGuidesEff(fPulserColor)
1541 / fQECam->GetPlexiglassQE();
1542
1543 const Float_t avphotrelvar = avpix.GetPheFFactorMethodRelVar()
1544 + qepix.GetDefaultQERelVar(fPulserColor) / fNumInnerFFactorMethodUsed
1545 + qepix.GetPMTCollectionEffRelVar()
1546 + qepix.GetLightGuidesEffRelVar(fPulserColor)
1547 + fQECam->GetPlexiglassQERelVar();
1548
1549 const UInt_t nareas = fGeom->GetNumAreas();
1550
1551 //
1552 // Set the results in the MCalibrationChargeCam
1553 //
1554 fCam->SetNumPhotonsFFactorMethod (avphotons);
1555 if (avphotrelvar > 0.)
1556 fCam->SetNumPhotonsFFactorMethodErr(TMath::Sqrt( avphotrelvar * avphotons * avphotons));
1557
1558 TArrayF lowlim (nareas);
1559 TArrayF upplim (nareas);
1560 TArrayD avffactorphotons (nareas);
1561 TArrayD avffactorphotvar (nareas);
1562 TArrayI numffactor (nareas);
1563
1564 const UInt_t npixels = fGeom->GetNumPixels();
1565
1566 MHCamera camffactor(*fGeom,"Camffactor","F-Factor in Camera");
1567
1568 for (UInt_t i=0; i<npixels; i++)
1569 {
1570
1571 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1572 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1573 MBadPixelsPix &bad = (*fBadPixels)[i];
1574
1575 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1576 continue;
1577
1578 const Float_t photons = avphotons / fGeom->GetPixRatio(i);
1579 const Float_t qe = pix.GetPheFFactorMethod() / photons ;
1580
1581 const Float_t qerelvar = avphotrelvar + pix.GetPheFFactorMethodRelVar();
1582
1583 qepix.SetQEFFactor ( qe , fPulserColor );
1584 qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor );
1585 qepix.SetFFactorMethodValid( kTRUE , fPulserColor );
1586
1587 if (!qepix.UpdateFFactorMethod())
1588 *fLog << warn << GetDescriptor()
1589 << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl;
1590
1591 //
1592 // The following pixels are those with deviating sigma, but otherwise OK,
1593 // probably those with stars during the pedestal run, but not the cal. run.
1594 //
1595 if (!pix.CalcMeanFFactor( photons , avphotrelvar ))
1596 {
1597 (*fBadPixels)[i].SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1598 (*fBadPixels)[i].SetUnsuitable ( MBadPixelsPix::kUnreliableRun );
1599 continue;
1600 }
1601
1602 const Int_t aidx = (*fGeom)[i].GetAidx();
1603 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1604
1605 camffactor.Fill(i,ffactor);
1606 camffactor.SetUsed(i);
1607
1608 avffactorphotons[aidx] += ffactor;
1609 avffactorphotvar[aidx] += ffactor*ffactor;
1610 numffactor[aidx]++;
1611 }
1612
1613 for (UInt_t i=0; i<nareas; i++)
1614 {
1615
1616 if (numffactor[i] == 0)
1617 {
1618 *fLog << warn << GetDescriptor() << ": No pixels with valid total F-Factor found "
1619 << "in area index: " << i << endl;
1620 continue;
1621 }
1622
1623 avffactorphotvar[i] = (avffactorphotvar[i] - avffactorphotons[i]*avffactorphotons[i]/numffactor[i]) / (numffactor[i]-1.);
1624 avffactorphotons[i] = avffactorphotons[i] / numffactor[i];
1625
1626 if (avffactorphotvar[i] < 0.)
1627 {
1628 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of total F-Factor found "
1629 << "in area index: " << i << endl;
1630 continue;
1631 }
1632
1633 lowlim [i] = 1.1; // Lowest known F-Factor of a PMT
1634 upplim [i] = avffactorphotons[i] + fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);
1635
1636 TArrayI area(1);
1637 area[0] = i;
1638
1639 TH1D *hist = camffactor.ProjectionS(TArrayI(),area,"_py",100);
1640 hist->Fit("gaus","Q");
1641 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1642 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1643 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1644
1645 if (IsDebug())
1646 camffactor.DrawClone();
1647
1648 if (ndf < 2)
1649 {
1650 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1651 << "in the camera with area index: " << i << endl;
1652 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1653 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1654 delete hist;
1655 continue;
1656 }
1657
1658 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1659
1660 if (prob < 0.001)
1661 {
1662 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1663 << "in the camera with area index: " << i << endl;
1664 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1665 << " is smaller than 0.001 " << endl;
1666 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1667 delete hist;
1668 continue;
1669 }
1670
1671 *fLog << inf << GetDescriptor() << ": Mean F-Factor "
1672 << "with area index #" << i << ": "
1673 << Form("%4.2f+-%4.2f",mean,sigma) << endl;
1674
1675 lowlim [i] = 1.1;
1676 upplim [i] = mean + fFFactorErrLimit*sigma;
1677
1678 delete hist;
1679 }
1680
1681 *fLog << endl;
1682
1683 for (UInt_t i=0; i<npixels; i++)
1684 {
1685
1686 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1687 MBadPixelsPix &bad = (*fBadPixels)[i];
1688
1689 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1690 continue;
1691
1692 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1693 const Int_t aidx = (*fGeom)[i].GetAidx();
1694
1695 if ( ffactor < lowlim[aidx] || ffactor > upplim[aidx] )
1696 {
1697 *fLog << warn << GetDescriptor() << ": Overall F-Factor "
1698 << Form("%5.2f",ffactor) << " out of range ["
1699 << Form("%5.2f,%5.2f",lowlim[aidx],upplim[aidx]) << "] pixel " << i << endl;
1700
1701 bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1702 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1703 }
1704 }
1705
1706 for (UInt_t i=0; i<npixels; i++)
1707 {
1708
1709 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1710 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1711 MBadPixelsPix &bad = (*fBadPixels)[i];
1712
1713 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1714 {
1715 qepix.SetFFactorMethodValid(kFALSE,fPulserColor);
1716 pix.SetFFactorMethodValid(kFALSE);
1717 pix.SetExcluded();
1718 continue;
1719 }
1720 }
1721}
1722
1723
1724// ------------------------------------------------------------------------
1725//
1726// Loop over pixels:
1727//
1728// - Continue, if not MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() and set:
1729// MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor)
1730//
1731// - Calculate the quantum efficiency with the formula:
1732//
1733// QE = Num.Phes / MCalibrationChargeBlindPix::GetFluxInsidePlexiglass()
1734// / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE()
1735//
1736// - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor );
1737// - Set Variance of QE in MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor );
1738// - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor)
1739//
1740// - Call MCalibrationQEPix::UpdateBlindPixelMethod()
1741//
1742void MCalibrationChargeCalc::FinalizeBlindPixelQECam()
1743{
1744
1745 const UInt_t npixels = fGeom->GetNumPixels();
1746
1747 //
1748 // Set the results in the MCalibrationChargeCam
1749 //
1750 if (fBlindPixel)
1751 {
1752 if (fBlindPixel->IsFluxInsidePlexiglassAvailable())
1753 {
1754
1755 const Float_t photons = fBlindPixel->GetFluxInsidePlexiglass() * (*fGeom)[0].GetA()
1756 / fQECam->GetPlexiglassQE();
1757 fCam->SetNumPhotonsBlindPixelMethod(photons);
1758
1759 const Float_t photrelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1760 + fQECam->GetPlexiglassQERelVar();
1761 if (photrelvar > 0.)
1762 fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1763 }
1764 }
1765 //
1766 // With the knowledge of the overall photon flux, calculate the
1767 // quantum efficiencies after the Blind Pixel and PIN Diode method
1768 //
1769 for (UInt_t i=0; i<npixels; i++)
1770 {
1771
1772 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1773
1774 if (!fBlindPixel)
1775 {
1776 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1777 continue;
1778 }
1779
1780 if (!fBlindPixel->IsFluxInsidePlexiglassAvailable())
1781 {
1782 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1783 continue;
1784 }
1785
1786 MBadPixelsPix &bad = (*fBadPixels)[i];
1787 if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1788 {
1789 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1790 continue;
1791 }
1792
1793 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1794 MGeomPix &geo = (*fGeom)[i];
1795
1796 const Float_t qe = pix.GetPheFFactorMethod()
1797 / fBlindPixel->GetFluxInsidePlexiglass()
1798 / geo.GetA()
1799 * fQECam->GetPlexiglassQE();
1800
1801 const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1802 + fQECam->GetPlexiglassQERelVar()
1803 + pix.GetPheFFactorMethodRelVar();
1804
1805 qepix.SetQEBlindPixel ( qe , fPulserColor );
1806 qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor );
1807 qepix.SetBlindPixelMethodValid( kTRUE , fPulserColor );
1808
1809 if (!qepix.UpdateBlindPixelMethod())
1810 *fLog << warn << GetDescriptor()
1811 << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl;
1812 }
1813}
1814
1815// ------------------------------------------------------------------------
1816//
1817// Loop over pixels:
1818//
1819// - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set:
1820// MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor)
1821//
1822// - Calculate the quantum efficiency with the formula:
1823//
1824// QE = Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA()
1825//
1826// - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor );
1827// - Set Variance of QE in MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor );
1828// - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor)
1829//
1830// - Call MCalibrationQEPix::UpdatePINDiodeMethod()
1831//
1832void MCalibrationChargeCalc::FinalizePINDiodeQECam()
1833{
1834
1835 const UInt_t npixels = fGeom->GetNumPixels();
1836
1837 //
1838 // With the knowledge of the overall photon flux, calculate the
1839 // quantum efficiencies after the PIN Diode method
1840 //
1841 for (UInt_t i=0; i<npixels; i++)
1842 {
1843
1844 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1845
1846 if (!fPINDiode)
1847 {
1848 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1849 continue;
1850 }
1851
1852 if (!fPINDiode->IsFluxOutsidePlexiglassAvailable())
1853 {
1854 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1855 continue;
1856 }
1857
1858 MBadPixelsPix &bad = (*fBadPixels)[i];
1859
1860 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1861 {
1862 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1863 continue;
1864 }
1865
1866 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1867 MGeomPix &geo = (*fGeom)[i];
1868
1869 const Float_t qe = pix.GetPheFFactorMethod()
1870 / fPINDiode->GetFluxOutsidePlexiglass()
1871 / geo.GetA();
1872
1873 const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar();
1874
1875 qepix.SetQEPINDiode ( qe , fPulserColor );
1876 qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor );
1877 qepix.SetPINDiodeMethodValid( kTRUE , fPulserColor );
1878
1879 if (!qepix.UpdatePINDiodeMethod())
1880 *fLog << warn << GetDescriptor()
1881 << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl;
1882 }
1883}
1884
1885// ------------------------------------------------------------------------
1886//
1887// Loop over pixels:
1888//
1889// - Call MCalibrationQEPix::UpdateCombinedMethod()
1890//
1891void MCalibrationChargeCalc::FinalizeCombinedQECam()
1892{
1893
1894 const UInt_t npixels = fGeom->GetNumPixels();
1895
1896 for (UInt_t i=0; i<npixels; i++)
1897 {
1898
1899 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1900 MBadPixelsPix &bad = (*fBadPixels)[i];
1901
1902 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1903 {
1904 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1905 continue;
1906 }
1907
1908 qepix.UpdateCombinedMethod();
1909 }
1910}
1911
1912// -----------------------------------------------------------------------------------------------
1913//
1914// - Print out statistics about BadPixels of type UnsuitableType_t
1915// - store numbers of bad pixels of each type in fCam
1916//
1917void MCalibrationChargeCalc::FinalizeUnsuitablePixels()
1918{
1919
1920 *fLog << inf << endl;
1921 *fLog << GetDescriptor() << ": Charge Calibration status:" << endl;
1922 *fLog << dec << setfill(' ');
1923
1924 const Int_t nareas = fGeom->GetNumAreas();
1925
1926 TArrayI counts(nareas);
1927
1928 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1929 {
1930 MBadPixelsPix &bad = (*fBadPixels)[i];
1931 if (!bad.IsBad())
1932 {
1933 const Int_t aidx = (*fGeom)[i].GetAidx();
1934 counts[aidx]++;
1935 }
1936 }
1937
1938 if (fGeom->InheritsFrom("MGeomCamMagic"))
1939 *fLog << " " << setw(7) << "Successfully calibrated Pixels: "
1940 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1941
1942 counts.Reset();
1943
1944 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1945 {
1946 MBadPixelsPix &bad = (*fBadPixels)[i];
1947 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1948 {
1949 const Int_t aidx = (*fGeom)[i].GetAidx();
1950 counts[aidx]++;
1951 }
1952 }
1953
1954 for (Int_t aidx=0; aidx<nareas; aidx++)
1955 fCam->SetNumUnsuitable(counts[aidx], aidx);
1956
1957 if (fGeom->InheritsFrom("MGeomCamMagic"))
1958 *fLog << " " << setw(7) << "Uncalibrated Pixels: "
1959 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1960
1961 counts.Reset();
1962
1963 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1964 {
1965 MBadPixelsPix &bad = (*fBadPixels)[i];
1966 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
1967 {
1968 const Int_t aidx = (*fGeom)[i].GetAidx();
1969 counts[aidx]++;
1970 }
1971 }
1972
1973 for (Int_t aidx=0; aidx<nareas; aidx++)
1974 fCam->SetNumUnreliable(counts[aidx], aidx);
1975
1976 *fLog << " " << setw(7) << "Unreliable Pixels: "
1977 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1978
1979}
1980
1981// -----------------------------------------------------------------------------------------------
1982//
1983// Print out statistics about BadPixels of type UncalibratedType_t
1984//
1985void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
1986{
1987
1988 UInt_t countinner = 0;
1989 UInt_t countouter = 0;
1990
1991 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1992 {
1993 MBadPixelsPix &bad = (*fBadPixels)[i];
1994 if (bad.IsUncalibrated(typ))
1995 {
1996 if (fGeom->GetPixRatio(i) == 1.)
1997 countinner++;
1998 else
1999 countouter++;
2000 }
2001 }
2002
2003 *fLog << " " << setw(7) << text
2004 << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
2005}
2006
2007// --------------------------------------------------------------------------
2008//
2009// Set the path for output file
2010//
2011void MCalibrationChargeCalc::SetOutputPath(TString path)
2012{
2013 fOutputPath = path;
2014 if (fOutputPath.EndsWith("/"))
2015 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
2016}
2017
2018// --------------------------------------------------------------------------
2019//
2020// Set the output file
2021//
2022void MCalibrationChargeCalc::SetOutputFile(TString file)
2023{
2024 fOutputFile = file;
2025}
2026
2027// --------------------------------------------------------------------------
2028//
2029// Get the output file
2030//
2031const char* MCalibrationChargeCalc::GetOutputFile()
2032{
2033 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
2034}
2035
2036// --------------------------------------------------------------------------
2037//
2038// Read the environment for the following data members:
2039// - fChargeLimit
2040// - fChargeErrLimit
2041// - fChargeRelErrLimit
2042// - fDebug
2043// - fFFactorErrLimit
2044// - fLambdaErrLimit
2045// - fLambdaCheckErrLimit
2046// - fPheErrLimit
2047//
2048Int_t MCalibrationChargeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
2049{
2050
2051 Bool_t rc = kFALSE;
2052 if (IsEnvDefined(env, prefix, "ChargeLimit", print))
2053 {
2054 SetChargeLimit(GetEnvValue(env, prefix, "ChargeLimit", fChargeLimit));
2055 rc = kTRUE;
2056 }
2057 if (IsEnvDefined(env, prefix, "ChargeErrLimit", print))
2058 {
2059 SetChargeErrLimit(GetEnvValue(env, prefix, "ChargeErrLimit", fChargeErrLimit));
2060 rc = kTRUE;
2061 }
2062 if (IsEnvDefined(env, prefix, "ChargeRelErrLimit", print))
2063 {
2064 SetChargeRelErrLimit(GetEnvValue(env, prefix, "ChargeRelErrLimit", fChargeRelErrLimit));
2065 rc = kTRUE;
2066 }
2067 if (IsEnvDefined(env, prefix, "Debug", print))
2068 {
2069 SetDebug(GetEnvValue(env, prefix, "Debug", IsDebug()));
2070 rc = kTRUE;
2071 }
2072 if (IsEnvDefined(env, prefix, "FFactorErrLimit", print))
2073 {
2074 SetFFactorErrLimit(GetEnvValue(env, prefix, "FFactorErrLimit", fFFactorErrLimit));
2075 rc = kTRUE;
2076 }
2077 if (IsEnvDefined(env, prefix, "LambdaErrLimit", print))
2078 {
2079 SetLambdaErrLimit(GetEnvValue(env, prefix, "LambdaErrLimit", fLambdaErrLimit));
2080 rc = kTRUE;
2081 }
2082 if (IsEnvDefined(env, prefix, "LambdaCheckLimit", print))
2083 {
2084 SetLambdaCheckLimit(GetEnvValue(env, prefix, "LambdaCheckLimit", fLambdaCheckLimit));
2085 rc = kTRUE;
2086 }
2087 if (IsEnvDefined(env, prefix, "PheErrLimit", print))
2088 {
2089 SetPheErrLimit(GetEnvValue(env, prefix, "PheErrLimit", fPheErrLimit));
2090 rc = kTRUE;
2091 }
2092 // void SetPulserColor(const MCalibrationCam::PulserColor_t col) { fPulserColor = col; }
2093
2094 return rc;
2095}
Note: See TracBrowser for help on using the repository browser.