Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3671) +++ trunk/MagicSoft/Mars/Changelog (revision 3672) @@ -154,9 +154,13 @@ - updated and enlarged documentation - * mcalib/MCalibrationPix.[h,cc] * mcalib/MHCalibrationChargePINDiode.cc * mcalib/MHCalibrationChargeCam.[h,cc] - rename Setter for fMeanVar from SetMeanErr() to SetMeanVar() - rename arrays of fAverageSigmaErr to AverageSigmaVar + + * mcalib/MCalibrationPix.[h,cc] + * mcalib/MCalibrationChargePix.[h,cc] + * mcalib/MCalibrationChargeBlindPix.[h,cc] + * mcalib/MCalibrationChargePINDiode.[h,cc] - add Getters to the relative variance Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.cc (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.cc (revision 3672) @@ -197,4 +197,90 @@ // -------------------------------------------------------------------------- // +// Return -1 if fLambdaVar is smaller than 0. +// Return -1 if fLambda is 0. +// Return fLambdaVar / (fLambda * fLambda ) +// +Float_t MCalibrationChargeBlindPix::GetLambdaRelVar() const +{ + if (fLambdaVar < 0.) + return -1.; + + if (fLambda == 0.) + return -1.; + + return fLambdaVar / fLambda / fLambda ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkBlindPixelQEGreenErr is smaller than 0. +// Return -1 if gkBlindPixelQEGreen is 0. +// Return gkBlindPixelQEGreenErr^2 / (gkBlindPixelQEGreen^2 ) +// +const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEGreenRelVar() const +{ + if (gkBlindPixelQEGreenErr < 0.) + return -1.; + + if (gkBlindPixelQEGreen == 0.) + return -1.; + + return gkBlindPixelQEGreenErr * gkBlindPixelQEGreenErr / gkBlindPixelQEGreen / gkBlindPixelQEGreen ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkBlindPixelQEBlueErr is smaller than 0. +// Return -1 if gkBlindPixelQEBlue is 0. +// Return gkBlindPixelQEBlueErr^2 / gkBlindPixelQEBlue^2 +// +const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEBlueRelVar() const +{ + if (gkBlindPixelQEBlueErr < 0.) + return -1.; + + if (gkBlindPixelQEBlue == 0.) + return -1.; + + return gkBlindPixelQEBlueErr * gkBlindPixelQEBlueErr / gkBlindPixelQEBlue / gkBlindPixelQEBlue ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkBlindPixelQEUVErr is smaller than 0. +// Return -1 if gkBlindPixelQEUV is 0. +// Return gkBlindPixelQEUVErr ^2 / gkBlindPixelQEUV^2 +// +const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEUVRelVar() const +{ + if (gkBlindPixelQEUVErr < 0.) + return -1.; + + if (gkBlindPixelQEUV == 0.) + return -1.; + + return gkBlindPixelQEUVErr * gkBlindPixelQEUVErr / gkBlindPixelQEUV / gkBlindPixelQEUV ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkBlindPixelQECT1Err is smaller than 0. +// Return -1 if gkBlindPixelQECT1 is 0. +// Return gkBlindPixelQECT1Err ^2 / gkBlindPixelQECT1^2 +// +const Float_t MCalibrationChargeBlindPix::GetBlindPixelQECT1RelVar() const +{ + if (gkBlindPixelQECT1Err < 0.) + return -1.; + + if (gkBlindPixelQECT1 == 0.) + return -1.; + + return gkBlindPixelQECT1Err * gkBlindPixelQECT1Err / gkBlindPixelQECT1 / gkBlindPixelQECT1 ; +} + + +// -------------------------------------------------------------------------- +// // Test bit kChargeFitValid // @@ -267,60 +353,28 @@ // The blind pixel has exactly 100 mm^2 area (with negligible error), // - fFluxInsidePlexiglass = fLambda*gkBlindPixelArea; - - // - // Start calculation of number of photons relative Variance - // - const Float_t lambdaRelVar = fLambdaVar / ( fLambda * fLambda ); - fFluxInsidePlexiglassVar = lambdaRelVar; - switch (fColor) { case kGREEN: - fFluxInsidePlexiglass /= gkBlindPixelQEGreen; - fFluxInsidePlexiglassVar += gkBlindPixelQEGreenErr * gkBlindPixelQEGreenErr - / ( gkBlindPixelQEGreen * gkBlindPixelQEGreen ); - - fFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttGreen); // correct for absorption + fFluxInsidePlexiglass = fLambda * gkBlindPixelQEGreen * TMath::Power(10,gkBlindPixelAttGreen); // attenuation has negligible error + fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEGreenRelVar(); break; case kBLUE: - fFluxInsidePlexiglass /= gkBlindPixelQEBlue; - fFluxInsidePlexiglassVar += gkBlindPixelQEBlueErr * gkBlindPixelQEBlueErr - / ( gkBlindPixelQEBlue * gkBlindPixelQEBlue ); - - fFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttBlue); // correct for absorption + fFluxInsidePlexiglass = fLambda * gkBlindPixelQEBlue * TMath::Power(10,gkBlindPixelAttBlue); // attenuation has negligible error + fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEBlueRelVar(); break; case kUV: - fFluxInsidePlexiglass /= gkBlindPixelQEUV; - fFluxInsidePlexiglassVar += gkBlindPixelQEUVErr* gkBlindPixelQEUVErr - / ( gkBlindPixelQEUV * gkBlindPixelQEUV ); - - fFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttUV); // correct for absorption + fFluxInsidePlexiglass = fLambda * gkBlindPixelQEUV * TMath::Power(10,gkBlindPixelAttUV); // attenuation has negligible error + fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEUVRelVar(); break; case kCT1: default: - fFluxInsidePlexiglass /= gkBlindPixelQECT1; - fFluxInsidePlexiglassVar += gkBlindPixelQECT1Err * gkBlindPixelQECT1Err - / ( gkBlindPixelQECT1 * gkBlindPixelQECT1 ); - - fFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttCT1); // correct for absorption + fFluxInsidePlexiglass = fLambda * gkBlindPixelQECT1 * TMath::Power(10,gkBlindPixelAttCT1); // attenuation has negligible error + fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQECT1RelVar(); break; } - - *fLog << inf << endl; - *fLog << inf << " Photon flux [ph/mm^2] inside Plexiglass: " - << fFluxInsidePlexiglass << endl; - - if (fFluxInsidePlexiglass < 0.) - return kFALSE; - - if (fFluxInsidePlexiglassVar < 0.) - return kFALSE; - - SetFluxInsidePlexiglassAvailable(kTRUE); // @@ -329,7 +383,15 @@ fFluxInsidePlexiglassVar *= fFluxInsidePlexiglass * fFluxInsidePlexiglass; - *fLog << inf << " Error on photon flux [ph/mm^2] inside Plexiglass: " - << GetFluxInsidePlexiglass() << endl; + if (fFluxInsidePlexiglass < 0.) + return kFALSE; + + if (fFluxInsidePlexiglassVar < 0.) + return kFALSE; + + SetFluxInsidePlexiglassAvailable(kTRUE); + *fLog << inf << endl; + *fLog << inf << " Photon flux [ph/mm^2] inside Plexiglass: " + << Form("%5.3f%s%5.3f",fFluxInsidePlexiglass," +- ",GetFluxInsidePlexiglassErr()) << endl; return kTRUE; Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.h (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.h (revision 3672) @@ -44,5 +44,10 @@ PulserColor_t fColor; -public: + const Float_t GetBlindPixelQEGreenRelVar() const; + const Float_t GetBlindPixelQEBlueRelVar () const; + const Float_t GetBlindPixelQEUVRelVar () const; + const Float_t GetBlindPixelQECT1RelVar () const; + + public: MCalibrationChargeBlindPix(const char *name=NULL, const char *title=NULL); @@ -52,17 +57,17 @@ // Setters - void SetColor ( const PulserColor_t color ) { fColor = color; } - void SetLambda ( const Float_t f ) { fLambda = f; } - void SetLambdaErr ( const Float_t f ) { fLambdaVar = f*f; } - void SetLambdaCheck ( const Float_t f ) { fLambdaCheck = f; } - void SetLambdaCheckErr ( const Float_t f ) { fLambdaCheckErr = f; } - void SetMu0 ( const Float_t f ) { fMu0 = f; } - void SetMu0Err ( const Float_t f ) { fMu0Err = f; } - void SetMu1 ( const Float_t f ) { fMu1 = f; } - void SetMu1Err ( const Float_t f ) { fMu1Err = f; } - void SetSigma0 ( const Float_t f ) { fSigma0 = f; } - void SetSigma0Err ( const Float_t f ) { fSigma0Err = f; } - void SetSigma1 ( const Float_t f ) { fSigma1 = f; } - void SetSigma1Err ( const Float_t f ) { fSigma1Err = f; } + void SetColor ( const PulserColor_t color ) { fColor = color; } + void SetLambda ( const Float_t f ) { fLambda = f; } + void SetLambdaVar ( const Float_t f ) { fLambdaVar = f; } + void SetLambdaCheck ( const Float_t f ) { fLambdaCheck = f; } + void SetLambdaCheckErr ( const Float_t f ) { fLambdaCheckErr = f; } + void SetMu0 ( const Float_t f ) { fMu0 = f; } + void SetMu0Err ( const Float_t f ) { fMu0Err = f; } + void SetMu1 ( const Float_t f ) { fMu1 = f; } + void SetMu1Err ( const Float_t f ) { fMu1Err = f; } + void SetSigma0 ( const Float_t f ) { fSigma0 = f; } + void SetSigma0Err ( const Float_t f ) { fSigma0Err = f; } + void SetSigma1 ( const Float_t f ) { fSigma1 = f; } + void SetSigma1Err ( const Float_t f ) { fSigma1Err = f; } void SetOscillating ( const Bool_t b=kTRUE); @@ -75,4 +80,5 @@ Float_t GetLambda() const { return fLambda; } Float_t GetLambdaErr() const; + Float_t GetLambdaRelVar() const; Float_t GetLambdaCheck() const { return fLambdaCheck; } Float_t GetLambdaCheckErr() const { return fLambdaCheckErr; } Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.cc (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.cc (revision 3672) @@ -193,4 +193,120 @@ } +// -------------------------------------------------------------------------- +// +// Return -1 if fNumPhotonsVar is smaller than 0. +// Return square root of fNumPhotonsVar +// +Float_t MCalibrationChargePINDiode::GetNumPhotonsErr() const +{ + if (fNumPhotonsVar < 0.) + return -1.; + + return TMath::Sqrt(fNumPhotonsVar); +} + +// -------------------------------------------------------------------------- +// +// Return -1 if fNumPhotonsVar is smaller than 0. +// Return -1 if fNumPhotons is 0. +// Return fNumPhotonsVar / (fNumPhotons^2 ) +// +Float_t MCalibrationChargePINDiode::GetNumPhotonsRelVar() const +{ + if (fNumPhotonsVar < 0.) + return -1.; + + if (fNumPhotons == 0.) + return -1.; + + return fNumPhotonsVar / fNumPhotons / fNumPhotons ; +} + + +// -------------------------------------------------------------------------- +// +// Return -1 if gkPINDiodeQEGreenErr is smaller than 0. +// Return -1 if gkPINDiodeQEGreen is 0. +// Return gkPINDiodeQEGreenErr^2 / (gkPINDiodeQEGreen^2 ) +// +const Float_t MCalibrationChargePINDiode::GetPINDiodeQEGreenRelVar() const +{ + if (gkPINDiodeQEGreenErr < 0.) + return -1.; + + if (gkPINDiodeQEGreen == 0.) + return -1.; + + return gkPINDiodeQEGreenErr * gkPINDiodeQEGreenErr / gkPINDiodeQEGreen / gkPINDiodeQEGreen ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkPINDiodeQEBlueErr is smaller than 0. +// Return -1 if gkPINDiodeQEBlue is 0. +// Return gkPINDiodeQEBlueErr^2 / gkPINDiodeQEBlue^2 +// +const Float_t MCalibrationChargePINDiode::GetPINDiodeQEBlueRelVar() const +{ + if (gkPINDiodeQEBlueErr < 0.) + return -1.; + + if (gkPINDiodeQEBlue == 0.) + return -1.; + + return gkPINDiodeQEBlueErr * gkPINDiodeQEBlueErr / gkPINDiodeQEBlue / gkPINDiodeQEBlue ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkPINDiodeQEUVErr is smaller than 0. +// Return -1 if gkPINDiodeQEUV is 0. +// Return gkPINDiodeQEUVErr ^2 / gkPINDiodeQEUV^2 +// +const Float_t MCalibrationChargePINDiode::GetPINDiodeQEUVRelVar() const +{ + if (gkPINDiodeQEUVErr < 0.) + return -1.; + + if (gkPINDiodeQEUV == 0.) + return -1.; + + return gkPINDiodeQEUVErr * gkPINDiodeQEUVErr / gkPINDiodeQEUV / gkPINDiodeQEUV ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkPINDiodeQECT1Err is smaller than 0. +// Return -1 if gkPINDiodeQECT1 is 0. +// Return gkPINDiodeQECT1Err ^2 / gkPINDiodeQECT1^2 +// +const Float_t MCalibrationChargePINDiode::GetPINDiodeQECT1RelVar() const +{ + if (gkPINDiodeQECT1Err < 0.) + return -1.; + + if (gkPINDiodeQECT1 == 0.) + return -1.; + + return gkPINDiodeQECT1Err * gkPINDiodeQECT1Err / gkPINDiodeQECT1 / gkPINDiodeQECT1 ; +} + +// -------------------------------------------------------------------------- +// +// Return -1 if gkSolidAngleRatioErr is smaller than 0. +// Return -1 if gkSolidAngleRatio is 0. +// Return gkSolidAngleRatioErr ^2 / gkSolidAngleRatio^2 +// +const Float_t MCalibrationChargePINDiode::GetSolidAngleRatioRelVar() const +{ + if (gkSolidAngleRatioErr < 0.) + return -1.; + + if (gkSolidAngleRatio == 0.) + return -1.; + + return gkSolidAngleRatioErr * gkSolidAngleRatioErr / gkSolidAngleRatio / gkSolidAngleRatio ; +} + // -------------------------------------------------------------------------- @@ -257,60 +373,41 @@ // Correct for the distance between camera and PIN Diode and for different areas. // - fFluxOutsidePlexiglass = fNumPhotons * gkSolidAngleRatio; - - // - // Start calculation of number of photons relative Variance (!!) - // - const Float_t numphotRelVar = fNumPhotonsVar / ( fNumPhotons * fNumPhotons ) ; - const Float_t solidangleRelVar = gkSolidAngleRatioErr* gkSolidAngleRatioErr - / ( gkSolidAngleRatio * gkSolidAngleRatio ); - Float_t fluxRelVar = numphotRelVar + solidangleRelVar ; - switch (fColor) { case kGREEN: - fFluxOutsidePlexiglass /= gkPINDiodeQEGreen; - fluxRelVar += gkPINDiodeQEGreenErr* gkPINDiodeQEGreenErr - / ( gkPINDiodeQEGreen * gkPINDiodeQEGreen ); + fFluxOutsidePlexiglass = fNumPhotons * gkSolidAngleRatio * gkPINDiodeQEGreen; + fFluxOutsidePlexiglassVar = GetNumPhotonsRelVar() + GetSolidAngleRatioRelVar() + GetPINDiodeQEGreenRelVar(); break; case kBLUE: - fFluxOutsidePlexiglass /= gkPINDiodeQEBlue; - fluxRelVar += gkPINDiodeQEBlueErr* gkPINDiodeQEBlueErr - / ( gkPINDiodeQEBlue * gkPINDiodeQEBlue ); + fFluxOutsidePlexiglass = fNumPhotons * gkSolidAngleRatio * gkPINDiodeQEBlue; + fFluxOutsidePlexiglassVar = GetNumPhotonsRelVar() + GetSolidAngleRatioRelVar() + GetPINDiodeQEBlueRelVar(); break; case kUV: - fFluxOutsidePlexiglass /= gkPINDiodeQEUV; - fluxRelVar += gkPINDiodeQEUVErr* gkPINDiodeQEUVErr - / ( gkPINDiodeQEUV * gkPINDiodeQEUV ); + fFluxOutsidePlexiglass = fNumPhotons * gkSolidAngleRatio * gkPINDiodeQEUV; + fFluxOutsidePlexiglassVar = GetNumPhotonsRelVar() + GetSolidAngleRatioRelVar() + GetPINDiodeQEUVRelVar(); break; case kCT1: default: - fFluxOutsidePlexiglass /= gkPINDiodeQECT1; - fluxRelVar += gkPINDiodeQECT1Err* gkPINDiodeQECT1Err - / ( gkPINDiodeQECT1 * gkPINDiodeQECT1 ); + fFluxOutsidePlexiglass = fNumPhotons * gkSolidAngleRatio * gkPINDiodeQECT1; + fFluxOutsidePlexiglassVar = GetNumPhotonsRelVar() + GetSolidAngleRatioRelVar() + GetPINDiodeQECT1RelVar(); break; } - - - *fLog << inf << endl; - *fLog << inf << " Mean Photon flux [ph/mm^2] outside Plexiglass: " - << fFluxOutsidePlexiglass << endl; - - if (fFluxOutsidePlexiglass < 0.) - return kFALSE; - - if (fluxRelVar < 0.) - return kFALSE; - - SetFluxOutsidePlexiglassAvailable(); // // Finish calculation of errors -> convert from relative variance to absolute variance // - fFluxOutsidePlexiglassVar *= fluxRelVar * fFluxOutsidePlexiglass * fFluxOutsidePlexiglass; - - *fLog << inf << " Error on Photon flux [ph/mm^2] outside Plexiglass: " - << GetFluxOutsidePlexiglassErr() << endl; + fFluxOutsidePlexiglassVar *= fFluxOutsidePlexiglass * fFluxOutsidePlexiglass; + + if (fFluxOutsidePlexiglass < 0.) + return kFALSE; + + if (fFluxOutsidePlexiglassVar < 0.) + return kFALSE; + + SetFluxOutsidePlexiglassAvailable(); + *fLog << inf << endl; + *fLog << inf << " Mean Photon flux [ph/mm^2] outside Plexiglass: " + << Form("%5.3f%s%5.3f",fFluxOutsidePlexiglass," +- ",GetFluxOutsidePlexiglassErr()) << endl; return kTRUE; Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.h (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePINDiode.h (revision 3672) @@ -47,5 +47,11 @@ kFluxOutsidePlexiglassAvailable }; // Possible bits to be set -public: + const Float_t GetPINDiodeQEGreenRelVar() const; + const Float_t GetPINDiodeQEBlueRelVar () const; + const Float_t GetPINDiodeQEUVRelVar () const; + const Float_t GetPINDiodeQECT1RelVar () const; + const Float_t GetSolidAngleRatioRelVar() const; + + public: MCalibrationChargePINDiode(const char *name=NULL, const char *title=NULL); @@ -75,12 +81,15 @@ // Pedestals - Float_t GetAbsTimeMean() const { return fAbsTimeMean; } - Float_t GetAbsTimeRms () const { return fAbsTimeRms; } - Float_t GetPed() const { return fPed; } - Float_t GetPedRms() const { return fPedRms; } + Float_t GetAbsTimeMean () const { return fAbsTimeMean; } + Float_t GetAbsTimeRms () const { return fAbsTimeRms; } + Float_t GetNumPhotons () const { return fNumPhotons; } + Float_t GetNumPhotonsErr() const; + Float_t GetNumPhotonsRelVar() const; + Float_t GetPed() const { return fPed; } + Float_t GetPedRms() const { return fPedRms; } - Bool_t IsChargeFitValid() const; - Bool_t IsTimeFitValid() const; - Bool_t IsOscillating() const; + Bool_t IsChargeFitValid() const; + Bool_t IsTimeFitValid() const; + Bool_t IsOscillating() const; Bool_t IsFluxOutsidePlexiglassAvailable() const; Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc (revision 3672) @@ -274,7 +274,58 @@ fPedVar = pederr*pederr; } - - -// -------------------------------------------------------------------------- + +// ------------------------------------------------------------------------------- +// +// Get the conversion Error Hi-Gain to Low-Gain: +// - If fConversionHiLoVar is smaller than 0 (i.e. has not yet been set), return -1. +// +Float_t MCalibrationChargePix::GetConversionHiLoErr() const +{ + if (fConversionHiLoVar < 0.) + return -1.; + + return TMath::Sqrt(fConversionHiLoVar); +} + +// -------------------------------------------------------------------------- +// +// Get the relative variance of the conversion factor between higain and logain: +// - If fConversionHiLo is 0, return -1. +// - If fConversionHiLoVar is smaller than 0, return -1. +// - Else returns: fConversionHiLoVar / fConversionHiLo^2 +// +const Float_t MCalibrationChargePix::GetConversionHiLoRelVar() const +{ + + if (fConversionHiLoVar < 0.) + return -1.; + + if (fConversionHiLo == 0.) + return -1.; + + return fConversionHiLoVar / (fConversionHiLo * fConversionHiLo); +} + + +// -------------------------------------------------------------------------- +// +// Get the relative variance of the conversion factor between higain and logain: +// - If gkFFactor is 0, return -1. +// - If gkFFactorErr is smaller than 0, return -1. +// - Else returns: gkFFactorErr^2 / gkFFactor*^2 +// +const Float_t MCalibrationChargePix::GetFFactorRelVar() const +{ + + if (gkFFactorErr < 0.) + return -1.; + + if (gkFFactor == 0.) + return -1.; + + return gkFFactorErr * gkFFactorErr / (gkFFactor * gkFFactor); +} + + // // Get the Error of the Mean pedestals: @@ -348,15 +399,11 @@ Float_t MCalibrationChargePix::GetConvertedLoGainMeanErr() const { - - if (fLoGainMeanVar < 0.) - return -1.; - - if (fLoGainMean < 0.) - return -1.; - - const Float_t chargeRelVar = fLoGainMeanVar - /( fLoGainMean * fLoGainMean ); - - return TMath::Sqrt(chargeRelVar+GetConversionHiLoRelVar()) * GetConvertedLoGainMean(); + + const Float_t logainrelvar = GetLoGainMeanRelVar(); + + if (logainrelvar < 0.) + return -1.; + + return TMath::Sqrt(logainrelvar + GetConversionHiLoRelVar()) * GetConvertedLoGainMean(); } @@ -445,14 +492,28 @@ // -------------------------------------------------------------------------- // -// 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 relative variance of the reduced Sigma: +// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1. +// - Calculate the relative variance of the reduced sigma squares with the formula: +// reduced sigma rel. variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare +// - Test bit kHiGainSaturation: +// If yes, returns the sum of the relative variances of the reduced sigma and fConversionHiLo +// Else returns the relative variance of the reduced sigma +// +Float_t MCalibrationChargePix::GetRSigmaRelVar() const +{ + + if (fRSigmaSquareVar < 0) + return -1; + + // + // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) + // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) + // + const Float_t rsigmaRelVar = 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare ); + + if (IsHiGainSaturation()) + return rsigmaRelVar + GetConversionHiLoRelVar(); + else + return rsigmaRelVar; } @@ -705,15 +766,12 @@ } - const Float_t mean = IsHiGainSaturation() ? fLoGainMean : fHiGainMean ; - const Float_t var = IsHiGainSaturation() ? fLoGainMeanVar : fHiGainMeanVar; - // // Square all variables in order to avoid applications of square root // - const Float_t meanSquare = mean * mean; - const Float_t meanSquareRelVar = 4.* var / meanSquare; + const Float_t meanSquare = GetMean() * GetMean(); + const Float_t meanSquareRelVar = 4.* GetMeanRelVar(); const Float_t ffactorsquare = gkFFactor * gkFFactor; - const Float_t ffactorsquareRelVar = 4.* gkFFactorErr * gkFFactorErr / ffactorsquare; + const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar(); const Float_t rsigmaSquareRelVar = fRSigmaSquareVar / fRSigmaSquare; @@ -804,7 +862,7 @@ const Float_t conversionSquareRelVar = 4.* GetConversionHiLoRelVar(); - const Float_t logainNsbSquare = higainNsbSquare / conversionSquare; - const Float_t logainNsbSquareVar = ( higainNsbSquareRelVar + conversionSquareRelVar ) - * logainNsbSquare * logainNsbSquare; + const Float_t logainNsbSquare = higainNsbSquare / conversionSquare; + const Float_t logainNsbSquareVar = ( higainNsbSquareRelVar + conversionSquareRelVar ) + * logainNsbSquare * logainNsbSquare; fLoGainPedRmsSquare = logainNsbSquare + elecRmsSquare; @@ -812,16 +870,2 @@ } -// -------------------------------------------------------------------------- -// -// Get the relative variance of the conversion factor between higain and logain: -// - If fConversionHiLo is smaller than 0, return -1. -// - Else returns: fConversionHiLoVar / (fConversionHiLo * fConversionHiLo) -// -const Float_t MCalibrationChargePix::GetConversionHiLoRelVar() const -{ - if (fConversionHiLo == 0.) - return 0.; - - return fConversionHiLoVar / (fConversionHiLo * fConversionHiLo); -} - Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.h (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.h (revision 3672) @@ -59,6 +59,6 @@ const Float_t GetConversionHiLoRelVar() const; - const Float_t GetChargeRelVar() const; - + const Float_t GetFFactorRelVar() const; + public: @@ -93,10 +93,10 @@ // Getters - Float_t GetAbsTimeMean () const { return fAbsTimeMean; } - Float_t GetAbsTimeRms () const { return fAbsTimeRms; } - Float_t GetConversionHiLo () const { return fConversionHiLo; } - Float_t GetConversionHiLoErr() const; + Float_t GetAbsTimeMean () const { return fAbsTimeMean; } + Float_t GetAbsTimeRms () const { return fAbsTimeRms; } + Float_t GetConversionHiLo () const { return fConversionHiLo; } + Float_t GetConversionHiLoErr () const; Float_t GetConversionBlindPixelMethodErr() const; - Float_t GetConversionFFactorMethodErr() const; + Float_t GetConversionFFactorMethodErr () const; Float_t GetConversionPINDiodeMethodErr() const; Float_t GetConversionCombinedMethodErr() const; @@ -121,4 +121,5 @@ Float_t GetRSigma () const; Float_t GetRSigmaErr () const; + Float_t GetRSigmaRelVar () const; Float_t GetTotalFFactorBlindPixelMethod() const { return fTotalFFactorBlindPixelMethod; } Float_t GetTotalFFactorBlindPixelMethodErr() const; @@ -136,5 +137,5 @@ // Miscellaneous - void CalcLoGainPedestal(const Float_t logainsamples); + void CalcLoGainPedestal(const Float_t logainsamples); Bool_t CalcReducedSigma(); Bool_t CalcFFactorMethod(); Index: trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3672) @@ -208,4 +208,18 @@ // -------------------------------------------------------------------------- // +// Get the Low Gain Mean Rel Variance +// +Float_t MCalibrationPix::GetLoGainMeanRelVar() const +{ + if (fLoGainMeanVar < 0.) + return -1.; + if (fLoGainMean == 0.) + return -1.; + + return fLoGainMeanVar / ( fLoGainMean * fLoGainMean); +} + +// -------------------------------------------------------------------------- +// // Get the Low Gain Sigma Error: Takes square root of fHiGainSigmaVar // Index: trunk/MagicSoft/Mars/mcalib/MCalibrationPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationPix.h (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationPix.h (revision 3672) @@ -73,10 +73,11 @@ Float_t GetHiGainNumPickup() const { return fHiGainNumPickup; } - virtual Float_t GetLoGainMean () const { return fLoGainMean ; } - virtual Float_t GetLoGainMeanErr () const; - virtual Float_t GetLoGainProb () const { return fLoGainProb ; } - virtual Float_t GetLoGainSigma () const { return fLoGainSigma ; } - virtual Float_t GetLoGainSigmaErr () const; - virtual Float_t GetLoGainNumPickup() const { return fLoGainNumPickup; } + virtual Float_t GetLoGainMean () const { return fLoGainMean ; } + virtual Float_t GetLoGainMeanErr () const; + virtual Float_t GetLoGainMeanRelVar() const; + virtual Float_t GetLoGainProb () const { return fLoGainProb ; } + virtual Float_t GetLoGainSigma () const { return fLoGainSigma ; } + virtual Float_t GetLoGainSigmaErr () const; + virtual Float_t GetLoGainNumPickup () const { return fLoGainNumPickup; } Float_t GetMean () const { return IsHiGainSaturation() ? GetLoGainMean() : GetHiGainMean() ; } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 3671) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 3672) @@ -441,17 +441,18 @@ fBlindPix->SetSinglePheFitOK(); - fBlindPix->SetLambda ( fLambda ); - fBlindPix->SetMu0 ( fMu0 ); - fBlindPix->SetMu0Err ( fMu0Err ); - fBlindPix->SetMu1 ( fMu1 ); - fBlindPix->SetMu1Err ( fMu1Err ); - fBlindPix->SetSigma0 ( fSigma0 ); - fBlindPix->SetSigma0Err ( fSigma0Err ); - fBlindPix->SetSigma1 ( fSigma1 ); - fBlindPix->SetSigma1Err ( fSigma1Err ); - fBlindPix->SetProb ( fProb ); - - fBlindPix->SetLambdaCheck ( fLambdaCheck ); - fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr ); + fBlindPix->SetLambda ( fLambda ); + fBlindPix->SetLambdaVar ( fLambdaErr*fLambdaErr ); + fBlindPix->SetMu0 ( fMu0 ); + fBlindPix->SetMu0Err ( fMu0Err ); + fBlindPix->SetMu1 ( fMu1 ); + fBlindPix->SetMu1Err ( fMu1Err ); + fBlindPix->SetSigma0 ( fSigma0 ); + fBlindPix->SetSigma0Err ( fSigma0Err ); + fBlindPix->SetSigma1 ( fSigma1 ); + fBlindPix->SetSigma1Err ( fSigma1Err ); + fBlindPix->SetProb ( fProb ); + + fBlindPix->SetLambdaCheck ( fLambdaCheck ); + fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr ); return kTRUE;