Change for a realistic time distribution of the raw signal

diff --git a/PHOS/AliPHOS.cxx b/PHOS/AliPHOS.cxx
index 42af88070d9feba0c6fbf86cf0193ba069620e5d..9eff25a38469336a04a9020ad364621ee40be6d9 100644 (file)
--- a/PHOS/AliPHOS.cxx
+++ b/PHOS/AliPHOS.cxx
@@ -45,8 +45,17 @@ class TFile;
 #include "AliPHOSSDigitizer.h"
 #include "AliPHOSDigit.h"
 #include "AliAltroBuffer.h"
 #include "AliPHOSSDigitizer.h"
 #include "AliPHOSDigit.h"
 #include "AliAltroBuffer.h"

+#include "AliLog.h"

 
 ClassImp(AliPHOS)
 
 ClassImp(AliPHOS)

+
+Double_t AliPHOS::fgCapa        = 1.;        // 1pF 
+Int_t    AliPHOS::fgOrder       = 2 ;
+Double_t AliPHOS::fgTimeMax     = 2.56E-5 ;  // each sample is over 100 ns fTimeMax/fTimeBins
+Double_t AliPHOS::fgTimePeak    = 4.1E-6 ;   // 4 micro seconds
+Double_t AliPHOS::fgTimeTrigger = 100E-9 ;      // 100ns, just for a reference
+ 
+

 //____________________________________________________________________________
   AliPHOS:: AliPHOS() : AliDetector()
 {
 //____________________________________________________________________________
   AliPHOS:: AliPHOS() : AliDetector()
 {


@@ -55,6 +64,7 @@ ClassImp(AliPHOS)
   fQATask = 0;
   fTreeQA = 0;
   fDebug  = 0; 
   fQATask = 0;
   fTreeQA = 0;
   fDebug  = 0; 

+

 }
 
 //____________________________________________________________________________
 }
 
 //____________________________________________________________________________


@@ -65,11 +75,11 @@ AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title)
   fQATask = 0 ; 
   fTreeQA = 0 ;
   fDebug  = 0 ; 
   fQATask = 0 ; 
   fTreeQA = 0 ;
   fDebug  = 0 ; 

-  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     = GetGeometry()->GetNModules() + 1 ;   
+    // offset added to the module id to distinguish high and low gain data

 }
 
 //____________________________________________________________________________
 }
 
 //____________________________________________________________________________


@@ -390,23 +400,24 @@ void AliPHOS::CreateMaterials()
 void AliPHOS::Digits2Raw()
 {
 // convert digits of the current event to raw data
 void AliPHOS::Digits2Raw()
 {
 // convert digits of the current event to raw data

+  
+  AliPHOSLoader * loader = dynamic_cast<AliPHOSLoader*>(fLoader) ; 

 
   // get the digits
 
   // get the digits

-  AliPHOSGetter * gime = AliPHOSGetter::Instance(AliRunLoader::GetGAliceName()) ; 
-  if (!gime) {
-    Error("Digits2Raw", "PHOS 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;
   }
 
 
   if (!digits) {
     Error("Digits2Raw", "no digits found !");
     return;
   }
 

+  // get the digitizer 
+  loader->LoadDigitizer();
+  AliPHOSDigitizer * digitizer = dynamic_cast<AliPHOSDigitizer *>(loader->Digitizer())  ; 
+  

   // get the geometry
   // get the geometry

-  AliPHOSGeometry* geom = gime->PHOSGeometry();
+  AliPHOSGeometry* geom = GetGeometry();

   if (!geom) {
     Error("Digits2Raw", "no geometry found !");
     return;
   if (!geom) {
     Error("Digits2Raw", "no geometry found !");
     return;


@@ -414,7 +425,7 @@ void AliPHOS::Digits2Raw()
 
   // some digitization constants
   const Int_t    kDDLOffset = 0x600; // assigned to PHOS
 
   // some digitization constants
   const Int_t    kDDLOffset = 0x600; // assigned to PHOS

-  const Int_t    kThreshold = 3; // skip digits below this threshold
+  const Int_t    kThreshold = 1; // skip digits below this threshold

 
   AliAltroBuffer* buffer = NULL;
   Int_t prevDDL = -1;
 
   AliAltroBuffer* buffer = NULL;
   Int_t prevDDL = -1;


@@ -423,7 +434,7 @@ void AliPHOS::Digits2Raw()
 
   // loop over digits (assume ordered digits)
   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
 
   // loop over digits (assume ordered digits)
   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {

-    AliPHOSDigit* digit = gime->Digit(iDigit);
+    AliPHOSDigit* digit = dynamic_cast<AliPHOSDigit *>(digits->At(iDigit)) ;

     if (digit->GetAmp() < kThreshold) 
       continue;
     Int_t relId[4];
     if (digit->GetAmp() < kThreshold) 
       continue;
     Int_t relId[4];


@@ -458,24 +469,26 @@ void AliPHOS::Digits2Raw()
     }
 
     // out of time range signal (?)
     }
 
     // out of time range signal (?)

-    if (digit->GetTimeR() > fTimeMax) {
+    if (digit->GetTimeR() > GetRawFormatTimeMax() ) {

       buffer->FillBuffer(digit->GetAmp());
       buffer->FillBuffer(digit->GetAmp());

-      buffer->FillBuffer(fkTimeBins);  // time bin
-      buffer->FillBuffer(3);          // bunch length
+      buffer->FillBuffer(GetRawFormatTimeBins() );  // time bin
+      buffer->FillBuffer(3);          // bunch length      

       buffer->WriteTrailer(3, relId[3], relId[2], module);  // trailer
       
       buffer->WriteTrailer(3, relId[3], relId[2], module);  // trailer
       

-    // simulate linear rise and gaussian decay of signal
+    // calculate the time response function

     } else {
     } else {

-      Bool_t highGain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; 
-      
+      Double_t energy = 0 ;  
+      if ( digit->GetId() <= geom->GetNModules() *  geom->GetNCristalsInModule())
+       energy = digit->GetAmp() * digitizer->GetEMCchannel() + digitizer->GetEMCpedestal() ; 
+      else 
+       energy = digit->GetAmp() * digitizer->GetCPVchannel() + digitizer->GetCPVpedestal() ;
+        
+      RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; 

       
       

-      // write low and eventually high gain channel

       buffer->WriteChannel(relId[3], relId[2], module, 
       buffer->WriteChannel(relId[3], relId[2], module, 

-                          fkTimeBins, adcValuesLow, kThreshold);
-      if (highGain) {
-       buffer->WriteChannel(relId[3], relId[2], module + fHighGainOffset, 
-                            fkTimeBins, adcValuesHigh, 1);
-      }
+                          GetRawFormatTimeBins(), adcValuesHigh, kThreshold);
+      buffer->WriteChannel(relId[3], relId[2], module + fLowGainOffset, 
+                          GetRawFormatTimeBins(), adcValuesLow, kThreshold);

     }
   }
   
     }
   }
   


@@ -486,7 +499,7 @@ void AliPHOS::Digits2Raw()
     delete buffer;
   }
   
     delete buffer;
   }
   

-  gime->PhosLoader()->UnloadDigits();
+  loader->UnloadDigits();

 }
 
 //____________________________________________________________________________
 }
 
 //____________________________________________________________________________


@@ -510,53 +523,70 @@ AliLoader* AliPHOS::MakeLoader(const char* topfoldername)
 }
 
 //__________________________________________________________________
 }
 
 //__________________________________________________________________

-Double_t AliPHOS::RawResponseFunction(Double_t *x, Double_t *par)
+Double_t AliPHOS::RawResponseFunction(Double_t *x, Double_t *par) 

 {
   // Shape of the electronics raw reponse:
 {
   // 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]) ;
-  
-  return signal < 0. ? 0. : signal ; 
+  Double_t signal ;
+  Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ; 
+
+  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 ;  
+}
+
+//__________________________________________________________________
+Double_t AliPHOS::RawResponseFunctionMax(Double_t charge, Double_t gain) 
+{
+  return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder) 
+     / ( fgCapa * TMath::Exp(fgOrder) ) );  
+

 }
 
 //__________________________________________________________________
 }
 
 //__________________________________________________________________

-Bool_t AliPHOS::RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const 
+Bool_t AliPHOS::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, 
 {
   // for a start time dtime and an amplitude damp given by digit, 

-  // calculates the raw sampled response
+  // calculates the raw sampled response AliPHOS::RawResponseFunction

 
   const Int_t kRawSignalOverflow = 0x3FF ; 
 
   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 ; 

 }
 
 //____________________________________________________________________________
 }
 
 //____________________________________________________________________________

diff --git a/PHOS/AliPHOS.h b/PHOS/AliPHOS.h
index 0903baaaa875390d9933e3e57f6adf648f48a752..acb66fb2566e93a15f1acac92e6b8e07cc4d8ae8 100644 (file)
--- a/PHOS/AliPHOS.h
+++ b/PHOS/AliPHOS.h
@@ -48,15 +48,24 @@ public:
   {return AliPHOSGeometry::GetInstance(GetTitle(),"") ;  }
   virtual void    Hits2SDigits();
   virtual Int_t   IsVersion(void) const = 0 ;  
   {return AliPHOSGeometry::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 ; }    
   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);
   AliPHOSQAChecker * QAChecker() {return fQATask;}  
   virtual AliLoader* MakeLoader(const char* topfoldername);
   AliPHOSQAChecker * QAChecker() {return fQATask;}  

-  static Double_t  RawResponseFunction(Double_t *x, Double_t *par) ; 

   virtual void    SetTreeAddress();   
   virtual TTree * TreeQA() const {return fTreeQA; } 
   virtual const TString Version() const {return TString(" ") ; } 
   virtual void    SetTreeAddress();   
   virtual TTree * TreeQA() const {return fTreeQA; } 
   virtual const TString Version() const {return TString(" ") ; } 


@@ -67,20 +76,23 @@ public:
 
 protected:
 
 
 protected:
 

-  Bool_t  RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const ; 
+    Bool_t   RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const ; 

 
 
 
 

-  AliPHOSQAChecker * fQATask ; //! PHOS checkers container
-  TTree * fTreeQA ;            // the QA tree that contains the alarms
-  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 
-
-  ClassDef(AliPHOS,3) // Photon Spectrometer Detector (base class)
-
+  AliPHOSQAChecker * fQATask ;          //! PHOS checkers container
+  TTree * fTreeQA ;                     // the QA tree that contains the alarms
+  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(AliPHOS,4) // Photon Spectrometer Detector (base class)

 } ;
 
 #endif // ALIPHOS_H
 } ;
 
 #endif // ALIPHOS_H

diff --git a/PHOS/AliPHOSGetter.cxx b/PHOS/AliPHOSGetter.cxx
index c8f9b6365e7ed30b47daac1c193af3eda36e47c0..9f94b708c39f7a1e9c16b588a51ef53919d0a711 100644 (file)
--- a/PHOS/AliPHOSGetter.cxx
+++ b/PHOS/AliPHOSGetter.cxx
@@ -46,6 +46,7 @@
 #include <TParticle.h>
 #include <TH1D.h>
 #include <TF1.h>
 #include <TParticle.h>
 #include <TH1D.h>
 #include <TF1.h>

+#include <TGraph.h>

 
 // --- Standard library ---
 
 
 // --- Standard library ---
 


@@ -61,6 +62,7 @@
 #include "AliStack.h"  
 #include "AliPHOSRawStream.h"
 #include "AliRawReaderFile.h"
 #include "AliStack.h"  
 #include "AliPHOSRawStream.h"
 #include "AliRawReaderFile.h"

+#include "AliLog.h"

 
 ClassImp(AliPHOSGetter)
   
 
 ClassImp(AliPHOSGetter)
   


@@ -519,6 +521,74 @@ Bool_t AliPHOSGetter::OpenESDFile()
   return rv ; 
 }
 
   return rv ; 
 }
 

+
+//____________________________________________________________________________ 
+void AliPHOSGetter::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) {
+    Int_t index ; 
+    for (index = 0; index < PHOS()->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 /= PHOS()->RawResponseFunctionMax(PHOS()->GetRawFormatLowCharge(), 
+                                               PHOS()->GetRawFormatLowGain()) ;
+    if ( timezero1 + PHOS()->GetRawFormatTimePeak() < PHOS()->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() - PHOS()->GetRawFormatTimePeak() - PHOS()->GetRawFormatTimeTrigger() ;
+    }
+  } else {
+    timezero1 = timezero2 = signalmax = timemax = 0. ;
+    signalF->FixParameter(0, PHOS()->GetRawFormatHighCharge()) ; 
+    signalF->FixParameter(1, PHOS()->GetRawFormatHighGain()) ; 
+    Int_t index ; 
+    for (index = 0; index < PHOS()->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 /= PHOS()->RawResponseFunctionMax(PHOS()->GetRawFormatHighCharge(), 
+                                               PHOS()->GetRawFormatHighGain()) ;;
+    if ( timezero1 + PHOS()->GetRawFormatTimePeak() < PHOS()->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() - PHOS()->GetRawFormatTimePeak() - PHOS()->GetRawFormatTimeTrigger() ;
+    }
+  }
+  
+  if (time == 0. && energy == 0.) 
+    amp = 0 ; 
+  else {
+  AliPHOSDigitizer * digitizer = Digitizer() ; 
+  amp = static_cast<Int_t>( (energy - digitizer->GetEMCpedestal()) / digitizer->GetEMCchannel() + 0.5 ) ; 
+  }
+}
+

 //____________________________________________________________________________ 
 Int_t AliPHOSGetter::ReadRaw(Int_t event)
 {
 //____________________________________________________________________________ 
 Int_t AliPHOSGetter::ReadRaw(Int_t event)
 {


@@ -527,80 +597,67 @@ Int_t AliPHOSGetter::ReadRaw(Int_t event)
 
   AliRawReaderFile rawReader(event) ; 
   AliPHOSRawStream in(&rawReader);
 
   AliRawReaderFile rawReader(event) ; 
   AliPHOSRawStream in(&rawReader);

-  
-  const Int_t kHighGainIdOffset = PHOSGeometry()->GetNModules()
-    * PHOSGeometry()->GetNPhi()
-    * PHOSGeometry()->GetNZ() 
-    * 2 ;

 
   Bool_t first = kTRUE ;
   
 
   Bool_t first = kTRUE ;
   

-  TH1D hLowGain("hLowGain", "Low Gain", 1000, 0., PHOS()->GetRawFormatTimeMax()) ; 
-  TH1D hHighGain("hHighGain", "High Gain", 1000, 0., PHOS()->GetRawFormatTimeMax()) ; 
-  
-  // fit half the gaussian decay rather than AliPHOS::RawResponseFunction because thiswould give a floating point
-  // exception during error calculation ... To solve... 
-  TF1 * gauss = new TF1("gauss", "gaus", 
-                        PHOS()->GetRawFormatTimePeak(), 
-                        PHOS()->GetRawFormatTimeMax() ) ;   
+  TGraph * gLowGain  = new TGraph(PHOS()->GetRawFormatTimeBins()) ; 
+  TGraph * gHighGain = new TGraph(PHOS()->GetRawFormatTimeBins()) ; 
+ 
+  TF1 * signalF = new TF1("signal", AliPHOS::RawResponseFunction, 0, PHOS()->GetRawFormatTimeMax(), 4);
+  signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero") ; 

 
   Int_t relId[4], id ;
 
   Int_t relId[4], id ;

-  Bool_t hgflag = kFALSE ; 
+  Bool_t lowGainFlag = kFALSE ; 

  
   TClonesArray * digits = Digits() ;
   digits->Clear() ; 
   Int_t idigit = 0 ; 
  
   TClonesArray * digits = Digits() ;
   digits->Clear() ; 
   Int_t idigit = 0 ; 

+  Int_t amp = 0 ; 
+  Double_t time = 0. ; 

   
   while ( in.Next() ) { // PHOS entries loop 
   
   while ( in.Next() ) { // PHOS entries loop 

-        

     if ( (in.IsNewRow() || in.IsNewColumn() || in.IsNewModule()) ) {
       if (!first) {
     if ( (in.IsNewRow() || in.IsNewColumn() || in.IsNewModule()) ) {
       if (!first) {

-       hLowGain.Fit(gauss, "QRON") ; 
-       Int_t ampL =  static_cast<Int_t>(gauss->Eval(gauss->GetParameter(2)) + 0.5) ; 
-       Double_t timeL  = PHOS()->GetRawFormatTimePeak() - gauss->GetParameter(2) ;
-       if (timeL < 0 ) // happens with noise 
-         timeL =  PHOS()->GetRawFormatTimePeak() ;  
-       if (hgflag) {
-         hHighGain.Fit(gauss, "QRON") ; 
-         Int_t ampH =  static_cast<Int_t>(gauss->Eval(gauss->GetParameter(2)) + 0.5) ; 
-         Double_t timeH  = PHOS()->GetRawFormatTimePeak() - gauss->GetParameter(2) ; 
-         if (timeH < 0 ) // happens with noise 
-           timeH =  PHOS()->GetRawFormatTimePeak() ;  
-         new((*digits)[idigit]) AliPHOSDigit( -1, id+kHighGainIdOffset, ampH, timeH) ;
-         idigit++ ; 
-       }
-       else {
-         new((*digits)[idigit]) AliPHOSDigit( -1, id, ampL, timeL) ;
+       
+       FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, amp, time) ; 
+       if (amp > 0) {
+         new((*digits)[idigit]) AliPHOSDigit( -1, id, amp, time) ;     

          idigit++ ; 
        }
       }
       first = kFALSE ; 
          idigit++ ; 
        }
       }
       first = kFALSE ; 

-      hLowGain.Reset() ; 
-      hHighGain.Reset() ; 

       relId[0] = in.GetModule() ;
       relId[0] = in.GetModule() ;

-      if ( relId[0] >= PHOS()->GetRawFormatHighGainOffset() ) { 
-       relId[0] -=  PHOS()->GetRawFormatHighGainOffset() ; 
-       hgflag = kTRUE ; 
+      if ( relId[0] >= PHOS()->GetRawFormatLowGainOffset() ) { 
+       relId[0] -=  PHOS()->GetRawFormatLowGainOffset() ; 
+       lowGainFlag = kTRUE ; 

       } else 
       } else 

-         hgflag = kFALSE ; 
+       lowGainFlag = kFALSE ; 

       relId[1] = 0 ; 
       relId[2] = in.GetRow() ; 
       relId[3] = in.GetColumn() ; 
       PHOSGeometry()->RelToAbsNumbering(relId, id) ;   
     }
       relId[1] = 0 ; 
       relId[2] = in.GetRow() ; 
       relId[3] = in.GetColumn() ; 
       PHOSGeometry()->RelToAbsNumbering(relId, id) ;   
     }

-    if (hgflag)
-      hHighGain.Fill(
-                    in.GetTime() * PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 
-                    in. GetSignal() * PHOS()->GetRawFormatHighGainFactor() ) ;
+    if (lowGainFlag)
+      gLowGain->SetPoint(in.GetTime(), 
+                    in.GetTime(), 
+                    in.GetSignal()) ;

     else 
     else 

-      hLowGain.Fill(
-                   in.GetTime() * PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 
-                   in. GetSignal() ) ;
+      gHighGain->SetPoint(in.GetTime(), 
+                        in.GetTime() * PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 
+                        in.GetSignal() ) ;
+

   } // PHOS entries loop
   } // PHOS entries loop

+
+  FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, amp, time) ; 
+  if (amp > 0 ) {
+    new((*digits)[idigit]) AliPHOSDigit( -1, id, amp, time) ;
+    idigit++ ; 
+  }
+
+  delete signalF ; 
+  delete gLowGain, gHighGain ; 

   
   

-  delete gauss ; 
-  
-  return Digits()->GetEntriesFast() ; 
+  return digits->GetEntriesFast() ; 

 }
 
 //____________________________________________________________________________ 
 }
 
 //____________________________________________________________________________ 

diff --git a/PHOS/AliPHOSGetter.h b/PHOS/AliPHOSGetter.h
index 420db469029c78b3e77b623a030d418bd1e77aac..fafad386b8a95f231fa224c3982f934fcc7f1fa5 100644 (file)
--- a/PHOS/AliPHOSGetter.h
+++ b/PHOS/AliPHOSGetter.h
@@ -20,6 +20,8 @@
 #include "TObject.h"  
 class TParticle ;
 class TTree ; 
 #include "TObject.h"  
 class TParticle ;
 class TTree ; 

+class TGraph ;
+class TF1 ;

 
 // --- Standard library ---
 
 
 // --- Standard library ---
 


@@ -205,7 +207,9 @@ private:
   Int_t ReadTreeE(Int_t event) ;    
   Bool_t OpenESDFile() ;
   void ReadPrimaries(void) ;
   Int_t ReadTreeE(Int_t event) ;    
   Bool_t OpenESDFile() ;
   void ReadPrimaries(void) ;

-  
+
+  void FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Int_t & amp, Double_t & time) ; 
+

 private:
   
   AliPHOSBeamTestEvent * fBTE ;           //! Header if BeamTest Event
 private:
   
   AliPHOSBeamTestEvent * fBTE ;           //! Header if BeamTest Event