Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3027) +++ /trunk/MagicSoft/Mars/Changelog (revision 3028) @@ -11,6 +11,10 @@ * mcalib/MCalibrationPix.[h,cc] + * mcalib/MCalibrationCam.cc - reshuffled HiLoGainConversion, calculation of F-Factor method, - calculation of total F-Factor of the readout. + + * mcalib/MCalibrationCam.cc + - include relative times to display with GetPixelContent Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3027) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3028) @@ -35,7 +35,8 @@ // - 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 +// with the Munich definition of the F-Factor, thus: +// F = Sigma(Out)/Mean(Out) * Mean(In)/Sigma(In) +// Mean F-Factor = 1.15 +// Error F-Factor = 0.02 // ///////////////////////////////////////////////////////////////////////////// @@ -50,13 +51,14 @@ using namespace std; -const Float_t MCalibrationPix::gkElectronicPedRms = 1.5; -const Float_t MCalibrationPix::gkErrElectronicPedRms = 0.3; -const Float_t MCalibrationPix::gkFFactor = 1.32; -const Float_t MCalibrationPix::gkFFactorError = 0.04; -const Float_t MCalibrationPix::gkChargeLimit = 3.; -const Float_t MCalibrationPix::gkChargeErrLimit = 0.; -const Float_t MCalibrationPix::gkChargeRelErrLimit = 1.; -const Float_t MCalibrationPix::gkTimeLimit = 1.5; -const Float_t MCalibrationPix::gkTimeErrLimit = 3.; +const Float_t MCalibrationPix::gkElectronicPedRms = 1.5; +const Float_t MCalibrationPix::gkErrElectronicPedRms = 0.3; +const Float_t MCalibrationPix::gkFFactor = 1.15; +const Float_t MCalibrationPix::gkFFactorError = 0.02; +const Float_t MCalibrationPix::gkChargeLimit = 3.; +const Float_t MCalibrationPix::gkChargeErrLimit = 0.; +const Float_t MCalibrationPix::gkChargeRelErrLimit = 1.; +const Float_t MCalibrationPix::gkTimeLimit = 1.5; +const Float_t MCalibrationPix::gkTimeErrLimit = 3.; +const Float_t MCalibrationPix::gkConvFFactorRelErrorLimit = 0.1; // -------------------------------------------------------------------------- @@ -115,12 +117,18 @@ fSigmaCharge = -1.; fErrSigmaCharge = -1.; - fRSigmaSquare = -1.; + fRSigmaCharge = -1.; + fErrRSigmaCharge = -1.; + fChargeProb = -1.; fPed = -1.; fPedRms = -1.; fErrPedRms = 0.; - fTime = -1.; - fErrTime = -1.; - fSigmaTime = -1.; + + fNumHiGainSamples = -1.; + fNumLoGainSamples = -1.; + + fMeanTimeOffset = -1.; + fMeanTimeOffsetError = -1.; + fTimingPrecision = -1.; fTimeProb = -1.; fTimeFirstHiGain = 0 ; @@ -134,13 +142,18 @@ fPheFFactorMethod = -1.; fPheFFactorMethodError = -1.; - fConversionFFactorMethod = -1.; - fConversionBlindPixelMethod = -1.; - fConversionPINDiodeMethod = -1.; - fConversionErrorFFactorMethod = -1.; - fConversionErrorBlindPixelMethod = -1.; - fConversionErrorPINDiodeMethod = -1.; - fConversionSigmaFFactorMethod = -1.; - fConversionSigmaBlindPixelMethod = -1.; - fConversionSigmaPINDiodeMethod = -1.; + + fMeanConversionFFactorMethod = -1.; + fMeanConversionBlindPixelMethod = -1.; + fMeanConversionPINDiodeMethod = -1.; + + fErrorConversionFFactorMethod = -1.; + fErrorConversionBlindPixelMethod = -1.; + fErrorConversionPINDiodeMethod = -1.; + + fSigmaConversionFFactorMethod = -1.; + fSigmaConversionBlindPixelMethod = -1.; + fSigmaConversionPINDiodeMethod = -1.; + + fFactorCalculated = kFALSE; } @@ -160,5 +173,6 @@ // Set the pedestals from outside // -void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms) +void MCalibrationPix::SetPedestal(const Float_t ped, const Float_t pedrms, + const Float_t higainsamp, const Float_t logainsamp ) { @@ -166,4 +180,7 @@ fPedRms = pedrms; + fNumHiGainSamples = higainsamp; + fNumLoGainSamples = logainsamp; + } @@ -174,7 +191,7 @@ void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig) { - fConversionFFactorMethod = c; - fConversionErrorFFactorMethod = err; - fConversionSigmaFFactorMethod = sig; + fMeanConversionFFactorMethod = c; + fErrorConversionFFactorMethod = err; + fSigmaConversionFFactorMethod = sig; } @@ -186,7 +203,7 @@ void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig) { - fConversionBlindPixelMethod = c; - fConversionErrorBlindPixelMethod = err; - fConversionSigmaBlindPixelMethod = sig; + fMeanConversionBlindPixelMethod = c; + fErrorConversionBlindPixelMethod = err; + fSigmaConversionBlindPixelMethod = sig; } @@ -197,7 +214,7 @@ void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig) { - fConversionPINDiodeMethod = c ; - fConversionErrorPINDiodeMethod = err; - fConversionSigmaPINDiodeMethod = sig; + fMeanConversionPINDiodeMethod = c ; + fErrorConversionPINDiodeMethod = err; + fSigmaConversionPINDiodeMethod = sig; } @@ -310,4 +327,85 @@ } +Float_t MCalibrationPix::GetPheFFactorMethod() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + + return fPheFFactorMethod; + +} + +Float_t MCalibrationPix::GetPheFFactorMethodError() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + + return fPheFFactorMethodError; + +} + + +Float_t MCalibrationPix::GetTotalFFactor() +{ + if (!fFactorCalculated) + CalcFFactorMethod(); + + return (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod); +} + +Float_t MCalibrationPix::GetTotalFFactorError() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + + const Float_t rsigmaChargeRelErrSquare = fErrRSigmaCharge * fErrRSigmaCharge + / (fRSigmaCharge * fRSigmaCharge) ; + const Float_t rChargeRelErrSquare = fErrCharge * fErrCharge + / (fCharge * fCharge) ; + const Float_t rPheRelErrSquare = fPheFFactorMethodError * fPheFFactorMethodError + / (fPheFFactorMethod * fPheFFactorMethod) ; + + return TMath::Sqrt(rsigmaChargeRelErrSquare+rChargeRelErrSquare+rPheRelErrSquare); +} + + +Float_t MCalibrationPix::GetMeanConversionFFactorMethod() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + + return fMeanConversionFFactorMethod; + +} + +Float_t MCalibrationPix::GetErrorConversionFFactorMethod() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + + return fErrorConversionFFactorMethod; + +} + +Float_t MCalibrationPix::GetSigmaConversionFFactorMethod() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + + return fSigmaConversionFFactorMethod; + +} + + +Float_t MCalibrationPix::GetTimingPrecisionError() const +{ + return fMeanTimeOffsetError/2.; +} @@ -337,6 +435,10 @@ } -Bool_t MCalibrationPix::IsFFactorMethodValid() const -{ +Bool_t MCalibrationPix::IsFFactorMethodValid() +{ + + if (!fFactorCalculated) + CalcFFactorMethod(); + return TESTBIT(fFlags, kFFactorMethodValid); } @@ -403,12 +505,10 @@ if (fHist->UseLoGain()) SetHiGainSaturation(); - + // // 4) Fit the Lo Gain histograms with a Gaussian // - if(fHist->FitCharge()) - { - SETBIT(fFlags,kFitted); - } + if (fHist->FitCharge()) + SETBIT(fFlags,kFitted); else { @@ -422,4 +522,5 @@ // // 6) Retrieve the results and store them in this class + // If fFitted is false, we get the means and RMS of the histogram!! // fCharge = fHist->GetChargeMean(); @@ -429,4 +530,10 @@ fChargeProb = fHist->GetChargeProb(); + // + // Calculate the conversion factors + // + if (TESTBIT(fFlags,kHiGainSaturation)) + ApplyLoGainConversion(); + if (CheckChargeFitValidity()) SETBIT(fFlags,kChargeFitValid); @@ -437,125 +544,157 @@ } - // - // 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.)) - { + return kTRUE; + +} + +// +// Calculate the number of photo-electrons after the F-Factor method +// Calculate the errors of the F-Factor method +// +Bool_t MCalibrationPix::CalcFFactorMethod() +{ + + if ( (fCharge == -1.) + || (fErrCharge < 0.) + || (fSigmaCharge < 0.) + || (fPedRms < 0.) ) + { + *fLog << warn << GetDescriptor() << "Cannot calculate the FFactor Method! " + << "Some of the needed parameters are not available "; + CLRBIT(fFlags,kFFactorMethodValid); + return kFALSE; + } + + // + // 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 ffactorsquare = gkFFactor * gkFFactor; + const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorError * gkFFactorError / ffactorsquare; + // + // Now the absolute error squares + // + const Float_t sigmaSquare = fSigmaCharge*fSigmaCharge; + const Float_t sigmaSquareErrSquare = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare; + + Float_t pedRmsSquare = fPedRms* fPedRms; + Float_t pedRmsSquareErrSquare = 4.*fErrPedRms*fErrPedRms * pedRmsSquare; + + if (!TESTBIT(fFlags,kHiGainSaturation)) + { /* HiGain */ + + pedRmsSquare *= fNumHiGainSamples; + pedRmsSquareErrSquare *= fNumHiGainSamples*fNumHiGainSamples; + } + else + { /* LoGain */ // - // Square all variables in order to avoid applications of square root + // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it + // from the HI GAIN (all calculation per slice up to now): // - // First the relative error squares + // We extract the pure NSB contribution: // - const Float_t chargeSquare = fCharge* fCharge; - const Float_t chargeSquareRelErrSquare = 4.*fErrCharge*fErrCharge / chargeSquare; - - const Float_t fFactorRelErrSquare = gkFFactorError * gkFFactorError / (gkFFactor * gkFFactor); + const Float_t elecRmsSquare = gkElectronicPedRms *gkElectronicPedRms; + const Float_t elecRmsSquareErrSquare = 4.*gkErrElectronicPedRms*gkErrElectronicPedRms * elecRmsSquare; + + Float_t nsbSquare = pedRmsSquare - elecRmsSquare; + Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare) + / (nsbSquare * nsbSquare) ; + + if (nsbSquare < 0.) + nsbSquare = 0.; + // - // Now the absolute error squares + // Now, we divide the NSB by the conversion factor and + // add it quadratically to the electronic noise // - const Float_t sigmaSquare = fSigmaCharge* fSigmaCharge; - const Float_t sigmaSquareErrSquare = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare; - - const Float_t elecRmsSquare = gkElectronicPedRms* gkElectronicPedRms; - const Float_t elecRmsSquareErrSquare = 4.*gkErrElectronicPedRms*gkErrElectronicPedRms * 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) */ + const Float_t conversionSquare = fConversionHiLo *fConversionHiLo; + const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoError*fConversionHiLoError/conversionSquare; + + const Float_t convertedNsbSquare = nsbSquare / conversionSquare; + const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare) + * convertedNsbSquare * convertedNsbSquare; + + pedRmsSquare = convertedNsbSquare + elecRmsSquare; + pedRmsSquareErrSquare = convertedNsbSquareErrSquare + elecRmsSquareErrSquare; // - // Calculate the reduced sigmas + // Now, correct for the number of used FADC slices in the LoGain: // - 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) ; - + pedRmsSquare *= fNumLoGainSamples; + pedRmsSquareErrSquare *= fNumLoGainSamples*fNumLoGainSamples; // - // Calculate the number of phe's from the F-Factor method - // (independent on Hi Gain or Lo Gain) + // Correct also for the conversion to Hi-Gain: // - fPheFFactorMethod = gkFFactor * 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 ( IsChargeFitValid() && - (fConversionFFactorMethod > 0.) && - (fConversionErrorFFactorMethod/fConversionFFactorMethod < 0.1) ) - SETBIT(fFlags,kFFactorMethodValid); - else - CLRBIT(fFlags,kFFactorMethodValid); - - } /* if ((fPed > 0.) && (fPedRms > 0.)) */ + pedRmsSquare *= fConversionHiLo*fConversionHiLo; + pedRmsSquareErrSquare *= fConversionHiLo*fConversionHiLo*fConversionHiLo*fConversionHiLo; + + } /* if (HiGainSaturation) */ + + // + // Calculate the reduced sigmas + // + const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare; + if (rsigmachargesquare <= 0.) + { + *fLog << warn + << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel " + << fPixId << endl; + CLRBIT(fFlags,kFFactorMethodValid); + fFactorCalculated = kTRUE; + return kFALSE; + } + + const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare) + / (rsigmachargesquare * rsigmachargesquare) ; + + fRSigmaCharge = TMath::Sqrt(rsigmachargesquare); + fErrRSigmaCharge = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare); + + + // + // Calculate the number of phe's from the F-Factor method + // (independent on Hi Gain or Lo Gain) + // + fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmachargesquare; + + const Float_t pheFFactorRelErrSquare = ffactorsquareRelErrSquare + + chargeSquareRelErrSquare + + rSigmaSquareRelErrSquare ; + + fPheFFactorMethodError = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod; + + const Float_t chargeRelErrSquare = fErrCharge*fErrCharge / (fCharge * fCharge); + + fMeanConversionFFactorMethod = fPheFFactorMethod / fCharge ; + fErrorConversionFFactorMethod = ( pheFFactorRelErrSquare + chargeRelErrSquare ) + * fMeanConversionFFactorMethod * fMeanConversionFFactorMethod; + + const Float_t convrelerror = fErrorConversionFFactorMethod/fMeanConversionFFactorMethod; + + if ( (fMeanConversionFFactorMethod > 0.) && (convrelerror < gkConvFFactorRelErrorLimit)) + SETBIT(fFlags,kFFactorMethodValid); + + fFactorCalculated = kTRUE; + + fSigmaConversionFFactorMethod = GetTotalFFactor()*TMath::Sqrt(fMeanConversionFFactorMethod); return kTRUE; - -} - -// -// The check return kTRUE if: -// -// 0) No value is nan +} + + +// +// The check returns kTRUE if: +// +// 0) Pixel has BIT fitted set: +// This means: +// a) No result is a nan +// b) The NDF is not smaller than fNDFLimit (5) +// c) The Probability is greater than gkProbLimit (default 0.001 == 99.9%) // 1) Pixel has a fitted charge greater than 3*PedRMS // 2) Pixel has a fit error greater than 0. @@ -567,24 +706,19 @@ { - if (TMath::IsNaN(fCharge) - || 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,kFitted)) + return kFALSE; if (TESTBIT(fFlags,kExcludeQualityCheck)) return kTRUE; - Float_t equivpedestal = GetPedRms(); - - if (TESTBIT(fFlags,kHiGainSaturation)) - equivpedestal /= fConversionHiLo; + Float_t pedestal; + + if (!TESTBIT(fFlags,kHiGainSaturation)) /* higain */ + pedestal = GetPedRms()*TMath::Sqrt(fNumHiGainSamples); + else /* logain */ + pedestal = GetPedRms()*TMath::Sqrt(fNumLoGainSamples); - if (fCharge < gkChargeLimit*equivpedestal) + + if (fCharge < gkChargeLimit*pedestal) { *fLog << warn << "WARNING: Fitted Charge is smaller than " @@ -614,5 +748,5 @@ } - if (fSigmaCharge < equivpedestal) + if (fSigmaCharge < pedestal) { *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " @@ -635,7 +769,7 @@ { - if (TMath::IsNaN(fTime) - || TMath::IsNaN(fErrTime) - || TMath::IsNaN(fSigmaTime) + if (TMath::IsNaN(fMeanTimeOffset) + || TMath::IsNaN(fMeanTimeOffsetError) + || TMath::IsNaN(fTimingPrecision) || TMath::IsNaN(fTimeProb)) { @@ -648,15 +782,15 @@ return kTRUE; - if (TMath::Abs(fTime) > gkTimeLimit) - { - *fLog << warn << "WARNING: Abs(Fitted Rel. Time) is greater than " - << gkTimeLimit << " in Pixel " << fPixId << endl; - return kFALSE; - } - - if (fErrTime > gkTimeErrLimit) - { - *fLog << warn << "WARNING: Error of Fitted Time is smaller than " - << gkTimeErrLimit << " in Pixel " << fPixId << endl; + if (TMath::Abs(fMeanTimeOffset) > gkTimeLimit) + { + *fLog << warn << "WARNING: Abs(Fitted Rel. Time) " << TMath::Abs(fMeanTimeOffset) + << " is greater than " << gkTimeLimit << " in Pixel " << fPixId << endl; + return kFALSE; + } + + if (fMeanTimeOffsetError > gkTimeErrLimit) + { + *fLog << warn << "WARNING: Error of Fitted Time " << fMeanTimeOffsetError + << " is smaller than " << gkTimeErrLimit << " in Pixel " << fPixId << endl; return kFALSE; } @@ -727,10 +861,6 @@ Bool_t MCalibrationPix::CheckOscillations() { - - return kTRUE; } - - void MCalibrationPix::ApplyLoGainConversion() @@ -773,8 +903,8 @@ } - fTime = fHist->GetRelTimeMean(); - fErrTime = fHist->GetRelTimeMeanErr(); - fSigmaTime = fHist->GetRelTimeSigma(); - fTimeProb = fHist->GetRelTimeProb(); + fMeanTimeOffset = fHist->GetRelTimeMean(); + fMeanTimeOffsetError = fHist->GetRelTimeMeanErr(); + fTimingPrecision = fHist->GetRelTimeSigma(); + fTimeProb = fHist->GetRelTimeProb(); fAbsTimeMean = fHist->GetAbsTimeMean();