Change for a realistic time distribution of the raw signal

[u/mrichter/AliRoot.git] / PHOS / AliPHOSGetter.cxx

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 <TGraph.h>
 
 // --- Standard library ---
 
@@ -61,6 +62,7 @@
 #include "AliStack.h"  
 #include "AliPHOSRawStream.h"
 #include "AliRawReaderFile.h"
+#include "AliLog.h"
 
 ClassImp(AliPHOSGetter)
   
@@ -519,6 +521,74 @@ Bool_t AliPHOSGetter::OpenESDFile()
   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)
 {
@@ -527,80 +597,67 @@ Int_t AliPHOSGetter::ReadRaw(Int_t event)
 
   AliRawReaderFile rawReader(event) ; 
   AliPHOSRawStream in(&rawReader);
-  
-  const Int_t kHighGainIdOffset = PHOSGeometry()->GetNModules()
-    * PHOSGeometry()->GetNPhi()
-    * PHOSGeometry()->GetNZ() 
-    * 2 ;
 
   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 ;
-  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. ; 
   
   while ( in.Next() ) { // PHOS entries loop 
-        
     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 ; 
-      hLowGain.Reset() ; 
-      hHighGain.Reset() ; 
       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 
-         hgflag = kFALSE ; 
+       lowGainFlag = kFALSE ; 
       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 
-      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
+
+  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() ; 
 }
 
 //____________________________________________________________________________