/* ======================================================================== *\ ! ! * ! * 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), fCharge(-1.), fErrCharge(-1.), fSigmaCharge(-1.), fErrSigmaCharge(-1.), fRSigmaSquare(-1.), fChargeProb(-1.), fPed(-1.), fPedRms(-1.), fErrPedRms(0.), fElectronicPedRms(1.5), fErrElectronicPedRms(0.3), fTime(-1.), fSigmaTime(-1.), fTimeChiSquare(-1.), fFactor(1.32), fFactorError(0.04), fPheFFactorMethod(-1.), fPheFFactorMethodError(-1.), fConversionFFactorMethod(-1.), fConversionBlindPixelMethod(-1.), fConversionPINDiodeMethod(-1.), fConversionErrorFFactorMethod(-1.), fConversionErrorBlindPixelMethod(-1.), fConversionErrorPINDiodeMethod(-1.), fConversionSigmaFFactorMethod(-1.), fConversionSigmaBlindPixelMethod(-1.), fConversionSigmaPINDiodeMethod(-1.), fHiGainSaturation(kFALSE), fFitValid(kFALSE), fFitted(kFALSE), fBlindPixelMethodValid(kFALSE), fFFactorMethodValid(kFALSE), fPINDiodeMethodValid(kFALSE) { 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 "); } MCalibrationPix::~MCalibrationPix() { delete fHist; } void MCalibrationPix::DefinePixId(Int_t i) { fPixId = i; fHist->ChangeHistId(i); } // ------------------------------------------------------------------------ // // Invalidate values // void MCalibrationPix::Clear(Option_t *o) { fHist->Reset(); } // -------------------------------------------------------------------------- // // 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 print out the fit results // 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 5*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 << "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->FitChargeLoGain()) { *fLog << warn << "Could not fit Lo Gain charges of pixel " << fPixId << endl; // // 5) In case of failure print out the fit results // fHist->PrintChargeFitResult(); } } else { // // 4) Fit the Hi Gain histograms with a Gaussian // if(!fHist->FitChargeHiGain()) { *fLog << warn << "Could not fit Hi Gain charges of pixel " << fPixId << endl; // // 5) In case of failure print out the fit results // fHist->PrintChargeFitResult(); } } // // 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 (fCharge <= 0.) { *fLog << warn << "Cannot apply calibration: Mean Fitted Charges are smaller than 0 in pixel " << fPixId << endl; return kFALSE; } if (fErrCharge > 0.) fFitted = kTRUE; if (CheckChargeFitValidity()) fFitValid = kTRUE; // // 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 (fHiGainSaturation) { // // 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 (fHiGainSaturation) */ // // Calculate the reduced sigmas // fRSigmaSquare = sigmaSquare - pedRmsSquare; if (fRSigmaSquare <= 0.) { *fLog << warn << "Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel " << fPixId << endl; if (fHiGainSaturation) 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 (fHiGainSaturation) 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) ) fFFactorMethodValid = kTRUE; } /* if ((fPed > 0.) && (fPedRms > 0.)) */ return kTRUE; } // // The check return kTRUE if: // // 1) Pixel has a fitted charge greater than 5*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 // Bool_t MCalibrationPix::CheckChargeFitValidity() { Float_t equivpedestal = GetPedRms(); if (fHiGainSaturation) equivpedestal /= fConversionHiLo; if (fCharge < 5.*equivpedestal) { *fLog << warn << "WARNING: Fitted Charge is smaller than 5 Pedestal RMS in Pixel " << fPixId << endl; return kFALSE; } if (fErrCharge < 0.) { *fLog << warn << "WARNING: Error of Fitted Charge is smaller than 0 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() { Float_t lowerrange; Float_t upperrange; if (fHiGainSaturation) { 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; } 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; } // -------------------------------------------------------------------------- // // Set the pedestals from outside // void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms) { fPed = ped; fPedRms = pedrms; } // -------------------------------------------------------------------------- // // 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 (fHiGainSaturation) { if(!fHist->FitTimeLoGain()) { *fLog << warn << "Could not fit Lo Gain times of pixel " << fPixId << endl; fHist->PrintTimeFitResult(); return kFALSE; } } // // Fit the High Gain // else { if(!fHist->FitTimeHiGain()) { *fLog << warn << "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()) fFitValid = kFALSE; return kTRUE; }