Cloned PHOS for the raw data format handling
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 28 Aug 2004 08:17:27 +0000 (08:17 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 28 Aug 2004 08:17:27 +0000 (08:17 +0000)
EMCAL/AliEMCAL.cxx
EMCAL/AliEMCAL.h
EMCAL/AliEMCALGetter.cxx
EMCAL/AliEMCALGetter.h
EMCAL/AliEMCALRawStream.h

index 1e95f60..74be313 100644 (file)
@@ -46,6 +46,12 @@ class TFile;
 #include "AliAltroBuffer.h"
 
 ClassImp(AliEMCAL)
+Double_t AliEMCAL::fgCapa        = 1.;        // 1pF 
+Int_t    AliEMCAL::fgOrder       = 2 ;
+Double_t AliEMCAL::fgTimeMax     = 2.56E-5 ;  // each sample is over 100 ns fTimeMax/fTimeBins
+Double_t AliEMCAL::fgTimePeak    = 4.1E-6 ;   // 4 micro seconds
+Double_t AliEMCAL::fgTimeTrigger = 100E-9 ;      // 100ns, just for a reference
 //____________________________________________________________________________
 AliEMCAL::AliEMCAL():AliDetector()
 {
@@ -58,11 +64,10 @@ AliEMCAL::AliEMCAL(const char* name, const char* title): AliDetector(name,title)
 {
   //   ctor : title is used to identify the layout
 
-  fTimeMax  = 1.28E-5 ; 
-  fTimePeak = 2.0E-6 ;
-  fTimeRes = 1.5E-6 ; 
-  fHighGainFactor = 40;
-  fHighGainOffset = 0x200 ; 
+  fHighCharge        = 8.2 ;          // adjusted for a high gain range of 5.12 GeV (10 bits)
+  fHighGain          = 6.64 ; 
+  fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits) 
+  fLowGainOffset     = 1 ;            // offset added to the module id to distinguish high and low gain data
 }
 
 //____________________________________________________________________________
@@ -75,6 +80,10 @@ AliEMCAL::~AliEMCAL()
 void AliEMCAL::Copy(AliEMCAL & emcal)
 {
   TObject::Copy(emcal) ; 
+  emcal.fHighCharge        = fHighCharge ;
+  emcal.fHighGain          = fHighGain ; 
+  emcal.fHighLowGainFactor = fHighLowGainFactor ;  
+  emcal.fLowGainOffset     = fLowGainOffset;   
 }
 
 //____________________________________________________________________________
@@ -172,24 +181,24 @@ void AliEMCAL::CreateMaterials()
 //____________________________________________________________________________
 void AliEMCAL::Digits2Raw()
 {
-// convert digits of the current event to raw data
+  // convert digits of the current event to raw data
+  AliEMCALLoader * loader = dynamic_cast<AliEMCALLoader*>(fLoader) ; 
 
   // get the digits
-  AliEMCALGetter * gime = AliEMCALGetter::Instance(AliRunLoader::GetGAliceName()) ; 
-  if (!gime) {
-    Error("Digits2Raw", "EMCAL Getter not instantiated") ;
-    return ; 
-  }
-  gime->Event(gime->EventNumber(), "D") ; 
-  TClonesArray* digits = gime->Digits() ;
+  loader->LoadDigits();
+  TClonesArray* digits = loader->Digits() ;
 
   if (!digits) {
     Error("Digits2Raw", "no digits found !");
     return;
   }
 
+  // get the digitizer 
+  loader->LoadDigitizer();
+  AliEMCALDigitizer * digitizer = dynamic_cast<AliEMCALDigitizer *>(loader->Digitizer())  ; 
+  
   // get the geometry
-  AliEMCALGeometry* geom = gime->EMCALGeometry();
+  AliEMCALGeometry* geom = GetGeometry();
   if (!geom) {
     Error("Digits2Raw", "no geometry found !");
     return;
@@ -197,17 +206,16 @@ void AliEMCAL::Digits2Raw()
 
   // some digitization constants
   const Int_t    kDDLOffset = 0x800;
-  const Int_t    kThreshold = 3;
+  const Int_t    kThreshold = 1;
   const Int_t    kChannelsperDDL = 897 ; 
-
   AliAltroBuffer* buffer = NULL;
   Int_t prevDDL = -1;
   Int_t adcValuesLow[fkTimeBins];
   Int_t adcValuesHigh[fkTimeBins];
-
+  
   // loop over digits (assume ordered digits)
   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
-    AliEMCALDigit* digit = gime->Digit(iDigit);
+    AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
     if (digit->GetAmp() < kThreshold) 
       continue;
     Int_t iDDL = digit->GetId() / kChannelsperDDL ;
@@ -233,32 +241,37 @@ void AliEMCAL::Digits2Raw()
     }
 
     // out of time range signal (?)
-    if (digit->GetTimeR() > fTimeMax) {
+    if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
       buffer->FillBuffer(digit->GetAmp());
-      buffer->FillBuffer(fkTimeBins);  // time bin
+      buffer->FillBuffer(GetRawFormatTimeBins() );  // time bin
       buffer->FillBuffer(3);          // bunch length
       buffer->WriteTrailer(3, idDDL, 0, 0);  // trailer
-      
-    // simulate linear rise and gaussian decay of signal
+
+      // calculate the time response function
     } else {
-      Bool_t highGain = kFALSE;
-
-      // write low and eventually high gain channel
-      buffer->WriteChannel(iDDL, 0, 0, 
-                          fkTimeBins, adcValuesLow, kThreshold);
-      if (highGain) {
-       buffer->WriteChannel(iDDL, 0, fHighGainOffset, 
-                            fkTimeBins, adcValuesHigh, 1);
-      }
+      Double_t energy = 0 ;  
+      energy = digit->GetAmp() * digitizer->GetECAchannel() + digitizer->GetECApedestal() ; 
+      
+      Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; 
+      
+      if (lowgain) 
+       buffer->WriteChannel(iDDL, 0, fLowGainOffset, 
+                            GetRawFormatTimeBins(), adcValuesLow, kThreshold);
+      else 
+       buffer->WriteChannel(iDDL, 0, 0, 
+                            GetRawFormatTimeBins(), adcValuesHigh, kThreshold);
+      
     }
   }
+  
   // write real header and close last file
   if (buffer) {
     buffer->Flush();
     buffer->WriteDataHeader(kFALSE, kFALSE);
     delete buffer;
   }
-  gime->EmcalLoader()->UnloadDigits();
+
+  loader->UnloadDigits();
 }
 
 //____________________________________________________________________________
@@ -285,50 +298,67 @@ AliLoader* AliEMCAL::MakeLoader(const char* topfoldername)
 Double_t AliEMCAL::RawResponseFunction(Double_t *x, Double_t *par)
 {
   // Shape of the electronics raw reponse:
-  // 1. the signal rises linearly from par[4] to par[1] to reach the maximu par[3]
-  // 2. the signal decays with a gaussian shape for par[4]+par[1] with a sigma of par[2]
-  
-  Float_t xx = x[0] ; 
-  Double_t signal = 0. ; 
+  // It is a semi-gaussian, 2nd order Gamma function of the general form
+  // v(t) = n**n * Q * A**n / C *(t/tp)**n * exp(-n * t/tp) with 
+  // tp : peaking time par[0]
+  // n  : order of the function
+  // C  : integrating capacitor in the preamplifier
+  // A  : open loop gain of the preamplifier
+  // Q  : the total APD charge to be measured Q = C * energy
   
-  if (xx < par[4] + par[1])    // signal is rising
-    signal = (gRandom->Rndm() + par[3]) * (xx - par[4]) / (par[1] - par[4]) ; 
-  else                                         // signal is decaying
-    signal = (gRandom->Rndm() + par[3]) * TMath::Gaus(xx, par[4] + par[1], par[2]) ;
+  Double_t signal ;
+  Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ; 
   
-  return signal < 0. ? 0. : signal ; 
+  if (xx < 0 || xx > fgTimeMax) 
+    signal = 0. ;  
+  else { 
+    Double_t fac = par[0] * TMath::Power(fgOrder, fgOrder) * TMath::Power(par[1], fgOrder) / fgCapa ; 
+    signal = fac * par[2] * TMath::Power(xx / fgTimePeak, fgOrder) * TMath::Exp(-fgOrder * (xx / fgTimePeak)) ; 
+  }
+  return signal ;  
 }
 
 //__________________________________________________________________
-Bool_t AliEMCAL::RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const 
+Double_t AliEMCAL::RawResponseFunctionMax(Double_t charge, Double_t gain) 
+{
+  return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder) 
+     / ( fgCapa * TMath::Exp(fgOrder) ) );  
+
+}
+//__________________________________________________________________
+Bool_t AliEMCAL::RawSampledResponse(
+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
+  // calculates the raw sampled response AliEMCAL::RawResponseFunction
 
   const Int_t kRawSignalOverflow = 0x3FF ; 
-  Bool_t highGain = kFALSE ; 
-
-  TF1 f1("signal", RawResponseFunction, 0, fkTimeBins, 5);
-
-  f1.SetParNames("Time Max", "Peaking time", "Decay width", "Amp max", "Start time") ; 
-  f1.SetParameter(0, fTimeMax) ; 
-  f1.SetParameter(1, fTimePeak) ; 
-  f1.SetParameter(2, fTimeRes) ; 
-  f1.SetParameter(3, damp) ; 
-  f1.SetParameter(4, dtime) ; 
-
-  for (Int_t iTime = 0; iTime < fkTimeBins; iTime++) {
-    Double_t time = iTime * fTimeMax/fkTimeBins;
-    Double_t signal = f1.Eval(time) ;     
-    adcL[iTime] =  static_cast<Int_t>(signal + 0.5) ;
-    if ( adcL[iTime] > kRawSignalOverflow)  // larger than 10 bits 
-      adcL[iTime] = kRawSignalOverflow ;
-    adcH[iTime] = static_cast<Int_t>(0.5 + (signal / fHighGainFactor)) ;
-    if (adcH[iTime] > 0) 
-      highGain = kTRUE;  
+  Bool_t lowGain = kFALSE ; 
+
+  TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
+
+  for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
+    signalF.SetParameter(0, GetRawFormatHighCharge() ) ; 
+    signalF.SetParameter(1, GetRawFormatHighGain() ) ; 
+    signalF.SetParameter(2, damp) ; 
+    signalF.SetParameter(3, dtime) ; 
+    Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
+    Double_t signal = signalF.Eval(time) ;     
+    if ( static_cast<Int_t>(signal+0.5) > kRawSignalOverflow ){  // larger than 10 bits 
+      signal = kRawSignalOverflow ;
+      lowGain = kTRUE ; 
+    }
+    adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
+
+    signalF.SetParameter(0, GetRawFormatLowCharge() ) ;     
+    signalF.SetParameter(1, GetRawFormatLowGain() ) ; 
+    signal = signalF.Eval(time) ;  
+    if ( static_cast<Int_t>(signal+0.5) > kRawSignalOverflow)  // larger than 10 bits 
+      signal = kRawSignalOverflow ;
+    adcL[iTime] = static_cast<Int_t>(0.5 + signal ) ; 
+
   }
-  return highGain ; 
+  return lowGain ; 
 }
 
 //____________________________________________________________________________
index 540942f..edf62e2 100644 (file)
@@ -45,35 +45,49 @@ class AliEMCAL : public AliDetector {
   {return AliEMCALGeometry::GetInstance(GetTitle(),"") ;  }   
   virtual void    Hits2SDigits();
   virtual Int_t   IsVersion(void) const = 0 ;   
-  Int_t GetRawFormatHighGainFactor() const { return fHighGainFactor ; }  
-  Int_t GetRawFormatHighGainOffset() const { return fHighGainOffset ; }  
+  // Raw Read Out
+  Double_t GetRawFormatCapa() const { return fgCapa ; }   
+  Double_t GetRawFormatHighCharge() const { return fHighCharge ; }  
+  Double_t GetRawFormatHighGain() const { return fHighGain ; }  
+  Double_t GetRawFormatHighLowGainFactor() const { return fHighLowGainFactor ; }  
+  Double_t GetRawFormatLowCharge() const { return ( fHighCharge *  fHighLowGainFactor ) ; }  
+  Double_t GetRawFormatLowGain() const { return ( fHighGain / fHighLowGainFactor ) ; }  
+  Int_t GetRawFormatLowGainOffset() const { return fLowGainOffset ; }  
+  Int_t GetRawFormatOrder() const { return fgOrder ; }   
   Int_t GetRawFormatTimeBins() const { return fkTimeBins ; }    
-  Double_t GetRawFormatTimeMax() const { return fTimeMax ; }   
-  Double_t GetRawFormatTimePeak() const { return fTimePeak ; }    
-  Double_t GetRawFormatTimeRes() const { return fTimeRes ; }   
+  Double_t GetRawFormatTimeMax() const { return fgTimeMax ; }   
+  Double_t GetRawFormatTimePeak() const { return fgTimePeak ; }    
+  Double_t GetRawFormatTimeTrigger() const { return fgTimeTrigger ; }   
+  static Double_t RawResponseFunction(Double_t *x, Double_t *par) ; 
+  static Double_t RawResponseFunctionMax(Double_t charge, Double_t gain) ;
+  //  
   virtual AliLoader* MakeLoader(const char* topfoldername);
-  static Double_t  RawResponseFunction(Double_t *x, Double_t *par) ; 
   virtual void SetTreeAddress() ;              
   virtual const TString Version() const {return TString(" ") ; }   
   AliEMCAL & operator = (const AliEMCAL & /*rvalue*/)  {
     Fatal("operator =", "not implemented") ;  return *this ; }
  
 protected:
+  
+  Bool_t   RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const ; 
 
- Bool_t  RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const ; 
 
   Int_t fBirkC0;    // constants for Birk's Law implementation
   Double_t fBirkC1; // constants for Birk's Law implementation
   Double_t fBirkC2; // constants for Birk's Law implementation
-  Int_t    fHighGainFactor ;   // High gain attenuation factor of the raw RO signal
-  Int_t    fHighGainOffset ;   // offset added to the module id to distinguish high and low gain data
-  static const Int_t fkTimeBins = 256 ;     // number of sampling bins of the raw RO signal  
-  Double_t fTimeMax ;          // maximum sampled time of the raw RO signal
-  Double_t fTimePeak ;         // peaking time of the raw RO signal
-  Double_t fTimeRes ;          // decay rime width of the raw RO signal 
+  static Double_t fgCapa ;              // capacitor of the preamplifier for the raw RO signal
+  Double_t fHighCharge ;                // high charge (to convert energy to charge) for the raw RO signal
+  Double_t fHighGain ;                  // high gain for the raw RO signal
+  Double_t fHighLowGainFactor ;         // high to low gain factor for the raw RO signal
+  Int_t    fLowGainOffset ;             // to separate high from low gain in the DDL
+  static Int_t fgOrder ;                // order of the gamma function for the RO signal
+  static const Int_t fkTimeBins = 256 ; // number of sampling bins of the raw RO signal  
+  static Double_t fgTimeMax ;           // maximum sampled time of the raw RO signal                             
+  static Double_t fgTimePeak ;          // peaking time of the raw RO signal                                    
+  static Double_t fgTimeTrigger ;       // time of the trigger for the RO signal 
 
-  ClassDef(AliEMCAL,8) // Electromagnetic calorimeter (base class)
-
-} ;
+  ClassDef(AliEMCAL,9) // Electromagnetic calorimeter (base class)
+    
+    } ;
 
 #endif // ALIEMCAL_H
index 9ef6fb7..0cfb84e 100644 (file)
 #include <TFile.h>
 #include <TROOT.h>
 #include <TSystem.h>
-#include <TH1D.h>
 #include <TF1.h>
+#include <TGraph.h>
+//#include <TCanvas.h>
+//#include <TFrame.h>
 
 // --- Standard library ---
 
@@ -457,6 +459,110 @@ void AliEMCALGetter::ReadPrimaries()
 }
 
 //____________________________________________________________________________ 
+void AliEMCALGetter::FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Int_t & amp, Double_t & time)
+{
+  // Fits the raw signal time distribution 
+
+  const Int_t kNoiseThreshold = 0 ;
+  Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
+  Double_t signal = 0., signalmax = 0. ;       
+  Double_t energy = time = 0. ; 
+
+  if (lowGainFlag) {
+    timezero1 = timezero2 = signalmax = timemax = 0. ;
+    signalF->FixParameter(0, EMCAL()->GetRawFormatLowCharge()) ; 
+    signalF->FixParameter(1, EMCAL()->GetRawFormatLowGain()) ; 
+    Int_t index ; 
+    for (index = 0; index < EMCAL()->GetRawFormatTimeBins(); index++) {
+      gLowGain->GetPoint(index, time, signal) ; 
+      if (signal > kNoiseThreshold && timezero1 == 0.) 
+       timezero1 = time ;
+      if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
+       timezero2 = time ; 
+      if (signal > signalmax) {
+       signalmax = signal ; 
+       timemax   = time ; 
+      }
+    }
+    signalmax /= EMCAL()->RawResponseFunctionMax(EMCAL()->GetRawFormatLowCharge(), 
+                                               EMCAL()->GetRawFormatLowGain()) ;
+    if ( timezero1 + EMCAL()->GetRawFormatTimePeak() < EMCAL()->GetRawFormatTimeMax() * 0.4 ) { // else its noise 
+      signalF->SetParameter(2, signalmax) ; 
+      signalF->SetParameter(3, timezero1) ;                
+      gLowGain->Fit(signalF, "QRON", "", 0., timezero2); //, "QRON") ; 
+      energy = signalF->GetParameter(2) ; 
+      time   = signalF->GetMaximumX() - EMCAL()->GetRawFormatTimePeak() - EMCAL()->GetRawFormatTimeTrigger() ;
+    }
+  } else {
+    timezero1 = timezero2 = signalmax = timemax = 0. ;
+    signalF->FixParameter(0, EMCAL()->GetRawFormatHighCharge()) ; 
+    signalF->FixParameter(1, EMCAL()->GetRawFormatHighGain()) ; 
+    Int_t index ; 
+    for (index = 0; index < EMCAL()->GetRawFormatTimeBins(); index++) {
+      gHighGain->GetPoint(index, time, signal) ;               
+      if (signal > kNoiseThreshold && timezero1 == 0.) 
+       timezero1 = time ;
+      if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
+       timezero2 = time ; 
+      if (signal > signalmax) {
+       signalmax = signal ;   
+       timemax   = time ; 
+      }
+    }
+    signalmax /= EMCAL()->RawResponseFunctionMax(EMCAL()->GetRawFormatHighCharge(), 
+                                               EMCAL()->GetRawFormatHighGain()) ;;
+    if ( timezero1 + EMCAL()->GetRawFormatTimePeak() < EMCAL()->GetRawFormatTimeMax() * 0.4 ) { // else its noise  
+      signalF->SetParameter(2, signalmax) ; 
+      signalF->SetParameter(3, timezero1) ;               
+      gHighGain->Fit(signalF, "QRON", "", 0., timezero2) ; 
+      energy = signalF->GetParameter(2) ; 
+      time   = signalF->GetMaximumX() - EMCAL()->GetRawFormatTimePeak() - EMCAL()->GetRawFormatTimeTrigger() ;
+    }
+  }
+  
+  if (time == 0. && energy == 0.) 
+    amp = 0 ; 
+  else {
+  AliEMCALDigitizer * digitizer = Digitizer() ; 
+  amp = static_cast<Int_t>( (energy - digitizer->GetECApedestal()) / digitizer->GetECAchannel() + 0.5 ) ; 
+  }
+  // dessin
+//   TCanvas * c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);
+//   c1->SetFillColor(42);
+//   c1->SetGrid();
+//   gLowGain->SetLineColor(2);
+//   gLowGain->SetLineWidth(4);
+//   gLowGain->SetMarkerColor(4);
+//   gLowGain->SetMarkerStyle(21);
+//   gLowGain->SetTitle("Lowgain");
+//   gLowGain->GetXaxis()->SetTitle("X title");
+//   gLowGain->GetYaxis()->SetTitle("Y title");
+//   gLowGain->Draw("ACP");
+  
+//   c1->Update();
+//   c1->GetFrame()->SetFillColor(21);
+//   c1->GetFrame()->SetBorderSize(12);
+//   c1->Modified();
+  
+//   TCanvas * c2 = new TCanvas("c2","A Simple Graph Example",200,10,700,500);
+//   c2->SetFillColor(42);
+//   c2->SetGrid();
+//   gHighGain->SetLineColor(2);
+//   gHighGain->SetLineWidth(4);
+//   gHighGain->SetMarkerColor(4);
+//   gHighGain->SetMarkerStyle(21);
+//   gHighGain->SetTitle("Highgain");
+//   gHighGain->GetXaxis()->SetTitle("X title");
+//   gHighGain->GetYaxis()->SetTitle("Y title");
+//     gHighGain->Draw("ACP");
+  
+//   c2->Update();
+//   c2->GetFrame()->SetFillColor(21);
+//   c2->GetFrame()->SetBorderSize(12);
+//   c2->Modified();
+}
+
+//____________________________________________________________________________ 
 Int_t AliEMCALGetter::ReadRaw(Int_t event)
 {
   // reads the raw format data, converts it into digits format and store digits in Digits()
@@ -465,73 +571,60 @@ Int_t AliEMCALGetter::ReadRaw(Int_t event)
   AliRawReaderFile rawReader(event) ; 
   AliEMCALRawStream in(&rawReader);
   
-  const Int_t kHighGainIdOffset = EMCALGeometry()->GetNTowers()
-    * EMCALGeometry()->GetNPhi()
-    * EMCALGeometry()->GetNZ() 
-    * 2 ;
-  
   Bool_t first = kTRUE ;
-  
-  TH1D hLowGain("hLowGain", "Low Gain", 1000, 0., EMCAL()->GetRawFormatTimeMax()) ; 
-  TH1D hHighGain("hHighGain", "High Gain", 1000, 0., EMCAL()->GetRawFormatTimeMax()) ; 
-  
-  // fit half the gaussian decay rather than AliEMCAL::RawResponseFunction because thiswould give a floating point
-  // exception during error calculation ... To solve... 
-  TF1 * gauss = new TF1("gauss", "gaus", 
-                       EMCAL()->GetRawFormatTimePeak(), 
-                       EMCAL()->GetRawFormatTimeMax() ) ;   
+  TF1 * signalF = new TF1("signal", AliEMCAL::RawResponseFunction, 0, EMCAL()->GetRawFormatTimeMax(), 4);
+  signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero") ; 
   
   Int_t id = -1;
-  Bool_t hgflag = kFALSE ; 
-  
+  Bool_t lowGainFlag = kFALSE ; 
+
   TClonesArray * digits = Digits() ;
   digits->Clear() ; 
   Int_t idigit = 0 ; 
-  
+  Int_t amp = 0 ; 
+  Double_t time = 0. ; 
+
+  TGraph * gLowGain = new TGraph(EMCAL()->GetRawFormatTimeBins()) ; 
+  TGraph * gHighGain= new TGraph(EMCAL()->GetRawFormatTimeBins()) ;  
+
   while ( in.Next() ) { // EMCAL entries loop 
-    
     if ( in.IsNewId() ) {
       if (!first) {
-       hLowGain.Fit(gauss, "QRON") ; 
-       Int_t ampL =  static_cast<Int_t>(gauss->Eval(gauss->GetParameter(2)) + 0.5) ; 
-       Double_t timeL  = EMCAL()->GetRawFormatTimePeak() - gauss->GetParameter(2) ;
-       if (timeL < 0 ) // happens with noise 
-         timeL =  EMCAL()->GetRawFormatTimePeak() ;  
-       if (hgflag) {
-         hHighGain.Fit(gauss, "QRON") ; 
-         Int_t ampH =  static_cast<Int_t>(gauss->Eval(gauss->GetParameter(2)) + 0.5) ; 
-         Double_t timeH  = EMCAL()->GetRawFormatTimePeak() - gauss->GetParameter(2) ; 
-         if (timeH < 0 ) // happens with noise 
-           timeH =  EMCAL()->GetRawFormatTimePeak() ;  
-         new((*digits)[idigit]) AliEMCALDigit( -1, -1, id+kHighGainIdOffset, ampH, timeH) ;
+       FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, amp, time) ; 
+       if (amp > 0) {
+         new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, amp, time) ;        
          idigit++ ; 
        }
-       else {
-         new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, ampL, timeL) ;
-         idigit++ ; 
-       }
-      }
+       Int_t index ; 
+       for (index = 0; index < EMCAL()->GetRawFormatTimeBins(); index++) {
+         gLowGain->SetPoint(index, index * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 0) ;  
+         gHighGain->SetPoint(index, index * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 0) ; 
+       } 
+      }  
       first = kFALSE ; 
-      hLowGain.Reset() ; 
-      hHighGain.Reset() ; 
       id = in.GetId() ; 
-      if (id > 9999 ) { // fixme 
-       hgflag = kTRUE ; 
-      } else 
-       hgflag = kFALSE ; 
+      if (in.GetModule() == EMCAL()->GetRawFormatLowGainOffset() ) 
+       lowGainFlag = kTRUE ; 
+      else 
+       lowGainFlag = kFALSE ; 
     }
-    if (hgflag)
-      hHighGain.Fill(
-                    in.GetTime() * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
-                    in. GetSignal() * EMCAL()->GetRawFormatHighGainFactor() ) ;
+    if (lowGainFlag)
+      gLowGain->SetPoint(in.GetTime(), 
+                        in.GetTime()* EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
+                        in.GetSignal()) ;
     else 
-      hLowGain.Fill(
-                   in.GetTime() * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
-                   in. GetSignal() ) ;
+      gHighGain->SetPoint(in.GetTime(), 
+                         in.GetTime() * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
+                         in.GetSignal() ) ;
+    
   } // EMCAL entries loop
-  
-  delete gauss ; 
-  
+  digits->Sort() ; 
+
+  delete signalF ; 
+  delete gLowGain, gHighGain ; 
+    
   return Digits()->GetEntriesFast() ; 
 }
   
index dddd0ae..3b89de7 100644 (file)
@@ -21,6 +21,8 @@
 #include "TClonesArray.h"
 class TParticle ;
 class TTree ; 
+class TGraph ;
+class TF1 ;
 
 // --- Standard library ---
 
@@ -207,6 +209,8 @@ private:
 
 
   void ReadPrimaries(void) ;
+  void FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Int_t & amp, Double_t & time) ; 
 
 private:
 
index d1dc27f..ac441c0 100644 (file)
@@ -25,6 +25,7 @@ public :
   AliEMCALRawStream(AliRawReader* rawReader);
   
   Int_t            GetId() const {return fPad;};
+  Int_t            GetModule() const {return fSector;}
   Int_t            GetPrevId() const {return fPrevPad;};
   Int_t            GetSignal() const {return fSignal;};
   Int_t            GetTime() const {return fTime;};