Context Navigation

← Previous Change
Next Change →

Changeset 3650 for trunk/MagicSoft/Mars/mcalib

Timestamp:

04/05/04 16:44:56 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mcalib

Files:

: 1 added
: 4 edited

MCalibrationChargePix.cc (modified) (35 diffs)
MCalibrationPedCam.cc (added)
MHCalibrationCam.cc (modified) (2 diffs)
MHCalibrationChargeCam.cc (modified) (2 diffs)
MHCalibrationRelTimeCam.cc (modified) (2 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc

-              r3645
+              r3650
 \* ======================================================================== */
 /////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// 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:
 //
 …
 //   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"
 …
 // 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)
 …
 // ------------------------------------------------------------------------
 //
+// Invalidate values
+// Sets:
+// - all flags to 0
+// - all variables to -1.
+//
+// Calls:
+// - MCalibrationPix::Clear()
 //
 void MCalibrationChargePix::Clear(Option_t *o)
 …
   SetCombinedMethodValid    ( kFALSE );
   fRSigma                           =  -1.;
   fRSigmaVar                        =  -1.;
+  fRSigmaSquare                     =  -1.;
+  fRSigmaSquareVar                  =  -1.;
   fPed                              =  -1.;
 …
   fPedVar                           =  -1.;
   fLoGainPedRms                     =  -1.;
   fLoGainPedRmsVar                  =  -1.;
+  fLoGainPedRmsSquare               =  -1.;
+  fLoGainPedRmsSquareVar            =  -1.;
   fAbsTimeMean                      =  -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;
+}
+// --------------------------------------------------------------------------
+//
+// 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)
 …
 // --------------------------------------------------------------------------
 //
+// 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)
 …
 // --------------------------------------------------------------------------
 //
 // Set the Excluded Bit from outside
+// Set the BlindPixelMethod Validity Bit from outside
 //
 void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b )
 …
 // --------------------------------------------------------------------------
 //
+// 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 )
 …
 // --------------------------------------------------------------------------
 //
 // Set the Excluded Bit from outside
+// Set the PINDiodeMethod Validity Bit from outside
 //
 void MCalibrationChargePix::SetPINDiodeMethodValid(const Bool_t b )
 …
 // --------------------------------------------------------------------------
 //
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
                                  /( 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
                                 /( 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// 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
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// Test bit kBlindPixelMethodValid
+//
 Bool_t MCalibrationChargePix::IsBlindPixelMethodValid() const
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// Test bit kFFactorMethodValid
+//
 Bool_t MCalibrationChargePix::IsFFactorMethodValid()   const
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// Test bit kPINDiodeMethodValid
+//
 Bool_t MCalibrationChargePix::IsPINDiodeMethodValid()  const
+{
 …
+}
+// --------------------------------------------------------------------------
+//
+// Test bit kCombinedMethodValid
+//
 Bool_t MCalibrationChargePix::IsCombinedMethodValid()  const
+{
 …
+}
 //
+// ----------------------------------------------------------------------------
 //
+// - 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
 …
       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);
 …
   // 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;
   //
 …
   // (independent on Hi Gain or Lo Gain)
   //
   fPheFFactorMethod = ffactorsquare * meanSquare / rsigmaSquare;
+  fPheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare;
   if (fPheFFactorMethod < fPheFFactorMethodLimit)
 …
+// ----------------------------------------------------------------------------
+//
+// - 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);
 …
   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.;
   //
 …
   // 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);
+}

trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc

-              r3642
+              r3650
   hist.Renorm();
   //
   // 4) Check for oscillations
 …
   hist.CreateFourierSpectrum();
   if (!hist.IsFourierSpectrumOK())
     bad.SetUncalibrated( osctyp );

trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc

-              r3639
+              r3650
+{
+  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");
 …
       return kFALSE;
+  }
-  const Int_t npixels  = fGeom->GetNumPixels();
-  const Int_t nsectors = fGeom->GetNumSectors();
-  const Int_t nareas   = fGeom->GetNumAreas();
   if (fHiGainArray->GetEntries()==0)

trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc

-              r3642
+              r3650
+{
   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)
 …
+{
   FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
+  FitHiGainArrays((*fCam),*fBadPixels,
                   MBadPixelsPix::kRelTimeNotFitted,
                   MBadPixelsPix::kRelTimeOscillating);
   FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
+  FitLoGainArrays((*fCam),*fBadPixels,
                   MBadPixelsPix::kRelTimeNotFitted,
                   MBadPixelsPix::kRelTimeOscillating);

Note: See TracChangeset for help on using the changeset viewer.

Download in other formats: