Modified to be used in trunk (Cynthia H.)

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

diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx
index ebeb4c357c8eb7316523ac1118c9c36ad32da030..afd87404a9e9a58d4eb50617efb1e2c0744c00c9 100644 (file)
--- a/EMCAL/AliEMCALRawUtils.cxx
+++ b/EMCAL/AliEMCALRawUtils.cxx
@@ -56,6 +56,7 @@
 #include "TSystem.h"
 
 #include "AliLog.h"
+#include "AliRun.h"
 #include "AliRunLoader.h"
 #include "AliCaloAltroMapping.h"
 #include "AliAltroBuffer.h"
@@ -68,7 +69,7 @@
 #include "AliEMCALGeometry.h"
 #include "AliEMCALDigitizer.h"
 #include "AliEMCALDigit.h"
-
+#include "AliEMCAL.h"
 
 ClassImp(AliEMCALRawUtils)
 
@@ -79,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),
@@ -101,12 +104,18 @@ AliEMCALRawUtils::AliEMCALRawUtils()
     fMapping[i] = (AliAltroMapping*)maps->At(i);
   }
 
-  fGeom = AliEMCALGeometry::GetInstance();
-  if(!fGeom) {
-    fGeom = AliEMCALGeometry::GetInstance("","");
-    if(!fGeom) AliFatal(Form("Could not get geometry!!"));
+  //To make sure we match with the geometry in a simulation file,
+  //let's try to get it first.  If not, take the default geometry
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) {
+    fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  } else {
+    AliInfo(Form("Using default geometry in raw reco"));
+    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
   }
 
+  if(!fGeom) AliFatal(Form("Could not get geometry!"));
+
 }
 
 //____________________________________________________________________________
@@ -267,13 +276,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
@@ -329,6 +331,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++;
@@ -336,9 +339,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);
     }
        
@@ -427,6 +432,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;
@@ -456,6 +463,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 
@@ -515,7 +536,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
@@ -529,11 +549,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 ;