Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2715) +++ trunk/MagicSoft/Mars/Changelog (revision 2716) @@ -12,5 +12,7 @@ * manalysis/MCalibrationPix.cc - introduced error calculation for the F-Factor method - + + * macros/calibration.C + - display now the errors in the F-Factor method correctly 2003/12/18: Abelardo Moralejo Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h (revision 2715) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h (revision 2716) @@ -99,6 +99,7 @@ UInt_t GetNumPixels() const { return fNumPixels; } - Bool_t IsPixelUsed(Int_t idx) const; - Bool_t IsPixelFitted(Int_t idx) const; + Bool_t IsPixelUsed(Int_t idx) const; + Bool_t IsPixelFitted(Int_t idx) const; + Bool_t IsPixelFitValid(Int_t idx) const; MCalibrationPix &operator[](Int_t i); @@ -114,7 +115,6 @@ void DrawPixelContent(Int_t num) const; - MCalibrationPix *GetCalibrationPix(Int_t idx) const; - MCalibrationBlindPix *GetBlindPixel() const { return fBlindPixel; } - MCalibrationPINDiode *GetPINDiode() const { return fPINDiode; } + MCalibrationBlindPix *GetBlindPixel() const { return fBlindPixel; } + MCalibrationPINDiode *GetPINDiode() const { return fPINDiode; } void SetColor(CalibrationColor_t color) { fColor = color; } Index: trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h (revision 2715) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h (revision 2716) @@ -17,5 +17,6 @@ // The conversion factor between High Gain and Low Gain -const Float_t gkConversionHiLo = 10.; +const Float_t gkConversionHiLo = 10.; +const Float_t gkConversionHiLoError = 2.5; // The penalty constant to produce overflow in the histogram @@ -40,9 +41,9 @@ // -// Area of Inner Pixel w.r.t. Blind Pixel (which is 1 cm²) +// Area of Inner Pixel w.r.t. Blind Pixel (which is 1 sq. cm) // // Hexagone of diagonal axis b = 3.5 cm // straight axis a = 3.0 cm +- 2% -// Area = sqrt(3)*a²/2 = 7.79 cm² +- 4% = 7.8 +- 0.3 cm² +// Area = sqrt(3)*a*a/2 = 7.79 sq.cm +- 4% = 7.8 +- 0.3 sq.cm // const Float_t gkCalibrationInnerPixelArea = 7.8; @@ -53,5 +54,5 @@ // Hexagone of diagonal axis b = 7.0 cm // straight axis a = 6.0 cm +- 1% -// Area = sqrt(3)*a²/2 = +// Area = sqrt(3)*a*a/2 = // const Float_t gkCalibrationOutervsInnerPixelArea = 4.00; @@ -68,5 +69,5 @@ // Hexagone of diagonal axis b = 3.5 cm // straight axis a = 3.0 cm +- 2% -// Area = sqrt(3)*a²/2 = 7.79 cm² +- 4% = 7.8 +- 0.3 cm² +// Area = sqrt(3)*a*a/2 = 7.79 sq.cm +- 4% = 7.8 +- 0.3 sq.cm // // Distance of PIN Diode to pulser D1: 1.5 +- 0.3 m @@ -74,7 +75,7 @@ // // -// A(Inner Pixel) D1² -// conversion C = -------------- * ----- = 0.054 -// A(PIN Diode) D2² +// A(Inner Pixel) D1*D1 +// conversion C = -------------- * ------ = 0.054 +// A(PIN Diode) D2*D2 // // Delta C / C = sqrt((Delta A(IP)/A(IP))² + 4 * ( (Delta D1/D1)² + (Delta D2/D2)² ) Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc (revision 2715) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc (revision 2716) @@ -32,4 +32,5 @@ ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationPix.h" +#include "MCalibrationConfig.h" #include "MLog.h" @@ -54,4 +55,7 @@ fPed(-1.), fPedRms(-1.), + fErrPedRms(0.), + fElectronicPedRms(1.5), + fErrElectronicPedRms(0.3), fTime(-1.), fSigmaTime(-1.), @@ -68,6 +72,5 @@ fConversionSigmaBlindPixelMethod(-1.), fConversionSigmaPINDiodeMethod(-1.), - fHiGainSaturation(kFALSE), - fElectronicPedRms(1.5) + fHiGainSaturation(kFALSE) { @@ -103,4 +106,127 @@ } + +Bool_t MCalibrationPix::IsFitted() +{ + return (fCharge > 0.) && (fErrCharge > 0.) ; +} + +// -------------------------------------------------------------------------- +// +// Return TRUE if: +// +// 1) Pixel has a fitted charge greater than 3*PedRMS +// 2) Pixel has a fit error greater than 0. +// 3) Pixel has a fit Probability greater than 0.0001 +// 4) Pixel has a charge sigma bigger than its Pedestal RMS +// 5) If FitTimes is used, +// the mean arrival time is at least 1.5 slices from the edge +// +Bool_t MCalibrationPix::IsFitValid() +{ + + if (fCharge < 3.*GetPedRms()) + return kFALSE; + + if (fErrCharge <= 0.) + return kFALSE; + + if (!fHist->IsFitOK()) + return kFALSE; + + if (fSigmaCharge <= fPedRms) + return kFALSE; + + if (fHist->GetTimeLowerFitRange() < fTime < fHist->GetTimeLowerFitRange()+1.) + return kFALSE; + + if (fHist->GetTimeUpperFitRange()-1. < fTime < fHist->GetTimeUpperFitRange()) + return kFALSE; + + return kTRUE; + +} + +// -------------------------------------------------------------------------- +// +// Return TRUE if: +// +// 1) Pixel is FitValid +// 2) Conversion Factor is bigger than 0. +// 3) The error of the conversion factor is smaller than 10% +// +Bool_t MCalibrationPix::IsBlindPixelMethodValid() +{ + + if (!IsFitValid()) + return kFALSE; + + if (fConversionBlindPixelMethod <= 0.) + return kFALSE; + + if (fConversionErrorBlindPixelMethod/fConversionBlindPixelMethod > 0.1) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Return TRUE if: +// +// 1) Pixel is FitValid +// 2) Conversion Factor is bigger than 0. +// 3) The error of the conversion factor is smaller than 10% +// +Bool_t MCalibrationPix::IsFFactorMethodValid() +{ + + if (!IsFitValid()) + return kFALSE; + + if (fConversionFFactorMethod <= 0.) + return kFALSE; + + if (fConversionErrorFFactorMethod/fConversionFFactorMethod > 0.1) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Return TRUE if: +// +// 1) Pixel is FitValid +// 2) Conversion Factor is bigger than 0. +// 3) The error of the conversion factor is smaller than 10% +// +Bool_t MCalibrationPix::IsPINDiodeMethodValid() +{ + + if (!IsFitValid()) + return kFALSE; + + if (fConversionPINDiodeMethod <= 0.) + return kFALSE; + + if (fConversionErrorPINDiodeMethod/fConversionPINDiodeMethod > 0.1) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// 1) Return if the charge distribution is already succesfully fitted +// 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 +// Bool_t MCalibrationPix::FitCharge() { @@ -123,5 +249,4 @@ *fLog << warn << "Could not fit Lo Gain charges of pixel " << fPixId << endl; fHist->PrintChargeFitResult(); - // return kFALSE; } } @@ -132,5 +257,4 @@ *fLog << warn << "Could not fit Hi Gain charges of pixel " << fPixId << endl; fHist->PrintChargeFitResult(); - // return kFALSE; } } @@ -143,39 +267,113 @@ fChargeProb = fHist->GetChargeProb(); + if (fCharge <= 0.) + { + *fLog << warn << "Cannot apply calibration: Mean Fitted Charges are smaller than 0 in pixel " + << fPixId << endl; + return kFALSE; + } + + + // + // Apply the F-Factor Method + // if ((fPed > 0.) && (fPedRms > 0.)) { - Float_t pedrmssquare = fPedRms*fPedRms; - Float_t sigmasquare = fSigmaCharge*fSigmaCharge; + // + // Square all variables in order to avoid applications of square root + // + const Float_t chargesquare = fCharge*fCharge; + const Float_t chargessquarerelerrsquare = 4.*fErrCharge*fErrCharge/chargesquare; + + const Float_t sigmasquare = fSigmaCharge*fSigmaCharge; + const Float_t sigmasquareerr = 2.*fErrSigmaCharge*fSigmaCharge; + + const Float_t pedrmssquare = fPedRms*fPedRms; + const Float_t pedrmssquareerr = 2.*fErrPedRms*fPedRms; + + const Float_t elecrmssquare = fElectronicPedRms*fElectronicPedRms; + const Float_t elecrmssquareerr = 2.*fErrElectronicPedRms*fElectronicPedRms; + + const Float_t conversionsquare = gkConversionHiLo *gkConversionHiLo; + const Float_t conversionsquareerr = 2.*gkConversionHiLoError*gkConversionHiLo; + + const Float_t ffactorrelerrsquare = fFactorError * fFactorError / fFactor / fFactor; + + Float_t rsigmasquarerelerrsquare = 0.; + Float_t pheffactorrelerrsquare = 0.; 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 nsbsquareerrsquare = pedrmssquareerr * pedrmssquareerr + + elecrmssquareerr * elecrmssquareerr; + + if (nsbsquare < 0.) + nsbsquare = 0.; - Float_t logainrmssquare = fElectronicPedRms*fElectronicPedRms; - Float_t nsbsquare = pedrmssquare - logainrmssquare; - // Float_t logainrms = fElectronicPedRms + (TMath::Sqrt(fPedRms*fPedRms - fElectronicPedRms*fElectronicPedRms)); - - if (nsbsquare > 0.) - logainrmssquare = nsbsquare/100. + logainrmssquare; - - fRSigmaSquare = sigmasquare - logainrmssquare; + // + // Now, we divide the NSB by the conversion factor and + // add it quadratically to the electronic noise + // + Float_t logainrmssquare = nsbsquare/conversionsquare + elecrmssquare; + Float_t logainrmssquareerrsquare = nsbsquareerrsquare/conversionsquare/conversionsquare + + elecrmssquareerr * elecrmssquareerr ; + // + // Calculate the reduced sigma with the new "Pedestal RMS" + // + fRSigmaSquare = sigmasquare - logainrmssquare; + rsigmasquarerelerrsquare = (sigmasquareerr * sigmasquareerr + logainrmssquareerrsquare) + / (fRSigmaSquare * fRSigmaSquare); if (fRSigmaSquare > 0.) - fPheFFactorMethod = fFactor*(fCharge*fCharge/100.) / fRSigmaSquare; - - } - else /* if (fHiGainSaturation) */ + { + fPheFFactorMethod = fFactor * chargesquare / fRSigmaSquare; + pheffactorrelerrsquare = ffactorrelerrsquare + + chargessquarerelerrsquare + + rsigmasquarerelerrsquare ; + } + + } /* if (fHiGainSaturation) */ + else { - fRSigmaSquare = sigmasquare - pedrmssquare; - fPheFFactorMethod = fFactor * fCharge*fCharge / fRSigmaSquare; + fRSigmaSquare = sigmasquare - pedrmssquare; + rsigmasquarerelerrsquare = ( sigmasquareerr * sigmasquareerr + + pedrmssquareerr * pedrmssquareerr) + / (fRSigmaSquare * fRSigmaSquare); + + fPheFFactorMethod = fFactor * chargesquare / fRSigmaSquare; + // + // We first calculate the squared relative error in order to save the + // TMath::Sqrt + // + pheffactorrelerrsquare = ( ffactorrelerrsquare + + chargessquarerelerrsquare + + rsigmasquarerelerrsquare ); + + + } /* if (fHiGainSaturation) */ + + + if (fPheFFactorMethod <= 0.) + { + *fLog << warn << "Cannot apply F-Factor calibration: Number of PhE smaller than 0 in pixel " + << fPixId << endl; + return kFALSE; } - - if (fCharge > 0.) - fConversionFFactorMethod = fPheFFactorMethod / fCharge ; - - else - *fLog << warn << "Cannot apply F-Factor method: Reduced Sigmas are smaller than 0 in pixel: " - << fPixId << endl; + fConversionFFactorMethod = fPheFFactorMethod / fCharge ; + fConversionErrorFFactorMethod = fConversionFFactorMethod * + ( pheffactorrelerrsquare + + (fErrCharge * fErrCharge / chargesquare ) ); + + fPheFFactorMethodError = TMath::Sqrt(pheffactorrelerrsquare) * fPheFFactorMethod; } /* if ((fPed > 0.) && (fPedRms > 0.)) */ @@ -186,4 +384,8 @@ +// -------------------------------------------------------------------------- +// +// Set the pedestals from outside +// void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms) { @@ -194,7 +396,21 @@ } +// -------------------------------------------------------------------------- +// +// 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) { @@ -206,4 +422,8 @@ } } + + // + // Fit the High Gain + // else { Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h (revision 2715) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h (revision 2716) @@ -23,4 +23,7 @@ Float_t fPed; // The mean pedestal (from MPedestalPix) Float_t fPedRms; // The pedestal RMS (from MPedestalPix) + Float_t fErrPedRms; // The error of the pedestal RMS (from MPedestalPix) + Float_t fElectronicPedRms; // The pure electronic component of the RMS + Float_t fErrElectronicPedRms; // The error of the pure electronic component of the RMS Float_t fTime; // The mean arrival time after the fit @@ -29,21 +32,21 @@ Float_t fFactor; // The laboratory F-factor + Float_t fFactorError; // The laboratory F-factor Error Float_t fPheFFactorMethod; // The number of Phe's calculated after the F-factor method + Float_t fPheFFactorMethodError; // The error on the number of Phe's calculated after the F-factor method - Float_t fConversionFFactorMethod; // The conversion factor to Ph's calculated after the F-factor method + Float_t fConversionFFactorMethod; // The conversion factor to Ph's calculated after the F-factor method Float_t fConversionBlindPixelMethod; // The conversion factor to Ph's calculated after the Blind Pixel method Float_t fConversionPINDiodeMethod; // The conversion factor to Ph's calculated after the PIN Diode method - Float_t fConversionErrorFFactorMethod; // The conversion factor to Phe's calculated after the F-factor method + Float_t fConversionErrorFFactorMethod; // The conversion factor to Phe's calculated after the F-factor method Float_t fConversionErrorBlindPixelMethod; // The conversion factor to Phe's calculated after the Blind Pixel method Float_t fConversionErrorPINDiodeMethod; // The conversion factor to Phe's calculated after the PIN Diode method - Float_t fConversionSigmaFFactorMethod; // The conversion factor to Ph's calculated after the F-factor method + Float_t fConversionSigmaFFactorMethod; // The conversion factor to Ph's calculated after the F-factor method Float_t fConversionSigmaBlindPixelMethod; // The conversion factor to Ph's calculated after the Blind Pixel method Float_t fConversionSigmaPINDiodeMethod; // The conversion factor to Phd's calculated after the PIN Diode method Bool_t fHiGainSaturation; // Is Lo-Gain used at all? - - Float_t fElectronicPedRms; MHCalibrationPixel *fHist; //! Pointer to the histograms performing the fits, etc. @@ -98,10 +101,12 @@ Float_t GetSigmaConversionBlindPixelMethod() const { return fConversionSigmaBlindPixelMethod ; } - Float_t GetMeanConversionFFactorMethod() const { return fConversionFFactorMethod ; } - Float_t GetErrorConversionFFactorMethod() const { return fConversionErrorFFactorMethod ; } - Float_t GetSigmaConversionFFactorMethod() const { return fConversionSigmaFFactorMethod ; } - Float_t GetPheFFactorMethod() const { return fPheFFactorMethod; } + Float_t GetMeanConversionFFactorMethod() const { return fConversionFFactorMethod ; } + Float_t GetErrorConversionFFactorMethod() const { return fConversionErrorFFactorMethod ; } + Float_t GetSigmaConversionFFactorMethod() const { return fConversionSigmaFFactorMethod ; } + + Float_t GetPheFFactorMethod() const { return fPheFFactorMethod; } + Float_t GetPheFFactorMethodError() const { return fPheFFactorMethodError; } - Float_t GetMeanConversionPINDiodeMethod() const { return fConversionPINDiodeMethod ; } + Float_t GetMeanConversionPINDiodeMethod() const { return fConversionPINDiodeMethod ; } Float_t GetErrorConversionPINDiodeMethod() const { return fConversionErrorPINDiodeMethod ; } Float_t GetSigmaConversionPINDiodeMethod() const { return fConversionSigmaPINDiodeMethod ; } @@ -117,5 +122,10 @@ Bool_t FillRChargevsTimeLoGain(Float_t rq, Int_t t) { return fHist->FillChargevsNLoGain(rq,t); } - Bool_t IsValid() const { return fCharge >= 3.*GetPedRms() || fErrCharge >= 0; } + Bool_t IsFitValid(); + Bool_t IsFitted(); + Bool_t IsBlindPixelMethodValid(); + Bool_t IsFFactorMethodValid(); + Bool_t IsPINDiodeMethodValid(); + Int_t GetPixId() const { return fPixId; } void DefinePixId(Int_t i);