fblanco and dsilverm - add proper pedestal, and pedestal rms, calculation
authordsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 May 2009 14:26:15 +0000 (14:26 +0000)
committerdsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 May 2009 14:26:15 +0000 (14:26 +0000)
EMCAL/AliCaloCalibPedestal.cxx
EMCAL/AliCaloCalibPedestal.h

index b176399..d29863c 100644 (file)
@@ -36,9 +36,9 @@
 
 #include "TH1.h"
 #include "TFile.h"
-
 #include <fstream>
 #include <iostream>
+#include <cmath>
 
 #include "AliCaloRawStream.h"
 
@@ -54,8 +54,8 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
   TObject(),
   fPedestalLowGain(),
   fPedestalHighGain(),
-  fSampleLowGain(),
-  fSampleHighGain(),
+  fPedestalRMSLowGain(),
+  fPedestalRMSHighGain(),
   fPeakMinusPedLowGain(),
   fPeakMinusPedHighGain(),
   fPedestalLowGainDiff(),
@@ -82,7 +82,10 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
   fRowMultiplier(0),
   fCaloString(),
   fMapping(NULL),
-  fRunNumber(-1)
+  fRunNumber(-1),
+  fSelectPedestalSamples(kTRUE), 
+  fFirstPedestalSample(0),
+  fLastPedestalSample(15)
 {
   //Default constructor. First we set the detector-type related constants.
   if (detectorType == kPhos) {
@@ -129,20 +132,20 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
                                         fColumns, 0.0, fColumns, 
                                         fRows, fRowMin, fRowMax,"s"));
     //All Samples, low gain
-    name = "hSamplelowgain";
+    name = "hPedestalRMSlowgain";
     name += i;
-    title = "All Samples, low gain, module ";
+    title = "Pedestal RMS, low gain, module ";
     title += i; 
-    fSampleLowGain.Add(new TProfile2D(name, title,
+    fPedestalRMSLowGain.Add(new TProfile2D(name, title,
                                        fColumns, 0.0, fColumns, 
                                        fRows, fRowMin, fRowMax,"s"));
   
     //All Samples, high gain
-    name = "hSamplehighgain";
+    name = "hPedestalRMShighgain";
     name += i;
-    title = "All Samples, high gain, module ";
+    title = "Pedestal RMS, high gain, module ";
     title += i; 
-    fSampleHighGain.Add(new TProfile2D(name, title,
+    fPedestalRMSHighGain.Add(new TProfile2D(name, title,
                                         fColumns, 0.0, fColumns, 
                                         fRows, fRowMin, fRowMax,"s"));
   
@@ -176,8 +179,8 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
   fPedestalLowGain.Compress();
   fPedestalHighGain.Compress();
-  fSampleLowGain.Compress();
-  fSampleHighGain.Compress();
+  fPedestalRMSLowGain.Compress();
+  fPedestalRMSHighGain.Compress();
   fPeakMinusPedLowGain.Compress();
   fPeakMinusPedHighGain.Compress();
   fDeadMap.Compress();
@@ -203,8 +206,8 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
   TObject(ped),
   fPedestalLowGain(),
   fPedestalHighGain(),
-  fSampleLowGain(),
-  fSampleHighGain(),
+  fPedestalRMSLowGain(),
+  fPedestalRMSHighGain(),
   fPeakMinusPedLowGain(),
   fPeakMinusPedHighGain(),
   fPedestalLowGainDiff(),
@@ -231,15 +234,18 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
   fRowMultiplier(ped.GetRowMultiplier()),
   fCaloString(ped.GetCaloString()),
   fMapping(NULL), //! note that we are not copying the map info
-  fRunNumber(ped.GetRunNumber())
+  fRunNumber(ped.GetRunNumber()),
+  fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
+  fFirstPedestalSample(ped.GetFirstPedestalSample()),
+  fLastPedestalSample(ped.GetLastPedestalSample())
 {
   // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
   //DS: this has not really been tested yet..
   for (int i = 0; i < fModules; i++) {
     fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
     fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
-    fSampleLowGain.Add( ped.GetSampleProfileLowGain(i) );
-    fSampleHighGain.Add( ped.GetSampleProfileHighGain(i) );
+    fPedestalRMSLowGain.Add( ped.GetPedRMSProfileLowGain(i) );
+    fPedestalRMSHighGain.Add( ped.GetPedRMSProfileHighGain(i) );
     fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
     fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
 
@@ -249,8 +255,8 @@ AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
   fPedestalLowGain.Compress();
   fPedestalHighGain.Compress();
-  fSampleLowGain.Compress();
-  fSampleHighGain.Compress();
+  fPedestalRMSLowGain.Compress();
+  fPedestalRMSHighGain.Compress();
   fPeakMinusPedLowGain.Compress();
   fPeakMinusPedHighGain.Compress();
   fDeadMap.Compress();
@@ -331,18 +337,52 @@ Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
 { 
   // Method to process=analyze one event in the data stream
   if (!in) return kFALSE; //Return right away if there's a null pointer
-  
   fNEvents++; // one more event
-  int sample, i = 0; //The sample temp, and the sample number in current event.
-  int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
+  
+  // indices for the reading
+  int sample = 0;
   int gain = 0;
+  int time = 0;
+  // counters
+  int i = 0; // the sample number in current event.
+  int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
+
+  // for the pedestal calculation
+  int sampleSum = 0; // sum of samples
+  int squaredSampleSum = 0; // sum of samples squared
+  int nSum = 0; // number of samples in sum
+  // calc. quantities
+  double mean = 0, squaredMean = 0, rms = 0;
   
   while (in->Next()) {
     sample = in->GetSignal(); //Get the adc signal
     if (sample < min) min = sample;
     if (sample > max) max = sample;
     i++;
+
+    // should we add it for the pedestal calculation?
+    time = in->GetTime();
+    if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
+        !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all 
+      sampleSum += sample;
+      squaredSampleSum += sample*sample;
+      nSum++;
+    }
+
     if ( i >= in->GetTimeLength()) {
+      // calculate pedesstal estimate: mean of possibly selected samples
+      if (nSum > 0) {
+       mean = sampleSum / (1.0 * nSum);
+       squaredMean = squaredSampleSum / (1.0 * nSum);
+       // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
+       rms = sqrt(squaredMean - mean*mean); 
+      }
+      else {
+       mean = 0;
+       squaredMean = 0;
+       rms  = 0;
+      }
+
       //If we're here then we're done with this tower
       gain = -1; // init to not valid value
       if ( in->IsLowGain() ) {
@@ -351,7 +391,7 @@ Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
       else if ( in->IsHighGain() ) {
        gain = 1;
       }
-
+      
       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
       if (arrayPos >= fModules) {
        //TODO: return an error message, if appopriate (perhaps if debug>0?)
@@ -367,18 +407,28 @@ Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
       //NOTE: coordinates are (column, row) for the profiles
       if (gain == 0) {
        //fill the low gain histograms
-       ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), min);
        ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
-       ((TProfile2D*)fSampleLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), sample);
+       if (nSum>0) { // only fill pedestal info in case it could be calculated
+         ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean); 
+         ((TProfile2D*)fPedestalRMSLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
+       }
       } 
-      else if (gain == 1) {//fill the high gain ones
-       ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), min);
+      else if (gain == 1) {
+       //fill the high gain ones
        ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
-       ((TProfile2D*)fSampleHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), sample);
+       if (nSum>0) { // only fill pedestal info in case it could be calculated
+         ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean); 
+         ((TProfile2D*)fPedestalRMSHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
+       }
       }//end if gain
 
+      // reset counters
       max = fgkSampleMin; min = fgkSampleMax;
       i = 0;
+      // also pedestal calc counters
+      sampleSum = 0; // sum of samples
+      squaredSampleSum = 0; // sum of samples squared
+      nSum = 0; // number of samples in sum
     
     }//End if end of tower
    
@@ -412,12 +462,13 @@ Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHi
     }
     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
       fPedestalHighGain[i]->Write();
+      Printf("save %d", i);
     }
-    if( ((TProfile2D *)fSampleLowGain[i])->GetEntries() || saveEmptyHistos) {
-      fSampleLowGain[i]->Write();
+    if( ((TProfile2D *)fPedestalRMSLowGain[i])->GetEntries() || saveEmptyHistos) {
+      fPedestalRMSLowGain[i]->Write();
     }
-    if( ((TProfile2D *)fSampleHighGain[i])->GetEntries() || saveEmptyHistos) {
-      fSampleHighGain[i]->Write();
+    if( ((TProfile2D *)fPedestalRMSHighGain[i])->GetEntries() || saveEmptyHistos) {
+      fPedestalRMSHighGain[i]->Write();
     }
   } 
   
index 1daf865..7f1ce7c 100644 (file)
@@ -62,22 +62,22 @@ class AliCaloCalibPedestal : public TObject {
   // Main profiles:
   TProfile2D * GetPedProfileLowGain(int i) const {return (TProfile2D*)fPedestalLowGain[i];};   // Return a pointer to the low-gain pedestal profile
   TProfile2D * GetPedProfileHighGain(int i) const {return (TProfile2D*)fPedestalHighGain[i];}; // Return a pointer to the high-gain pedestal profile
-  TProfile2D * GetSampleProfileLowGain(int i) const {return (TProfile2D*)fSampleLowGain[i];};  // Return a pointer to the low-gain profile of all samples
-  TProfile2D * GetSampleProfileHighGain(int i) const {return (TProfile2D*)fSampleHighGain[i];};        // Return a pointer to the high-gain profile of all samples
-  TProfile2D * GetPeakProfileLowGain(int i) const {return (TProfile2D*)fPeakMinusPedLowGain[i];};      // Return a pointer to the low-gain pedestal profile
-  TProfile2D * GetPeakProfileHighGain(int i) const {return (TProfile2D*)fPeakMinusPedHighGain[i];};    // Return a pointer to the high-gain pedestal profile
+  TProfile2D * GetPedRMSProfileLowGain(int i) const {return (TProfile2D*)fPedestalRMSLowGain[i];};     // Return a pointer to the low-gain rms profile 
+  TProfile2D * GetPedRMSProfileHighGain(int i) const {return (TProfile2D*)fPedestalRMSHighGain[i];};   // Return a pointer to the high-gain rms profile 
+  TProfile2D * GetPeakProfileLowGain(int i) const {return (TProfile2D*)fPeakMinusPedLowGain[i];};      // Return a pointer to the low-gain peak-pedestal profile
+  TProfile2D * GetPeakProfileHighGain(int i) const {return (TProfile2D*)fPeakMinusPedHighGain[i];};    // Return a pointer to the high-gain peak-pedestal profile
   
   // Differences to references:
   TProfile2D * GetPedProfileLowGainDiff(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPedestalLowGainDiff[i];};    // Return a pointer to the low-gain pedestal profile difference
   TProfile2D * GetPedProfileHighGainDiff(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPedestalHighGainDiff[i];};  // Return a pointer to the high-gain pedestal profile difference
-  TProfile2D * GetPeakProfileLowGainDiff(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedLowGainDiff[i];};       // Return a pointer to the low-gain pedestal profile difference
-  TProfile2D * GetPeakProfileHighGainDiff(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedHighGainDiff[i];};     // Return a pointer to the high-gain pedestal profile difference
+  TProfile2D * GetPeakProfileLowGainDiff(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedLowGainDiff[i];};       // Return a pointer to the low-gain peak-pedestal profile difference
+  TProfile2D * GetPeakProfileHighGainDiff(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedHighGainDiff[i];};     // Return a pointer to the high-gain peak-pedestal profile difference
   
   // Ratio to references:
   TProfile2D * GetPedProfileLowGainRatio(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPedestalLowGainRatio[i];};  // Return a pointer to the low-gain pedestal profile ratio
   TProfile2D * GetPedProfileHighGainRatio(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPedestalHighGainRatio[i];};        // Return a pointer to the high-gain pedestal profile ratio
-  TProfile2D * GetPeakProfileLowGainRatio(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedLowGainRatio[i];};     // Return a pointer to the low-gain pedestal profile ratio
-  TProfile2D * GetPeakProfileHighGainRatio(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedHighGainRatio[i];};   // Return a pointer to the high-gain pedestal profile ratio
+  TProfile2D * GetPeakProfileLowGainRatio(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedLowGainRatio[i];};     // Return a pointer to the low-gain peak-pedestal profile ratio
+  TProfile2D * GetPeakProfileHighGainRatio(int i){ValidateComparisonProfiles(); return (TProfile2D*)fPeakMinusPedHighGainRatio[i];};   // Return a pointer to the high-gain peak-pedestal profile ratio
   
   TH2D * GetDeadMap(int i) const {return (TH2D*)fDeadMap[i];};
 
@@ -97,6 +97,14 @@ class AliCaloCalibPedestal : public TObject {
   int GetRunNumber() const {return fRunNumber;};
   int GetRefRunNumber() const {if (fReference) return fReference->GetRunNumber(); else return -1;};
 
+  // Possibility to select only some samples for the pedestal calculation
+  void SetSelectPedestalSamples(Bool_t flag = kFALSE) {fSelectPedestalSamples = flag;} // select to to use only some range of samples for pedestal calc.
+  Bool_t GetSelectPedestalSamples() const {return fSelectPedestalSamples;} // select to to use only some range of samples for pedestal calc.
+  void SetFirstPedestalSample(int i) {fFirstPedestalSample = i;} // first sample to use
+  void SetLastPedestalSample(int i) {fLastPedestalSample = i;} // last sample to use
+  int GetFirstPedestalSample() const {return fFirstPedestalSample;}; // first sample to use
+  int GetLastPedestalSample() const {return fLastPedestalSample;}; // last sample to use
+
   // Basic counters
   int GetNEvents() const {return fNEvents;};
   int GetNChanFills() const {return fNChanFills;};
@@ -118,7 +126,6 @@ class AliCaloCalibPedestal : public TObject {
   AliCaloCalibPedestal * GetReference() const {return fReference;}; //Get the reference object. Needed for debug, will probably be removed later
   void ComputeDeadTowers(int threshold = 5, const char * deadMapFile = 0);//Computes the dead tower values
   
-  
   //Saving functions
   Bool_t SaveHistograms(TString fileName, Bool_t saveEmptyHistos = kFALSE); //Saves the histograms to a .root file
   
@@ -131,8 +138,8 @@ class AliCaloCalibPedestal : public TObject {
   //since we have only around 12 objects (maximum) in the array anyway.
   TObjArray fPedestalLowGain; // pedestal info for low gain
   TObjArray fPedestalHighGain; // pedestal info for high gain
-  TObjArray fSampleLowGain; // all sample info for low gain
-  TObjArray fSampleHighGain; // all sample info for high gain
+  TObjArray fPedestalRMSLowGain; // pedestal rms info for low gain
+  TObjArray fPedestalRMSHighGain; // pedestal rms info for high gain
   TObjArray fPeakMinusPedLowGain; // (peak-pedestal) info for low gain
   TObjArray fPeakMinusPedHighGain; // (peak-pedestal) info for high gain
   
@@ -171,6 +178,9 @@ class AliCaloCalibPedestal : public TObject {
   TString fCaloString; // id for which detector type we have 
   AliCaloAltroMapping **fMapping;    //! Altro Mapping object
   int fRunNumber; //The run number. Needs to be set by the user.
+  Bool_t fSelectPedestalSamples; // select to to use only some range of samples for pedestal calc.
+  int fFirstPedestalSample; // first sample to use
+  int fLastPedestalSample; // last sample to use
 
   //Constants needed by the class
   static const int fgkSampleMax = 1023; // highest possible sample value (10-bit = 0x3ff)
@@ -184,7 +194,7 @@ class AliCaloCalibPedestal : public TObject {
   static const int fgkEmCalCols = 48; // number of columns per module for EMCAL
   static const int fgkEmCalModules = 12; // number of modules for EMCAL
   
-  ClassDef(AliCaloCalibPedestal,2)
+  ClassDef(AliCaloCalibPedestal,3)
 
 };