Changeset 2716 for trunk

Timestamp:

12/18/03 17:26:36 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• manalysis/MCalibrationCam.h (modified) (2 diffs)
• manalysis/MCalibrationConfig.h (modified) (5 diffs)
• manalysis/MCalibrationPix.cc (modified) (10 diffs)
• manalysis/MCalibrationPix.h (modified) (4 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -12 +12 @@
    * 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

trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h

@@ -99 +99 @@
   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 +115 @@
   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; }

trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h

@@ -17 +17 @@
 
 // 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 +41 @@
 
 //
-// 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 +54 @@
 // 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 +69 @@
 // 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 +75 @@
 //
 //
-//                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)² )

trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc

@@ -32 +32 @@
 /////////////////////////////////////////////////////////////////////////////
 #include "MCalibrationPix.h"
+#include "MCalibrationConfig.h"
 
 #include "MLog.h"
@@ -54 +55 @@
       fPed(-1.),
       fPedRms(-1.),
+      fErrPedRms(0.),
+      fElectronicPedRms(1.5),
+      fErrElectronicPedRms(0.3),
       fTime(-1.),
       fSigmaTime(-1.),
@@ -68 +72 @@
       fConversionSigmaBlindPixelMethod(-1.),
       fConversionSigmaPINDiodeMethod(-1.),
-      fHiGainSaturation(kFALSE),
-      fElectronicPedRms(1.5)
+      fHiGainSaturation(kFALSE)
 {
 
@@ -103 +106 @@
 }
 
+
+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 +249 @@
           *fLog << warn << "Could not fit Lo Gain charges of pixel " << fPixId << endl;
           fHist->PrintChargeFitResult();
-          //      return kFALSE;
         }
     }
@@ -132 +257 @@
           *fLog << warn << "Could not fit Hi Gain charges of pixel " << fPixId << endl;
           fHist->PrintChargeFitResult();
-          //      return kFALSE;
         }
     }
@@ -143 +267 @@
   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 +384 @@
 
 
+// --------------------------------------------------------------------------
+//
+// Set the pedestals from outside
+//
 void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms)
 {
@@ -194 +396 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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 +422 @@
         }
     }
+
+  //
+  // Fit the High Gain
+  //
   else
     {

trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h

@@ -23 +23 @@
   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 +32 @@
   
   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 +101 @@
   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 +122 @@
   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);