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

Last change on this file since 4231 was 4228, checked in by gaug, 21 years ago
*** empty log message ***
File size: 55.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// - FinalizePINDiode()
52// - FinalizeFFactorQECam()
53// - FinalizeBlindPixelQECam()
54// - FinalizePINDiodeQECam()
55//
56// Input Containers:
57// MCalibrationChargeCam
58// MCalibrationChargeBlindPix
59// MCalibrationChargePINDiode
60// MCalibrationQECam
61// MPedestalCam
62// MBadPixelsCam
63// MGeomCam
64// MTime
65//
66// Output Containers:
67// MCalibrationChargeCam
68// MCalibrationChargeBlindPix
69// MCalibrationChargePINDiode
70// MCalibrationQECam
71// MBadPixelsCam
72//
73//
74// Preliminary description of the calibration in photons (email from 12/02/04)
75//
76// Why calibrating in photons:
77// ===========================
78//
79// At the Barcelona meeting in 2002, we decided to calibrate the camera in
80// photons. This for the following reasons:
81//
82// * The physical quantity arriving at the camera are photons. This is
83// the direct physical information from the air shower. The photons
84// have a flux and a spectrum.
85//
86// * The photon fluxes depend mostly on the shower energy (with
87// corrections deriving from the observation conditions), while the photon
88// spectra depend mostly on the observation conditions: zenith angle,
89// quality of the air, also the impact parameter of the shower.
90//
91// * The photomultiplier, in turn, has different response properties
92// (quantum efficiencies) for photons of different colour. (Moreover,
93// different pixels have slightly different quantum efficiencies).
94// The resulting number of photo-electrons is then amplified (linearly)
95// with respect to the photo-electron flux.
96//
97// * In the ideal case, one would like to disentagle the effects
98// of the observation conditions from the primary particle energy (which
99// one likes to measure). To do so, one needs:
100//
101// 1) A reliable calibration relating the FADC counts to the photo-electron
102// flux -> This is accomplished with the F-Factor method.
103//
104// 2) A reliable calibration of the wavelength-dependent quantum efficiency
105// -> This is accomplished with the combination of the three methods,
106// together with QE-measurements performed by David in order to do
107// the interpolation.
108//
109// 3) A reliable calibration of the observation conditions. This means:
110// - Tracing the atmospheric conditions -> LIDAR
111// - Tracing the observation zenith angle -> Drive System
112//
113// 4) Some knowlegde about the impact parameter:
114// - This is the only part which cannot be accomplished well with a
115// single telescope. We would thus need to convolute the spectrum
116// over the distribution of impact parameters.
117//
118//
119// How an ideal calibration would look like:
120// =========================================
121//
122// We know from the combined PIN-Diode and Blind-Pixel Method the response of
123// each pixel to well-measured light fluxes in three representative
124// wavelengths (green, blue, UV). We also know the response to these light
125// fluxes in photo-electrons. Thus, we can derive:
126//
127// - conversion factors to photo-electrons
128// - conversion factors to photons in three wavelengths.
129//
130// Together with David's measurements and some MC-simulation, we should be
131// able to derive tables for typical Cherenkov-photon spectra - convoluted
132// with the impact parameters and depending on the athmospheric conditions
133// and the zenith angle (the "outer parameters").
134//
135// From these tables we can create "calibration tables" containing some
136// effective quantum efficiency depending on these outer parameters and which
137// are different for each pixel.
138//
139// In an ideal MCalibrate, one would thus have to convert first the FADC
140// slices to Photo-electrons and then, depending on the outer parameters,
141// look up the effective quantum efficiency and get the mean number of
142// photons which is then used for the further analysis.
143//
144// How the (first) MAGIC calibration should look like:
145// ===================================================
146//
147// For the moment, we have only one reliable calibration method, although
148// with very large systematic errors. This is the F-Factor method. Knowing
149// that the light is uniform over the whole camera (which I would not at all
150// guarantee in the case of the CT1 pulser), one could in principle already
151// perform a relative calibration of the quantum efficiencies in the UV.
152// However, the spread in QE at UV is about 10-15% (according to the plot
153// that Abelardo sent around last time. The spread in photo-electrons is 15%
154// for the inner pixels, but much larger (40%) for the outer ones.
155//
156// I'm not sure if we can already say that we have measured the relative
157// difference in quantum efficiency for the inner pixels and produce a first
158// QE-table for each pixel. To so, I would rather check in other wavelengths
159// (which we can do in about one-two weeks when the optical transmission of
160// the calibration trigger is installed).
161//
162// Thus, for the moment being, I would join Thomas proposal to calibrate in
163// photo-electrons and apply one stupid average quantum efficiency for all
164// pixels. This keeping in mind that we will have much preciser information
165// in about one to two weeks.
166//
167//
168// What MCalibrate should calculate and what should be stored:
169// ===========================================================
170//
171// It is clear that in the end, MCerPhotEvt will store photons.
172// MCalibrationCam stores the conversionfactors to photo-electrons and also
173// some tables of how to apply the conversion to photons, given the outer
174// parameters. This is not yet implemented and not even discussed.
175//
176// To start, I would suggest that we define the "average quantum efficiency"
177// (maybe something like 25+-3%) and apply them equally to all
178// photo-electrons. Later, this average factor can be easily replaced by a
179// pixel-dependent factor and later by a (pixel-dependent) table.
180//
181//
182//
183//////////////////////////////////////////////////////////////////////////////
184#include "MCalibrationChargeCalc.h"
185
186#include <TSystem.h>
187#include <TH1.h>
188
189#include "MLog.h"
190#include "MLogManip.h"
191
192#include "MParList.h"
193
194#include "MRawRunHeader.h"
195#include "MRawEvtPixelIter.h"
196
197#include "MGeomCam.h"
198#include "MGeomPix.h"
199
200#include "MPedestalCam.h"
201#include "MPedestalPix.h"
202
203#include "MCalibrationChargeCam.h"
204#include "MCalibrationChargePix.h"
205#include "MCalibrationChargePINDiode.h"
206#include "MCalibrationChargeBlindPix.h"
207
208#include "MExtractedSignalCam.h"
209#include "MExtractedSignalPix.h"
210#include "MExtractedSignalBlindPixel.h"
211#include "MExtractedSignalPINDiode.h"
212
213#include "MBadPixelsCam.h"
214#include "MBadPixelsPix.h"
215
216#include "MCalibrationQECam.h"
217#include "MCalibrationQEPix.h"
218
219#include "MCalibrationCam.h"
220
221ClassImp(MCalibrationChargeCalc);
222
223using namespace std;
224
225const Float_t MCalibrationChargeCalc::fgChargeLimit = 2.5;
226const Float_t MCalibrationChargeCalc::fgChargeErrLimit = 0.;
227const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit = 1.;
228const Float_t MCalibrationChargeCalc::fgLambdaErrLimit = 0.2;
229const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit = 0.5;
230const Float_t MCalibrationChargeCalc::fgPheErrLimit = 3.5;
231const Float_t MCalibrationChargeCalc::fgFFactorErrLimit = 3.;
232// --------------------------------------------------------------------------
233//
234// Default constructor.
235//
236// Sets all pointers to NULL
237//
238// Calls AddToBranchList for:
239// - MRawEvtData.fHiGainPixId
240// - MRawEvtData.fLoGainPixId
241// - MRawEvtData.fHiGainFadcSamples
242// - MRawEvtData.fLoGainFadcSamples
243//
244// Initializes:
245// - fChargeLimit to fgChargeLimit
246// - fChargeErrLimit to fgChargeErrLimit
247// - fChargeRelErrLimit to fgChargeRelErrLimit
248// - fFFactorErrLimit to fgFFactorErrLimit
249// - fLambdaCheckLimit to fgLambdaCheckLimit
250// - fLambdaErrLimit to fgLambdaErrLimit
251// - fPheErrLimit to fgPheErrLimit
252// - fPulserColor to MCalibrationCam::kCT1
253// - fOutputPath to "."
254// - fOutputFile to "ChargeCalibStat.txt"
255//
256// Calls:
257// - Clear()
258//
259MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
260 : fBadPixels(NULL), fCam(NULL), fBlindPixel(NULL), fPINDiode(NULL),
261 fQECam(NULL), fGeom(NULL), fPedestals(NULL), fEvtTime(NULL)
262{
263
264 fName = name ? name : "MCalibrationChargeCalc";
265 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
266
267 AddToBranchList("MRawEvtData.fHiGainPixId");
268 AddToBranchList("MRawEvtData.fLoGainPixId");
269 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
270 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
271
272 SetChargeLimit();
273 SetChargeErrLimit();
274 SetChargeRelErrLimit();
275 SetFFactorErrLimit();
276 SetLambdaCheckLimit();
277 SetLambdaErrLimit();
278 SetPheErrLimit();
279 SetPulserColor(MCalibrationCam::kNONE);
280 SetOutputPath();
281 SetOutputFile();
282
283 Clear();
284
285}
286
287// --------------------------------------------------------------------------
288//
289// Sets:
290// - all variables to 0.,
291// - all flags to kFALSE
292//
293void MCalibrationChargeCalc::Clear(const Option_t *o)
294{
295
296 fNumHiGainSamples = 0.;
297 fNumLoGainSamples = 0.;
298 fSqrtHiGainSamples = 0.;
299 fSqrtLoGainSamples = 0.;
300 SkipHiLoGainCalibration( kFALSE );
301}
302
303
304// -----------------------------------------------------------------------------------
305//
306// The following container are searched for and execution aborted if not in MParList:
307// - MPedestalCam
308//
309// The following containers are searched and created if they were not found:
310//
311// - MCalibrationQECam
312// - MBadPixelsCam
313//
314// The following output containers are only searched, but not created. If they
315// cannot be found, the corresponding calibration part is only skipped.
316//
317// - MTime
318//
319Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
320{
321
322 //
323 // Containers that have to be there.
324 //
325 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
326 if (!fPedestals)
327 {
328 *fLog << err << "MPedestalCam not found... aborting" << endl;
329 return kFALSE;
330 }
331
332 //
333 // Containers that are created in case that they are not there.
334 //
335 fQECam = (MCalibrationQECam*)pList->FindCreateObj("MCalibrationQECam");
336 if (!fQECam)
337 {
338 *fLog << err << "Cannot find nor create MCalibrationQECam... aborting" << endl;
339 return kFALSE;
340 }
341
342 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
343 if (!fBadPixels)
344 {
345 *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
346 return kFALSE;
347 }
348
349
350 fEvtTime = (MTime*)pList->FindObject("MTime");
351
352 //
353 // Check the pulser colour --> FIXME: this solution is only valid until the arrival of the DM's
354 //
355 if (fPulserColor == MCalibrationCam::kNONE)
356 {
357 *fLog << endl;
358 *fLog << err << GetDescriptor()
359 << ": No Pulser colour has been chosen. Since the installation of the IFAE pulser box,"
360 << " you HAVE to provide the LEDs colour, otherwise there is no calibration. " << endl;
361 *fLog << "See e.g. the macro calibration.C " << endl;
362 return kFALSE;
363 }
364
365 return kTRUE;
366}
367
368
369// --------------------------------------------------------------------------
370//
371// Search for the following input containers and abort if not existing:
372// - MGeomCam
373// - MCalibrationChargeCam
374//
375// Search for the following input containers and give a warning if not existing:
376// - MCalibrationChargeBlindPix
377// - MCalibrationChargePINDiode
378//
379// It retrieves the following variables from MCalibrationChargeCam:
380//
381// - fNumHiGainSamples
382// - fNumLoGainSamples
383//
384// It defines the PixId of every pixel in:
385//
386// - MCalibrationChargeCam
387// - MCalibrationQECam
388//
389// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
390//
391// - MCalibrationChargePix
392// - MCalibrationQEPix
393//
394// Sets the pulser colour and tests if it has not changed w.r.t. fPulserColor in:
395//
396// - MCalibrationChargeCam
397// - MCalibrationChargeBlindPix (if existing)
398// - MCalibrationChargePINDiode (if existing)
399//
400Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
401{
402
403 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
404 if (!fGeom)
405 {
406 *fLog << err << "No MGeomCam found... aborting." << endl;
407 return kFALSE;
408 }
409
410 fCam = (MCalibrationChargeCam*)pList->FindObject("MCalibrationChargeCam");
411 if (!fCam)
412 {
413 *fLog << err << "Cannot find MCalibrationChargeCam... aborting" << endl;
414 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationChargeCam before..." << endl;
415 return kFALSE;
416 }
417
418 //
419 // Optional Containers
420 //
421 fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindObject("MCalibrationChargeBlindPix");
422 if (!fBlindPixel)
423 {
424 *fLog << endl;
425 *fLog << warn << GetDescriptor()
426 << ": MCalibrationChargeBlindPix not found... no Blind Pixel method! " << endl;
427 }
428
429 fPINDiode = (MCalibrationChargePINDiode*)pList->FindObject("MCalibrationChargePINDiode");
430 if (!fPINDiode)
431 {
432 *fLog << endl;
433 *fLog << warn << GetDescriptor()
434 << ": MCalibrationChargePINDiode not found... no PIN Diode method! " << endl;
435 }
436
437
438 //
439 // Initialize the pulser colours
440 //
441 if (fCam->GetPulserColor() == MCalibrationCam::kNONE)
442 {
443 fCam->SetPulserColor( fPulserColor );
444
445 if (fBlindPixel)
446 fBlindPixel->SetColor( fPulserColor );
447
448 if (fPINDiode)
449 fPINDiode->SetColor( fPulserColor );
450 }
451
452 if (fPulserColor != fCam->GetPulserColor())
453 {
454 *fLog << err << GetDescriptor()
455 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeCam" << endl;
456 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
457 return kFALSE;
458 }
459
460 if (fBlindPixel)
461 if (fPulserColor != fBlindPixel->GetColor())
462 {
463 *fLog << err << GetDescriptor()
464 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeBlindPix." << endl;
465 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
466 return kFALSE;
467 }
468
469 if (fPINDiode)
470 if (fPulserColor != fPINDiode->GetColor())
471 {
472 *fLog << err << GetDescriptor()
473 << ": Pulser colour has changed w.r.t. last file in MCalibrationChargePINDiode." << endl;
474 *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
475 return kFALSE;
476 }
477
478
479 fNumHiGainSamples = fCam->GetNumHiGainFADCSlices();
480 fNumLoGainSamples = fCam->GetNumLoGainFADCSlices();
481 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
482 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
483
484 UInt_t npixels = fGeom->GetNumPixels();
485
486 for (UInt_t i=0; i<npixels; i++)
487 {
488
489 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam) [i];
490 MCalibrationQEPix &pqe = (MCalibrationQEPix&) (*fQECam)[i];
491 MBadPixelsPix &bad = (*fBadPixels)[i];
492
493 pix.SetPixId(i);
494 pqe.SetPixId(i);
495
496 if (bad.IsBad())
497 {
498 pix.SetExcluded();
499 pqe.SetExcluded();
500 continue;
501 }
502
503 }
504
505 return kTRUE;
506}
507
508// ----------------------------------------------------------------------------------
509//
510// Nothing to be done in Process, but have a look at MHCalibrationChargeCam, instead
511//
512Int_t MCalibrationChargeCalc::Process()
513{
514 return kTRUE;
515}
516
517// -----------------------------------------------------------------------
518//
519// Return if number of executions is null.
520//
521// First loop over pixels, average areas and sectors, call:
522// - FinalizePedestals()
523// - FinalizeCharges()
524// for every entry. Count number of valid pixels in loop and return kFALSE
525// if there are none (the "Michele check").
526//
527// Call FinalizeBadPixels()
528//
529// Call FinalizeFFactorMethod() (second and third loop over pixels and areas)
530//
531// Call FinalizeBlindPixel()
532// Call FinalizePINDiode()
533//
534// Call FinalizeFFactorQECam() (fourth loop over pixels and areas)
535// Call FinalizeBlindPixelQECam() (fifth loop over pixels and areas)
536// Call FinalizePINDiodeQECam() (sixth loop over pixels and areas)
537//
538// Call FinalizeUnsuitablePixels()
539//
540// Call MParContainer::SetReadyToSave() for fCam, fQECam, fBadPixels and
541// fBlindPixel and fPINDiode if they exist
542//
543// Print out some statistics
544//
545Int_t MCalibrationChargeCalc::PostProcess()
546{
547
548 if (GetNumExecutions()==0)
549 return kFALSE;
550
551 if (fPINDiode)
552 if (!fPINDiode->IsValid())
553 {
554 *fLog << warn << GetDescriptor()
555 << ": MCalibrationChargePINDiode is declared not valid... no PIN Diode method! " << endl;
556 fPINDiode = NULL;
557 }
558
559 if (fBlindPixel)
560 if (!fBlindPixel->IsValid())
561 {
562 *fLog << warn << GetDescriptor()
563 << ": MCalibrationChargeBlindPix is declared not valid... no Blind Pixel method! " << endl;
564 fBlindPixel = NULL;
565 }
566
567 //
568 // First loop over pixels, call FinalizePedestals and FinalizeCharges
569 //
570 Int_t nvalid = 0;
571
572 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
573 {
574
575 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[pixid];
576 //
577 // Check if the pixel has been excluded from the fits
578 //
579 if (pix.IsExcluded())
580 continue;
581
582 MPedestalPix &ped = (*fPedestals)[pixid];
583 MBadPixelsPix &bad = (*fBadPixels)[pixid];
584 const Int_t aidx = (*fGeom)[pixid].GetAidx();
585
586 FinalizePedestals(ped,pix,aidx);
587
588 if (FinalizeCharges(pix,bad))
589 nvalid++;
590 }
591
592 //
593 // The Michele check ...
594 //
595 if (nvalid == 0)
596 {
597 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
598 << "Did you forget to fill the histograms "
599 << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
600 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
601 << "instead of a calibration run " << endl;
602 return kFALSE;
603 }
604
605 for (UInt_t aidx=0; aidx<fGeom->GetNumAreas(); aidx++)
606 {
607
608 const MPedestalPix &ped = fPedestals->GetAverageArea(aidx);
609 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
610
611 FinalizePedestals(ped,pix,aidx);
612 FinalizeCharges(pix, fCam->GetAverageBadArea(aidx));
613 }
614
615 for (UInt_t sector=0; sector<fGeom->GetNumSectors(); sector++)
616 {
617
618 const MPedestalPix &ped = fPedestals->GetAverageSector(sector);
619
620 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
621 FinalizePedestals(ped,pix, 0);
622 FinalizeCharges(pix, fCam->GetAverageBadSector(sector));
623 }
624
625 //
626 // Finalize Bad Pixels
627 //
628 FinalizeBadPixels();
629
630 //
631 // Finalize F-Factor method
632 //
633 if (!FinalizeFFactorMethod())
634 {
635 *fLog << warn << "Could not calculate the photons flux from the F-Factor method " << endl;
636 fCam->SetFFactorMethodValid(kFALSE);
637 return kFALSE;
638 }
639 else
640 fCam->SetFFactorMethodValid(kTRUE);
641
642 //
643 // Finalize Blind Pixel
644 //
645 if (FinalizeBlindPixel())
646 fQECam->SetBlindPixelMethodValid(kTRUE);
647 else
648 fQECam->SetBlindPixelMethodValid(kFALSE);
649
650 //
651 // Finalize PIN Diode
652 //
653 if (FinalizePINDiode())
654 fQECam->SetPINDiodeMethodValid(kTRUE);
655 else
656 fQECam->SetPINDiodeMethodValid(kFALSE);
657
658 //
659 // Finalize QE Cam
660 //
661 FinalizeFFactorQECam();
662 FinalizeBlindPixelQECam();
663 FinalizePINDiodeQECam();
664
665 //
666 // Re-direct the output to an ascii-file from now on:
667 //
668 MLog asciilog;
669 asciilog.SetOutputFile(GetOutputFile(),kTRUE);
670 SetLogStream(&asciilog);
671 //
672 // Finalize calibration statistics
673 //
674 FinalizeUnsuitablePixels();
675
676 fCam ->SetReadyToSave();
677 fQECam ->SetReadyToSave();
678 fBadPixels->SetReadyToSave();
679
680 if (fBlindPixel)
681 fBlindPixel->SetReadyToSave();
682 if (fPINDiode)
683 fPINDiode->SetReadyToSave();
684
685 *fLog << inf << endl;
686 *fLog << GetDescriptor() << ": Fatal errors statistics:" << endl;
687
688 PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal,
689 Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS: "));
690 PrintUncalibrated(MBadPixelsPix::kChargeRelErrNotValid,
691 Form("%s%2.1f%s","Signal Error bigger than ",fChargeRelErrLimit," times Mean Signal: "));
692 PrintUncalibrated(MBadPixelsPix::kChargeSigmaNotValid,
693 "Signal Sigma smaller than Pedestal RMS: ");
694 PrintUncalibrated(MBadPixelsPix::kLoGainSaturation,
695 "Pixels with Low Gain Saturation: ");
696 PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin,
697 Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s): "));
698 PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins,
699 Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s): "));
700 PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes,
701 "Pixels with deviating number of phes: ");
702
703 *fLog << inf << endl;
704 *fLog << GetDescriptor() << ": Unreliable errors statistics:" << endl;
705
706 PrintUncalibrated(MBadPixelsPix::kHiGainOscillating,
707 "Pixels with changing Hi Gain signal over time: ");
708 PrintUncalibrated(MBadPixelsPix::kLoGainOscillating,
709 "Pixels with changing Lo Gain signal over time: ");
710 PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted,
711 "Pixels with unsuccesful Gauss fit to the Hi Gain: ");
712 PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted,
713 "Pixels with unsuccesful Gauss fit to the Lo Gain: ");
714
715 SetLogStream(&gLog);
716
717 return kTRUE;
718}
719
720// ----------------------------------------------------------------------------------
721//
722// Retrieves pedestal and pedestal RMS from MPedestalPix
723// Retrieves total entries from MPedestalCam
724// Sets pedestal*fNumHiGainSamples and pedestal*fNumLoGainSamples in MCalibrationChargePix
725// Sets pedRMS *fSqrtHiGainSamples and pedRMS *fSqrtLoGainSamples in MCalibrationChargePix
726//
727// If the flag MCalibrationPix::IsHiGainSaturation() is set, call also:
728// - MCalibrationChargePix::CalcLoGainPedestal()
729//
730void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal, const Int_t aidx)
731{
732
733 //
734 // get the pedestals
735 //
736 const Float_t pedes = ped.GetPedestal();
737 const Float_t prms = ped.GetPedestalRms();
738 const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
739
740 //
741 // set them in the calibration camera
742 //
743 if (cal.IsHiGainSaturation())
744 {
745 cal.SetPedestal(pedes* fNumLoGainSamples,
746 prms * fSqrtLoGainSamples,
747 prms * fNumLoGainSamples / num);
748 cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples, aidx);
749 }
750 else
751 {
752 cal.SetPedestal(pedes* fNumHiGainSamples,
753 prms * fSqrtHiGainSamples,
754 prms * fNumHiGainSamples / num);
755 }
756
757}
758
759// ----------------------------------------------------------------------------------------------------
760//
761// Check fit results validity. Bad Pixels flags are set if:
762//
763// 1) Pixel has a mean smaller than fChargeLimit*PedRMS ( Flag: MBadPixelsPix::kChargeIsPedestal)
764// 2) Pixel has a mean error smaller than fChargeErrLimit ( Flag: MBadPixelsPix::kChargeErrNotValid)
765// 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error
766// ( Flag: MBadPixelsPix::kChargeRelErrNotValid)
767// 4) Pixel has a sigma bigger than its Pedestal RMS ( Flag: MBadPixelsPix::kChargeSigmaNotValid )
768//
769// Further returns if flags: MBadPixelsPix::kUnsuitableRun is set
770//
771// Calls MCalibrationChargePix::CalcReducedSigma() and sets flag: MBadPixelsPix::kChargeIsPedestal
772// if not succesful.
773//
774// Calls MCalibrationChargePix::CalcFFactorMethod() and sets flag: MBadPixelsPix::kDeviatingNumPhes)
775// if not succesful.
776//
777Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad)
778{
779
780 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
781 return kFALSE;
782
783 if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
784 {
785 *fLog << warn << GetDescriptor() << ": Fitted Charge: " << cal.GetMean() << " is smaller than "
786 << fChargeLimit << " Pedestal RMS: " << cal.GetPedRms() << " in Pixel " << cal.GetPixId() << endl;
787 bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
788 }
789
790 if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr())
791 {
792 *fLog << warn << GetDescriptor() << ": Fitted Charge: " << cal.GetMean() << " is smaller than "
793 << fChargeRelErrLimit << "* its error: " << cal.GetMeanErr()
794 << " in Pixel " << cal.GetPixId() << endl;
795 bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
796 }
797
798 if (cal.GetSigma() < cal.GetPedRms())
799 {
800 *fLog << warn << GetDescriptor() << ": Sigma of Fitted Charge: " << cal.GetSigma()
801 << " smaller than Pedestal RMS: " << cal.GetPedRms() << " in Pixel " << cal.GetPixId() << endl;
802 bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
803 }
804
805 if (!cal.CalcReducedSigma())
806 {
807 *fLog << warn << GetDescriptor()
808 << ": Could not calculate reduced sigmas of pixel: " << cal.GetPixId() << endl;
809 bad.SetUncalibrated(MBadPixelsPix::kChargeIsPedestal);
810 return kFALSE;
811 }
812
813 if (!cal.CalcFFactorMethod())
814 {
815 *fLog << warn << GetDescriptor()
816 << ": Could not calculate F-Factor of pixel: " << cal.GetPixId() << endl;
817 bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
818 return kFALSE;
819 }
820
821 return kTRUE;
822}
823
824
825
826// -----------------------------------------------------------------------------------
827//
828// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
829// - MBadPixelsPix::kChargeIsPedestal
830// - MBadPixelsPix::kChargeRelErrNotValid
831// - MBadPixelsPix::kChargeSigmaNotValid
832// - MBadPixelsPix::kMeanTimeInFirstBin
833// - MBadPixelsPix::kMeanTimeInLast2Bins
834// - MBadPixelsPix::kDeviatingNumPhes
835//
836// - Call MCalibrationPix::SetExcluded() for the bad pixels
837//
838void MCalibrationChargeCalc::FinalizeBadPixels()
839{
840
841 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
842 {
843
844 MBadPixelsPix &bad = (*fBadPixels)[i];
845 MCalibrationPix &pix = (*fCam)[i];
846
847 if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal))
848 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
849
850 if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid ))
851 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
852
853 if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ))
854 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
855
856 if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ))
857 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
858
859 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ))
860 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
861
862 if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ))
863 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
864
865 if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes ))
866 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
867
868 if (bad.IsUnsuitable( MBadPixelsPix::kUnsuitableRun ))
869 pix.SetExcluded();
870
871 }
872}
873
874// ------------------------------------------------------------------------
875//
876//
877// First loop: Calculate a mean and mean RMS of photo-electrons per area index
878// Include only pixels which are not MBadPixelsPix::kUnsuitableRun or
879// MBadPixelsPix::kUnreliableRun (see FinalizeBadPixels()) and set
880// MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case.
881//
882// Second loop: Get weighted mean number of photo-electrons and its RMS including
883// only pixels with flag MCalibrationChargePix::IsFFactorMethodValid()
884// and further exclude those deviating by more than fPheErrLimit mean
885// sigmas from the mean (obtained in first loop). Set
886// MBadPixelsPix::kDeviatingNumPhes if excluded.
887//
888// Set weighted mean and variance of photo-electrons per area index in:
889// average area pixels of MCalibrationChargeCam (obtained from:
890// MCalibrationChargeCam::GetAverageArea() )
891//
892// Set weighted mean and variance of photo-electrons per sector in:
893// average sector pixels of MCalibrationChargeCam (obtained from:
894// MCalibrationChargeCam::GetAverageSector() )
895//
896Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod()
897{
898
899 const UInt_t npixels = fGeom->GetNumPixels();
900 const UInt_t nareas = fGeom->GetNumAreas();
901 const UInt_t nsectors = fGeom->GetNumSectors();
902
903 Float_t lowlim [nareas];
904 Float_t upplim [nareas];
905 Float_t areavars [nareas];
906 Float_t areaweights [nareas], sectorweights [nsectors];
907 Float_t areaphes [nareas], sectorphes [nsectors];
908 Int_t numareavalid[nareas], numsectorvalid[nsectors];
909
910 memset(lowlim ,0, nareas * sizeof(Float_t));
911 memset(upplim ,0, nareas * sizeof(Float_t));
912 memset(areaphes ,0, nareas * sizeof(Float_t));
913 memset(areavars ,0, nareas * sizeof(Float_t));
914 memset(areaweights ,0, nareas * sizeof(Float_t));
915 memset(numareavalid ,0, nareas * sizeof(Int_t ));
916 memset(sectorweights ,0, nsectors * sizeof(Float_t));
917 memset(sectorphes ,0, nsectors * sizeof(Float_t));
918 memset(numsectorvalid,0, nsectors * sizeof(Int_t ));
919
920 //
921 // First loop: Get mean number of photo-electrons and the RMS
922 // The loop is only to recognize later pixels with very deviating numbers
923 //
924 for (UInt_t i=0; i<npixels; i++)
925 {
926
927 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam) [i];
928 MBadPixelsPix &bad = (*fBadPixels)[i];
929
930 if (!pix.IsFFactorMethodValid())
931 continue;
932
933 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
934 {
935 pix.SetFFactorMethodValid(kFALSE);
936 continue;
937 }
938
939 // if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
940 // continue;
941
942 const Float_t nphe = pix.GetPheFFactorMethod();
943 const Float_t nvar = pix.GetPheFFactorMethod()*pix.GetPheFFactorMethod();
944 const Int_t aidx = (*fGeom)[i].GetAidx();
945
946 areaphes [aidx] += nphe;
947 areavars [aidx] += nvar;
948 numareavalid[aidx] ++;
949 }
950
951 for (UInt_t i=0; i<nareas; i++)
952 {
953 if (numareavalid[i] == 0)
954 {
955 *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
956 << "in area index: " << i << endl;
957 continue;
958 }
959
960 areaphes[i] = areaphes[i] / numareavalid[i];
961 areavars[i] = (areavars[i] - areaphes[i]*areaphes[i]/numareavalid[i]) / (numareavalid[i]-1.);
962
963 if (areavars[i] > 0.)
964 areavars[i] = TMath::Sqrt(areavars[i]);
965 else
966 {
967 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of photo-electrons found "
968 << "in area index: " << i << endl;
969 continue;
970 }
971
972 lowlim [i] = areaphes[i] - fPheErrLimit*TMath::Sqrt(areavars[i]);
973 upplim [i] = areaphes[i] + fPheErrLimit*TMath::Sqrt(areavars[i]);
974 }
975
976 memset(numareavalid,0,nareas*sizeof(Int_t));
977 memset(areaphes ,0,nareas*sizeof(Int_t));
978 memset(areavars ,0,nareas*sizeof(Int_t));
979
980 //
981 // Second loop: Get mean number of photo-electrons and its RMS excluding
982 // pixels deviating by more than fPheErrLimit sigma.
983 // Set the conversion factor FADC counts to photo-electrons
984 //
985 for (UInt_t i=0; i<npixels; i++)
986 {
987
988 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
989
990 if (!pix.IsFFactorMethodValid())
991 continue;
992
993 const Float_t nvar = pix.GetPheFFactorMethodVar();
994
995 if (nvar <= 0.)
996 {
997 pix.SetFFactorMethodValid(kFALSE);
998 continue;
999 }
1000
1001 MBadPixelsPix &bad = (*fBadPixels)[i];
1002
1003 const Int_t aidx = (*fGeom)[i].GetAidx();
1004 const Int_t sector = (*fGeom)[i].GetSector();
1005 const Float_t nphe = pix.GetPheFFactorMethod();
1006
1007 if ( nphe < lowlim[aidx] || nphe > upplim[aidx] )
1008 {
1009 *fLog << warn << GetDescriptor() << ": Deviating number of photo-electrons: "
1010 << Form("%4.2f",nphe) << " out of accepted limits: ["
1011 << Form("%4.2f%s%4.2f",lowlim[aidx],",",upplim[aidx]) << "] in pixel " << i << endl;
1012 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1013 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1014 pix.SetExcluded();
1015 continue;
1016 }
1017
1018 // const Float_t weight = 1./nvar;
1019 // areaweights [aidx] += weight;
1020 // areaphes [aidx] += weight*nphe;
1021 areaweights [aidx] += nphe*nphe;
1022 areaphes [aidx] += nphe;
1023 numareavalid [aidx] ++;
1024 // sectorweights [sector] += weight;
1025 // sectorphes [sector] += weight*nphe;
1026 sectorweights [sector] += nphe*nphe;
1027 sectorphes [sector] += nphe;
1028 numsectorvalid[sector] ++;
1029 }
1030
1031 for (UInt_t aidx=0; aidx<nareas; aidx++)
1032 {
1033
1034 MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);
1035
1036 areaweights[aidx] -= areaphes[aidx]*areaphes[aidx]/numareavalid[aidx];
1037 areaphes[aidx] /= numareavalid[aidx];
1038 areaweights[aidx] /= numareavalid[aidx]-1.;
1039
1040 if (areaweights[aidx] <= 0. || areaphes[aidx] <= 0.)
1041 {
1042 *fLog << warn << " Mean number of phe's from area index " << aidx << " cannot be calculated: "
1043 << " Mean number of phes: " << areaphes[aidx]
1044 << " Variance: " << areaweights[aidx] << endl;
1045 apix.SetFFactorMethodValid(kFALSE);
1046 continue;
1047 }
1048
1049 *fLog << inf << "Replacing number photo-electrons of average area idx " << aidx << ": "
1050 << Form("%5.3f%s%5.3f",apix.GetPheFFactorMethod()," +- ",apix.GetPheFFactorMethodErr()) << endl;
1051 *fLog << inf << " by average number of photo-electrons from area idx " << aidx << ": "
1052 << Form("%5.3f%s%5.3f",areaphes[aidx]," +- ",TMath::Sqrt(areaweights[aidx])) << endl;
1053
1054 // apix.SetPheFFactorMethod ( areaphes[aidx]/ areaweights[aidx] );
1055 // apix.SetPheFFactorMethodVar( 1. / areaweights[aidx] );
1056 apix.SetPheFFactorMethod ( areaphes[aidx] );
1057 apix.SetPheFFactorMethodVar( areaweights[aidx] / numareavalid[aidx] );
1058 apix.SetFFactorMethodValid ( kTRUE );
1059
1060 }
1061
1062 for (UInt_t sector=0; sector<nsectors; sector++)
1063 {
1064
1065 sectorweights[sector] -= sectorphes[sector]*sectorphes[sector]/numsectorvalid[sector];
1066 sectorphes[sector] /= numsectorvalid[sector];
1067 sectorweights[sector] /= numsectorvalid[sector]-1.;
1068
1069 MCalibrationChargePix &spix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
1070
1071 if (sectorweights[sector] <= 0. || sectorphes[sector] <= 0.)
1072 {
1073 *fLog << warn << " Mean number of phe's from sector " << sector << " cannot be calculated: "
1074 << " Mean number of phes: " << sectorphes[sector]
1075 << " Variance: " << sectorweights[sector] << endl;
1076 spix.SetFFactorMethodValid(kFALSE);
1077 continue;
1078 }
1079
1080 *fLog << inf << "Replacing number photo-electrons of average sector " << sector << ": "
1081 << Form("%5.3f%s%5.3f",spix.GetPheFFactorMethod()," +- ",spix.GetPheFFactorMethodErr()) << endl;
1082 *fLog << inf << " by average number photo-electrons from sector " << sector << ": "
1083 << Form("%5.3f%s%5.3f",sectorphes[sector]," +- ",TMath::Sqrt(sectorweights[sector])) << endl;
1084
1085 // spix.SetPheFFactorMethod ( sectorphes[sector]/ sectorweights[sector] );
1086 // spix.SetPheFFactorMethodVar( 1. / sectorweights[sector] );
1087 spix.SetPheFFactorMethod ( sectorphes[sector] );
1088 spix.SetPheFFactorMethodVar( sectorweights[sector] / numsectorvalid[sector]);
1089 spix.SetFFactorMethodValid ( kTRUE );
1090
1091 }
1092
1093 return kTRUE;
1094}
1095
1096
1097// ------------------------------------------------------------------------
1098//
1099// Returns kFALSE if pointer to MCalibrationChargeBlindPix is NULL
1100//
1101// The check returns kFALSE if:
1102//
1103// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
1104// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
1105//
1106// Calls:
1107// - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
1108//
1109Bool_t MCalibrationChargeCalc::FinalizeBlindPixel()
1110{
1111
1112 if (!fBlindPixel)
1113 return kFALSE;
1114
1115 const Float_t lambda = fBlindPixel->GetLambda();
1116 const Float_t lambdaerr = fBlindPixel->GetLambdaErr();
1117 const Float_t lambdacheck = fBlindPixel->GetLambdaCheck();
1118
1119 if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
1120 {
1121 *fLog << warn << GetDescriptor()
1122 << Form("%s%4.2f%s%4.2f%s%4.2f%s",": Lambda: ",lambda," and Lambda-Check: ",
1123 lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel ")
1124 << endl;
1125 return kFALSE;
1126 }
1127
1128 if (lambdaerr > fLambdaErrLimit)
1129 {
1130 *fLog << warn << GetDescriptor()
1131 << Form("%s%4.2f%s%4.2f%s",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
1132 fLambdaErrLimit," in Blind Pixel ") << endl;
1133 return kFALSE;
1134 }
1135
1136 if (!fBlindPixel->CalcFluxInsidePlexiglass())
1137 {
1138 *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
1139 << "will skip Blind Pixel Calibration " << endl;
1140 return kFALSE;
1141 }
1142
1143 return kTRUE;
1144}
1145
1146// ------------------------------------------------------------------------
1147//
1148// Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL
1149//
1150// The check returns kFALSE if:
1151//
1152// 1) PINDiode has a fitted charge smaller than fChargeLimit*PedRMS
1153// 2) PINDiode has a fit error smaller than fChargeErrLimit
1154// 3) PINDiode has a fitted charge smaller its fChargeRelErrLimit times its charge error
1155// 4) PINDiode has a charge sigma smaller than its Pedestal RMS
1156//
1157// Calls:
1158// - MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass()
1159//
1160Bool_t MCalibrationChargeCalc::FinalizePINDiode()
1161{
1162
1163 if (!fPINDiode)
1164 return kFALSE;
1165
1166 if (fPINDiode->GetMean() < fChargeLimit*fPINDiode->GetPedRms())
1167 {
1168 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1169 << fChargeLimit << " Pedestal RMS in PINDiode " << endl;
1170 return kFALSE;
1171 }
1172
1173 if (fPINDiode->GetMeanErr() < fChargeErrLimit)
1174 {
1175 *fLog << warn << GetDescriptor() << ": Error of Fitted Charge is smaller than "
1176 << fChargeErrLimit << " in PINDiode " << endl;
1177 return kFALSE;
1178 }
1179
1180 if (fPINDiode->GetMean() < fChargeRelErrLimit*fPINDiode->GetMeanErr())
1181 {
1182 *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
1183 << fChargeRelErrLimit << "* its error in PINDiode " << endl;
1184 return kFALSE;
1185 }
1186
1187 if (fPINDiode->GetSigma() < fPINDiode->GetPedRms())
1188 {
1189 *fLog << warn << GetDescriptor()
1190 << ": Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl;
1191 return kFALSE;
1192 }
1193
1194
1195 if (!fPINDiode->CalcFluxOutsidePlexiglass())
1196 {
1197 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
1198 << "will skip PIN Diode Calibration " << endl;
1199 return kFALSE;
1200 }
1201
1202 return kTRUE;
1203}
1204
1205// ------------------------------------------------------------------------
1206//
1207// Calculate the average number of photons outside the plexiglass with the
1208// formula:
1209//
1210// av.Num.photons(area index) = av.Num.Phes(area index)
1211// / MCalibrationQEPix::GetDefaultQE(fPulserColor)
1212// / MCalibrationQECam::GetPlexiglassQE()
1213//
1214// Calculate the variance on the average number of photons.
1215//
1216// Loop over pixels:
1217//
1218// - Continue, if not MCalibrationChargePix::IsFFactorMethodValid() and set:
1219// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor)
1220//
1221// - Call MCalibrationChargePix::CalcMeanFFactor(av.Num.photons) and set:
1222// MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) if not succesful
1223//
1224// - Calculate the quantum efficiency with the formula:
1225//
1226// QE = ( Num.Phes / av.Num.photons ) * MGeomCam::GetPixRatio()
1227//
1228// - Set QE in MCalibrationQEPix::SetQEFFactor ( QE, fPulserColor );
1229// - Set Variance of QE in MCalibrationQEPix::SetQEFFactorVar ( Variance, fPulserColor );
1230// - Set bit MCalibrationQEPix::SetFFactorMethodValid(kTRUE,fPulserColor)
1231//
1232// - Call MCalibrationQEPix::UpdateFFactorMethod()
1233//
1234void MCalibrationChargeCalc::FinalizeFFactorQECam()
1235{
1236
1237 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCam->GetAverageArea(0);
1238 MCalibrationQEPix &qepix = (MCalibrationQEPix&) fQECam->GetAverageArea(0);
1239
1240 const Float_t avphotons = avpix.GetPheFFactorMethod()
1241 / qepix.GetDefaultQE(fPulserColor)
1242 / fQECam->GetPlexiglassQE();
1243
1244 const Float_t avphotrelvar = avpix.GetPheFFactorMethodRelVar()
1245 + qepix.GetDefaultQERelVar(fPulserColor)
1246 + fQECam->GetPlexiglassQERelVar();
1247
1248 const UInt_t nareas = fGeom->GetNumAreas();
1249
1250 //
1251 // Set the results in the MCalibrationChargeCam
1252 //
1253 fCam->SetNumPhotonsFFactorMethod (avphotons);
1254 if (avphotrelvar > 0.)
1255 fCam->SetNumPhotonsFFactorMethodErr(TMath::Sqrt( avphotrelvar * avphotons * avphotons));
1256
1257 Float_t lowlim [nareas];
1258 Float_t upplim [nareas];
1259 Float_t avffactorphotons [nareas];
1260 Float_t avffactorphotvar [nareas];
1261 Int_t numffactor [nareas];
1262
1263 memset(lowlim ,0, nareas * sizeof(Float_t));
1264 memset(upplim ,0, nareas * sizeof(Float_t));
1265 memset(avffactorphotons,0, nareas * sizeof(Float_t));
1266 memset(avffactorphotvar,0, nareas * sizeof(Float_t));
1267 memset(numffactor ,0, nareas * sizeof(Int_t));
1268
1269 const UInt_t npixels = fGeom->GetNumPixels();
1270
1271 for (UInt_t i=0; i<npixels; i++)
1272 {
1273
1274 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1275 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1276
1277 if (!pix.IsFFactorMethodValid())
1278 {
1279 qepix.SetExcluded();
1280 qepix.SetFFactorMethodValid(kFALSE,fPulserColor);
1281 continue;
1282 }
1283
1284 const Float_t photons = avphotons / fGeom->GetPixRatio(i);
1285 const Float_t qe = pix.GetPheFFactorMethod() / photons ;
1286
1287 if (!pix.CalcMeanFFactor( photons , avphotrelvar ))
1288 {
1289 pix.SetFFactorMethodValid(kFALSE);
1290 qepix.SetFFactorMethodValid(kFALSE, fPulserColor);
1291 (*fBadPixels)[i].SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
1292 }
1293
1294 const Float_t qerelvar = avphotrelvar + pix.GetPheFFactorMethodRelVar();
1295
1296 qepix.SetQEFFactor ( qe , fPulserColor );
1297 qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor );
1298 qepix.SetFFactorMethodValid( kTRUE , fPulserColor );
1299
1300 if (!qepix.UpdateFFactorMethod())
1301 *fLog << warn << GetDescriptor()
1302 << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl;
1303
1304 const Int_t aidx = (*fGeom)[i].GetAidx();
1305
1306 avffactorphotons[aidx] += pix.GetMeanFFactorFADC2Phot();
1307 avffactorphotvar[aidx] += pix.GetMeanFFactorFADC2PhotVar();
1308 numffactor[aidx]++;
1309 }
1310
1311 for (UInt_t i=0; i<nareas; i++)
1312 {
1313 avffactorphotons[i] /= numffactor[i];
1314 avffactorphotvar[i] /= numffactor[i];
1315 lowlim [i] = avffactorphotons[i] - fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);
1316 upplim [i] = avffactorphotons[i] + fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);
1317 }
1318
1319 for (UInt_t i=0; i<npixels; i++)
1320 {
1321
1322 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1323
1324 if (!pix.IsFFactorMethodValid())
1325 continue;
1326
1327 const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
1328 MBadPixelsPix &bad = (*fBadPixels)[i];
1329
1330 const Int_t aidx = (*fGeom)[i].GetAidx();
1331
1332 if ( ffactor < lowlim[aidx] || ffactor > upplim[aidx] )
1333 {
1334 *fLog << warn << GetDescriptor() << ": Deviating F-Factor: "
1335 << Form("%4.2f",ffactor) << " out of accepted limits: ["
1336 << Form("%4.2f%s%4.2f",lowlim[aidx],",",upplim[aidx]) << "] in pixel " << i << endl;
1337 bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
1338 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
1339 pix.SetExcluded();
1340 }
1341 }
1342}
1343
1344
1345// ------------------------------------------------------------------------
1346//
1347// Loop over pixels:
1348//
1349// - Continue, if not MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() and set:
1350// MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor)
1351//
1352// - Calculate the quantum efficiency with the formula:
1353//
1354// QE = Num.Phes / MCalibrationChargeBlindPix::GetFluxInsidePlexiglass()
1355// / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE()
1356//
1357// - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor );
1358// - Set Variance of QE in MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor );
1359// - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor)
1360//
1361// - Call MCalibrationQEPix::UpdateBlindPixelMethod()
1362//
1363void MCalibrationChargeCalc::FinalizeBlindPixelQECam()
1364{
1365
1366 const UInt_t npixels = fGeom->GetNumPixels();
1367
1368 //
1369 // Set the results in the MCalibrationChargeCam
1370 //
1371 if (fBlindPixel)
1372 {
1373 if (fBlindPixel->IsFluxInsidePlexiglassAvailable())
1374 {
1375
1376 const Float_t photons = fBlindPixel->GetFluxInsidePlexiglass() * (*fGeom)[0].GetA()
1377 / fQECam->GetPlexiglassQE();
1378 fCam->SetNumPhotonsBlindPixelMethod(photons);
1379
1380 const Float_t photrelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1381 + fQECam->GetPlexiglassQERelVar();
1382 if (photrelvar > 0.)
1383 fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));
1384 }
1385 }
1386 //
1387 // With the knowledge of the overall photon flux, calculate the
1388 // quantum efficiencies after the Blind Pixel and PIN Diode method
1389 //
1390 for (UInt_t i=0; i<npixels; i++)
1391 {
1392
1393 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1394
1395 if (!fBlindPixel)
1396 {
1397 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1398 continue;
1399 }
1400
1401 if (!fBlindPixel->IsFluxInsidePlexiglassAvailable())
1402 {
1403 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1404 continue;
1405 }
1406
1407 MBadPixelsPix &bad = (*fBadPixels)[i];
1408
1409 if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1410 {
1411 qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
1412 continue;
1413 }
1414
1415 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1416 MGeomPix &geo = (*fGeom)[i];
1417
1418 const Float_t qe = pix.GetPheFFactorMethod()
1419 / fBlindPixel->GetFluxInsidePlexiglass()
1420 / geo.GetA()
1421 * fQECam->GetPlexiglassQE();
1422
1423 const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar()
1424 + fQECam->GetPlexiglassQERelVar()
1425 + pix.GetPheFFactorMethodRelVar();
1426
1427 qepix.SetQEBlindPixel ( qe , fPulserColor );
1428 qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor );
1429 qepix.SetBlindPixelMethodValid( kTRUE , fPulserColor );
1430
1431 if (!qepix.UpdateBlindPixelMethod())
1432 *fLog << warn << GetDescriptor()
1433 << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl;
1434 }
1435}
1436
1437// ------------------------------------------------------------------------
1438//
1439// Loop over pixels:
1440//
1441// - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set:
1442// MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor)
1443//
1444// - Calculate the quantum efficiency with the formula:
1445//
1446// QE = Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA()
1447//
1448// - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor );
1449// - Set Variance of QE in MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor );
1450// - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor)
1451//
1452// - Call MCalibrationQEPix::UpdatePINDiodeMethod()
1453//
1454void MCalibrationChargeCalc::FinalizePINDiodeQECam()
1455{
1456
1457 const UInt_t npixels = fGeom->GetNumPixels();
1458
1459 //
1460 // With the knowledge of the overall photon flux, calculate the
1461 // quantum efficiencies after the PIN Diode method
1462 //
1463 for (UInt_t i=0; i<npixels; i++)
1464 {
1465
1466 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam)[i];
1467
1468 if (!fPINDiode)
1469 {
1470 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1471 continue;
1472 }
1473
1474 if (!fPINDiode->IsFluxOutsidePlexiglassAvailable())
1475 {
1476 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1477 continue;
1478 }
1479
1480 MBadPixelsPix &bad = (*fBadPixels)[i];
1481
1482 if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
1483 {
1484 qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
1485 continue;
1486 }
1487
1488 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
1489 MGeomPix &geo = (*fGeom)[i];
1490
1491 const Float_t qe = pix.GetPheFFactorMethod()
1492 / fPINDiode->GetFluxOutsidePlexiglass()
1493 / geo.GetA();
1494
1495 const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar();
1496
1497 qepix.SetQEPINDiode ( qe , fPulserColor );
1498 qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor );
1499 qepix.SetPINDiodeMethodValid( kTRUE , fPulserColor );
1500
1501 if (!qepix.UpdatePINDiodeMethod())
1502 *fLog << warn << GetDescriptor()
1503 << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl;
1504 }
1505}
1506
1507// -----------------------------------------------------------------------------------------------
1508//
1509// - Print out statistics about BadPixels of type UnsuitableType_t
1510// - store numbers of bad pixels of each type in fCam
1511//
1512void MCalibrationChargeCalc::FinalizeUnsuitablePixels()
1513{
1514
1515 *fLog << inf << endl;
1516 *fLog << GetDescriptor() << ": Charge Calibration status:" << endl;
1517 *fLog << dec << setfill(' ');
1518
1519 const Int_t nareas = fGeom->GetNumAreas();
1520
1521 Int_t counts[nareas];
1522 memset(counts,0,nareas*sizeof(Int_t));
1523
1524 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1525 {
1526 MBadPixelsPix &bad = (*fBadPixels)[i];
1527 if (!bad.IsBad())
1528 {
1529 const Int_t aidx = (*fGeom)[i].GetAidx();
1530 counts[aidx]++;
1531 }
1532 }
1533
1534 if (fGeom->InheritsFrom("MGeomCamMagic"))
1535 *fLog << " " << setw(7) << "Successfully calibrated Pixels: "
1536 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1537
1538 memset(counts,0,nareas*sizeof(Int_t));
1539
1540 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1541 {
1542 MBadPixelsPix &bad = (*fBadPixels)[i];
1543 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1544 {
1545 const Int_t aidx = (*fGeom)[i].GetAidx();
1546 counts[aidx]++;
1547 }
1548 }
1549
1550 for (Int_t aidx=0; aidx<nareas; aidx++)
1551 fCam->SetNumUnsuitable(counts[aidx], aidx);
1552
1553 if (fGeom->InheritsFrom("MGeomCamMagic"))
1554 *fLog << " " << setw(7) << "Uncalibrated Pixels: "
1555 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1556
1557 memset(counts,0,nareas*sizeof(Int_t));
1558
1559 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1560 {
1561 MBadPixelsPix &bad = (*fBadPixels)[i];
1562 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
1563 {
1564 const Int_t aidx = (*fGeom)[i].GetAidx();
1565 counts[aidx]++;
1566 }
1567 }
1568
1569 for (Int_t aidx=0; aidx<nareas; aidx++)
1570 fCam->SetNumUnreliable(counts[aidx], aidx);
1571
1572 *fLog << " " << setw(7) << "Unreliable Pixels: "
1573 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
1574
1575}
1576
1577// -----------------------------------------------------------------------------------------------
1578//
1579// Print out statistics about BadPixels of type UncalibratedType_t
1580//
1581void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
1582{
1583
1584 UInt_t countinner = 0;
1585 UInt_t countouter = 0;
1586
1587 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
1588 {
1589 MBadPixelsPix &bad = (*fBadPixels)[i];
1590 if (bad.IsUncalibrated(typ))
1591 {
1592 if (fGeom->GetPixRatio(i) == 1.)
1593 countinner++;
1594 else
1595 countouter++;
1596 }
1597 }
1598
1599 *fLog << " " << setw(7) << text
1600 << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
1601}
1602
1603// --------------------------------------------------------------------------
1604//
1605// Set the path for output file
1606//
1607void MCalibrationChargeCalc::SetOutputPath(TString path)
1608{
1609 fOutputPath = path;
1610 if (fOutputPath.EndsWith("/"))
1611 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
1612}
1613
1614void MCalibrationChargeCalc::SetOutputFile(TString file)
1615{
1616 fOutputFile = file;
1617}
1618
1619// --------------------------------------------------------------------------
1620//
1621// Get the output file
1622//
1623const char* MCalibrationChargeCalc::GetOutputFile()
1624{
1625 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
1626}
Note: See TracBrowser for help on using the repository browser.