/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Markus Gaug 11/2003 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MCalibrationChargeCam // // Hold the whole Calibration results of the camera: // // 1) MCalibrationChargeCam initializes a TClonesArray whose elements are // pointers to MCalibrationChargePix Containers // 2) It initializes a pointer to an MCalibrationBlindPix container // 3) It initializes a pointer to an MCalibrationPINDiode container // // // The calculated values (types of GetPixelContent) are: // // -------------------------------------------------------------------------- // // The types are as follows: // // Fitted values: // ============== // // 0: Fitted Charge // 1: Error of fitted Charge // 2: Sigma of fitted Charge // 3: Error of Sigma of fitted Charge // // Useful variables derived from the fit results: // ============================================= // // 4: Returned probability of Gauss fit to Charge distribution // 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2) // 6: Error Reduced Sigma of fitted Charge // 7: Reduced Sigma per Charge // 8: Error of Reduced Sigma per Charge // // Results of the different calibration methods: // ============================================= // // 9: Number of Photo-electrons obtained with the F-Factor method // 10: Error on Number of Photo-electrons obtained with the F-Factor method // 11: Mean conversion factor obtained with the F-Factor method // 12: Error on the mean conversion factor obtained with the F-Factor method // 13: Overall F-Factor of the readout obtained with the F-Factor method // 14: Error on Overall F-Factor of the readout obtained with the F-Factor method // 15: Pixels with valid calibration by the F-Factor-Method // 16: Mean conversion factor obtained with the Blind Pixel method // 17: Error on the mean conversion factor obtained with the Blind Pixel method // 18: Overall F-Factor of the readout obtained with the Blind Pixel method // 19: Error on Overall F-Factor of the readout obtained with the Blind Pixel method // 20: Pixels with valid calibration by the Blind Pixel-Method // 21: Mean conversion factor obtained with the PIN Diode method // 22: Error on the mean conversion factor obtained with the PIN Diode method // 23: Overall F-Factor of the readout obtained with the PIN Diode method // 24: Error on Overall F-Factor of the readout obtained with the PIN Diode method // 25: Pixels with valid calibration by the PIN Diode-Method // // Localized defects: // ================== // // 26: Excluded Pixels // 27: Number of probable pickup events in the Hi Gain // 28: Number of probable pickup events in the Lo Gain // // Other classifications of pixels: // ================================ // // 29: Pixels with saturated Hi-Gain // // Used Pedestals: // =============== // // 30: Mean Pedestal over the entire range of signal extraction // 31: Error on the Mean Pedestal over the entire range of signal extraction // 32: Pedestal RMS over the entire range of signal extraction // 33: Error on the Pedestal RMS over the entire range of signal extraction // // Calculated absolute arrival times (very low precision!): // ======================================================== // // 34: Absolute Arrival time of the signal // 35: RMS of the Absolute Arrival time of the signal // ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationChargeCam.h" #include "MCalibrationCam.h" #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MGeomCam.h" #include "MGeomPix.h" #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" #include "MCalibrationChargePix.h" #include "MCalibrationChargeBlindPix.h" #include "MCalibrationChargePINDiode.h" ClassImp(MCalibrationChargeCam); using namespace std; const Float_t MCalibrationChargeCam::fgAverageQE = 0.18; const Float_t MCalibrationChargeCam::fgAverageQEErr = 0.02; const Float_t MCalibrationChargeCam::fgConvFFactorRelErrLimit = 0.35; const Float_t MCalibrationChargeCam::fgPheFFactorRelErrLimit = 5.; // -------------------------------------------------------------------------- // // Default constructor. // // Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry // Later, a call to MCalibrationChargeCam::InitSize(Int_t size) has to be performed // // Creates an MCalibrationBlindPix container // MCalibrationChargeCam::MCalibrationChargeCam(const char *name, const char *title) : fOffsets(NULL), fSlopes(NULL), fOffvsSlope(NULL) { fName = name ? name : "MCalibrationChargeCam"; fTitle = title ? title : "Storage container for the Calibration Information in the camera"; fPixels = new TClonesArray("MCalibrationChargePix",1); fAverageAreas = new TClonesArray("MCalibrationChargePix",1); fAverageSectors = new TClonesArray("MCalibrationChargePix",1); fMeanFluxPhesInnerPixel = 0.; fMeanFluxPhesInnerPixelVar = 0.; fMeanFluxPhesOuterPixel = 0.; fMeanFluxPhesOuterPixelVar = 0.; CLRBIT(fFlags,kBlindPixelMethodValid); CLRBIT(fFlags,kFFactorMethodValid); CLRBIT(fFlags,kPINDiodeMethodValid); SetAverageQE(); SetConvFFactorRelErrLimit(); SetPheFFactorRelErrLimit(); } // -------------------------------------------------------------------------- // // Delete the TClonesArray of MCalibrationPix containers // Delete the MCalibrationPINDiode and the MCalibrationBlindPix // // Delete the histograms if they exist // MCalibrationChargeCam::~MCalibrationChargeCam() { if (fOffsets) delete fOffsets; if (fSlopes) delete fSlopes; if (fOffvsSlope) delete fOffvsSlope; } // -------------------------------------- // void MCalibrationChargeCam::Clear(Option_t *o) { fMeanFluxPhesInnerPixel = 0.; fMeanFluxPhesInnerPixelVar = 0.; fMeanFluxPhesOuterPixel = 0.; fMeanFluxPhesOuterPixelVar = 0.; CLRBIT(fFlags,kBlindPixelMethodValid); CLRBIT(fFlags,kFFactorMethodValid); CLRBIT(fFlags,kPINDiodeMethodValid); MCalibrationCam::Clear(); return; } void MCalibrationChargeCam::SetFFactorMethodValid(const Bool_t b) { b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid); } void MCalibrationChargeCam::SetBlindPixelMethodValid(const Bool_t b) { b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); } void MCalibrationChargeCam::SetPINDiodeMethodValid(const Bool_t b) { b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); } Float_t MCalibrationChargeCam::GetMeanFluxPhesInnerPixelErr() const { if (fMeanFluxPhesInnerPixelVar <= 0.) return -1.; return TMath::Sqrt(fMeanFluxPhesInnerPixelVar); } Float_t MCalibrationChargeCam::GetMeanFluxPhesOuterPixelErr() const { if (fMeanFluxPhesOuterPixelVar <= 0.) return -1.; return TMath::Sqrt(fMeanFluxPhesOuterPixelVar); } Float_t MCalibrationChargeCam::GetMeanFluxPhotonsInnerPixelErr() const { if (fMeanFluxPhotonsInnerPixelVar <= 0.) return -1.; return TMath::Sqrt(fMeanFluxPhotonsInnerPixelVar); } Float_t MCalibrationChargeCam::GetMeanFluxPhotonsOuterPixelErr() const { if (fMeanFluxPhotonsOuterPixelVar <= 0.) return -1.; return TMath::Sqrt(fMeanFluxPhotonsOuterPixelVar); } Bool_t MCalibrationChargeCam::IsBlindPixelMethodValid() const { return TESTBIT(fFlags,kBlindPixelMethodValid); } Bool_t MCalibrationChargeCam::IsPINDiodeMethodValid() const { return TESTBIT(fFlags,kPINDiodeMethodValid); } // -------------------------------------------------------------------------- // // Print first the well fitted pixels // and then the ones which are not FitValid // void MCalibrationChargeCam::Print(Option_t *o) const { *fLog << all << GetDescriptor() << ":" << endl; int id = 0; *fLog << all << "Calibrated pixels:" << endl; *fLog << all << endl; TIter Next(fPixels); MCalibrationChargePix *pix; while ((pix=(MCalibrationChargePix*)Next())) { if (!pix->IsExcluded()) { *fLog << all << "Pix " << pix->GetPixId() << ": Ped. Rms: " << pix->GetPedRms() << " +- " << pix->GetPedRmsErr() << " Mean signal: " << pix->GetMean() << " +- " << pix->GetSigma() << " Reduced Sigma: " << pix->GetRSigma() << " Nr Phe's: " << pix->GetPheFFactorMethod() << " Saturated? :" << pix->IsHiGainSaturation() << endl; id++; } } *fLog << all << id << " pixels" << endl; id = 0; *fLog << all << endl; *fLog << all << "Excluded pixels:" << endl; *fLog << all << endl; id = 0; TIter Next4(fPixels); while ((pix=(MCalibrationChargePix*)Next4())) { if (pix->IsExcluded()) { *fLog << all << pix->GetPixId() << endl; id++; } } *fLog << all << id << " Excluded pixels " << endl; *fLog << endl; TIter Next5(fAverageAreas); while ((pix=(MCalibrationChargePix*)Next5())) { *fLog << all << "Average Area:" << " Ped. Rms: " << pix->GetPedRms() << " +- " << pix->GetPedRmsErr() << " Mean signal: " << pix->GetMean() << " +- " << pix->GetMeanErr() << " Sigma signal: " << pix->GetSigma() << " +- "<< pix->GetSigmaErr() << " Reduced Sigma: " << pix->GetRSigma() << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl; } TIter Next6(fAverageSectors); while ((pix=(MCalibrationChargePix*)Next5())) { *fLog << all << "Average Sector:" << " Ped. Rms: " << pix->GetPedRms() << " +- " << pix->GetPedRmsErr() << " Mean signal: " << pix->GetMean() << " +- " << pix->GetMeanErr() << " Sigma signal: " << pix->GetSigma() << " +- "<< pix->GetSigmaErr() << " Reduced Sigma: " << pix->GetRSigma() << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl; } } // -------------------------------------------------------------------------- // // The types are as follows: // // Fitted values: // ============== // // 0: Fitted Charge // 1: Error of fitted Charge // 2: Sigma of fitted Charge // 3: Error of Sigma of fitted Charge // // Useful variables derived from the fit results: // ============================================= // // 4: Returned probability of Gauss fit to Charge distribution // 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2) // 6: Error Reduced Sigma of fitted Charge // 7: Reduced Sigma per Charge // 8: Error of Reduced Sigma per Charge // // Useful variables derived from the fit results: // ============================================= // // 4: Returned probability of Gauss fit to Charge distribution // 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2) // 6: Error Reduced Sigma of fitted Charge // 7: Reduced Sigma per Charge // 8: Error of Reduced Sigma per Charge // // Results of the different calibration methods: // ============================================= // // 9: Number of Photo-electrons obtained with the F-Factor method // 10: Error on Number of Photo-electrons obtained with the F-Factor method // 11: Mean conversion factor obtained with the F-Factor method // 12: Error on the mean conversion factor obtained with the F-Factor method // 13: Overall F-Factor of the readout obtained with the F-Factor method // 14: Error on Overall F-Factor of the readout obtained with the F-Factor method // 15: Pixels with valid calibration by the F-Factor-Method // 16: Mean conversion factor obtained with the Blind Pixel method // 17: Error on the mean conversion factor obtained with the Blind Pixel method // 18: Overall F-Factor of the readout obtained with the Blind Pixel method // 19: Error on Overall F-Factor of the readout obtained with the Blind Pixel method // 20: Pixels with valid calibration by the Blind Pixel-Method // 21: Mean conversion factor obtained with the PIN Diode method // 22: Error on the mean conversion factor obtained with the PIN Diode method // 23: Overall F-Factor of the readout obtained with the PIN Diode method // 24: Error on Overall F-Factor of the readout obtained with the PIN Diode method // 25: Pixels with valid calibration by the PIN Diode-Method // // Localized defects: // ================== // // 26: Excluded Pixels // 27: Number of probable pickup events in the Hi Gain // 28: Number of probable pickup events in the Lo Gain // // Other classifications of pixels: // ================================ // // 29: Pixels with saturated Hi-Gain // // Used Pedestals: // =============== // // 30: Mean Pedestal over the entire range of signal extraction // 31: Error on the Mean Pedestal over the entire range of signal extraction // 32: Pedestal RMS over the entire range of signal extraction // 33: Error on the Pedestal RMS over the entire range of signal extraction // // Calculated absolute arrival times (very low precision!): // ======================================================== // // 34: Absolute Arrival time of the signal // 35: RMS of the Absolute Arrival time of the signal // Bool_t MCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const { if (idx > GetSize()) return kFALSE; Float_t area = cam[idx].GetA(); if (area == 0) return kFALSE; MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[idx]; switch (type) { case 0: if (pix.IsExcluded()) return kFALSE; val = pix.GetMean(); break; case 1: if (pix.IsExcluded()) return kFALSE; val = pix.GetMeanErr(); break; case 2: if (pix.IsExcluded()) return kFALSE; val = pix.GetSigma(); break; case 3: if (pix.IsExcluded()) return kFALSE; val = pix.GetSigmaErr(); break; case 4: if (pix.IsExcluded()) return kFALSE; val = pix.GetProb(); break; case 5: if (pix.IsExcluded()) return kFALSE; if (pix.GetRSigma() == -1.) return kFALSE; val = pix.GetRSigma(); break; case 6: if (pix.IsExcluded()) return kFALSE; if (pix.GetRSigma() == -1.) return kFALSE; val = pix.GetRSigmaErr(); break; case 7: if (pix.IsExcluded()) return kFALSE; if (pix.GetRSigma() == -1.) return kFALSE; if (pix.GetMean() == 0.) return kFALSE; val = pix.GetRSigma() / pix.GetMean(); break; case 8: if (pix.IsExcluded()) return kFALSE; if (pix.GetRSigma() <= 0.) return kFALSE; if (pix.GetMean() <= 0.) return kFALSE; if (pix.GetRSigmaErr() <= 0.) return kFALSE; if (pix.GetMeanErr() <= 0.) return kFALSE; // relative error Rsigma square val = pix.GetRSigmaErr()* pix.GetRSigmaErr() / (pix.GetRSigma() * pix.GetRSigma() ); // relative error square val += pix.GetMeanErr() * pix.GetMeanErr() / (pix.GetMean() * pix.GetMean() ); // calculate relative error out of squares val = TMath::Sqrt(val) ; // multiply with value to get absolute error val *= pix.GetRSigma() / pix.GetMean(); break; case 9: if (pix.IsExcluded() || !pix.IsFFactorMethodValid()) return kFALSE; val = pix.GetPheFFactorMethod(); break; case 10: if (pix.IsExcluded() || !pix.IsFFactorMethodValid()) return kFALSE; val = pix.GetPheFFactorMethodErr(); break; case 11: if (pix.IsExcluded() || !pix.IsFFactorMethodValid()) return kFALSE; val = pix.GetMeanConversionFFactorMethod(); break; case 12: if (pix.IsExcluded() || !pix.IsFFactorMethodValid()) return kFALSE; val = pix.GetConversionFFactorMethodErr(); break; case 13: if (pix.IsExcluded() || !pix.IsFFactorMethodValid()) return kFALSE; val = pix.GetTotalFFactorFFactorMethod(); break; case 14: if (pix.IsExcluded() || !pix.IsFFactorMethodValid()) return kFALSE; val = pix.GetTotalFFactorFFactorMethodErr(); break; case 15: if (pix.IsExcluded()) return kFALSE; if (pix.IsFFactorMethodValid()) val = 1; else return kFALSE; break; case 16: if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid()) return kFALSE; val = pix.GetMeanConversionBlindPixelMethod(); break; case 17: if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid()) return kFALSE; val = pix.GetConversionBlindPixelMethodErr(); break; case 18: if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid()) return kFALSE; val = pix.GetTotalFFactorBlindPixelMethod(); break; case 19: if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid()) return kFALSE; val = pix.GetTotalFFactorBlindPixelMethodErr(); break; case 20: if (pix.IsExcluded()) return kFALSE; if (pix.IsBlindPixelMethodValid()) val = 1; else return kFALSE; break; case 21: if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid()) return kFALSE; val = pix.GetMeanConversionPINDiodeMethod(); break; case 22: if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid()) return kFALSE; val = pix.GetConversionPINDiodeMethodErr(); break; case 23: if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid()) return kFALSE; val = pix.GetTotalFFactorPINDiodeMethod(); break; case 24: if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid()) return kFALSE; val = pix.GetTotalFFactorPINDiodeMethodErr(); break; case 25: if (pix.IsExcluded()) return kFALSE; if (pix.IsPINDiodeMethodValid()) val = 1; else return kFALSE; break; case 26: if (pix.IsExcluded()) val = 1.; else return kFALSE; break; case 27: if (pix.IsExcluded()) return kFALSE; val = pix.GetHiGainNumPickup(); break; case 28: if (pix.IsExcluded()) return kFALSE; val = pix.GetLoGainNumPickup(); break; case 29: if (pix.IsExcluded()) return kFALSE; val = pix.IsHiGainSaturation(); break; case 30: if (pix.IsExcluded()) return kFALSE; val = pix.GetPed(); break; case 31: if (pix.IsExcluded()) return kFALSE; val = pix.GetPedErr(); break; case 32: if (pix.IsExcluded()) return kFALSE; val = pix.GetPedRms(); break; case 33: if (pix.IsExcluded()) return kFALSE; val = pix.GetPedErr()/2.; break; case 34: if (pix.IsExcluded()) return kFALSE; val = pix.GetAbsTimeMean(); break; case 35: if (pix.IsExcluded()) return kFALSE; val = pix.GetAbsTimeRms(); break; default: return kFALSE; } return val!=-1.; } // -------------------------------------------------------------------------- // // What MHCamera needs in order to draw an individual pixel in the camera // void MCalibrationChargeCam::DrawPixelContent(Int_t idx) const { MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[idx]; pix.DrawClone(); } // // Calculate the weighted mean of the phe's of all inner and outer pixels, respectively. // Bad pixels are excluded from the calculation. Two loops are performed to exclude pixels // which are fPheFFactorRelLimit sigmas from the mean. // Bool_t MCalibrationChargeCam::CalcMeanFluxPhotonsFFactorMethod(const MGeomCam &geom, MBadPixelsCam &bad) { const Float_t avQERelVar = fAverageQEVar / (fAverageQE * fAverageQE ); Float_t sumweightsinner = 0.; Float_t sumphesinner = 0.; Float_t sumweightsouter = 0.; Float_t sumphesouter = 0.; Int_t validinner = 0; Int_t validouter = 0; TIter Next(fPixels); MCalibrationChargePix *pix; while ((pix=(MCalibrationChargePix*)Next())) { if (!pix->IsFFactorMethodValid()) continue; const Int_t idx = pix->GetPixId(); if(!bad[idx].IsCalibrationResultOK()) continue; const Float_t nphe = pix->GetPheFFactorMethod(); const Float_t npheerr = pix->GetPheFFactorMethodErr(); const Float_t ratio = geom.GetPixRatio(idx); if (npheerr > 0.) { const Float_t weight = nphe*nphe/npheerr/npheerr; // // first the inner pixels: // if (ratio == 1.) { sumweightsinner += weight; sumphesinner += weight*nphe; validinner++; } else { // // now the outers // sumweightsouter += weight; sumphesouter += weight*nphe; validouter++; } } /* if npheerr != 0 */ } /* while ((pix=(MCalibrationChargePix*)Next())) */ if (sumweightsinner <= 0. || sumphesinner <= 0.) { *fLog << warn << " Mean number of phe's from inner pixels cannot be calculated: " << " Sum of weights: " << sumweightsinner << " Sum of weighted phes: " << sumphesinner << endl; return kFALSE; } else { fMeanFluxPhesInnerPixel = sumphesinner/sumweightsinner; fMeanFluxPhesInnerPixelVar = (1./sumweightsinner)*fMeanFluxPhesInnerPixel*fMeanFluxPhesInnerPixel; } if (sumweightsouter <= 0. || sumphesouter <= 0.) { *fLog << warn << " Mean number of phe's from outer pixels cannot be calculated: " << " Sum of weights or sum of weighted phes is 0. " << endl; } else { fMeanFluxPhesOuterPixel = sumphesouter/sumweightsouter; fMeanFluxPhesOuterPixelVar = (1./sumweightsouter)*fMeanFluxPhesOuterPixel*fMeanFluxPhesOuterPixel; } Float_t meanFluxPhotonsRelVar = fMeanFluxPhesInnerPixelVar / (fMeanFluxPhesInnerPixel * fMeanFluxPhesInnerPixel); fMeanFluxPhotonsInnerPixel = fMeanFluxPhesInnerPixel/fAverageQE; fMeanFluxPhotonsInnerPixelVar = (meanFluxPhotonsRelVar + avQERelVar) * fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel; fMeanFluxPhotonsOuterPixel = 4. *fMeanFluxPhotonsInnerPixel; fMeanFluxPhotonsOuterPixelVar = 16.*fMeanFluxPhotonsInnerPixelVar; *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): " << fMeanFluxPhesInnerPixel << " +- " << GetMeanFluxPhesInnerPixelErr() << endl; *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): " << fMeanFluxPhotonsInnerPixel << " +- " << GetMeanFluxPhotonsInnerPixelErr() << endl; // // Here starts the second loop discarting pixels out of the range: // const Float_t innervar = (Float_t)validinner*fPheFFactorRelVarLimit*fMeanFluxPhesInnerPixelVar; const Float_t outervar = (Float_t)validouter*fPheFFactorRelVarLimit*fMeanFluxPhesOuterPixelVar; Float_t innererr; Float_t outererr; if (innervar > 0.) innererr = TMath::Sqrt(innervar); else { *fLog << warn << GetDescriptor() << ": Variance of mean number of photo-electrons for inner pixels " << " is smaller than 0. " << endl; SetFFactorMethodValid(kFALSE); return kFALSE; } if (outervar > 0.) outererr = TMath::Sqrt(outervar); else { *fLog << warn << GetDescriptor() << ": Variance of mean number of photo-electrons for outer pixels " << " is smaller than 0. " << endl; SetFFactorMethodValid(kFALSE); return kFALSE; } const Float_t lowerpheinnerlimit = fMeanFluxPhesInnerPixel - innererr; const Float_t upperpheinnerlimit = fMeanFluxPhesInnerPixel + innererr; const Float_t lowerpheouterlimit = fMeanFluxPhesOuterPixel - outererr; const Float_t upperpheouterlimit = fMeanFluxPhesOuterPixel + outererr; sumweightsinner = 0.; sumphesinner = 0.; sumweightsouter = 0.; sumphesouter = 0.; TIter Next2(fPixels); MCalibrationChargePix *pix2; while ((pix2=(MCalibrationChargePix*)Next2())) { if (!pix2->IsFFactorMethodValid()) continue; const Int_t idx = pix2->GetPixId(); if(!bad[idx].IsCalibrationResultOK()) continue; const Float_t nphe = pix2->GetPheFFactorMethod(); const Float_t npheerr = pix2->GetPheFFactorMethodErr(); const Float_t ratio = geom.GetPixRatio(idx); if (npheerr > 0.) { const Float_t weight = nphe*nphe/npheerr/npheerr; // // first the inner pixels: // if (ratio == 1.) { if (nphe < lowerpheinnerlimit || nphe > upperpheinnerlimit) { bad[idx].SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes); bad[idx].SetUnsuitable( MBadPixelsPix::kUnreliableRun); continue; } sumweightsinner += weight; sumphesinner += weight*nphe; } else { // // now the outers // if (nphe < lowerpheouterlimit || nphe > upperpheouterlimit) { bad[idx].SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes); bad[idx].SetUnsuitable( MBadPixelsPix::kUnreliableRun); continue; } sumweightsouter += weight; sumphesouter += weight*nphe; } } /* if npheerr != 0 */ } /* while ((pix2=(MCalibrationChargePix*)Next2())) */ if (sumweightsinner <= 0. || sumphesinner <= 0.) { *fLog << warn << " Mean number of phe's from inner pixels cannot be calculated: " << " Sum of weights: " << sumweightsinner << " Sum of weighted phes: " << sumphesinner << endl; return kFALSE; } else { fMeanFluxPhesInnerPixel = sumphesinner/sumweightsinner; fMeanFluxPhesInnerPixelVar = (1./sumweightsinner)*fMeanFluxPhesInnerPixel*fMeanFluxPhesInnerPixel; } if (sumweightsouter <= 0. || sumphesouter <= 0.) { *fLog << warn << " Mean number of phe's from outer pixels cannot be calculated: " << " Sum of weights or sum of weighted phes is 0. " << endl; } else { fMeanFluxPhesOuterPixel = sumphesouter/sumweightsouter; fMeanFluxPhesOuterPixelVar = (1./sumweightsouter)*fMeanFluxPhesOuterPixel*fMeanFluxPhesOuterPixel; } meanFluxPhotonsRelVar = fMeanFluxPhesInnerPixelVar / (fMeanFluxPhesInnerPixel * fMeanFluxPhesInnerPixel); fMeanFluxPhotonsInnerPixel = fMeanFluxPhesInnerPixel/fAverageQE; fMeanFluxPhotonsInnerPixelVar = (meanFluxPhotonsRelVar + avQERelVar) * fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel; fMeanFluxPhotonsOuterPixel = 4. *fMeanFluxPhotonsInnerPixel; fMeanFluxPhotonsOuterPixelVar = 16.*fMeanFluxPhotonsInnerPixelVar; *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): " << fMeanFluxPhesInnerPixel << " +- " << GetMeanFluxPhesInnerPixelErr() << endl; *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): " << fMeanFluxPhotonsInnerPixel << " +- " << GetMeanFluxPhotonsInnerPixelErr() << endl; return kTRUE; } void MCalibrationChargeCam::ApplyFFactorCalibration(const MGeomCam &geom, const MBadPixelsCam &bad) { const Float_t meanphotRelVar = fMeanFluxPhotonsInnerPixelVar /( fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel ); TIter Next(fPixels); MCalibrationChargePix *pix; while ((pix=(MCalibrationChargePix*)Next())) { const Int_t idx = pix->GetPixId(); if (!(bad[idx].IsCalibrationResultOK())) { pix->SetFFactorMethodValid(kFALSE); continue; } if (!(pix->IsFFactorMethodValid())) continue; const Float_t ratio = geom.GetPixRatio(idx); // // Calculate the conversion factor between PHOTONS and FADC counts // // Nphot = Nphe / avQE // conv = Nphot / FADC counts // Float_t conv; if (ratio == 1.) conv = fMeanFluxPhotonsInnerPixel / pix->GetMean(); else conv = fMeanFluxPhotonsOuterPixel / pix->GetMean(); if (conv <= 0.) { pix->SetFFactorMethodValid(kFALSE); continue; } const Float_t chargeRelVar = pix->GetMeanErr() * pix->GetMeanErr() / ( pix->GetMean() * pix->GetMean()); const Float_t rsigmaRelVar = pix->GetRSigmaErr() * pix->GetRSigmaErr() / (pix->GetRSigma() * pix->GetRSigma()) ; const Float_t convrelvar = meanphotRelVar + chargeRelVar; if (convrelvar > fConvFFactorRelVarLimit) { *fLog << warn << GetDescriptor() << ": Conversion Factor F-Factor Method Rel. Variance: " << convrelvar << " above limit of: " << fConvFFactorRelVarLimit << " in pixel: " << idx << endl; pix->SetFFactorMethodValid(kFALSE); continue; } // // Calculate the Total F-Factor of the camera (in photons) // const Float_t totalFFactor = (pix->GetRSigma()/pix->GetMean()) *TMath::Sqrt(fMeanFluxPhotonsInnerPixel); // // Calculate the error of the Total F-Factor of the camera ( in photons ) // const Float_t totalFFactorVar = rsigmaRelVar + chargeRelVar + meanphotRelVar; if (convrelvar > 0. && conv > 0.) pix->SetConversionFFactorMethod(conv, TMath::Sqrt(convrelvar)*conv, totalFFactor*TMath::Sqrt(conv)); pix->SetTotalFFactorFFactorMethod( totalFFactor ); if (totalFFactorVar > 0.) pix->SetTotalFFactorFFactorMethodErr(TMath::Sqrt(totalFFactorVar)); pix->SetFFactorMethodValid(); } } void MCalibrationChargeCam::ApplyBlindPixelCalibration(const MGeomCam &geom, const MBadPixelsCam &bad, const MCalibrationChargeBlindPix &blindpix) { Float_t flux = blindpix.GetFluxInsidePlexiglass(); Float_t fluxerr = blindpix.GetFluxInsidePlexiglassErr(); TIter Next(fPixels); MCalibrationChargePix *pix; while ((pix=(MCalibrationChargePix*)Next())) { const Int_t idx = pix->GetPixId(); if (!(bad[idx].IsCalibrationResultOK())) { pix->SetBlindPixelMethodValid(kFALSE); continue; } const Float_t charge = pix->GetMean(); const Float_t area = geom[idx].GetA(); const Float_t chargeerr = pix->GetMeanErr(); const Float_t nphot = flux * area; const Float_t nphoterr = fluxerr * area; const Float_t conversion = nphot/charge; Float_t conversionerr; conversionerr = nphoterr/charge * nphoterr/charge ; conversionerr += chargeerr/charge * chargeerr/charge * conversion*conversion; conversionerr = TMath::Sqrt(conversionerr); const Float_t conversionsigma = 0.; pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma); if (conversionerr/conversion < 0.1) pix->SetBlindPixelMethodValid(); } } void MCalibrationChargeCam::ApplyPINDiodeCalibration(const MGeomCam &geom, const MBadPixelsCam &bad, const MCalibrationChargePINDiode &pindiode) { Float_t flux = pindiode.GetFluxOutsidePlexiglass(); Float_t fluxerr = pindiode.GetFluxOutsidePlexiglassErr(); TIter Next(fPixels); MCalibrationChargePix *pix; while ((pix=(MCalibrationChargePix*)Next())) { const Int_t idx = pix->GetPixId(); if (!(bad[idx].IsCalibrationResultOK())) { pix->SetPINDiodeMethodValid(kFALSE); continue; } const Float_t charge = pix->GetMean(); const Float_t area = geom[idx].GetA(); const Float_t chargeerr = pix->GetMeanErr(); const Float_t nphot = flux * area; const Float_t nphoterr = fluxerr * area; const Float_t conversion = nphot/charge; Float_t conversionerr; conversionerr = nphoterr/charge * nphoterr/charge ; conversionerr += chargeerr/charge * chargeerr/charge * conversion*conversion; if (conversionerr > 0.) conversionerr = TMath::Sqrt(conversionerr); const Float_t conversionsigma = 0.; pix->SetConversionPINDiodeMethod(conversion, conversionerr, conversionsigma); if (conversionerr/conversion < 0.1) pix->SetPINDiodeMethodValid(); } } Bool_t MCalibrationChargeCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx]; mean = pix.GetMeanConversionBlindPixelMethod(); err = pix.GetConversionBlindPixelMethodErr(); sigma = pix.GetSigmaConversionBlindPixelMethod(); return kTRUE; } Bool_t MCalibrationChargeCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx]; Float_t conv = pix.GetMeanConversionFFactorMethod(); if (conv < 0.) return kFALSE; mean = conv; err = pix.GetConversionFFactorMethodErr(); sigma = pix.GetSigmaConversionFFactorMethod(); return kTRUE; } //----------------------------------------------------------------------------------- // // Calculates the conversion factor between the integral of FADCs slices // (as defined in the signal extractor MExtractSignal.cc) // and the number of photons reaching the plexiglass for one Inner Pixel // Bool_t MCalibrationChargeCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx]; mean = pix.GetMeanConversionPINDiodeMethod(); err = pix.GetConversionPINDiodeMethodErr(); sigma = pix.GetSigmaConversionPINDiodeMethod(); return kFALSE; } //----------------------------------------------------------------------------------- // // Calculates the best combination of the three used methods possible // between the integral of FADCs slices // (as defined in the signal extractor MExtractSignal.cc) // and the number of photons reaching one Inner Pixel. // The procedure is not yet defined. // Bool_t MCalibrationChargeCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { return kFALSE; } /* void MCalibrationChargeCam::DrawHiLoFits() { if (!fOffsets) fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.); if (!fSlopes) fSlopes = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.); if (!fOffvsSlope) fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.); TIter Next(fPixels); MCalibrationPix *pix; MHCalibrationPixel *hist; while ((pix=(MCalibrationPix*)Next())) { hist = pix->GetHist(); hist->FitHiGainvsLoGain(); fOffsets->Fill(hist->GetOffset(),1.); fSlopes->Fill(hist->GetSlope(),1.); fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.); } TCanvas *c1 = new TCanvas(); c1->Divide(1,3); c1->cd(1); fOffsets->Draw(); gPad->Modified(); gPad->Update(); c1->cd(2); fSlopes->Draw(); gPad->Modified(); gPad->Update(); c1->cd(3); fOffvsSlope->Draw("col1"); gPad->Modified(); gPad->Update(); } */