incorporated raw signal fit parameters into AliEMCALRecParam

diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx
index ebf9031f9fde1f79c99d7976d118d16a937dbc8b..ebeb4c357c8eb7316523ac1118c9c36ad32da030 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.cxx
+++ b/EMCAL/AliEMCALRawUtils.cxx
@@ -73,9 +73,7 @@
 ClassImp(AliEMCALRawUtils)
 
 // Signal shape parameters
-Int_t    AliEMCALRawUtils::fgOrder         = 2 ;      // Order of gamma function 
 Double_t AliEMCALRawUtils::fgTimeBinWidth  = 100E-9 ; // each sample is 100 ns
-Double_t AliEMCALRawUtils::fgTau         = 235E-9 ;   // 235 ns (from CERN testbeam; not very accurate)
 Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 musec
 
 // some digitization constants
@@ -83,10 +81,17 @@ Int_t    AliEMCALRawUtils::fgThreshold = 1;
 Int_t    AliEMCALRawUtils::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
 
 AliEMCALRawUtils::AliEMCALRawUtils()
-  : fHighLowGainFactor(0.), fGeom(0), 
-    fOption("")
+  : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
+    fNPedSamples(0), fGeom(0), fOption("")
 {
+
+  //These are default parameters.  
+  //Can be re-set from without with setter functions
   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits) 
+  fOrder = 2;                         // order of gamma fn
+  fTau = 2.35;                        // in units of timebin, from CERN 2007 testbeam
+  fNoiseThreshold = 3;
+  fNPedSamples = 5;
 
   //Get Mapping RCU files from the AliEMCALRecParam                                 
   const TObjArray* maps = AliEMCALRecParam::GetMappings();
@@ -96,7 +101,6 @@ AliEMCALRawUtils::AliEMCALRawUtils()
     fMapping[i] = (AliAltroMapping*)maps->At(i);
   }
 
-
   fGeom = AliEMCALGeometry::GetInstance();
   if(!fGeom) {
     fGeom = AliEMCALGeometry::GetInstance("","");
@@ -109,6 +113,10 @@ AliEMCALRawUtils::AliEMCALRawUtils()
 AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU)
   : TObject(),
     fHighLowGainFactor(rawU.fHighLowGainFactor), 
+    fOrder(rawU.fOrder),
+    fTau(rawU.fTau),
+    fNoiseThreshold(rawU.fNoiseThreshold),
+    fNPedSamples(rawU.fNPedSamples),
     fGeom(rawU.fGeom), 
     fOption(rawU.fOption)
 {
@@ -124,6 +132,10 @@ AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU)
 
   if(this != &rawU) {
     fHighLowGainFactor = rawU.fHighLowGainFactor;
+    fOrder = rawU.fOrder;
+    fTau = rawU.fTau;
+    fNoiseThreshold = rawU.fNoiseThreshold;
+    fNPedSamples = rawU.fNPedSamples;
     fGeom = rawU.fGeom;
     fOption = rawU.fOption;
     fMapping[0] = rawU.fMapping[0];
@@ -268,11 +280,11 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
   //given raw signal being fit
 
   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5);
-  signalF->SetParameters(10.,0.,2.35,2.,5.); //set all defaults once, just to be safe
+  signalF->SetParameters(10.,0.,fTau,fOrder,5.); //set all defaults once, just to be safe
   signalF->SetParNames("amp","t0","tau","N","ped");
-  signalF->SetParameter(2,2.35); // tau in units of time bin
+  signalF->SetParameter(2,fTau); // tau in units of time bin
   signalF->SetParLimits(2,2,-1);
-  signalF->SetParameter(3,2); // order
+  signalF->SetParameter(3,fOrder); // order
   signalF->SetParLimits(3,2,-1);
   
   Int_t id =  -1;
@@ -335,7 +347,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
       gSig->SetPoint(index, index, 0) ;  
     } 
     // Reset starting parameters for fit function
-    signalF->SetParameters(10.,0.,2.35,2.,5.); //reset all defaults just to be safe
+    signalF->SetParameters(10.,0.,fTau,fOrder,5.); //reset all defaults just to be safe
 
   }; // EMCAL entries loop
   
@@ -369,14 +381,14 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain
     new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ;  
   }
   else { // a digit already exists, check range 
-         // (use high gain if signal < 800, otherwise low gain)
+         // (use high gain if signal < cut value, otherwise low gain)
     if (lowGain) { // new digit is low gain
-      if (digit->GetAmp() > 800) {  // use if stored digit is out of range
+      if (digit->GetAmp() > fgkOverflowCut) {  // use if stored digit is out of range
        digit->SetAmp(Int_t(fHighLowGainFactor * amp));
        digit->SetTime(time);
       }
     }
-    else if (amp < 800) { // new digit is high gain; use if not out of range
+    else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range
       digit->SetAmp(amp);
       digit->SetTime(time);
     }
@@ -388,13 +400,11 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
 {
   // Fits the raw signal time distribution; from AliEMCALGetter 
 
-  const Int_t kNoiseThreshold = 5;
-  const Int_t kNPedSamples = 5;
   amp = time = 0. ; 
   Double_t ped = 0;
   Int_t nPed = 0;
 
-  for (Int_t index = 0; index < kNPedSamples; index++) {
+  for (Int_t index = 0; index < fNPedSamples; index++) {
     Double_t ttime, signal;
     gSig->GetPoint(index, ttime, signal) ; 
     if (signal > 0) {
@@ -418,14 +428,14 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
   Int_t tmin_after_sig = gSig->GetN();
   Int_t n_ped_after_sig = 0;
 
-  for (Int_t i=kNPedSamples; i < gSig->GetN(); i++) {
+  for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) {
     Double_t ttime, signal;
     gSig->GetPoint(i, ttime, signal) ; 
     if (!max_found && signal > max) {
       i_max = i;
       max = signal;
     }
-    else if ( max > ped + kNoiseThreshold ) {
+    else if ( max > ped + fNoiseThreshold ) {
       max_found = 1;
       min_after_sig = signal;
       tmin_after_sig = i;
@@ -439,7 +449,7 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
         max_fit = tmin_after_sig;
         break;
       }
-      if ( signal < ped + kNoiseThreshold)
+      if ( signal < ped + fNoiseThreshold)
         n_ped_after_sig++;
       if (n_ped_after_sig >= 5) {  // include 5 pedestal bins after peak
         max_fit = i;
@@ -448,7 +458,7 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
     }
   }
 
-  if ( max - ped > kNoiseThreshold ) { // else its noise 
+  if ( max - ped > fNoiseThreshold ) { // else its noise 
     AliDebug(2,Form("Fitting max %d ped %d", max, ped));
     signalF->SetRange(0,max_fit);
 
@@ -505,7 +515,6 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
   // for a start time dtime and an amplitude damp given by digit, 
   // calculates the raw sampled response AliEMCAL::RawResponseFunction
 
-  const Int_t kRawSignalOverflow = 0x3FF ; 
   const Int_t pedVal = 32;
   Bool_t lowGain = kFALSE ; 
 
@@ -518,24 +527,24 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 5);
   signalF.SetParameter(0, damp) ; 
   signalF.SetParameter(1, dtime + fgTimeTrigger) ; 
-  signalF.SetParameter(2, fgTau) ; 
-  signalF.SetParameter(3, fgOrder);
+  signalF.SetParameter(2, fTau) ; 
+  signalF.SetParameter(3, fOrder);
   signalF.SetParameter(4, pedVal);
 
   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
     Double_t time = iTime * GetRawFormatTimeBinWidth() ;
     Double_t signal = signalF.Eval(time) ;     
     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
-    if ( adcH[iTime] > kRawSignalOverflow ){  // larger than 10 bits 
-      adcH[iTime] = kRawSignalOverflow ;
+    if ( adcH[iTime] > fgkRawSignalOverflow ){  // larger than 10 bits 
+      adcH[iTime] = fgkRawSignalOverflow ;
       lowGain = kTRUE ; 
     }
 
     signal /= fHighLowGainFactor;
 
     adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
-    if ( adcL[iTime] > kRawSignalOverflow)  // larger than 10 bits 
-      adcL[iTime] = kRawSignalOverflow ;
+    if ( adcL[iTime] > fgkRawSignalOverflow)  // larger than 10 bits 
+      adcL[iTime] = fgkRawSignalOverflow ;
   }
   return lowGain ; 
 }

diff --git a/EMCAL/AliEMCALRawUtils.h b/EMCAL/AliEMCALRawUtils.h
index 1a13a34934fffb857b125a3ad05535e68c86d07c..cdca0d0d5b79f6fb1efadc72ab8a8e963808d57d 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.h
+++ b/EMCAL/AliEMCALRawUtils.h
@@ -53,12 +53,21 @@ class AliEMCALRawUtils : public TObject {
   void AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time);
 
   // Signal shape parameters
-  Double_t GetRawFormatHighLowGainFactor() const {return fHighLowGainFactor ;}
-  static Int_t GetRawFormatOrder() { return fgOrder ; }   
+  Double_t GetRawFormatHighLowGainFactor() const { return fHighLowGainFactor ;}
+  Int_t GetRawFormatOrder()                const { return fOrder ; }   
+  Double_t GetRawFormatTau()               const { return fTau ; }    
+  Int_t GetNoiseThreshold()                const { return fNoiseThreshold; }
+  Int_t GetNPedSamples()                   const { return fNPedSamples; }
+
+  void SetRawFormatHighLowGainFactor(Double_t val) {fHighLowGainFactor=val;}
+  void SetRawFormatOrder(Int_t val)                {fOrder=val; }   
+  void SetRawFormatTau(Double_t val)               {fTau=val; }    
+  void SetNoiseThreshold(Int_t val)                {fNoiseThreshold=val; }
+  void SetNPedSamples(Int_t val)                   {fNPedSamples=val; }
+  
   static Int_t GetRawFormatTimeBins() { return fgkTimeBins ; }    
   static Double_t GetRawFormatTimeMax() { return fgkTimeBins*fgTimeBinWidth; }   
   static Double_t GetRawFormatTimeBinWidth() { return fgTimeBinWidth; }   
-  Double_t GetRawFormatTau() const { return fgTau ; }    
   Double_t GetRawFormatTimeTrigger() const { return fgTimeTrigger ; }
   Int_t GetRawFormatThreshold() const { return fgThreshold ; }       
   Int_t GetRawFormatDDLPerSuperModule() const { return fgDDLPerSuperModule ; } 
@@ -71,15 +80,19 @@ class AliEMCALRawUtils : public TObject {
   static Double_t RawResponseFunction(Double_t *x, Double_t *par); 
   Bool_t   RawSampledResponse(Double_t dtime, Double_t damp, Int_t * adcH, Int_t * adcL) const;  
 
-  ClassDef(AliEMCALRawUtils,1)
 
  private:
   Double_t fHighLowGainFactor ;         // high to low gain factor for the raw RO signal
-  static Int_t fgOrder ;                // order of the gamma function for the RO signal
-  static Double_t fgTau ;                // tau parameter of gamma function for the RO signal
-  static Double_t fgTimeTrigger ;       // time of the trigger for the RO signal 
+  Int_t fOrder ;                        // order of the gamma function for the RO signal
+  Double_t fTau ;                       // tau parameter of gamma function for the RO signal
+  Int_t fNoiseThreshold;                // threshold to consider signal or noise
+  Int_t fNPedSamples;                   // number of samples to use in pedestal calculation
 
+  static const Int_t fgkOverflowCut = 950;  // cut to discriminate overflowed channels
   static const Int_t fgkTimeBins = 256 ; // number of sampling bins of the raw RO signal  
+  static const Int_t fgkRawSignalOverflow = 0x3FF; // maximum signal (10 bits)
+
+  static Double_t fgTimeTrigger ;       // time of the trigger for the RO signal 
   static Double_t fgTimeBinWidth;       // maximum sampled time of the raw RO signal                             
   static Int_t fgThreshold;             // threshold
   static Int_t fgDDLPerSuperModule;     // number of DDL per SuperModule
@@ -88,6 +101,8 @@ class AliEMCALRawUtils : public TObject {
   AliAltroMapping*  fMapping[2];   //only two for now
 
   TString fOption;                      //! option passed from Reconstructor
+
+  ClassDef(AliEMCALRawUtils,2)          // utilities for raw signal fitting
 };
 
 #endif

diff --git a/EMCAL/AliEMCALRecParam.cxx b/EMCAL/AliEMCALRecParam.cxx
index 1e19be3b9f1796512c9aefdcbab6318f425c15a0..66c106fce2b869752f7fbed1ff0bca7ec4dad00a 100644 (file)
--- a/EMCAL/AliEMCALRecParam.cxx
+++ b/EMCAL/AliEMCALRecParam.cxx
@@ -35,7 +35,8 @@ TObjArray* AliEMCALRecParam::fgkMaps =0; //ALTRO mappings
 AliEMCALRecParam::AliEMCALRecParam():
   fClusteringThreshold(0.5),fW0(4.5),fMinECut(0.45), //clustering
   fTrkCutX(6.0), fTrkCutY(6.0), fTrkCutZ(6.0),  fTrkCutR(10.0),//track matching
-  fTrkCutAlphaMin(-50.0), fTrkCutAlphaMax(50.0), fTrkCutAngle(10000.0) //track matching
+  fTrkCutAlphaMin(-50.0), fTrkCutAlphaMax(50.0), fTrkCutAngle(10000.0), //track matching
+  fHighLowGainFactor(16.0), fOrderParameter(2), fTau(2.35), fNoiseThreshold(3), fNPedSamples(5) //raw signal
 {
   // default reco values
 
@@ -147,10 +148,6 @@ AliEMCALRecParam::AliEMCALRecParam():
   fPiZero10to60[5][0] =  0.002942;
   fPiZero10to60[5][1] = -3.976e-05;
 
-  //track matching
-
-
-
 }
 
 //-----------------------------------------------------------------------------
@@ -203,6 +200,9 @@ void AliEMCALRecParam::Print(Option_t *) const
 
   printf("\n");
 
+  AliInfo(Form("Raw signal parameters: \n gain factor=%f, order=%d, tau=%f, noise threshold=%d, nped samples=%d \n",
+              fHighLowGainFactor,fOrderParameter,fTau,fNoiseThreshold,fNPedSamples));
+
 }
 
 //-----------------------------------------------------------------------------                           

diff --git a/EMCAL/AliEMCALRecParam.h b/EMCAL/AliEMCALRecParam.h
index 1c8d0cf5abaaaa01cad7554eae97428b0f4e00f7..8d21df24c6e076c2a529ab5ef55ace3f21e4a801 100644 (file)
--- a/EMCAL/AliEMCALRecParam.h
+++ b/EMCAL/AliEMCALRecParam.h
@@ -8,7 +8,10 @@
 //-----------------------------------------------------------------------------
 // Container of EMCAL reconstruction parameters
 // The purpose of this object is to store it to OCDB
-// and retrieve it in AliEMCALClusterizerv1, AliEMCALPID and AliEMCALTracker
+// and retrieve it in AliEMCALClusterizerv1, AliEMCALPID,
+// AliEMCALTracker and use it to configure AliEMCALRawUtils
+// 
+// 
 // Author: Yuri Kharlov
 //-----------------------------------------------------------------------------
 
@@ -60,6 +63,20 @@ public:
   Double_t GetTrkCutAlphaMax() const {return fTrkCutAlphaMax;}
   Double_t GetTrkCutAngle() const    {return fTrkCutAngle;}
 
+  //Raw signal fitting (Jenn)
+  /* raw signal setters */
+  void SetHighLowGainFactor(Double_t value) {fHighLowGainFactor = value;}
+  void SetOrderParameter(Int_t value)       {fOrderParameter = value;}
+  void SetTau(Double_t value)               {fTau = value;}
+  void SetNoiseThreshold(Int_t value)       {fNoiseThreshold = value;}
+  void SetNPedSamples(Int_t value)          {fNPedSamples = value;}
+  /* raw signal getters */
+  Double_t GetHighLowGainFactor() const {return fHighLowGainFactor;}
+  Int_t    GetOrderParameter()    const {return fOrderParameter;}
+  Double_t GetTau()               const {return fTau;}
+  Int_t    GetNoiseThreshold()    const {return fNoiseThreshold;}
+  Int_t    GetNPedSamples()       const {return fNPedSamples;}
+ 
   virtual void Print(Option_t * option="") const ;
 
   static const  TObjArray* GetMappings();
@@ -85,9 +102,16 @@ private:
   Double_t  fTrkCutAlphaMax;       // cut on 'alpha' parameter for track matching (min)
   Double_t  fTrkCutAngle;          // cut on relative angle between different track points for track matching
 
+  //Raw signal fitting parameters (Jenn)
+  Double_t fHighLowGainFactor;     //gain factor to convert between high and low gain
+  Int_t    fOrderParameter;        //order parameter for raw signal fit
+  Double_t fTau;                   //decay constant for raw signal fit
+  Int_t    fNoiseThreshold;        //threshold to consider signal or noise
+  Int_t    fNPedSamples;           //number of time samples to use in pedestal calculation
+
   static TObjArray* fgkMaps;       // ALTRO mappings for RCU0..RCUX
 
-  ClassDef(AliEMCALRecParam,3)   // Reconstruction parameters
+  ClassDef(AliEMCALRecParam,4)   // Reconstruction parameters
 
 } ;
 

diff --git a/EMCAL/AliEMCALReconstructor.cxx b/EMCAL/AliEMCALReconstructor.cxx
index 6f65e85b5adf8045df840922989f5fd25d25150a..5ac1668309c92966b750da6198da4efee8cd5ebc 100644 (file)
--- a/EMCAL/AliEMCALReconstructor.cxx
+++ b/EMCAL/AliEMCALReconstructor.cxx
@@ -139,6 +139,13 @@ void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
 
   //must be done here because, in constructor, option is not yet known
   fgRawUtils->SetOption(GetOption());
+
+  fgRawUtils->SetRawFormatHighLowGainFactor(fgkRecParam->GetHighLowGainFactor());
+  fgRawUtils->SetRawFormatOrder(fgkRecParam->GetOrderParameter());
+  fgRawUtils->SetRawFormatTau(fgkRecParam->GetTau());
+  fgRawUtils->SetNoiseThreshold(fgkRecParam->GetNoiseThreshold());
+  fgRawUtils->SetNPedSamples(fgkRecParam->GetNPedSamples());
+
   fgRawUtils->Raw2Digits(rawReader,digitsArr);
 
   digitsTree->Fill();

diff --git a/EMCAL/Config/RecParam/Run0_999999_v0_s0.root b/EMCAL/Config/RecParam/Run0_999999_v0_s0.root
index 08f27cc6ac05eef689ccb02ca92d046cfa2d4206..22b1ad86aaa18cdcf5755a57b891c2789f6239b6 100644 (file)
Binary files a/EMCAL/Config/RecParam/Run0_999999_v0_s0.root and b/EMCAL/Config/RecParam/Run0_999999_v0_s0.root differ

diff --git a/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C b/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C
index e913be3c0dd2766bab987d93724738471f0a866a..4267c1567f73ad9d9a77e4e0193b2323d8cebf09 100644 (file)
--- a/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C
+++ b/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C
@@ -176,6 +176,13 @@ void SetRecParam()
   
   recParamDB->SetPiZero10to60(5,0,  0.002942);
   recParamDB->SetPiZero10to60(5,1, -3.976e-05);
+
+  // raw signal fitting
+  recParamDB->SetHighLowGainFactor(16.);
+  recParamDB->SetOrderParameter(2);
+  recParamDB->SetTau(2.35);
+  recParamDB->SetNoiseThreshold(3);
+  recParamDB->SetNPedSamples(5);
   
   // Store calibration data into database