/* ======================================================================== *\ ! ! * ! * 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-2001 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MCalibrationPix // // // // This is the storage container to hold informations about the pedestal // // (offset) value of one Pixel (PMT). // // // ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationPix.h" #include "MCalibrationConfig.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MCalibrationPix); using namespace std; // -------------------------------------------------------------------------- // // Default Constructor: // // The following values are initialized to meaningful values: // // - The Electronic Rms to 1.5 per FADC slice // - The uncertainty about the Electronic RMS to 0.3 per slice // - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17. // We use here the Square of the Munich definition, thus: // Mean F-Factor = 1.15*1.15 = 1.32 // Error F-Factor = 2.*0.02 = 0.04 // MCalibrationPix::MCalibrationPix(const char *name, const char *title) : fPixId(-1), fElectronicPedRms(1.5), fErrElectronicPedRms(0.3), fFactor(1.32), fFactorError(0.04), fChargeLimit(3.), fChargeErrLimit(0.), fChargeRelErrLimit(1.), fFlags(0) { fName = name ? name : "MCalibrationPixel"; fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results"; // // At the moment, we don't have a database, yet, // so we get it from the configuration file // fConversionHiLo = gkConversionHiLo; fConversionHiLoError = gkConversionHiLoError; fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel "); if (!fHist) *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl; Clear(); } MCalibrationPix::~MCalibrationPix() { delete fHist; } // ------------------------------------------------------------------------ // // Invalidate values // void MCalibrationPix::Clear(Option_t *o) { fHist->Reset(); CLRBIT(fFlags, kHiGainSaturation); CLRBIT(fFlags, kExcluded); CLRBIT(fFlags, kFitValid); CLRBIT(fFlags, kFitted); CLRBIT(fFlags, kBlindPixelMethodValid); CLRBIT(fFlags, kFFactorMethodValid); CLRBIT(fFlags, kPINDiodeMethodValid); fCharge = -1.; fErrCharge = -1.; fSigmaCharge = -1.; fErrSigmaCharge = -1.; fRSigmaSquare = -1.; fChargeProb = -1.; fPed = -1.; fPedRms = -1.; fErrPedRms = 0.; fTime = -1.; fSigmaTime = -1.; fTimeChiSquare = -1.; fPheFFactorMethod = -1.; fPheFFactorMethodError = -1.; fConversionFFactorMethod = -1.; fConversionBlindPixelMethod = -1.; fConversionPINDiodeMethod = -1.; fConversionErrorFFactorMethod = -1.; fConversionErrorBlindPixelMethod = -1.; fConversionErrorPINDiodeMethod = -1.; fConversionSigmaFFactorMethod = -1.; fConversionSigmaBlindPixelMethod = -1.; fConversionSigmaPINDiodeMethod = -1.; } void MCalibrationPix::DefinePixId(Int_t i) { fPixId = i; fHist->ChangeHistId(i); } // -------------------------------------------------------------------------- // // Set the pedestals from outside // void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms) { fPed = ped; fPedRms = pedrms; } // -------------------------------------------------------------------------- // // Set the conversion factors from outside (only for MC) // void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig) { fConversionFFactorMethod = c; fConversionErrorFFactorMethod = err; fConversionSigmaFFactorMethod = sig; } // -------------------------------------------------------------------------- // // Set the conversion factors from outside (only for MC) // void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig) { fConversionBlindPixelMethod = c; fConversionErrorBlindPixelMethod = err; fConversionSigmaBlindPixelMethod = sig; } // -------------------------------------------------------------------------- // // Set the conversion factors from outside (only for MC) // void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig) { fConversionPINDiodeMethod = c ; fConversionErrorPINDiodeMethod = err; fConversionSigmaPINDiodeMethod = sig; } // -------------------------------------------------------------------------- // // Set the Hi Gain Saturation Bit from outside (only for MC) // void MCalibrationPix::SetHiGainSaturation(Bool_t b) { if (b) { SETBIT(fFlags, kHiGainSaturation); fHist->SetUseLoGain(1); } else { CLRBIT(fFlags, kHiGainSaturation); fHist->SetUseLoGain(0); } } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetExcluded(Bool_t b ) { b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetExcludeQualityCheck(Bool_t b ) { b ? SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetFitValid(Bool_t b ) { b ? SETBIT(fFlags, kFitValid) : CLRBIT(fFlags, kFitValid); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetFitted(Bool_t b ) { b ? SETBIT(fFlags, kFitted) : CLRBIT(fFlags, kFitted); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetBlindPixelMethodValid(Bool_t b ) { b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetFFactorMethodValid(Bool_t b ) { b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationPix::SetPINDiodeMethodValid(Bool_t b ) { b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); } Bool_t MCalibrationPix::IsExcluded() const { return TESTBIT(fFlags,kExcluded); } Bool_t MCalibrationPix::IsFitValid() const { return TESTBIT(fFlags, kFitValid); } Bool_t MCalibrationPix::IsFitted() const { return TESTBIT(fFlags, kFitted); } Bool_t MCalibrationPix::IsBlindPixelMethodValid() const { return TESTBIT(fFlags, kBlindPixelMethodValid); } Bool_t MCalibrationPix::IsFFactorMethodValid() const { return TESTBIT(fFlags, kFFactorMethodValid); } Bool_t MCalibrationPix::IsPINDiodeMethodValid() const { return TESTBIT(fFlags, kPINDiodeMethodValid); } // -------------------------------------------------------------------------- // // 1) Return if the charge distribution is already succesfully fitted // or if the histogram is empty // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid // possible remaining cosmics to spoil the fit. // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram // 4) Fit the histograms with a Gaussian // 5) In case of failure set the bit kFitted to false // 6) Retrieve the results and store them in this class // 7) Calculate the number of photo-electrons after the F-Factor method // 8) Calculate the errors of the F-Factor method // // The fits are declared valid (fFitValid = kTRUE), if: // // 1) Pixel has a fitted charge greater than 3*PedRMS // 2) Pixel has a fit error greater than 0. // 3) Pixel has a fit Probability greater than 0.0001 // 4) Pixel has a charge sigma bigger than its Pedestal RMS // 5) If FitTimes is used, // the mean arrival time is at least 1.0 slices from the used edge slices // (this stage is only performed in the times fit) // // If the histogram is empty, all values are set to -1. // // The conversion factor after the F-Factor method is declared valid, if: // // 1) fFitValid is kTRUE // 2) Conversion Factor is bigger than 0. // 3) The error of the conversion factor is smaller than 10% // Bool_t MCalibrationPix::FitCharge() { // // 1) Return if the charge distribution is already succesfully fitted // or if the histogram is empty // if (fHist->IsFitOK() || fHist->IsEmpty()) return kTRUE; // // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid // possible remaining cosmics to spoil the fit. // // if (fPed && fPedRms) // fHist->SetLowerFitRange(1.5*fPedRms); // else // *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl; // // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram // if (fHist->UseLoGain()) SetHiGainSaturation(); // // 4) Fit the Lo Gain histograms with a Gaussian // if(fHist->FitCharge()) { SETBIT(fFlags,kFitted); } else { *fLog << warn << "WARNING: Could not fit charges of pixel " << fPixId << endl; // // 5) In case of failure set the bit kFitted to false // CLRBIT(fFlags,kFitted); } // // 6) Retrieve the results and store them in this class // fCharge = fHist->GetChargeMean(); fErrCharge = fHist->GetChargeMeanErr(); fSigmaCharge = fHist->GetChargeSigma(); fErrSigmaCharge = fHist->GetChargeSigmaErr(); fChargeProb = fHist->GetChargeProb(); if (CheckChargeFitValidity()) SETBIT(fFlags,kFitValid); else { CLRBIT(fFlags,kFitValid); return kFALSE; } // // 7) Calculate the number of photo-electrons after the F-Factor method // 8) Calculate the errors of the F-Factor method // if ((fPed > 0.) && (fPedRms > 0.)) { // // Square all variables in order to avoid applications of square root // // First the relative error squares // const Float_t chargeSquare = fCharge* fCharge; const Float_t chargeSquareRelErrSquare = 4.*fErrCharge*fErrCharge / chargeSquare; const Float_t fFactorRelErrSquare = fFactorError * fFactorError / (fFactor * fFactor); // // Now the absolute error squares // const Float_t sigmaSquare = fSigmaCharge* fSigmaCharge; const Float_t sigmaSquareErrSquare = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare; const Float_t elecRmsSquare = fElectronicPedRms* fElectronicPedRms; const Float_t elecRmsSquareErrSquare = 4.*fErrElectronicPedRms*fErrElectronicPedRms * elecRmsSquare; Float_t pedRmsSquare = fPedRms* fPedRms; Float_t pedRmsSquareErrSquare = 4.*fErrPedRms*fErrPedRms * pedRmsSquare; if (TESTBIT(fFlags,kHiGainSaturation)) { // // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it // from the Hi Gain: // // We extract the pure NSB contribution: // Float_t nsbSquare = pedRmsSquare - elecRmsSquare; Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare) / (nsbSquare * nsbSquare) ; if (nsbSquare < 0.) nsbSquare = 0.; // // Now, we divide the NSB by the conversion factor and // add it quadratically to the electronic noise // const Float_t conversionSquare = fConversionHiLo *fConversionHiLo; const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoError*fConversionHiLoError/conversionSquare; // // Calculate the new "Pedestal RMS" // const Float_t convertedNsbSquare = nsbSquare / conversionSquare; const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare) * convertedNsbSquare * convertedNsbSquare; pedRmsSquare = convertedNsbSquare + elecRmsSquare; pedRmsSquareErrSquare = convertedNsbSquareErrSquare + elecRmsSquareErrSquare; } /* if (kHiGainSaturation) */ // // Calculate the reduced sigmas // fRSigmaSquare = sigmaSquare - pedRmsSquare; if (fRSigmaSquare <= 0.) { *fLog << warn << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel " << fPixId << endl; if (TESTBIT(fFlags,kHiGainSaturation)) ApplyLoGainConversion(); return kFALSE; } const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare) / (fRSigmaSquare * fRSigmaSquare) ; // // Calculate the number of phe's from the F-Factor method // (independent on Hi Gain or Lo Gain) // fPheFFactorMethod = fFactor * chargeSquare / fRSigmaSquare; const Float_t pheFFactorRelErrSquare = fFactorRelErrSquare + chargeSquareRelErrSquare + rSigmaSquareRelErrSquare ; fPheFFactorMethodError = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod; // // Calculate the conversion factors // if (TESTBIT(fFlags,kHiGainSaturation)) ApplyLoGainConversion(); const Float_t chargeRelErrSquare = fErrCharge*fErrCharge / (fCharge * fCharge); fConversionFFactorMethod = fPheFFactorMethod / fCharge ; fConversionErrorFFactorMethod = ( pheFFactorRelErrSquare + chargeRelErrSquare ) * fConversionFFactorMethod * fConversionFFactorMethod; if ( IsFitValid() && (fConversionFFactorMethod > 0.) && (fConversionErrorFFactorMethod/fConversionFFactorMethod < 0.1) ) SETBIT(fFlags,kFFactorMethodValid); else CLRBIT(fFlags,kFFactorMethodValid); } /* if ((fPed > 0.) && (fPedRms > 0.)) */ return kTRUE; } // // The check return kTRUE if: // // 0) No value is nan // 1) Pixel has a fitted charge greater than 3*PedRMS // 2) Pixel has a fit error greater than 0. // 3) Pixel has a fitted charge greater its charge error // 4) Pixel has a fit Probability greater than 0.0001 // 5) Pixel has a charge sigma bigger than its Pedestal RMS // Bool_t MCalibrationPix::CheckChargeFitValidity() { if (TMath::IsNaN(fCharge) || TMath::IsNaN(fErrCharge) || TMath::IsNaN(fErrCharge) || TMath::IsNaN(fSigmaCharge) || TMath::IsNaN(fErrSigmaCharge) || TMath::IsNaN(fChargeProb)) { *fLog << warn << "WARNING: Some of the fit values are NAN in Pixel " << fPixId << endl; return kFALSE; } if (TESTBIT(fFlags,kExcludeQualityCheck)) return kTRUE; Float_t equivpedestal = GetPedRms(); if (TESTBIT(fFlags,kHiGainSaturation)) equivpedestal /= fConversionHiLo; if (fCharge < fChargeLimit*equivpedestal) { *fLog << warn << "WARNING: Fitted Charge is smaller than " << fChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl; return kFALSE; } if (fErrCharge < fChargeErrLimit) { *fLog << warn << "WARNING: Error of Fitted Charge is smaller than " << fChargeErrLimit << " in Pixel " << fPixId << endl; return kFALSE; } if (fCharge < fChargeRelErrLimit*fErrCharge) { *fLog << warn << "WARNING: Fitted Charge is smaller than " << fChargeRelErrLimit << "* its error in Pixel " << fPixId << endl; return kFALSE; } if (!fHist->IsFitOK()) { *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel " << fPixId << endl; return kFALSE; } if (fSigmaCharge < equivpedestal) { *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " << fPixId << endl; return kFALSE; } return kTRUE; } // // The check returns kTRUE if: // // The mean arrival time is at least 1.0 slices from the used edge slices // Bool_t MCalibrationPix::CheckTimeFitValidity() { if (TESTBIT(fFlags,kExcludeQualityCheck)) return kTRUE; Float_t lowerrange; Float_t upperrange; if (TESTBIT(fFlags,kHiGainSaturation)) { lowerrange = (Float_t)fHist->GetTimeLowerFitRangeLoGain()+1.; upperrange = (Float_t)fHist->GetTimeUpperFitRangeLoGain()+1.; } else { lowerrange = (Float_t)fHist->GetTimeLowerFitRangeHiGain()+1.; upperrange = (Float_t)fHist->GetTimeUpperFitRangeHiGain()+1.; } if (fTime < lowerrange) { *fLog << warn << "WARNING: Mean Fitted Time inside or smaller than first used FADC slice in Pixel " << fPixId << " time: " << fTime << " Range: " << lowerrange << endl; return kFALSE; } if (fTime > upperrange) { *fLog << warn << "WARNING: Mean Fitted Time inside or greater than last used FADC slice in Pixel " << fPixId << " time: " << fTime << " Range: " << upperrange << endl; return kFALSE; } return kTRUE; } // // The check returns kTRUE if: // // // Bool_t MCalibrationPix::CheckOscillations() { return kTRUE; } void MCalibrationPix::ApplyLoGainConversion() { const Float_t chargeRelErrSquare = fErrCharge*fErrCharge /( fCharge * fCharge); const Float_t sigmaRelErrSquare = fErrSigmaCharge*fErrSigmaCharge /( fSigmaCharge * fSigmaCharge); const Float_t conversionRelErrSquare = fConversionHiLoError*fConversionHiLoError /(fConversionHiLo * fConversionHiLo); fCharge *= fConversionHiLo; fErrCharge = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge; fSigmaCharge *= fConversionHiLo; fErrSigmaCharge = TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge; } // -------------------------------------------------------------------------- // // 1) Fit the arrival times // 2) Retrieve the results // 3) Note that because of the low number of bins, the NDf is sometimes 0, so // Root does not give a reasonable Probability, the Chisquare is more significant // // This fit has to be done AFTER the Charges fit, // otherwise only the Hi Gain will be fitted, even if there are no entries // // Bool_t MCalibrationPix::FitTime() { // // Fit the Low Gain // if (TESTBIT(fFlags,kHiGainSaturation)) { if(!fHist->FitTimeLoGain()) { *fLog << warn << "WARNING: Could not fit Lo Gain times of pixel " << fPixId << endl; // fHist->PrintTimeFitResult(); return kFALSE; } } // // Fit the High Gain // else { if(!fHist->FitTimeHiGain()) { *fLog << warn << "WARNING: Could not fit Hi Gain times of pixel " << fPixId << endl; // fHist->PrintTimeFitResult(); return kFALSE; } } fTime = fHist->GetTimeMean(); fSigmaTime = fHist->GetTimeSigma(); fTimeChiSquare = fHist->GetTimeChiSquare(); fTimeProb = fHist->GetTimeProb(); if (CheckTimeFitValidity()) SETBIT(fFlags,kFitValid); else CLRBIT(fFlags,kFitValid); return kTRUE; }