/* ======================================================================== *\ ! ! * ! * 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 02/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MCalibrationChargePix // // // // Storage container to hold informations about the calibration values // // values of one Pixel (PMT). // // // // The following values are initialized to meaningful values: // // - The Electronic Rms to 1.5 per FADC slice // - The uncertainty about the Electronic RMS to 0.3 per slice // - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17. // 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. // // // 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 // ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationChargePix.h" #include "MLog.h" #include "MLogManip.h" #include "MBadPixelsPix.h" ClassImp(MCalibrationChargePix); using namespace std; const Float_t MCalibrationChargePix::gkElectronicPedRms = 1.5; const Float_t MCalibrationChargePix::gkElectronicPedRmsErr = 0.3; const Float_t MCalibrationChargePix::gkFFactor = 1.15; const Float_t MCalibrationChargePix::gkFFactorErr = 0.02; const Float_t MCalibrationChargePix::gkConversionHiLo = 10.; const Float_t MCalibrationChargePix::gkConversionHiLoErr = 2.5; const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit = 5.; // -------------------------------------------------------------------------- // // Default Constructor: // MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title) : fPixId(-1), fFlags(0) { fName = name ? name : "MCalibrationChargePix"; fTitle = title ? title : "Container of the fit results of MHCalibrationChargePixs "; Clear(); // // At the moment, we don't have a database, yet, // so we get it from the configuration file // SetConversionHiLo(); SetConversionHiLoErr(); SetPheFFactorMethodLimit(); } // ------------------------------------------------------------------------ // // Invalidate values // void MCalibrationChargePix::Clear(Option_t *o) { SetHiGainSaturation ( kFALSE ); SetLoGainSaturation ( kFALSE ); SetExcluded ( kFALSE ); SetBlindPixelMethodValid ( kFALSE ); SetFFactorMethodValid ( kFALSE ); SetPINDiodeMethodValid ( kFALSE ); SetCombinedMethodValid ( kFALSE ); fHiGainMeanCharge = -1.; fHiGainMeanChargeVar = -1.; fHiGainSigmaCharge = -1.; fHiGainSigmaChargeVar = -1.; fHiGainChargeProb = -1.; fLoGainMeanCharge = -1.; fLoGainMeanChargeVar = -1.; fLoGainSigmaCharge = -1.; fLoGainSigmaChargeVar = -1.; fLoGainChargeProb = -1.; fRSigmaCharge = -1.; fRSigmaChargeVar = -1.; fHiGainNumPickup = -1; fLoGainNumPickup = -1; fPed = -1.; fPedRms = -1.; fPedVar = -1.; fLoGainPedRms = -1.; fLoGainPedRmsVar = -1.; fAbsTimeMean = -1.; fAbsTimeRms = -1.; fPheFFactorMethod = -1.; fPheFFactorMethodVar = -1.; fMeanConversionFFactorMethod = -1.; fMeanConversionBlindPixelMethod = -1.; fMeanConversionPINDiodeMethod = -1.; fMeanConversionCombinedMethod = -1.; fConversionFFactorMethodVar = -1.; fConversionBlindPixelMethodVar = -1.; fConversionPINDiodeMethodVar = -1.; fConversionCombinedMethodVar = -1.; fSigmaConversionFFactorMethod = -1.; fSigmaConversionBlindPixelMethod = -1.; fSigmaConversionPINDiodeMethod = -1.; fSigmaConversionCombinedMethod = -1.; fTotalFFactorFFactorMethod = -1.; fTotalFFactorBlindPixelMethod = -1.; fTotalFFactorPINDiodeMethod = -1.; fTotalFFactorCombinedMethod = -1.; } // -------------------------------------------------------------------------- // // 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; } void MCalibrationChargePix::SetMeanCharge( const Float_t f ) { if (IsHiGainSaturation()) fLoGainMeanCharge = f; else fHiGainMeanCharge = f; } void MCalibrationChargePix::SetMeanChargeErr( const Float_t f ) { if (IsHiGainSaturation()) fLoGainMeanChargeVar = f*f; else fHiGainMeanChargeVar = f*f; } void MCalibrationChargePix::SetSigmaCharge( const Float_t f ) { if (IsHiGainSaturation()) fLoGainSigmaCharge = f; else fHiGainSigmaCharge = f; } void MCalibrationChargePix::SetSigmaChargeErr( const Float_t f ) { if (IsHiGainSaturation()) fLoGainSigmaChargeVar = f*f; else fHiGainSigmaChargeVar = f*f; } // -------------------------------------------------------------------------- // // Set the conversion factors from outside (only for MC) // void MCalibrationChargePix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig) { fMeanConversionFFactorMethod = c; fConversionFFactorMethodVar = err*err; fSigmaConversionFFactorMethod = sig; } // -------------------------------------------------------------------------- // // 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) // void MCalibrationChargePix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig) { fMeanConversionPINDiodeMethod = c ; fConversionPINDiodeMethodVar = err*err; fSigmaConversionPINDiodeMethod = sig; } // -------------------------------------------------------------------------- // // Set the Hi Gain Saturation Bit from outside // void MCalibrationChargePix::SetHiGainSaturation(Bool_t b) { b ? SETBIT(fFlags, kHiGainSaturation) : CLRBIT(fFlags, kHiGainSaturation); } // -------------------------------------------------------------------------- // // Set the Lo Gain Saturation Bit from outside // void MCalibrationChargePix::SetLoGainSaturation(Bool_t b) { b ? SETBIT(fFlags, kLoGainSaturation) : CLRBIT(fFlags, kLoGainSaturation); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationChargePix::SetExcluded(Bool_t b ) { b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b ) { b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationChargePix::SetFFactorMethodValid(const Bool_t b ) { b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationChargePix::SetPINDiodeMethodValid(const Bool_t b ) { b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); } // -------------------------------------------------------------------------- // // Set the Excluded Bit from outside // void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b ) { b ? SETBIT(fFlags, kCombinedMethodValid) : CLRBIT(fFlags, kCombinedMethodValid); } Float_t MCalibrationChargePix::GetPedRms() const { return IsHiGainSaturation() ? fLoGainPedRms : fPedRms; } Float_t MCalibrationChargePix::GetPedRmsErr() const { return IsHiGainSaturation() ? TMath::Sqrt(fLoGainPedRmsVar) : TMath::Sqrt(fPedVar)/2.; } Float_t MCalibrationChargePix::GetPedErr() const { return TMath::Sqrt(fPedVar); } Float_t MCalibrationChargePix::GetMeanCharge() const { return IsHiGainSaturation() ? GetLoGainMeanCharge() : GetHiGainMeanCharge() ; } Float_t MCalibrationChargePix::GetMeanChargeErr() const { return IsHiGainSaturation() ? GetLoGainMeanChargeErr() : GetHiGainMeanChargeErr() ; } Float_t MCalibrationChargePix::GetChargeProb() const { return IsHiGainSaturation() ? fLoGainChargeProb : fHiGainChargeProb ; } Float_t MCalibrationChargePix::GetSigmaCharge() const { return IsHiGainSaturation() ? GetLoGainSigmaCharge() : GetHiGainSigmaCharge() ; } Float_t MCalibrationChargePix::GetSigmaChargeErr() const { return IsHiGainSaturation() ? GetLoGainSigmaChargeErr() : GetHiGainSigmaChargeErr() ; } Float_t MCalibrationChargePix::GetHiGainMeanChargeErr() const { return TMath::Sqrt(fHiGainMeanChargeVar); } Float_t MCalibrationChargePix::GetLoGainMeanCharge() const { return fLoGainMeanCharge * fConversionHiLo; } Float_t MCalibrationChargePix::GetLoGainMeanChargeErr() const { const Float_t chargeRelVar = fLoGainMeanChargeVar /( fLoGainMeanCharge * fLoGainMeanCharge ); const Float_t conversionRelVar = fConversionHiLoVar /( fConversionHiLo * fConversionHiLo ); return TMath::Sqrt(chargeRelVar+conversionRelVar) * GetLoGainMeanCharge(); } Float_t MCalibrationChargePix::GetLoGainSigmaCharge() const { return fLoGainSigmaCharge * fConversionHiLo; } Float_t MCalibrationChargePix::GetLoGainSigmaChargeErr() const { const Float_t sigmaRelVar = fLoGainSigmaChargeVar /( fLoGainSigmaCharge * fLoGainSigmaCharge ); const Float_t conversionRelVar = fConversionHiLoVar /( fConversionHiLo * fConversionHiLo ); return TMath::Sqrt(sigmaRelVar+conversionRelVar) * GetLoGainSigmaCharge(); } Float_t MCalibrationChargePix::GetHiGainSigmaChargeErr() const { return TMath::Sqrt(fHiGainSigmaChargeVar); } Float_t MCalibrationChargePix::GetRSigmaCharge() const { return IsHiGainSaturation() ? fRSigmaCharge*fConversionHiLo : fRSigmaCharge ; } Float_t MCalibrationChargePix::GetRSigmaChargeErr() const { if (IsHiGainSaturation()) { const Float_t rsigmaRelVar = fRSigmaChargeVar /( fRSigmaCharge * fRSigmaCharge ); const Float_t conversionRelVar = fConversionHiLoVar /( fConversionHiLo * fConversionHiLo ); return TMath::Sqrt(rsigmaRelVar+conversionRelVar) * GetRSigmaCharge(); } else return TMath::Sqrt(fRSigmaChargeVar); } Float_t MCalibrationChargePix::GetConversionHiLoErr() const { if (fConversionHiLoVar < 0.) return -1.; return TMath::Sqrt(fConversionHiLoVar); } Float_t MCalibrationChargePix::GetPheFFactorMethodErr() const { if (fPheFFactorMethodVar < 0.) return -1.; return TMath::Sqrt(fPheFFactorMethodVar); } Float_t MCalibrationChargePix::GetConversionCombinedMethodErr() const { if (fConversionCombinedMethodVar < 0.) return -1.; return TMath::Sqrt(fConversionCombinedMethodVar); } Float_t MCalibrationChargePix::GetConversionPINDiodeMethodErr() const { if (fConversionPINDiodeMethodVar < 0.) return -1.; return TMath::Sqrt(fConversionPINDiodeMethodVar); } Float_t MCalibrationChargePix::GetConversionBlindPixelMethodErr() const { if (fConversionBlindPixelMethodVar < 0.) return -1.; return TMath::Sqrt(fConversionBlindPixelMethodVar); } Float_t MCalibrationChargePix::GetConversionFFactorMethodErr() const { if (fConversionFFactorMethodVar < 0.) return -1.; return TMath::Sqrt(fConversionFFactorMethodVar); } Float_t MCalibrationChargePix::GetTotalFFactorCombinedMethodErr() const { if (fTotalFFactorCombinedMethodVar < 0.) return -1.; return TMath::Sqrt(fTotalFFactorCombinedMethodVar); } Float_t MCalibrationChargePix::GetTotalFFactorPINDiodeMethodErr() const { if (fTotalFFactorPINDiodeMethodVar < 0.) return -1.; return TMath::Sqrt(fTotalFFactorPINDiodeMethodVar); } Float_t MCalibrationChargePix::GetTotalFFactorBlindPixelMethodErr() const { if (fTotalFFactorBlindPixelMethodVar < 0.) return -1.; return TMath::Sqrt(fTotalFFactorBlindPixelMethodVar); } Float_t MCalibrationChargePix::GetTotalFFactorFFactorMethodErr() const { if (fTotalFFactorFFactorMethodVar < 0.) return -1.; return TMath::Sqrt(fTotalFFactorFFactorMethodVar); } Bool_t MCalibrationChargePix::IsExcluded() const { return TESTBIT(fFlags,kExcluded); } Bool_t MCalibrationChargePix::IsHiGainSaturation() const { return TESTBIT(fFlags,kHiGainSaturation); } Bool_t MCalibrationChargePix::IsLoGainSaturation() const { return TESTBIT(fFlags,kLoGainSaturation); } Bool_t MCalibrationChargePix::IsBlindPixelMethodValid() const { return TESTBIT(fFlags, kBlindPixelMethodValid); } Bool_t MCalibrationChargePix::IsFFactorMethodValid() const { return TESTBIT(fFlags, kFFactorMethodValid); } Bool_t MCalibrationChargePix::IsPINDiodeMethodValid() const { return TESTBIT(fFlags, kPINDiodeMethodValid); } Bool_t MCalibrationChargePix::IsCombinedMethodValid() const { return TESTBIT(fFlags, kCombinedMethodValid); } // // // Bool_t MCalibrationChargePix::CalcReducedSigma() { const Float_t sigmacharge = IsHiGainSaturation() ? fLoGainSigmaCharge : fHiGainSigmaCharge ; const Float_t sigmachargevar = IsHiGainSaturation() ? fLoGainSigmaChargeVar : fHiGainSigmaChargeVar; const Float_t sigmaSquare = sigmacharge * sigmacharge; const Float_t sigmaSquareVar = 4.* sigmachargevar * 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. } // // Calculate the reduced sigmas // const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare; if (rsigmachargesquare <= 0.) { *fLog << warn << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel " << fPixId << endl; return kFALSE; } fRSigmaCharge = TMath::Sqrt(rsigmachargesquare); fRSigmaChargeVar = 0.25 * (sigmaSquareVar + pedRmsSquareVar) / rsigmachargesquare; return kTRUE; } // // Calculate the number of photo-electrons after the F-Factor method // Calculate the errors of the F-Factor method // Bool_t MCalibrationChargePix::CalcFFactorMethod() { if (fRSigmaCharge < 0.) { SetFFactorMethodValid(kFALSE); return kFALSE; } const Float_t charge = IsHiGainSaturation() ? fLoGainMeanCharge : fHiGainMeanCharge ; const Float_t chargevar = IsHiGainSaturation() ? fLoGainMeanChargeVar : fHiGainMeanChargeVar; // // Square all variables in order to avoid applications of square root // // First the relative error squares // const Float_t chargeSquare = charge * charge; const Float_t chargeSquareRelVar = 4.* chargevar/ chargeSquare; const Float_t ffactorsquare = gkFFactor * gkFFactor; const Float_t ffactorsquareRelVar = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare; const Float_t rsigmaSquare = fRSigmaCharge * fRSigmaCharge; const Float_t rsigmaSquareRelVar = 4.* fRSigmaChargeVar / rsigmaSquare; // // Calculate the number of phe's from the F-Factor method // (independent on Hi Gain or Lo Gain) // fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmaSquare; if (fPheFFactorMethod < fPheFFactorMethodLimit) { SetFFactorMethodValid(kFALSE); return kFALSE; } // // Calculate the Error of Nphe // fPheFFactorMethodVar = (ffactorsquareRelVar + chargeSquareRelVar + rsigmaSquareRelVar) * fPheFFactorMethod * fPheFFactorMethod; SetFFactorMethodValid(kTRUE); return kTRUE; } void MCalibrationChargePix::CalcLoGainPedestal(Float_t logainsamples) { fElectronicPedRms = gkElectronicPedRms * TMath::Sqrt(logainsamples); fElectronicPedRmsVar = gkElectronicPedRmsErr * gkElectronicPedRmsErr * logainsamples; Float_t pedRmsSquare = fPedRms * fPedRms; Float_t pedRmsSquareVar = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2. // // 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): // // We extract the pure NSB contribution: // const Float_t elecRmsSquare = fElectronicPedRms * fElectronicPedRms; const Float_t elecRmsSquareVar = 4.*fElectronicPedRmsVar * elecRmsSquare; Float_t nsbSquare = pedRmsSquare - elecRmsSquare; Float_t nsbSquareRelVar = (pedRmsSquareVar + elecRmsSquareVar) / (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 convertedNsbSquare = nsbSquare / conversionSquare; const Float_t convertedNsbSquareVar = nsbSquareRelVar * convertedNsbSquare * convertedNsbSquare; pedRmsSquare = convertedNsbSquare + elecRmsSquare; pedRmsSquareVar = convertedNsbSquareVar + elecRmsSquareVar; fLoGainPedRms = TMath::Sqrt(pedRmsSquare); fLoGainPedRmsVar = 0.25 * pedRmsSquareVar / pedRmsSquare; }