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

Last change on this file since 4781 was 4781, checked in by gaug, 20 years ago
*** empty log message ***
File size: 69.6 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 Int_t num = fPedestals->GetTotalEntries();
792
793 //
794 // RMS error set by PedCalcFromLoGain, 0 in case MPedCalcPedRun was used.
795 //
796 const Float_t prmserr = num>0 ? prms/TMath::Sqrt(2.*num) : ped.GetPedestalRmsError();
797
798 //
799 // set them in the calibration camera
800 //
801 if (cal.IsHiGainSaturation())
802 {
803 cal.SetPedestal(pedes * fNumLoGainSamples,
804 prms * fSqrtLoGainSamples,
805 prmserr * fSqrtLoGainSamples);
806 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples, aidx);
807 }
808 else
809 {
810 cal.SetPedestal(pedes * fNumHiGainSamples,
811 prms * fSqrtHiGainSamples,
812 prmserr * fSqrtHiGainSamples);
813 }
814
815}
816
817// ----------------------------------------------------------------------------------------------------
818//
819// Check fit results validity. Bad Pixels flags are set if:
820//
821// 1) Pixel has a mean smaller than fChargeLimit*PedRMS ( Flag: MBadPixelsPix::kChargeIsPedestal)
822// 2) Pixel has a mean error smaller than fChargeErrLimit ( Flag: MBadPixelsPix::kChargeErrNotValid)
823// 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error
824// ( Flag: MBadPixelsPix::kChargeRelErrNotValid)
825// 4) Pixel has a sigma bigger than its Pedestal RMS ( Flag: MBadPixelsPix::kChargeSigmaNotValid )
826//
827// Further returns if flags: MBadPixelsPix::kUnsuitableRun is set
828//
829// Calls MCalibrationChargePix::CalcReducedSigma() and sets flag: MBadPixelsPix::kChargeIsPedestal
830// and returns kFALSE if not succesful.
831//
832// Calls MCalibrationChargePix::CalcFFactor() and sets flag: MBadPixelsPix::kDeviatingNumPhes)
833// and returns kFALSE if not succesful.
834//
835// Calls MCalibrationChargePix::CalcConvFFactor()and sets flag: MBadPixelsPix::kDeviatingNumPhes)
836// and returns kFALSE if not succesful.
837//
838Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what)
839{
840
841 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
842 return kFALSE;
843
844 if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
845 {
846 *fLog << warn << GetDescriptor()
847 << Form(": Fitted Charge: %5.2f is smaller than %2.1f",cal.GetMean(),fChargeLimit)
848 << Form(" Pedestal RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
849 bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
850 }
851
852 if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr())
853 {
854 *fLog << warn << GetDescriptor()
855 << Form(": Fitted Charge: %4.2f is smaller than %2.1f",cal.GetMean(),fChargeRelErrLimit)
856 << Form(" times its error: %4.2f in %s%4i",cal.GetMeanErr(),what,cal.GetPixId()) << endl;
857 bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
858 }
859
860 if (cal.GetSigma() < cal.GetPedRms())
861 {
862 *fLog << warn << GetDescriptor()
863 << Form(": Sigma of Fitted Charge: %6.2f is smaller than",cal.GetSigma())
864 << Form(" Ped. RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
865 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
866 return kFALSE;
867 }
868
869 if (!cal.CalcReducedSigma())
870 {
871 *fLog << warn << GetDescriptor()
872 << Form(": Could not calculate the reduced sigma in %s: ",what)
873 << Form(" %4i",cal.GetPixId())
874 << endl;
875 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
876 return kFALSE;
877 }
878
879 if (!cal.CalcFFactor())
880 {
881 *fLog << warn << GetDescriptor()
882 << Form(": Could not calculate the F-Factor in %s: ",what)
883 << Form(" %4i",cal.GetPixId())
884 << endl;
885 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
886 return kFALSE;
887 }
888
889 if (!cal.CalcConvFFactor())
890 {
891 *fLog << warn << GetDescriptor()
892 << Form(": Could not calculate the Conv. FADC counts to Phes in %s: ",what)
893 << Form(" %4i",cal.GetPixId())
894 << endl;
895 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
896 return kFALSE;
897 }
898
899 return kTRUE;
900}
901
902
903
904// -----------------------------------------------------------------------------------
905//
906// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
907// - MBadPixelsPix::kChargeIsPedestal
908// - MBadPixelsPix::kChargeRelErrNotValid
909// - MBadPixelsPix::kMeanTimeInFirstBin
910// - MBadPixelsPix::kMeanTimeInLast2Bins
911// - MBadPixelsPix::kDeviatingNumPhes
912//
913// - Call MCalibrationPix::SetExcluded() for the bad pixels
914//
915// Sets pixel to MBadPixelsPix::kUnreliableRun, if one of the following flags is set:
916// - MBadPixelsPix::kChargeSigmaNotValid
917//
918void MCalibrationChargeCalc::FinalizeBadPixels()
919{
920
921 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
922 {
923
924 MBadPixelsPix &bad = (*fBadPixels)[i];
925 MCalibrationPix &pix = (*fCam)[i];
926
927 if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal))
928 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
929
930 if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid ))
931 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
932
933 if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ))
934 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
935
936 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ))
937 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
938
939 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ))
940 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
941
942 if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes ))
943 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
944
945 if (bad.IsUnsuitable( MBadPixelsPix::kUnsuitableRun ))
946 pix.SetExcluded();
947
948 if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ))
949 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
950 }
951}
952
953// ------------------------------------------------------------------------
954//
955//
956// First loop: Calculate a mean and mean RMS of photo-electrons per area index
957// Include only pixels which are not MBadPixelsPix::kUnsuitableRun nor
958// MBadPixelsPix::kChargeSigmaNotValid (see FinalizeBadPixels()) and set
959// MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case.
960//
961// Second loop: Get mean number of photo-electrons and its RMS including
962// only pixels with flag MCalibrationChargePix::IsFFactorMethodValid()
963// and further exclude those deviating by more than fPheErrLimit mean
964// sigmas from the mean (obtained in first loop). Set
965// MBadPixelsPix::kDeviatingNumPhes if excluded.
966//
967// For the suitable pixels with flag MBadPixelsPix::kChargeSigmaNotValid
968// set the number of photo-electrons as the mean number of photo-electrons
969// calculated in that area index.
970//
971// Set weighted mean and variance of photo-electrons per area index in:
972// average area pixels of MCalibrationChargeCam (obtained from:
973// MCalibrationChargeCam::GetAverageArea() )
974//
975// Set weighted mean and variance of photo-electrons per sector in:
976// average sector pixels of MCalibrationChargeCam (obtained from:
977// MCalibrationChargeCam::GetAverageSector() )
978//
979//
980// Third loop: Set mean number of photo-electrons and its RMS in the pixels
981// only excluded as: MBadPixelsPix::kChargeSigmaNotValid
982//
983Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod()
984{
985
986 const UInt_t npixels = fGeom->GetNumPixels();
987 const UInt_t nareas = fGeom->GetNumAreas();
988 const UInt_t nsectors = fGeom->GetNumSectors();
989
990 Float_t lowlim [nareas];
991 Float_t upplim [nareas];
992 Double_t areavars [nareas];
993 Double_t areaweights[nareas], sectorweights [nsectors];
994 Double_t areaphes [nareas], sectorphes [nsectors];
995 Int_t numareavalid[nareas], numsectorvalid[nsectors];
996
997 memset(lowlim ,0, nareas * sizeof(Float_t));
998 memset(upplim ,0, nareas * sizeof(Float_t));
999 memset(areaphes ,0, nareas * sizeof(Double_t));
1000 memset(areavars ,0, nareas * sizeof(Double_t));
1001 memset(areaweights ,0, nareas * sizeof(Double_t));
1002 memset(numareavalid ,0, nareas * sizeof(Int_t ));
1003 memset(sectorweights ,0, nsectors * sizeof(Double_t));
1004 memset(sectorphes ,0, nsectors * sizeof(Double_t));
1005 memset(numsectorvalid,0, nsectors * sizeof(Int_t ));
1006
1007 //
1008 // First loop: Get mean number of photo-electrons and the RMS
1009 // The loop is only to recognize later pixels with very deviating numbers
1010 //
1011 MHCamera camphes(*fGeom,"Camphes","Phes in Camera");
1012
1013 for (UInt_t i=0; i<npixels; i++)
1014 {
1015
1016 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam) [i];
1017 MBadPixelsPix &bad = (*fBadPixels)[i];
1018
1019 if (!pix.IsFFactorMethodValid())
1020 continue;
1021
1022 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1023 {
1024 pix.SetFFactorMethodValid(kFALSE);
1025 continue;
1026 }
1027
1028 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1029 continue;
1030
1031 const Float_t nphe = pix.GetPheFFactorMethod();
1032 const Int_t aidx = (*fGeom)[i].GetAidx();
1033
1034 camphes.Fill(i,nphe);
1035 camphes.SetUsed(i);
1036
1037 areaphes [aidx] += nphe;
1038 areavars [aidx] += nphe*nphe;
1039 numareavalid[aidx] ++;
1040 }
1041
1042 for (UInt_t i=0; i<nareas; i++)
1043 {
1044 if (numareavalid[i] == 0)
1045 {
1046 *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
1047 << "in area index: " << i << endl;
1048 continue;
1049 }
1050
1051 if (numareavalid[i] == 1)
1052 areavars[i] = 0.;
1053 else if (numareavalid[i] == 0)
1054 {
1055 areaphes[i] = -1.;
1056 areaweights[i] = -1.;
1057 }
1058 else
1059 {
1060 areavars[i] = (areavars[i] - areaphes[i]*areaphes[i]/numareavalid[i]) / (numareavalid[i]-1);
1061 areaphes[i] = areaphes[i] / numareavalid[i];
1062 }
1063
1064 if (areavars[i] < 0.)
1065 {
1066 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of photo-electrons found "
1067 << "in area index: " << i << endl;
1068 continue;
1069 }
1070
1071 lowlim [i] = areaphes[i] - fPheErrLimit*TMath::Sqrt(areavars[i]);
1072 upplim [i] = areaphes[i] + fPheErrLimit*TMath::Sqrt(areavars[i]);
1073
1074 TArrayI area(1);
1075 area[0] = i;
1076
1077 TH1D *hist = camphes.ProjectionS(TArrayI(),area,"_py",100);
1078 hist->Fit("gaus","Q");
1079 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1080 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1081 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1082
1083 if (IsDebug())
1084 camphes.DrawClone();
1085
1086 if (ndf < 2)
1087 {
1088 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1089 << "in the camera with area index: " << i << endl;
1090 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1091 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1092 delete hist;
1093 continue;
1094 }
1095
1096 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1097
1098 if (prob < 0.001)
1099 {
1100 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons "
1101 << "in the camera with area index: " << i << endl;
1102 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1103 << " is smaller than 0.001 " << endl;
1104 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1105 delete hist;
1106 continue;
1107 }
1108
1109 *fLog << inf << GetDescriptor() << ": Mean number of photo-electrons "
1110 << "with area idx " << i << ": "
1111 << Form("%7.2f+-%6.2f",mean,sigma) << endl;
1112
1113 lowlim [i] = mean - fPheErrLimit*sigma;
1114 upplim [i] = mean + fPheErrLimit*sigma;
1115
1116 delete hist;
1117 }
1118
1119 *fLog << endl;
1120
1121 memset(numareavalid,0,nareas*sizeof(Int_t));
1122 memset(areaphes ,0,nareas*sizeof(Double_t));
1123 memset(areavars ,0,nareas*sizeof(Double_t));
1124
1125 //
1126 // Second loop: Get mean number of photo-electrons and its RMS excluding
1127 // pixels deviating by more than fPheErrLimit sigma.
1128 // Set the conversion factor FADC counts to photo-electrons
1129 //
1130 for (UInt_t i=0; i<npixels; i++)
1131 {
1132
1133 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1134
1135 if (!pix.IsFFactorMethodValid())
1136 continue;
1137
1138 MBadPixelsPix &bad = (*fBadPixels)[i];
1139
1140 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1141 continue;
1142
1143 const Float_t nvar = pix.GetPheFFactorMethodVar();
1144
1145 if (nvar <= 0.)
1146 {
1147 pix.SetFFactorMethodValid(kFALSE);
1148 continue;
1149 }
1150
1151 const Int_t aidx = (*fGeom)[i].GetAidx();
1152 const Int_t sector = (*fGeom)[i].GetSector();
1153 const Float_t area = (*fGeom)[i].GetA();
1154 const Float_t nphe = pix.GetPheFFactorMethod();
1155
1156 if ( nphe < lowlim[aidx] || nphe > upplim[aidx] )
1157 {
1158 *fLog << warn << GetDescriptor() << ": Number of phes: "
1159 << Form("%7.2f out of %3.1f sigma limit: ",nphe,fPheErrLimit)
1160 << Form("[%7.2f,%7.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
1161 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1162 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1163 pix.SetFFactorMethodValid(kFALSE);
1164 continue;
1165 }
1166
1167 areaweights [aidx] += nphe*nphe;
1168 areaphes [aidx] += nphe;
1169 numareavalid [aidx] ++;
1170
1171 if (aidx == 0)
1172 fNumInnerFFactorMethodUsed++;
1173
1174 sectorweights [sector] += nphe*nphe/area/area;
1175 sectorphes [sector] += nphe/area;
1176 numsectorvalid[sector] ++;
1177 }
1178
1179 *fLog << endl;
1180
1181 for (UInt_t aidx=0; aidx<nareas; aidx++)
1182 {
1183
1184 MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
1185
1186 if (numareavalid[aidx] == 1)
1187 areaweights[aidx] = 0.;
1188 else if (numareavalid[aidx] == 0)
1189 {
1190 areaphes[aidx] = -1.;
1191 areaweights[aidx] = -1.;
1192 }
1193 else
1194 {
1195 areaweights[aidx] = (areaweights[aidx] - areaphes[aidx]*areaphes[aidx]/numareavalid[aidx])
1196 / (numareavalid[aidx]-1);
1197 areaphes[aidx] /= numareavalid[aidx];
1198 }
1199
1200 if (areaweights[aidx] < 0. || areaphes[aidx] <= 0.)
1201 {
1202 *fLog << warn << GetDescriptor()
1203 << ": Mean number phes from area index " << aidx << " could not be calculated: "
1204 << " Mean: " << areaphes[aidx]
1205 << " Variance: " << areaweights[aidx] << endl;
1206 apix.SetFFactorMethodValid(kFALSE);
1207 continue;
1208 }
1209
1210 *fLog << inf << GetDescriptor()
1211 << ": Average total number phes in area idx " << aidx << ": "
1212 << Form("%7.2f%s%6.2f",areaphes[aidx]," +- ",TMath::Sqrt(areaweights[aidx])) << endl;
1213
1214 apix.SetPheFFactorMethod ( areaphes[aidx] );
1215 apix.SetPheFFactorMethodVar( areaweights[aidx] / numareavalid[aidx] );
1216 apix.SetFFactorMethodValid ( kTRUE );
1217
1218 }
1219
1220 *fLog << endl;
1221
1222 for (UInt_t sector=0; sector<nsectors; sector++)
1223 {
1224
1225 if (numsectorvalid[sector] == 1)
1226 sectorweights[sector] = 0.;
1227 else if (numsectorvalid[sector] == 0)
1228 {
1229 sectorphes[sector] = -1.;
1230 sectorweights[sector] = -1.;
1231 }
1232 else
1233 {
1234 sectorweights[sector] = (sectorweights[sector]
1235 - sectorphes[sector]*sectorphes[sector]/numsectorvalid[sector]
1236 )
1237 / (numsectorvalid[sector]-1.);
1238 sectorphes[sector] /= numsectorvalid[sector];
1239 }
1240
1241 MCalibrationChargePix &spix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
1242
1243 if (sectorweights[sector] < 0. || sectorphes[sector] <= 0.)
1244 {
1245 *fLog << warn << GetDescriptor()
1246 <<": Mean number phes per area for sector " << sector << " could not be calculated: "
1247 << " Mean: " << sectorphes[sector]
1248 << " Variance: " << sectorweights[sector] << endl;
1249 spix.SetFFactorMethodValid(kFALSE);
1250 continue;
1251 }
1252
1253 *fLog << inf << GetDescriptor()
1254 << ": Average number phes per area in sector " << sector << ": "
1255 << Form("%5.2f+-%4.2f [phe/mm^2]",sectorphes[sector],TMath::Sqrt(sectorweights[sector]))
1256 << endl;
1257
1258 spix.SetPheFFactorMethod ( sectorphes[sector] );
1259 spix.SetPheFFactorMethodVar( sectorweights[sector] / numsectorvalid[sector]);
1260 spix.SetFFactorMethodValid ( kTRUE );
1261
1262 }
1263
1264 //
1265 // Third loop: Set mean number of photo-electrons and its RMS in the pixels
1266 // only excluded as: MBadPixelsPix::kChargeSigmaNotValid
1267 //
1268 for (UInt_t i=0; i<npixels; i++)
1269 {
1270
1271 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1272 MBadPixelsPix &bad = (*fBadPixels)[i];
1273
1274 if (!pix.IsFFactorMethodValid())
1275 continue;
1276
1277 const Int_t aidx = (*fGeom)[i].GetAidx();
1278
1279 if (bad.IsUncalibrated(MBadPixelsPix::kChargeSigmaNotValid))
1280 {
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 Float_t lowlim [nareas];
1559 Float_t upplim [nareas];
1560 Double_t avffactorphotons [nareas];
1561 Double_t avffactorphotvar [nareas];
1562 Int_t numffactor [nareas];
1563
1564 memset(lowlim ,0, nareas * sizeof(Float_t));
1565 memset(upplim ,0, nareas * sizeof(Float_t));
1566 memset(avffactorphotons,0, nareas * sizeof(Double_t));
1567 memset(avffactorphotvar,0, nareas * sizeof(Double_t));
1568 memset(numffactor ,0, nareas * sizeof(Int_t));
1569
1570 const UInt_t npixels = fGeom->GetNumPixels();
1571
1572 MHCamera camffactor(*fGeom,"Camffactor","F-Factor in Camera");
1573
1574 for (UInt_t i=0; i<npixels; i++)
1575 {
1576
1577 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1578 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1579 MBadPixelsPix &bad = (*fBadPixels)[i];
1580
1581 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1582 continue;
1583
1584 const Float_t photons = avphotons / fGeom->GetPixRatio(i);
1585 const Float_t qe = pix.GetPheFFactorMethod() / photons ;
1586
1587 const Float_t qerelvar = avphotrelvar + pix.GetPheFFactorMethodRelVar();
1588
1589 qepix.SetQEFFactor ( qe , fPulserColor );
1590 qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor );
1591 qepix.SetFFactorMethodValid( kTRUE , fPulserColor );
1592
1593 if (!qepix.UpdateFFactorMethod())
1594 *fLog << warn << GetDescriptor()
1595 << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl;
1596
1597 //
1598 // The following pixels are those with deviating sigma, but otherwise OK,
1599 // probably those with stars during the pedestal run, but not the cal. run.
1600 //
1601 if (!pix.CalcMeanFFactor( photons , avphotrelvar ))
1602 {
1603 (*fBadPixels)[i].SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1604 (*fBadPixels)[i].SetUnsuitable ( MBadPixelsPix::kUnreliableRun );
1605 continue;
1606 }
1607
1608 const Int_t aidx = (*fGeom)[i].GetAidx();
1609 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1610
1611 camffactor.Fill(i,ffactor);
1612 camffactor.SetUsed(i);
1613
1614 avffactorphotons[aidx] += ffactor;
1615 avffactorphotvar[aidx] += ffactor*ffactor;
1616 numffactor[aidx]++;
1617 }
1618
1619 for (UInt_t i=0; i<nareas; i++)
1620 {
1621
1622 if (numffactor[i] == 0)
1623 {
1624 *fLog << warn << GetDescriptor() << ": No pixels with valid total F-Factor found "
1625 << "in area index: " << i << endl;
1626 continue;
1627 }
1628
1629 avffactorphotvar[i] = (avffactorphotvar[i] - avffactorphotons[i]*avffactorphotons[i]/numffactor[i]) / (numffactor[i]-1.);
1630 avffactorphotons[i] = avffactorphotons[i] / numffactor[i];
1631
1632 if (avffactorphotvar[i] < 0.)
1633 {
1634 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of total F-Factor found "
1635 << "in area index: " << i << endl;
1636 continue;
1637 }
1638
1639 lowlim [i] = 1.1; // Lowest known F-Factor of a PMT
1640 upplim [i] = avffactorphotons[i] + fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);
1641
1642 TArrayI area(1);
1643 area[0] = i;
1644
1645 TH1D *hist = camffactor.ProjectionS(TArrayI(),area,"_py",100);
1646 hist->Fit("gaus","Q");
1647 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
1648 const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
1649 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
1650
1651 if (IsDebug())
1652 camffactor.DrawClone();
1653
1654 if (ndf < 2)
1655 {
1656 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1657 << "in the camera with area index: " << i << endl;
1658 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
1659 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1660 delete hist;
1661 continue;
1662 }
1663
1664 const Double_t prob = hist->GetFunction("gaus")->GetProb();
1665
1666 if (prob < 0.001)
1667 {
1668 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor "
1669 << "in the camera with area index: " << i << endl;
1670 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
1671 << " is smaller than 0.001 " << endl;
1672 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
1673 delete hist;
1674 continue;
1675 }
1676
1677 *fLog << inf << GetDescriptor() << ": Mean F-Factor "
1678 << "with area index #" << i << ": "
1679 << Form("%4.2f+-%4.2f",mean,sigma) << endl;
1680
1681 lowlim [i] = 1.1;
1682 upplim [i] = mean + fFFactorErrLimit*sigma;
1683
1684 delete hist;
1685 }
1686
1687 *fLog << endl;
1688
1689 for (UInt_t i=0; i<npixels; i++)
1690 {
1691
1692 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1693 MBadPixelsPix &bad = (*fBadPixels)[i];
1694
1695 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1696 continue;
1697
1698 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1699 const Int_t aidx = (*fGeom)[i].GetAidx();
1700
1701 if ( ffactor < lowlim[aidx] || ffactor > upplim[aidx] )
1702 {
1703 *fLog << warn << GetDescriptor() << ": Overall F-Factor "
1704 << Form("%5.2f",ffactor) << " out of range ["
1705 << Form("%5.2f,%5.2f",lowlim[aidx],upplim[aidx]) << "] pixel " << i << endl;
1706
1707 bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1708 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1709 }
1710 }
1711
1712 for (UInt_t i=0; i<npixels; i++)
1713 {
1714
1715 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1716 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1717 MBadPixelsPix &bad = (*fBadPixels)[i];
1718
1719 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1720 {
1721 qepix.SetFFactorMethodValid(kFALSE,fPulserColor);
1722 pix.SetFFactorMethodValid(kFALSE);
1723 pix.SetExcluded();
1724 continue;
1725 }
1726 }
1727}
1728
1729
1730// ------------------------------------------------------------------------
1731//
1732// Loop over pixels:
1733//
1734// - Continue, if not MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() and set:
1735// MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor)
1736//
1737// - Calculate the quantum efficiency with the formula:
1738//
1739// QE = Num.Phes / MCalibrationChargeBlindPix::GetFluxInsidePlexiglass()
1740// / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE()
1741//
1742// - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor );
1743// - Set Variance of QE in MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor );
1744// - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor)
1745//
1746// - Call MCalibrationQEPix::UpdateBlindPixelMethod()
1747//
1748void MCalibrationChargeCalc::FinalizeBlindPixelQECam()
1749{
1750
1751 const UInt_t npixels = fGeom->GetNumPixels();
1752
1753 //
1754 // Set the results in the MCalibrationChargeCam
1755 //
1756 if (fBlindPixel)
1757 {
1758 if (fBlindPixel->IsFluxInsidePlexiglassAvailable())
1759 {
1760
1761 const Float_t photons = fBlindPixel->GetFluxInsidePlexiglass() * (*fGeom)[0].GetA()
1762 / fQECam->GetPlexiglassQE();
1763 fCam->SetNumPhotonsBlindPixelMethod(photons);
1764
1765 const Float_t photrelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1766 + fQECam->GetPlexiglassQERelVar();
1767 if (photrelvar > 0.)
1768 fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1769 }
1770 }
1771 //
1772 // With the knowledge of the overall photon flux, calculate the
1773 // quantum efficiencies after the Blind Pixel and PIN Diode method
1774 //
1775 for (UInt_t i=0; i<npixels; i++)
1776 {
1777
1778 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1779
1780 if (!fBlindPixel)
1781 {
1782 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1783 continue;
1784 }
1785
1786 if (!fBlindPixel->IsFluxInsidePlexiglassAvailable())
1787 {
1788 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1789 continue;
1790 }
1791
1792 MBadPixelsPix &bad = (*fBadPixels)[i];
1793 if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1794 {
1795 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1796 continue;
1797 }
1798
1799 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1800 MGeomPix &geo = (*fGeom)[i];
1801
1802 const Float_t qe = pix.GetPheFFactorMethod()
1803 / fBlindPixel->GetFluxInsidePlexiglass()
1804 / geo.GetA()
1805 * fQECam->GetPlexiglassQE();
1806
1807 const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1808 + fQECam->GetPlexiglassQERelVar()
1809 + pix.GetPheFFactorMethodRelVar();
1810
1811 qepix.SetQEBlindPixel ( qe , fPulserColor );
1812 qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor );
1813 qepix.SetBlindPixelMethodValid( kTRUE , fPulserColor );
1814
1815 if (!qepix.UpdateBlindPixelMethod())
1816 *fLog << warn << GetDescriptor()
1817 << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl;
1818 }
1819}
1820
1821// ------------------------------------------------------------------------
1822//
1823// Loop over pixels:
1824//
1825// - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set:
1826// MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor)
1827//
1828// - Calculate the quantum efficiency with the formula:
1829//
1830// QE = Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA()
1831//
1832// - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor );
1833// - Set Variance of QE in MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor );
1834// - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor)
1835//
1836// - Call MCalibrationQEPix::UpdatePINDiodeMethod()
1837//
1838void MCalibrationChargeCalc::FinalizePINDiodeQECam()
1839{
1840
1841 const UInt_t npixels = fGeom->GetNumPixels();
1842
1843 //
1844 // With the knowledge of the overall photon flux, calculate the
1845 // quantum efficiencies after the PIN Diode method
1846 //
1847 for (UInt_t i=0; i<npixels; i++)
1848 {
1849
1850 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1851
1852 if (!fPINDiode)
1853 {
1854 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1855 continue;
1856 }
1857
1858 if (!fPINDiode->IsFluxOutsidePlexiglassAvailable())
1859 {
1860 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1861 continue;
1862 }
1863
1864 MBadPixelsPix &bad = (*fBadPixels)[i];
1865
1866 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1867 {
1868 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1869 continue;
1870 }
1871
1872 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1873 MGeomPix &geo = (*fGeom)[i];
1874
1875 const Float_t qe = pix.GetPheFFactorMethod()
1876 / fPINDiode->GetFluxOutsidePlexiglass()
1877 / geo.GetA();
1878
1879 const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar();
1880
1881 qepix.SetQEPINDiode ( qe , fPulserColor );
1882 qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor );
1883 qepix.SetPINDiodeMethodValid( kTRUE , fPulserColor );
1884
1885 if (!qepix.UpdatePINDiodeMethod())
1886 *fLog << warn << GetDescriptor()
1887 << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl;
1888 }
1889}
1890
1891// -----------------------------------------------------------------------------------------------
1892//
1893// - Print out statistics about BadPixels of type UnsuitableType_t
1894// - store numbers of bad pixels of each type in fCam
1895//
1896void MCalibrationChargeCalc::FinalizeUnsuitablePixels()
1897{
1898
1899 *fLog << inf << endl;
1900 *fLog << GetDescriptor() << ": Charge Calibration status:" << endl;
1901 *fLog << dec << setfill(' ');
1902
1903 const Int_t nareas = fGeom->GetNumAreas();
1904
1905 Int_t counts[nareas];
1906 memset(counts,0,nareas*sizeof(Int_t));
1907
1908 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1909 {
1910 MBadPixelsPix &bad = (*fBadPixels)[i];
1911 if (!bad.IsBad())
1912 {
1913 const Int_t aidx = (*fGeom)[i].GetAidx();
1914 counts[aidx]++;
1915 }
1916 }
1917
1918 if (fGeom->InheritsFrom("MGeomCamMagic"))
1919 *fLog << " " << setw(7) << "Successfully calibrated Pixels: "
1920 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1921
1922 memset(counts,0,nareas*sizeof(Int_t));
1923
1924 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1925 {
1926 MBadPixelsPix &bad = (*fBadPixels)[i];
1927 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1928 {
1929 const Int_t aidx = (*fGeom)[i].GetAidx();
1930 counts[aidx]++;
1931 }
1932 }
1933
1934 for (Int_t aidx=0; aidx<nareas; aidx++)
1935 fCam->SetNumUnsuitable(counts[aidx], aidx);
1936
1937 if (fGeom->InheritsFrom("MGeomCamMagic"))
1938 *fLog << " " << setw(7) << "Uncalibrated Pixels: "
1939 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1940
1941 memset(counts,0,nareas*sizeof(Int_t));
1942
1943 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1944 {
1945 MBadPixelsPix &bad = (*fBadPixels)[i];
1946 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
1947 {
1948 const Int_t aidx = (*fGeom)[i].GetAidx();
1949 counts[aidx]++;
1950 }
1951 }
1952
1953 for (Int_t aidx=0; aidx<nareas; aidx++)
1954 fCam->SetNumUnreliable(counts[aidx], aidx);
1955
1956 *fLog << " " << setw(7) << "Unreliable Pixels: "
1957 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1958
1959}
1960
1961// -----------------------------------------------------------------------------------------------
1962//
1963// Print out statistics about BadPixels of type UncalibratedType_t
1964//
1965void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
1966{
1967
1968 UInt_t countinner = 0;
1969 UInt_t countouter = 0;
1970
1971 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1972 {
1973 MBadPixelsPix &bad = (*fBadPixels)[i];
1974 if (bad.IsUncalibrated(typ))
1975 {
1976 if (fGeom->GetPixRatio(i) == 1.)
1977 countinner++;
1978 else
1979 countouter++;
1980 }
1981 }
1982
1983 *fLog << " " << setw(7) << text
1984 << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
1985}
1986
1987// --------------------------------------------------------------------------
1988//
1989// Set the path for output file
1990//
1991void MCalibrationChargeCalc::SetOutputPath(TString path)
1992{
1993 fOutputPath = path;
1994 if (fOutputPath.EndsWith("/"))
1995 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
1996}
1997
1998// --------------------------------------------------------------------------
1999//
2000// Set the output file
2001//
2002void MCalibrationChargeCalc::SetOutputFile(TString file)
2003{
2004 fOutputFile = file;
2005}
2006
2007// --------------------------------------------------------------------------
2008//
2009// Get the output file
2010//
2011const char* MCalibrationChargeCalc::GetOutputFile()
2012{
2013 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
2014}
2015
2016// --------------------------------------------------------------------------
2017//
2018// Read the environment for the following data members:
2019// - fChargeLimit
2020// - fChargeErrLimit
2021// - fChargeRelErrLimit
2022// - fDebug
2023// - fFFactorErrLimit
2024// - fLambdaErrLimit
2025// - fLambdaCheckErrLimit
2026// - fPheErrLimit
2027//
2028Int_t MCalibrationChargeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
2029{
2030
2031 Bool_t rc = kFALSE;
2032 if (IsEnvDefined(env, prefix, "ChargeLimit", print))
2033 {
2034 SetChargeLimit(GetEnvValue(env, prefix, "ChargeLimit", fChargeLimit));
2035 rc = kTRUE;
2036 }
2037 if (IsEnvDefined(env, prefix, "ChargeErrLimit", print))
2038 {
2039 SetChargeErrLimit(GetEnvValue(env, prefix, "ChargeErrLimit", fChargeErrLimit));
2040 rc = kTRUE;
2041 }
2042 if (IsEnvDefined(env, prefix, "ChargeRelErrLimit", print))
2043 {
2044 SetChargeRelErrLimit(GetEnvValue(env, prefix, "ChargeRelErrLimit", fChargeRelErrLimit));
2045 rc = kTRUE;
2046 }
2047 if (IsEnvDefined(env, prefix, "Debug", print))
2048 {
2049 SetDebug(GetEnvValue(env, prefix, "Debug", IsDebug()));
2050 rc = kTRUE;
2051 }
2052 if (IsEnvDefined(env, prefix, "FFactorErrLimit", print))
2053 {
2054 SetFFactorErrLimit(GetEnvValue(env, prefix, "FFactorErrLimit", fFFactorErrLimit));
2055 rc = kTRUE;
2056 }
2057 if (IsEnvDefined(env, prefix, "LambdaErrLimit", print))
2058 {
2059 SetLambdaErrLimit(GetEnvValue(env, prefix, "LambdaErrLimit", fLambdaErrLimit));
2060 rc = kTRUE;
2061 }
2062 if (IsEnvDefined(env, prefix, "LambdaCheckLimit", print))
2063 {
2064 SetLambdaCheckLimit(GetEnvValue(env, prefix, "LambdaCheckLimit", fLambdaCheckLimit));
2065 rc = kTRUE;
2066 }
2067 if (IsEnvDefined(env, prefix, "PheErrLimit", print))
2068 {
2069 SetPheErrLimit(GetEnvValue(env, prefix, "PheErrLimit", fPheErrLimit));
2070 rc = kTRUE;
2071 }
2072 // void SetPulserColor(const MCalibrationCam::PulserColor_t col) { fPulserColor = col; }
2073
2074 return rc;
2075}
Note: See TracBrowser for help on using the repository browser.