added cutoff parameter in z direction to ignore clusters of large z (Gaute)

[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRawUtils.cxx

diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx
index 349148f6cbb5c137ef323f2d48b6a86f63cdc58d..1eb5826c9df1349d160a156bb7959d85b4633fa1 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.cxx
+++ b/EMCAL/AliEMCALRawUtils.cxx
@@ -69,7 +69,7 @@
 #include "AliEMCALGeometry.h"
 #include "AliEMCALDigitizer.h"
 #include "AliEMCALDigit.h"
-
+#include "AliEMCAL.h"
 
 ClassImp(AliEMCALRawUtils)
 
@@ -80,6 +80,8 @@ Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ;   // 15 time bins ~ 1.5 muse
 // some digitization constants
 Int_t    AliEMCALRawUtils::fgThreshold = 1;
 Int_t    AliEMCALRawUtils::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
+Int_t    AliEMCALRawUtils::fgPedestalValue = 32;      // pedestal value for digits2raw
+Double_t AliEMCALRawUtils::fgFEENoise = 3.;            // 3 ADC channels of noise (sampled)
 
 AliEMCALRawUtils::AliEMCALRawUtils()
   : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
@@ -116,6 +118,38 @@ AliEMCALRawUtils::AliEMCALRawUtils()
 
 }
 
+//____________________________________________________________________________
+AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry)
+  : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
+    fNPedSamples(0), fGeom(pGeometry), fOption("")
+{
+  //
+  // Initialize with the given geometry - constructor required by HLT
+  // HLT does not use/support AliRunLoader(s) instances
+  // This is a minimum intervention solution
+  // Comment by MPloskon@lbl.gov
+  //
+
+  //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();
+  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
+
+  for(Int_t i = 0; i < 2; i++) {
+    fMapping[i] = (AliAltroMapping*)maps->At(i);
+  }
+
+  if(!fGeom) AliFatal(Form("Could not get geometry!"));
+
+}
+
 //____________________________________________________________________________
 AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU)
   : TObject(),
@@ -274,13 +308,6 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
   // Select EMCAL DDL's;
   reader->Select("EMCAL");
 
-  TString option = GetOption();
-  if (option.Contains("OldRCUFormat"))
-    in.SetOldRCUFormat(kTRUE); // Needed for testbeam data
-  else
-    in.SetOldRCUFormat(kFALSE);
-
-
   //Updated fitting routine from 2007 beam test takes into account
   //possibility of two peaks in data and selects first one for fitting
   //Also sets some of the starting parameters based on the shape of the
@@ -336,6 +363,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
       row = in.GetRow();
       
       gSig->SetPoint(in.GetTime(), in.GetTime(), in.GetSignal()) ;
+
       if (in.GetTime() > maxTime)
         maxTime = in.GetTime();
       iTime++;
@@ -343,9 +371,11 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr)
 
     FitRaw(gSig, signalF, amp, time) ; 
     
-    if (amp > 0) {
+    if (amp > 0 && amp < 10000) {  //check both high and low end of
+                                  //result, 10000 is somewhat arbitrary
       AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
       //cout << "col " << col-40 << " row " << row-8 << " lowGain " << lowGain << " amp " << amp << endl;
+
       AddDigit(digitsArr, id, lowGain, (Int_t)amp, time);
     }
        
@@ -434,6 +464,8 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
   Float_t min_after_sig = 9999;
   Int_t tmin_after_sig = gSig->GetN();
   Int_t n_ped_after_sig = 0;
+  Int_t plateau_width = 0;
+  Int_t plateau_start = 9999;
 
   for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) {
     Double_t ttime, signal;
@@ -463,6 +495,20 @@ void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_
         break;
       }
     }
+    //Add check on plateau
+    if (signal >= fgkRawSignalOverflow - fNoiseThreshold) {
+      if(plateau_width == 0) plateau_start = i;
+      plateau_width++;
+    }
+  }
+
+  if(plateau_width > 0) {
+    for(int j = 0; j < plateau_width-2; j++) {
+      //Note, have to remove the same point N times because after each
+      //remove, the positions of all subsequent points have shifted down
+      //We leave the first and last as anchor points
+      gSig->RemovePoint(plateau_start+1);
+    }
   }
 
   if ( max - ped > fNoiseThreshold ) { // else its noise 
@@ -522,7 +568,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 pedVal = 32;
   Bool_t lowGain = kFALSE ; 
 
   // A:   par[0]   // Amplitude = peak value
@@ -536,11 +581,17 @@ const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
   signalF.SetParameter(1, dtime + fgTimeTrigger) ; 
   signalF.SetParameter(2, fTau) ; 
   signalF.SetParameter(3, fOrder);
-  signalF.SetParameter(4, pedVal);
+  signalF.SetParameter(4, fgPedestalValue);
 
   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
     Double_t time = iTime * GetRawFormatTimeBinWidth() ;
     Double_t signal = signalF.Eval(time) ;     
+
+    //According to Terry Awes, 13-Apr-2008
+    //add gaussian noise in quadrature to each sample
+    //Double_t noise = gRandom->Gaus(0.,fgFEENoise);
+    //signal = sqrt(signal*signal + noise*noise);
+
     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
     if ( adcH[iTime] > fgkRawSignalOverflow ){  // larger than 10 bits 
       adcH[iTime] = fgkRawSignalOverflow ;