Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3649) +++ /trunk/MagicSoft/Mars/Changelog (revision 3650) @@ -19,6 +19,15 @@ -*-*- END OF LINE -*-*- + 2004/04/05: Markus Gaug + + * mcalib/MCalibrationPedCam.[h,cc] + * mcalib/Makefile + * mcalib/CalibLinkDef.h + - new class to store information obtained from MHPedestalCam + + 2004/04/02: Markus Gaug + * mcalib/MCalibrationChargePix.[h,cc] * mcalib/MCalibrationPix.[h,cc] * mcalib/MCalibrationCam.[h,cc] Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc (revision 3649) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc (revision 3650) @@ -22,10 +22,9 @@ \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// -// // -// MCalibrationChargePix // -// // -// Storage container to hold informations about the calibration values // -// values of one Pixel (PMT). // -// // +// +// MCalibrationChargePix +// +// Storage container of the calibrated Charge of one pixel. +// // The following values are initialized to meaningful values: // @@ -35,54 +34,38 @@ // 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 -// -// - Average QE: (email David Paneque, 14.2.04): -// -// The conversion factor that comes purely from QE folded to a Cherenkov -// spectrum has to be multiplied by: -// * Plexiglass window -->> 0.96 X 0.96 -// * PMT photoelectron collection efficiency -->> 0.9 -// * Light guides efficiency -->> 0.94 -// -// Concerning the light guides effiency estimation... Daniel Ferenc -// is preparing some work (simulations) to estimate it. Yet so far, he has -// been busy with other stuff, and this work is still UNfinished. -// -// The estimation I did comes from: -// 1) Reflectivity of light guide walls is 85 % (aluminum) -// 2) At ZERO degree light incidence, 37% of the light hits such walls -// (0.15X37%= 5.6% of light lost) -// 3) When increasing the light incidence angle, more and more light hits -// the walls. -// -// However, the loses due to larger amount of photons hitting the walls is more -// or less counteracted by the fact that more and more photon trajectories cross -// the PMT photocathode twice, increasing the effective sensitivity of the PMT. -// -// Jurgen Gebauer did some quick measurements about this issue. I attach a -// plot. You can see that the angular dependence is (more or less) in agreement -// with a CosTheta function (below 20-25 degrees), -// which is the variation of teh entrance window cross section. So, in -// first approximation, no loses when increasing light incidence angle; -// and therefore, the factor 0.94. -// -// So, summarizing... I would propose the following conversion factors -// (while working with CT1 cal box) in order to get the final number of photons -// from the detected measured size in ADC counts. -// -// Nph = ADC * FmethodConversionFactor / ConvPhe-PhFactor -// -// FmethodConversionFactor ; measured for individual pmts -// -// ConvPhe-PhFactor = 0.98 * 0.23 * 0.90 * 0.94 * 0.96 * 0.96 = 0.18 -// -// I would not apply any smearing of this factor (which we have in nature), -// since we might be applying it to PMTs in the totally wrong direction. -// +// Mean F-Factor (gkFFactor) = 1.15 +// Error F-Factor (gkFFactorErr) = 0.02 +// +// The following variables are calculated inside this class: +// - fLoGainPedRmsSquare and LoGainPedRmsSquareVar (see CalcLoGainPedestal()) +// - fRSigmaSquare and fRSigmaSquareVar (see CalcReducedSigma() ) +// - fPheFFactorMethod and fPheFFactorMethodVar see CalcFFactorMethod() ) +// +// The following variables are set by MHCalibrationChargeCam: +// - fAbsTimeMean and fAbsTimeRms +// - all variables in MCalibrationPix +// +// The following variables are set by MCalibrationChargeCalc: +// - fPed, fPedVar and fPedRms +// - fMeanConversionFFactorMethod, fMeanConversionBlindPixelMethod, +// fMeanConversionPINDiodeMethod and fMeanConversionCombinedMethod +// - fConversionFFactorMethodVar, fConversionBlindPixelMethodVar +// fConversionPINDiodeMethodVar and fConversionCombinedMethodVar +// - fSigmaConversionFFactorMethodm, fSigmaConversionBlindPixelMethod +// fSigmaConversionPINDiodeMethod and fSigmaConversionCombinedMethod +// - fTotalFFactorFFactorMethod, fTotalFFactorBlindPixelMethod +// fTotalFFactorPINDiodeMethod and fTotalFFactorCombinedMethod +// - fTotalFFactorFFactorMethodVar, fTotalFFactorBlindPixelMethodVar, +// fTotalFFactorPINDiodeMethodVar and fTotalFFactorCombinedMethodVar +// +// The following variables are not yet implemented: +// - fConversionHiLo and fConversionHiLoVar (now set fixed to 10. +- 2.5) // // Error of all variables are calculated by error-propagation. Note that internally, // all error variables contain Variances in order to save the CPU-intensive square rooting // +// See also: MCalibrationChargeCam, MCalibrationChargeCalc, +// MHCalibrationChargeCam, MHCalibrationChargePix +// ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationChargePix.h" @@ -109,4 +92,13 @@ // Default Constructor: // +// Sets: +// - fCalibFlags to 0 +// - fConversionHiLo to fgConversionHiLo +// - fConversionHiLoErr to fgConversionHiLoErr +// - PheFFactorMethodLimit to fgPheFFactorMethodLimit +// +// Calls: +// - Clear() +// MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title) : fCalibFlags(0) @@ -131,5 +123,10 @@ // ------------------------------------------------------------------------ // -// Invalidate values +// Sets: +// - all flags to 0 +// - all variables to -1. +// +// Calls: +// - MCalibrationPix::Clear() // void MCalibrationChargePix::Clear(Option_t *o) @@ -141,6 +138,6 @@ SetCombinedMethodValid ( kFALSE ); - fRSigma = -1.; - fRSigmaVar = -1.; + fRSigmaSquare = -1.; + fRSigmaSquareVar = -1.; fPed = -1.; @@ -148,6 +145,6 @@ fPedVar = -1.; - fLoGainPedRms = -1.; - fLoGainPedRmsVar = -1.; + fLoGainPedRmsSquare = -1.; + fLoGainPedRmsSquareVar = -1.; fAbsTimeMean = -1.; @@ -183,19 +180,29 @@ // -------------------------------------------------------------------------- // -// Set the pedestals from outside -// -void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr) -{ - - fPed = ped; - fPedRms = pedrms; - fPedVar = pederr*pederr; -} - - - -// -------------------------------------------------------------------------- -// -// Set the conversion factors from outside (only for MC) +// Set the conversion factors Blind Pixel Method from outside (only for MC) +// +void MCalibrationChargePix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig) +{ + fMeanConversionBlindPixelMethod = c; + fConversionBlindPixelMethodVar = err*err; + fSigmaConversionBlindPixelMethod = sig; +} + + +// -------------------------------------------------------------------------- +// +// Set the conversion factors Combined Method from outside (only for MC) +// +void MCalibrationChargePix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig) +{ + fMeanConversionCombinedMethod = c; + fConversionCombinedMethodVar = err*err; + fSigmaConversionCombinedMethod = sig; +} + + +// -------------------------------------------------------------------------- +// +// Set the conversion factors F-Factor Method from outside (only for MC) // void MCalibrationChargePix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig) @@ -208,28 +215,5 @@ // -------------------------------------------------------------------------- // -// Set the conversion factors from outside (only for MC) -// -void MCalibrationChargePix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig) -{ - fMeanConversionCombinedMethod = c; - fConversionCombinedMethodVar = err*err; - fSigmaConversionCombinedMethod = sig; -} - - -// -------------------------------------------------------------------------- -// -// Set the conversion factors from outside (only for MC) -// -void MCalibrationChargePix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig) -{ - fMeanConversionBlindPixelMethod = c; - fConversionBlindPixelMethodVar = err*err; - fSigmaConversionBlindPixelMethod = sig; -} - -// -------------------------------------------------------------------------- -// -// Set the conversion factors from outside (only for MC) +// Set the conversion factors PIN Diode Method from outside (only for MC) // void MCalibrationChargePix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig) @@ -242,5 +226,5 @@ // -------------------------------------------------------------------------- // -// Set the Excluded Bit from outside +// Set the BlindPixelMethod Validity Bit from outside // void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b ) @@ -251,5 +235,14 @@ // -------------------------------------------------------------------------- // -// Set the Excluded Bit from outside +// Set the CombinedMethod Validity Bit from outside +// +void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b ) +{ + b ? SETBIT(fCalibFlags, kCombinedMethodValid) : CLRBIT(fCalibFlags, kCombinedMethodValid); +} + +// -------------------------------------------------------------------------- +// +// Set the FFactorMethod Validity Bit from outside // void MCalibrationChargePix::SetFFactorMethodValid(const Bool_t b ) @@ -260,5 +253,5 @@ // -------------------------------------------------------------------------- // -// Set the Excluded Bit from outside +// Set the PINDiodeMethod Validity Bit from outside // void MCalibrationChargePix::SetPINDiodeMethodValid(const Bool_t b ) @@ -269,28 +262,70 @@ // -------------------------------------------------------------------------- // -// Set the Excluded Bit from outside -// -void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b ) -{ - b ? SETBIT(fCalibFlags, kCombinedMethodValid) : CLRBIT(fCalibFlags, kCombinedMethodValid); -} - - +// Set the pedestals from outside +// +void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr) +{ + + fPed = ped; + fPedRms = pedrms; + fPedVar = pederr*pederr; +} + + +// -------------------------------------------------------------------------- +// +// Get the Error of the Mean pedestals: +// Returns square root of fPedVar +// +Float_t MCalibrationChargePix::GetPedErr() const +{ + return TMath::Sqrt(fPedVar); +} + +// -------------------------------------------------------------------------- +// +// Get the pedestals RMS: +// - Test bit kHiGainSaturation: +// If yes, return square root of fLoGainPedRmsSquare (if greater than 0, otherwise -1.), +// If no, return fPedRms +// Float_t MCalibrationChargePix::GetPedRms() const { - return IsHiGainSaturation() ? fLoGainPedRms : fPedRms; -} - + + if (IsHiGainSaturation()) + if (fLoGainPedRmsSquare < 0.) + return -1.; + else + return TMath::Sqrt(fLoGainPedRmsSquare); + + return fPedRms; +} + +// -------------------------------------------------------------------------- +// +// Get the Error of the pedestals RMS: +// - Test bit kHiGainSaturation: +// If yes, return kLoGainPedRms, else fPedRms +// Float_t MCalibrationChargePix::GetPedRmsErr() const { - return IsHiGainSaturation() ? TMath::Sqrt(fLoGainPedRmsVar) : TMath::Sqrt(fPedVar)/2.; -} - -Float_t MCalibrationChargePix::GetPedErr() const -{ - return TMath::Sqrt(fPedVar); -} - - + if (IsHiGainSaturation()) + if (fLoGainPedRmsSquareVar < 0.) + return -1.; + else + return TMath::Sqrt(0.25*fLoGainPedRmsSquareVar/fLoGainPedRmsSquare); + else + if (fPedVar < 0.) + return -1.; + else + return TMath::Sqrt(fPedVar)/2.; +} + + +// -------------------------------------------------------------------------- +// +// Get the Low Gain Mean: +// Returns fLoGainMean multiplied with fConversionHiLo +// Float_t MCalibrationChargePix::GetLoGainMean() const { @@ -298,4 +333,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the Error of the Low Gain Mean: +// Returns the quadratic sum of the relative low Gain Mean error and the +// the relative conversion High-to-Low error, mulitplied with GetLoGainMean() +// Float_t MCalibrationChargePix::GetLoGainMeanErr() const { @@ -304,10 +345,12 @@ /( fLoGainMean * fLoGainMean ); - const Float_t conversionRelVar = fConversionHiLoVar - /( fConversionHiLo * fConversionHiLo ); - - return TMath::Sqrt(chargeRelVar+conversionRelVar) * GetLoGainMean(); -} - + return TMath::Sqrt(chargeRelVar+GetConversionHiLoRelVar()) * GetLoGainMean(); +} + +// -------------------------------------------------------------------------- +// +// Get the Low Gain Sigma: +// Returns fLoGainSigma multiplied with fConversionHiLo +// Float_t MCalibrationChargePix::GetLoGainSigma() const { @@ -315,4 +358,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the Error of the Low Gain Sigma: +// Returns the quadratic sum of the relative low Gain Sigma error and the +// the relative conversion High-to-Low error, mulitplied with GetLoGainSigma() +// Float_t MCalibrationChargePix::GetLoGainSigmaErr() const { @@ -321,37 +370,75 @@ /( fLoGainSigma * fLoGainSigma ); - const Float_t conversionRelVar = fConversionHiLoVar - /( fConversionHiLo * fConversionHiLo ); - - return TMath::Sqrt(sigmaRelVar+conversionRelVar) * GetLoGainSigma(); -} - + return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetLoGainSigma(); +} + +// -------------------------------------------------------------------------- +// +// Get the reduced Sigma: +// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1. +// - Test bit kHiGainSaturation: +// If yes, return square root of fRSigmaSquare, multiplied with fConversionHiLo, +// If no , return square root of fRSigmaSquare +// Float_t MCalibrationChargePix::GetRSigma() const { - return IsHiGainSaturation() ? fRSigma*fConversionHiLo : fRSigma ; + if (fRSigmaSquare < 0) + return -1; + + const Float_t rsigma = TMath::Sqrt(fRSigmaSquare); + + return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ; } +// -------------------------------------------------------------------------- +// +// Get the error of the reduced Sigma: +// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Calculate the absolute variance of the reduced sigma with the formula: +// sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare +// - Test bit kHiGainSaturation: +// If yes, returns the quadratic sum of the relative reduced Sigma error and the +// the relative conversion High-to-Low error, mulitplied with GetRSigma() +// Else returns the square root of rel. (0.25*fRSigmaSquareVar / fRSigmaSquare) +// Float_t MCalibrationChargePix::GetRSigmaErr() const { - if (IsHiGainSaturation()) - { - const Float_t rsigmaRelVar = fRSigmaVar - /( fRSigma * fRSigma ); - const Float_t conversionRelVar = fConversionHiLoVar - /( fConversionHiLo * fConversionHiLo ); - return TMath::Sqrt(rsigmaRelVar+conversionRelVar) * GetRSigma(); - } + + if (fRSigmaSquareVar < 0) + return -1; + + // + // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) + // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) + // + const Float_t rsigmaVar = 0.25 * fRSigmaSquareVar / fRSigmaSquare; + + if (IsHiGainSaturation()) + return TMath::Sqrt(rsigmaVar/fRSigmaSquare + GetConversionHiLoRelVar()) * GetRSigma(); else - return TMath::Sqrt(fRSigmaVar); - -} - + return TMath::Sqrt(rsigmaVar); + +} + +// -------------------------------------------------------------------------- +// +// Get the conversion Error Hi-Gain to Low-Gain: +// - If fConversionHiLoVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fConversionHiLoVar +// Float_t MCalibrationChargePix::GetConversionHiLoErr() const { if (fConversionHiLoVar < 0.) return -1.; + return TMath::Sqrt(fConversionHiLoVar); } +// -------------------------------------------------------------------------- +// +// Get the error on the number of photo-electrons (F-Factor Method): +// - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fPheFFactorMethodVar +// Float_t MCalibrationChargePix::GetPheFFactorMethodErr() const { @@ -361,4 +448,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the mean conversion factor (Combined Method): +// - If fConversionCombinedMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fConversionCombinedMethodVar +// Float_t MCalibrationChargePix::GetConversionCombinedMethodErr() const { @@ -368,4 +461,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the mean conversion factor (PINDiode Method): +// - If fConversionPINDiodeMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fConversionPINDiodeMethodVar +// Float_t MCalibrationChargePix::GetConversionPINDiodeMethodErr() const { @@ -375,4 +474,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the mean conversion factor (BlindPixel  Method): +// - If fConversionBlindPixelMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fConversionBlindPixelMethodVar +// Float_t MCalibrationChargePix::GetConversionBlindPixelMethodErr() const { @@ -382,4 +487,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the mean conversion factor (FFactor  Method): +// - If fConversionFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fConversionFFactorMethodVar +// Float_t MCalibrationChargePix::GetConversionFFactorMethodErr() const { @@ -389,4 +500,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the total F-Factor (Combined  Method): +// - If fTotalFFactorCombinedMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fTotalFFactorCombinedMethodVar +// Float_t MCalibrationChargePix::GetTotalFFactorCombinedMethodErr() const { @@ -396,4 +513,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the total F-Factor (PIN Diode  Method): +// - If fTotalPINDiodeCombinedMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fTotalPINDiodeCombinedMethodVar +// Float_t MCalibrationChargePix::GetTotalFFactorPINDiodeMethodErr() const { @@ -403,4 +526,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the total F-Factor (Blind Pixel  Method): +// - If fTotalBlindPixelCombinedMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fTotalBlindPixelCombinedMethodVar +// Float_t MCalibrationChargePix::GetTotalFFactorBlindPixelMethodErr() const { @@ -410,4 +539,10 @@ } +// -------------------------------------------------------------------------- +// +// Get the error on the total F-Factor (F-Factor  Method): +// - If fTotalFFactorCombinedMethodVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Else returns the square root of fTotalFFactorCombinedMethodVar +// Float_t MCalibrationChargePix::GetTotalFFactorFFactorMethodErr() const { @@ -417,4 +552,8 @@ } +// -------------------------------------------------------------------------- +// +// Test bit kBlindPixelMethodValid +// Bool_t MCalibrationChargePix::IsBlindPixelMethodValid() const { @@ -422,4 +561,8 @@ } +// -------------------------------------------------------------------------- +// +// Test bit kFFactorMethodValid +// Bool_t MCalibrationChargePix::IsFFactorMethodValid() const { @@ -427,4 +570,8 @@ } +// -------------------------------------------------------------------------- +// +// Test bit kPINDiodeMethodValid +// Bool_t MCalibrationChargePix::IsPINDiodeMethodValid() const { @@ -432,4 +579,8 @@ } +// -------------------------------------------------------------------------- +// +// Test bit kCombinedMethodValid +// Bool_t MCalibrationChargePix::IsCombinedMethodValid() const { @@ -437,34 +588,46 @@ } -// +// ---------------------------------------------------------------------------- // +// - If fSigma is smaller than 0 (i.e. has not yet been set), return kFALSE +// - If fPedRms is smaller than 0 (i.e. has not yet been set), return kFALSE +// +// Calculate the reduced sigma: +// - Test bit IsHiGainSaturation() for the Sigma: +// If yes, take fLoGainSigma and fLoGainSigmaVar +// If no , take fHiGainSigma and fHiGainSigmaVar +// +// - Test bit IsHiGainSaturation() for the pedRMS: +// If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar +// If no , take fPedRms and fPedRmsVar +// +// - Calculate the reduced sigma with the formula: +// fRSigmaSquare = Sigma*Sigma - pedRMS*pedRMS +// +// - Calculate the variance of the reduced sigma with the formula: +// fRSigmaSquareVar = 4.* (sigmaVar*Sigma*Sigma + pedRmsVar*pedRMS*pedRMS) // Bool_t MCalibrationChargePix::CalcReducedSigma() { + + if (GetSigma() < 0.) + return kFALSE; + + if (GetPedRms() < 0.) + return kFALSE; const Float_t sigma = IsHiGainSaturation() ? fLoGainSigma : fHiGainSigma ; const Float_t sigmavar = IsHiGainSaturation() ? fLoGainSigmaVar : fHiGainSigmaVar; - - const Float_t sigmaSquare = sigma * sigma; - const Float_t sigmaSquareVar = 4.* sigmavar * sigmaSquare; - - Float_t pedRmsSquare ; - Float_t pedRmsSquareVar; - - if (IsHiGainSaturation()) - { - pedRmsSquare = fLoGainPedRms * fLoGainPedRms; - pedRmsSquareVar = 4.* fLoGainPedRmsVar * pedRmsSquare; - } - else - { - pedRmsSquare = fPedRms * fPedRms; - pedRmsSquareVar = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2. - } + const Float_t pedRmsSquare = IsHiGainSaturation() ? fLoGainPedRmsSquare : fPedRms*fPedRms; + const Float_t pedRmsSquareVar = IsHiGainSaturation() ? fLoGainPedRmsSquareVar : 0.25*fPedVar*pedRmsSquare; + + const Float_t sigmaSquare = sigma * sigma; + const Float_t sigmaSquareVar = sigmavar * sigmaSquare; + // // Calculate the reduced sigmas // - const Float_t rsigmasquare = sigmaSquare - pedRmsSquare; - if (rsigmasquare <= 0.) + fRSigmaSquare = sigmaSquare - pedRmsSquare; + if (fRSigmaSquare <= 0.) { *fLog << warn @@ -473,20 +636,40 @@ return kFALSE; } - - - fRSigma = TMath::Sqrt(rsigmasquare); - fRSigmaVar = 0.25 * (sigmaSquareVar + pedRmsSquareVar) / rsigmasquare; + + fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar); return kTRUE; } -// -// Calculate the number of photo-electrons after the F-Factor method -// Calculate the errors of the F-Factor method +// ------------------------------------------------------------------ +// +// If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), +// set kFFactorMethodValid to kFALSE and return kFALSE +// +// Calculate the number of photo-electrons with the F-Factor method: +// - Test bit IsHiGainSaturation() for the Mean Sum of FADC slices: +// If yes, take fLoGainMean and fLoGainMeanVar +// If no , take fHiGainMean and fHiGainMeanVar +// +// - Test bit IsHiGainSaturation() for the pedRMS: +// If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar +// If no , take fPedRms and fPedRmsVar +// +// - Calculate the number of photo-electrons with the formula: +// fPheFFactorMethod = gkFFactor*gkFFactor * mean * mean / fRSigmaSquare +// +// - Calculate the Variance on the photo-electrons with the formula: +// fPheFFactorMethodVar = ( 4. * gkFFactorErr * gkFFactorErr / ( gkFFactor * gkFFactor ) +// + 4. * fMeanVar / ( mean * mean ) +// + 4. * fRSigmaVar / ( fRSigma * fRSigma ) +// ) * fPheFFactor * fPheFFactor +// - if fPheFFactorMethod is less than fPheFFactorMethodLimit, +// Set kFFactorMethodValid to kFALSE and return kFALSE +// else: Set kFFactorMethodValid to kTRUE and return kTRUE // Bool_t MCalibrationChargePix::CalcFFactorMethod() { - if (fRSigma < 0.) + if (fRSigmaSquare < 0.) { SetFFactorMethodValid(kFALSE); @@ -500,14 +683,11 @@ // Square all variables in order to avoid applications of square root // - // First the relative error squares - // - const Float_t meanSquare = mean * mean; - const Float_t meanSquareRelVar = 4.* var/ meanSquare; - - const Float_t ffactorsquare = gkFFactor * gkFFactor; - const Float_t ffactorsquareRelVar = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare; - - const Float_t rsigmaSquare = fRSigma * fRSigma; - const Float_t rsigmaSquareRelVar = 4.* fRSigmaVar / rsigmaSquare; + const Float_t meanSquare = mean * mean; + const Float_t meanSquareRelVar = 4.* var / meanSquare; + + const Float_t ffactorsquare = gkFFactor * gkFFactor; + const Float_t ffactorsquareRelVar = 4.* gkFFactorErr * gkFFactorErr / ffactorsquare; + + const Float_t rsigmaSquareRelVar = fRSigmaSquareVar / fRSigmaSquare; // @@ -515,5 +695,5 @@ // (independent on Hi Gain or Lo Gain) // - fPheFFactorMethod = ffactorsquare * meanSquare / rsigmaSquare; + fPheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare; if (fPheFFactorMethod < fPheFFactorMethodLimit) @@ -534,6 +714,36 @@ +// ---------------------------------------------------------------------------- +// +// - If fPed is smaller than 0 (i.e. has not yet been set), return. +// - If fPedVar is smaller than 0 (i.e. has not yet been set), return. +// +// Calculate the electronic pedestal RMS with the formula: +// - elec. pedestal = gkElectronicPedRms * sqrt(logainsamples) +// +// Calculate the LONS ped. RMS contribution in the high-gain +// from the high gain Pedestal RMS with the formula: +// - HiGain NSB square = fPedRms * fPedRms - elec.ped.* elec.ped. +// - Var(HiGain NSB square) = fPedVar * fPedRms * fPedRms + 4.*elecPedRmsVar * elec.ped.* elec.ped. +// +// If PedRMS (LONS,lowgain) square is smaller than 0., set it to zero. (but not the error!) +// +// Convert the LONS ped. RMS contribution to the low-gain with the formula: +// - LoGain NSB square = - HiGain NSB square / (fConversionHiLo*fConversionHiLo) +// - Var(LoGain NSB square) = ( Var(HiGain NSB square) / (HiGain NSB square * HiGain NSB square) +// + GetConversionHiLoRelVar() ) +// * LoGain NSB square * LoGain NSB square +// +// - Low Gain Ped RMS Square = LoGain NSB square + elec.ped. square +// Var (Low Gain Ped RMS Square) = Var(LoGain NSB square) + Var(elec.ped. square) +// void MCalibrationChargePix::CalcLoGainPedestal(Float_t logainsamples) { + + if (fPedRms < 0.) + return; + + if (fPedVar < 0.) + return; const Float_t elecPedRms = gkElectronicPedRms * TMath::Sqrt(logainsamples); @@ -552,10 +762,10 @@ const Float_t elecRmsSquareVar = 4.*elecPedRmsVar * elecRmsSquare; - Float_t nsbSquare = pedRmsSquare - elecRmsSquare; - Float_t nsbSquareRelVar = (pedRmsSquareVar + elecRmsSquareVar) - / (nsbSquare * nsbSquare) ; - - if (nsbSquare < 0.) - nsbSquare = 0.; + Float_t higainNsbSquare = pedRmsSquare - elecRmsSquare; + Float_t higainNsbSquareRelVar = (pedRmsSquareVar + elecRmsSquareVar) + / (higainNsbSquare * higainNsbSquare) ; + + if (higainNsbSquare < 0.) + higainNsbSquare = 0.; // @@ -563,17 +773,21 @@ // add it quadratically to the electronic noise // - const Float_t conversionSquare = fConversionHiLo * fConversionHiLo; - const Float_t convertedNsbSquare = nsbSquare / conversionSquare; - const Float_t convertedNsbSquareVar = nsbSquareRelVar - * convertedNsbSquare * convertedNsbSquare; + const Float_t conversionSquare = fConversionHiLo * fConversionHiLo; + const Float_t conversionSquareRelVar = 4.* GetConversionHiLoRelVar(); + + const Float_t logainNsbSquare = higainNsbSquare / conversionSquare; + const Float_t logainNsbSquareVar = ( higainNsbSquareRelVar + conversionSquareRelVar ) + * logainNsbSquare * logainNsbSquare; - pedRmsSquare = convertedNsbSquare + elecRmsSquare; - pedRmsSquareVar = convertedNsbSquareVar + elecRmsSquareVar; - - fLoGainPedRms = TMath::Sqrt(pedRmsSquare); - fLoGainPedRmsVar = 0.25 * pedRmsSquareVar / pedRmsSquare; - + fLoGainPedRmsSquare = logainNsbSquare + elecRmsSquare; + fLoGainPedRmsSquareVar = logainNsbSquareVar + elecRmsSquareVar; } +const Float_t MCalibrationChargePix::GetConversionHiLoRelVar() const +{ + if (fConversionHiLo == 0.) + return 0.; + + return fConversionHiLoVar / (fConversionHiLo * fConversionHiLo); +} - Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationPedCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationPedCam.cc (revision 3650) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationPedCam.cc (revision 3650) @@ -0,0 +1,139 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MCalibrationPedCam +// +// Hold the pedestal Calibration results obtained from MHPedestalCam of the camera: +// +// 1) MCalibrationPedCam initializes a TClonesArray whose elements are +// pointers to MCalibrationPix Containers +// +// The calculated values (types of GetPixelContent) are: +// +// Fitted values: +// ============== +// +// 0: Fitted Pedestal +// 1: Error of fitted Pedestal +// 2: Sigma of fitted Pedestal +// 3: Error of Sigma of fitted Pedestal +// +// Useful variables derived from the fit results: +// ============================================= +// +// 4: Returned probability of Gauss fit to Pedestal distribution +// +///////////////////////////////////////////////////////////////////////////// +#include "MCalibrationPedCam.h" +#include "MCalibrationCam.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MCalibrationPix.h" + +ClassImp(MCalibrationPedCam); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry +// Later, a call to MCalibrationPedCam::InitSize(Int_t size) has to be performed +// +MCalibrationPedCam::MCalibrationPedCam(const char *name, const char *title) +{ + fName = name ? name : "MCalibrationPedCam"; + fTitle = title ? title : "Storage container for the Pedestal Calibration Results in the camera"; + + fPixels = new TClonesArray("MCalibrationPix",1); + fAverageAreas = new TClonesArray("MCalibrationPix",1); + fAverageSectors = new TClonesArray("MCalibrationPix",1); + +} + +// -------------------------------------------------------------------------- +// +// The calculated values (types of GetPixelContent) are: +// +// Fitted values: +// ============== +// +// 0: Fitted Pedestal +// 1: Error of fitted Pedestal +// 2: Sigma of fitted Pedestal +// 3: Error of Sigma of fitted Pedestal +// +// Useful variables derived from the fit results: +// ============================================= +// +// 4: Returned probability of Gauss fit to Pedestal distribution +// +Bool_t MCalibrationPedCam::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; + + if ((*this)[idx].IsExcluded()) + return kFALSE; + + switch (type) + { + case 0: + val = (*this)[idx].GetMean(); + break; + case 1: + val = (*this)[idx].GetMeanErr(); + break; + case 2: + val = (*this)[idx].GetSigma(); + break; + case 3: + val = (*this)[idx].GetSigmaErr(); + break; + case 4: + val = (*this)[idx].GetProb(); + break; + default: + return kFALSE; + } + + return val!=-1.; + +} Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc (revision 3649) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc (revision 3650) @@ -703,4 +703,5 @@ hist.Renorm(); + // // 4) Check for oscillations @@ -708,4 +709,5 @@ hist.CreateFourierSpectrum(); + if (!hist.IsFourierSpectrumOK()) bad.SetUncalibrated( osctyp ); Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc (revision 3649) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc (revision 3650) @@ -216,7 +216,24 @@ { - fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationChargeCam"); + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + fCam = (MCalibrationCam*)pList->FindObject("MCalibrationChargeCam"); if (!fCam) - return kFALSE; + { + fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam")); + if (!fCam) + { + gLog << err << "Cannot find nor create MCalibrationChargeCam ... abort." << endl; + return kFALSE; + } + else + { + fCam->InitSize(npixels); + fCam->InitAverageAreas(nareas); + fCam->InitAverageSectors(nsectors); + } + } MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); @@ -226,8 +243,4 @@ return kFALSE; } - - const Int_t npixels = fGeom->GetNumPixels(); - const Int_t nsectors = fGeom->GetNumSectors(); - const Int_t nareas = fGeom->GetNumAreas(); if (fHiGainArray->GetEntries()==0) Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 3649) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 3650) @@ -153,8 +153,21 @@ { - fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationRelTimeCam"); + fCam = (MCalibrationCam*)pList->FindObject("MCalibrationRelTimeCam"); if (!fCam) - return kFALSE; - + { + fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationRelTimeCam")); + if (!fCam) + { + gLog << err << "Cannot find nor create MCalibrationRelTimeCam ... abort." << endl; + return kFALSE; + } + else + { + fCam->InitSize(npixels); + fCam->InitAverageAreas(nareas); + fCam->InitAverageSectors(nsectors); + } + } + MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam"); if (!signal) @@ -376,9 +389,9 @@ { - FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, + FitHiGainArrays((*fCam),*fBadPixels, MBadPixelsPix::kRelTimeNotFitted, MBadPixelsPix::kRelTimeOscillating); - FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, + FitLoGainArrays((*fCam),*fBadPixels, MBadPixelsPix::kRelTimeNotFitted, MBadPixelsPix::kRelTimeOscillating);