First implementation of digitizer for AD
authormbroz <Michal.Broz@cern.ch>
Fri, 11 Jul 2014 20:41:26 +0000 (22:41 +0200)
committermbroz <Michal.Broz@cern.ch>
Fri, 11 Jul 2014 20:41:26 +0000 (22:41 +0200)
13 files changed:
AD/AliAD.cxx
AD/AliAD.h
AD/AliADCalibData.cxx
AD/AliADCalibData.h
AD/AliADConst.h [new file with mode: 0644]
AD/AliADDigitizer.cxx
AD/AliADDigitizer.h
AD/AliADReconstructor.cxx
AD/AliADSDigit.h
AD/AliADdigit.cxx
AD/AliADdigit.h
AD/AliADv1.cxx
AD/AliADv1.h

index ecac3bd..a29a69a 100644 (file)
@@ -254,6 +254,37 @@ AliDigitizer* AliAD::CreateDigitizer(AliDigitizationInput* digInput) const
 }
 
 //_____________________________________________________________________________
+void AliAD::Hits2Digits(){
+  //
+  // Converts hits to digits
+  //
+  // Creates the AD digitizer 
+  AliADDigitizer* dig = new AliADDigitizer(this,AliADDigitizer::kHits2Digits);
+
+  // Creates the digits
+  dig->Digitize("");
+
+  // deletes the digitizer
+  delete dig;
+}
+
+//_____________________________________________________________________________
+void AliAD::Hits2SDigits(){
+  //
+  // Converts hits to summable digits
+  //
+  // Creates the AD digitizer 
+  AliADDigitizer* dig = new AliADDigitizer(this,AliADDigitizer::kHits2SDigits);
+
+  // Creates the sdigits
+  dig->Digitize("");
+
+  // deletes the digitizer
+  delete dig;
+}
+
+
+//_____________________________________________________________________________
 
 void AliAD::Digits2Raw()
 {
index 16d2c87..f13aa84 100644 (file)
@@ -30,7 +30,9 @@ public:
   AliDigitizer*  CreateDigitizer(AliDigitizationInput* digInput) const;
   virtual AliTriggerDetector* CreateTriggerDetector() const { return new AliADTrigger();}
   
-  virtual      void     Digits2Raw();
+  virtual    void      Hits2Digits();
+  virtual    void      Hits2SDigits();
+  virtual    void      Digits2Raw();
   virtual    Bool_t     Raw2SDigits(AliRawReader*);
   virtual    void      SetADAToInstalled(Bool_t b){fSetADAToInstalled = b;}
   virtual    void      SetADCToInstalled(Bool_t b){fSetADCToInstalled = b;}
index ed515fa..78c8079 100644 (file)
 
 /* $Id: AliADCalibData.cxx,                                            */
 
+#include <TMath.h>
+#include <TObjString.h>
+#include <TMap.h>
+#include <TH1F.h>
+#include <TH2F.h>
 
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
 #include "AliADCalibData.h"
-#include "TList.h"
-#include "TCanvas.h"
+#include "AliADConst.h"
+#include "AliLog.h"
 
 ClassImp(AliADCalibData)
 
 
 //________________________________________________________________
-AliADCalibData::AliADCalibData()
+AliADCalibData::AliADCalibData():
+  fLightYields(NULL),
+  fPMGainsA(NULL),
+  fPMGainsB(NULL)
 {
-       for (Int_t imod = 0; imod < 60; imod++)
-       {
-               fEfficiencies[imod]=0.;
-               fRates[imod]=0.;
-               fModulesActivity[imod]=0.;
-       } 
-}
+  // default constructor
+  
+    for(int t=0; t<16; t++) {
+        fMeanHV[t]      = 100.0;
+        fWidthHV[t]     = 0.0; 
+       fTimeOffset[t]  = 5.0;
+        fTimeGain[t]    = 1.0;
+       fDeadChannel[t]= kFALSE;
+       fDiscriThr[t]  = 2.5;
+    }
+    for(int t=0; t<32; t++) {
+        fPedestal[t]    = 0.0;     
+        fSigma[t]       = 0.0;        
+        fADCmean[t]     = 0.0;      
+        fADCsigma[t]    = 0.0;
+    }
+    for(int i=0; i<kNCIUBoards ;i++) {
+       fTimeResolution[i]  = 25./256.;     // Default time resolution
+       fWidthResolution[i] = 25./64.;     // Default time width resolution
+       fMatchWindow[i] = 4;
+       fSearchWindow[i] = 16;
+       fTriggerCountOffset[i] = 3247;
+       fRollOver[i] = 3563;
+    }
 
+}
 //________________________________________________________________
 void AliADCalibData::Reset()
 {
@@ -42,41 +70,77 @@ void AliADCalibData::Reset()
 
 //________________________________________________________________
 AliADCalibData::AliADCalibData(const char* name) :
-  TNamed()
+  fLightYields(NULL),
+  fPMGainsA(NULL),
+  fPMGainsB(NULL)
 {
-  TString namst = "Calib_";
-  namst += name;
-  SetName(namst.Data());
-  SetTitle(namst.Data());
-
+  // Constructor
+   TString namst = "Calib_";
+   namst += name;
+   SetName(namst.Data());
+   SetTitle(namst.Data());
+   for(int t=0; t<16; t++) {
+       fMeanHV[t]      = 100.0;
+       fWidthHV[t]     = 0.0; 
+       fTimeOffset[t]  = 5.0;
+       fTimeGain[t]    = 1.0;
+       fDeadChannel[t]= kFALSE;
+       fDiscriThr[t]  = 2.5;
+    }
+   for(int t=0; t<32; t++) {
+       fPedestal[t]    = 0.0;     
+       fSigma[t]       = 0.0;        
+       fADCmean[t]     = 0.0;      
+       fADCsigma[t]    = 0.0;
+   }
+   for(int i=0; i<kNCIUBoards ;i++) {
+       fTimeResolution[i]  = 25./256.;    // Default time resolution in ns / channel
+       fWidthResolution[i] = 25./64.;     // Default time width resolution in ns / channel
+       fMatchWindow[i] = 4;
+       fSearchWindow[i] = 16;
+       fTriggerCountOffset[i] = 3247;
+       fRollOver[i] = 3563;
+   }
 }
 
 //________________________________________________________________
 AliADCalibData::AliADCalibData(const AliADCalibData& calibda) :
-  TNamed(calibda)
+  TNamed(calibda),
+  fLightYields(NULL),
+  fPMGainsA(NULL),
+  fPMGainsB(NULL)
 {
 // copy constructor
 
   SetName(calibda.GetName());
   SetTitle(calibda.GetName());
   
-  // there are 60 modules. Note that number of first module is 1 (one)
-  for(int t=0; t<60; t++) 
-  {
-       fEfficiencies[t] =calibda.GetEfficiency(t);
-       fRates[t] = calibda.GetRate(t);
-       fModulesActivity[t] = calibda.GetModuleActivity(t);
+  for(int t=0; t<32; t++) { 
+      fPedestal[t] = calibda.GetPedestal(t);
+      fSigma[t]    = calibda.GetSigma(t);
+      fADCmean[t]  = calibda.GetADCmean(t);
+      fADCsigma[t] = calibda.GetADCsigma(t); }
+      
+  for(int t=0; t<16; t++) { 
+      fMeanHV[t]       = calibda.GetMeanHV(t);
+      fWidthHV[t]      = calibda.GetWidthHV(t);        
+      fTimeOffset[t]   = calibda.GetTimeOffset(t);
+      fTimeGain[t]     = calibda.GetTimeGain(t); 
+      fDeadChannel[t]  = calibda.IsChannelDead(t);
+      fDiscriThr[t]    = calibda.GetDiscriThr(t);
+  }  
+  
+  for(int i=0; i<kNCIUBoards ;i++) {
+      fTimeResolution[i]  = calibda.GetTimeResolution(i);
+      fWidthResolution[i] = calibda.GetWidthResolution(i);       
+      fMatchWindow[i] = calibda.GetMatchWindow(i);
+      fSearchWindow[i] = calibda.GetSearchWindow(i);
+      fTriggerCountOffset[i] = calibda.GetTriggerCountOffset(i);
+      fRollOver[i] = calibda.GetRollOver(i);
   }
+  
 }
-//_______________________________________________________________
-void AliADCalibData::Draw(Option_t *)
-{
-
-  //fHits->Draw();
 
-}
 //________________________________________________________________
 AliADCalibData &AliADCalibData::operator =(const AliADCalibData& calibda)
 {
@@ -84,54 +148,148 @@ AliADCalibData &AliADCalibData::operator =(const AliADCalibData& calibda)
 
   SetName(calibda.GetName());
   SetTitle(calibda.GetName());
-  // there are 60 modules. Note that number of first module is 1 (one)
-  for(int t=0; t<60; t++) 
-  {
-       fEfficiencies[t] =calibda.GetEfficiency(t);
-       fRates[t] = calibda.GetRate(t);
-       fModulesActivity[t] = calibda.GetModuleActivity(t);
+  
+  for(int t=0; t<32; t++) {
+      fPedestal[t] = calibda.GetPedestal(t);
+      fSigma[t]    = calibda.GetSigma(t);
+      fADCmean[t]  = calibda.GetADCmean(t);
+      fADCsigma[t] = calibda.GetADCsigma(t); }
+      
+  for(int t=0; t<16; t++) {
+      fMeanHV[t]       = calibda.GetMeanHV(t);
+      fWidthHV[t]      = calibda.GetWidthHV(t);        
+      fTimeOffset[t]   = calibda.GetTimeOffset(t);
+      fTimeGain[t]     = calibda.GetTimeGain(t); 
+      fDeadChannel[t]  = calibda.IsChannelDead(t);
+      fDiscriThr[t]    = calibda.GetDiscriThr(t);
+  }   
+  for(int i=0; i<kNCIUBoards ;i++) {
+      fTimeResolution[i]  = calibda.GetTimeResolution(i);
+      fWidthResolution[i] = calibda.GetWidthResolution(i);       
+      fMatchWindow[i] = calibda.GetMatchWindow(i);
+      fSearchWindow[i] = calibda.GetSearchWindow(i);
+      fTriggerCountOffset[i] = calibda.GetTriggerCountOffset(i);
+      fRollOver[i] = calibda.GetRollOver(i);
   }
+   
   return *this;
+  
 }
-//_______________________________________________________________
-/*void AliADCalibData::AddHisto(TH1D *fHist)
-{
-    
 
+//________________________________________________________________
+AliADCalibData::~AliADCalibData()
+{
+  // destructor
+  if (fLightYields)
+    delete [] fLightYields;
+  if (fPMGainsA)
+    delete [] fPMGainsA;
+  if (fPMGainsB)
+    delete [] fPMGainsB;
+}
 
- = (TH1D*)fHist->Clone("hnew");
+//________________________________________________________________
+Int_t AliADCalibData::GetBoardNumber(Int_t channel)
+{
+  // Get FEE board number
+  // from offline channel index
+  if (channel >= 0 && channel < 16) return (channel);
+  //if (channel >=8 && channel < 16) return (channel / 2);
 
-     
-   
+  AliErrorClass(Form("Wrong channel index: %d",channel));
+  return -1;
 }
-*/
 
 //________________________________________________________________
-AliADCalibData::~AliADCalibData()
+Float_t AliADCalibData::GetLightYields(Int_t channel)
 {
-  
-}
+  // Get the light yield efficiency
+  // for a given channel
+  if (!fLightYields) InitLightYields();
+
+  if (channel >= 0 && channel < 64) {
+    return fLightYields[channel];
+  }
 
-                                                                                   
+  AliError(Form("Wrong channel index: %d",channel));
+  return 0;
+}
 
 //________________________________________________________________
-void AliADCalibData::SetEfficiencies(Float_t* Eff)
+void  AliADCalibData::InitLightYields()
 {
-  // there are 60 modules. Note that number of first module is 1 (one)
-  if(Eff) for(int t=0; t<60; t++) fEfficiencies[t] = Eff[t];
-  else for(int t=0; t<60; t++) fEfficiencies[t] = 0.0;
+  // Initialize the light yield factors
+  // Read from a separate OCDB entry
+  if (fLightYields) return;
+
+  AliCDBEntry *entry = AliCDBManager::Instance()->Get("VZERO/Calib/LightYields");
+  if (!entry) AliFatal("AD light yields are not found in OCDB !");
+  TH1F *yields = (TH1F*)entry->GetObject();
+
+  fLightYields = new Float_t[16];
+  for(Int_t i = 0 ; i < 16; ++i) {
+    fLightYields[i] = yields->GetBinContent(i+1);
+  }
 }
 
-void AliADCalibData::SetRates(Float_t* Rt)
+//________________________________________________________________
+void  AliADCalibData::InitPMGains()
 {
-   if(Rt) for (int t=0;t<60; t++) fRates[t] = Rt[t];
-else for (int t=0;t<60; t++) fRates[t] = 0.0;
+  // Initialize the PM gain factors
+  // Read from a separate OCDB entry
+  if (fPMGainsA) return;
+
+  AliCDBEntry *entry = AliCDBManager::Instance()->Get("VZERO/Calib/PMGains");
+  if (!entry) AliFatal("VZERO PM gains are not found in OCDB !");
+  TH2F *gains = (TH2F*)entry->GetObject();
+
+  fPMGainsA = new Float_t[16];
+  fPMGainsB = new Float_t[16];
+  for(Int_t i = 0 ; i < 16; ++i) {
+    fPMGainsA[i] = gains->GetBinContent(i+1,1);
+    fPMGainsB[i] = gains->GetBinContent(i+1,2);
+  }
 }
 
-void AliADCalibData::SetModulesActivity(Float_t* Mac)
+//________________________________________________________________
+Float_t AliADCalibData::GetGain(Int_t channel)
 {
-       if(Mac) for (int t=0;t<60;t++) fModulesActivity[t] = Mac[t];
-       else for (int t=0;t<60;t++) fModulesActivity[t] = 0.0;
+  // Computes the PM gains
+  // Argument passed is the PM number (aliroot numbering)
+  if (!fPMGainsA) InitPMGains();
+
+  // High Voltage retrieval from Calibration Data Base:  
+  Float_t hv = fMeanHV[channel];
+  Float_t gain = 0;
+  if (hv>0)
+    gain = TMath::Exp(fPMGainsA[channel]+fPMGainsB[channel]*TMath::Log(hv));
+  return gain;
 }
 
+//________________________________________________________________
+Float_t AliADCalibData::GetCalibDiscriThr(Int_t channel, Bool_t scaled)
+{
+  // The method returns actual TDC discri threshold
+  // extracted from the data.
+  //
+  // In case scaled flag is set the threshold is scaled
+  // so that to get 4.0 for a FEE threshold of 4.0.
+  // In this way we avoid a change in the slewing correction
+  // for the entire 2010 p-p data.
+  //
+  // The method is to be moved to OCDB object.
+
+  Float_t thr = GetDiscriThr(channel);
+
+  Float_t calThr = 0;
+  if (thr <= 1.) 
+    calThr = 3.1;
+  else if (thr >= 2.)
+    calThr = (3.1+1.15*thr-1.7);
+  else
+    calThr = (3.1-0.3*thr+0.3*thr*thr);
+
+  if (scaled) calThr *= 4./(3.1+1.15*4.-1.7);
+
+  return calThr;
+}
index 7645813..e0aab0d 100644 (file)
@@ -5,7 +5,7 @@
  * See cxx source for full Copyright notice                               */
 
 #include "TNamed.h"
-#include "TH1D.h"
+#include "AliADConst.h"
 class AliADCalibData: public TNamed {
 
  public:
@@ -17,26 +17,82 @@ class AliADCalibData: public TNamed {
   virtual ~AliADCalibData();
   void Reset();
 
-  Float_t* GetEfficiencies() const { return (float*)fEfficiencies; }
-  Float_t  GetEfficiency(Int_t i) const { return fEfficiencies[i-1];}
-  Float_t* GetRates() const {return (float*)fRates;}
-  Float_t GetRate(Int_t i) const {return fRates[i-1];}
-  Float_t* GetModulesActivity() const {return (float*)fModulesActivity;}
-  Float_t GetModuleActivity(Int_t i) const {return fModulesActivity[i-1];}
- // TList*  GetHistos()const {return Hist;} 
-  void SetRates(Float_t* Rt);
-  void SetRate(Float_t rate, Int_t mod){fRates[mod-1]=rate;}
-  void SetEfficiencies(Float_t* Eff);
-  void SetEfficiency(Float_t eff, Int_t mod) {fEfficiencies[mod-1]=eff;}
-  void Draw(Option_t *option="");
-  void SetModulesActivity(Float_t* Mac);
-  void SetModuleActivity(Float_t mac,Int_t mod){fModulesActivity[mod-1]=mac;}
+  Float_t  GetPedestal(Int_t channel)   const {return fPedestal[channel];}
+  Float_t* GetPedestal()   const {return (float*)fPedestal;}
+  Float_t  GetSigma(Int_t channel)   const {return fSigma[channel];}
+  Float_t* GetSigma()   const {return (float*)fSigma;}
+  Float_t  GetADCmean(Int_t channel)   const {return fADCmean[channel];}
+  Float_t* GetADCmean()   const {return (float*)fADCmean;}
+  Float_t  GetADCsigma(Int_t channel)  const {return fADCsigma[channel];}
+  Float_t* GetADCsigma()   const {return (float*)fADCsigma;}
+  Float_t  GetMeanHV(Int_t channel)    const {return fMeanHV[channel];}
+  Float_t* GetMeanHV()   const {return (float*)fMeanHV;} 
+  Float_t  GetWidthHV(Int_t channel)   const {return fWidthHV[channel];}
+  Float_t* GetWidthHV()   const {return (float*)fWidthHV;}
+  Bool_t   IsChannelDead(Int_t channel)        const {return fDeadChannel[channel];}
+  Bool_t*  GetDeadMap()   const {return (bool*)fDeadChannel;} 
+   
+  Float_t  GetGain(Int_t channel);
+  Float_t  GetTimeOffset(Int_t channel)        const {return fTimeOffset[channel];}
+  Float_t* GetTimeOffset()   const {return (float*)fTimeOffset;}
+  Float_t  GetTimeGain(Int_t channel)  const {return fTimeGain[channel];}
+  Float_t* GetTimeGain()   const {return (float*)fTimeGain;}
+
+  Float_t* GetTimeResolution() const {return (Float_t*) fTimeResolution;};
+  Float_t  GetTimeResolution(Int_t board ) const  {return ((board>=0 && board<kNCIUBoards)?fTimeResolution[board]:0);};
+
+  Float_t* GetWidthResolution() const {return (Float_t*) fWidthResolution;};
+  Float_t  GetWidthResolution(Int_t board ) const  {return ((board>=0 && board<kNCIUBoards)?fWidthResolution[board]:0);};
+
+  const UInt_t*  GetMatchWindow() const { return fMatchWindow; }
+  UInt_t   GetMatchWindow(Int_t board) const { return ((board>=0 && board<kNCIUBoards)?fMatchWindow[board]:0); }
+  const UInt_t*  GetSearchWindow() const { return fSearchWindow; }
+  UInt_t   GetSearchWindow(Int_t board) const { return ((board>=0 && board<kNCIUBoards)?fSearchWindow[board]:0); }
+  const UInt_t*  GetTriggerCountOffset() const { return fTriggerCountOffset; }
+  UInt_t   GetTriggerCountOffset(Int_t board) const { return ((board>=0 && board<kNCIUBoards)?fTriggerCountOffset[board]:0); }
+  const UInt_t*  GetRollOver() const { return fRollOver; }
+  UInt_t   GetRollOver(Int_t board) const { return ((board>=0 && board<kNCIUBoards)?fRollOver[board]:0); }
+
+  Float_t  GetDiscriThr(Int_t channel) const {return fDiscriThr[channel];}
+  Float_t* GetDiscriThr()   const {return (Float_t*)fDiscriThr;}
+
+  static Int_t GetBoardNumber(Int_t channel);
+  
+  Float_t  GetLightYields(Int_t channel);
+
+  Float_t *GetPMGainsA() const { return fPMGainsA; }
+  Float_t *GetPMGainsB() const { return fPMGainsB; }
+
+  Float_t  GetCalibDiscriThr(Int_t channel, Bool_t scaled);
 
 
  protected:
-  Float_t fEfficiencies[60];
-  Float_t fRates[60];
-  Float_t fModulesActivity[60];
+  void     InitLightYields();
+  void     InitPMGains();
+
+  Float_t  fPedestal[32];     // Mean pedestal values
+  Float_t  fSigma[32];        // Sigmas of pedestal peaks
+  Float_t  fADCmean[32];      // ADC mean values
+  Float_t  fADCsigma[32];     // ADC sigma values
+  Float_t  fMeanHV[16];        // Mean PMT HV needed to compute MIP value
+  Float_t  fWidthHV[16];       // Width of the PMT HV
+  
+  Float_t  fTimeOffset[16];    // Time offsets of the TDC
+  Float_t  fTimeGain[16];      // Gain factors of the TDC
+  Bool_t   fDeadChannel[16];   // List of dead channels
+  Float_t  fTimeResolution[kNCIUBoards]; // Time Resolution of the TDC (ns / channel)
+  Float_t  fWidthResolution[kNCIUBoards]; // Time Width Resolution of the TDC (ns / channel)
+
+  UInt_t   fMatchWindow[kNCIUBoards]; // HPTDC matching window (25ns units)
+  UInt_t   fSearchWindow[kNCIUBoards];// HPTDC search window (25ns units)
+  UInt_t   fTriggerCountOffset[kNCIUBoards]; // HPTDC trigger count offset (25ns units)
+  UInt_t   fRollOver[kNCIUBoards]; // HPTDC roll-over (25ns units)
+
+  Float_t  fDiscriThr[16];     // Discriminator thresholds
+
+  Float_t *fLightYields;       //! Light Yields channel by channel (read from separate OCDB entry)
+  Float_t *fPMGainsA;          //! PM gain factors channel by channel (read from separate OCDB entry)
+  Float_t *fPMGainsB;          //! PM gain factors channel by channel (read from separate OCDB entry)
 
   ClassDef(AliADCalibData,1)    // AD Calibration data
 };
diff --git a/AD/AliADConst.h b/AD/AliADConst.h
new file mode 100644 (file)
index 0000000..922e448
--- /dev/null
@@ -0,0 +1,20 @@
+#ifndef ALIVZEROCONST_H
+#define ALIVZEROCONST_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. */
+
+const Float_t kIntTimeRes = 0.39; // intrinsic time resolution of the scintillator
+const Float_t kV0CDelayCables = 8.1; // delay cables on the C side (in ns)
+const Float_t kV0Offset = -1338.6; // general V0 offset between the TDCs and the trigger
+const Int_t   kNClocks = 21; // Number of ADC clocks that are read out
+const Float_t kChargePerADC = 0.6e-12; // Charge per ADC
+const Int_t   kMinTDCWidth = 13; // minimum signal width measured by TDC
+const Int_t   kMaxTDCWidth = 128; // maximum signal width measured by TDC
+const Float_t kPMRespTime = 6.0; // PM response time (corresponds to 1.9 ns rise time)
+const Float_t kPMTransparency = 0.25; // Transparency of the first dynode of the PM
+const Float_t kPMNbOfSecElec = 6.0;   // Number of secondary electrons emitted from first dynode (per ph.e.)
+const Float_t kPhotoCathodeEfficiency = 0.18; // Photocathode efficiency
+const Int_t   kNCIUBoards = 16; //Number of CIU boards
+
+#endif
+
index 945e2d5..2e3ce2d 100644 (file)
@@ -38,6 +38,7 @@
 #include "AliRun.h"
 #include "AliAD.h"
 #include "AliADhit.h"
+#include "AliADConst.h"
 #include "AliRunLoader.h"
 #include "AliLoader.h"
 #include "AliGRPObject.h"
 #include "AliCDBManager.h"
 #include "AliCDBStorage.h"
 #include "AliCDBEntry.h"
-//#include "AliADCalibData.h"
+#include "AliADCalibData.h"
 #include "AliCTPTimeParams.h"
 #include "AliLHCClockPhase.h"
 #include "AliADdigit.h"
 #include "AliADDigitizer.h"
-//#include "AliADSDigit.h"
+#include "AliADSDigit.h"
 
 ClassImp(AliADDigitizer)
 
-//____________________________________________________________________________ 
-AliADDigitizer::AliADDigitizer()
-  :AliDigitizer(),
-   fNdigits(0),
-   fDigits(0)
+//____________________________________________________________________________
+ AliADDigitizer::AliADDigitizer()
+                   :AliDigitizer(),
+                    fCalibData(GetCalibData()),
+                    fNdigits(0),
+                    fDigits(0),
+                    fSignalShape(NULL),
+                    fPMResponse(NULL),
+                    fSinglePhESpectrum(NULL),
+                   fEvenOrOdd(kFALSE),
+                   fTask(kHits2Digits),
+                   fAD(NULL)
 {
   // default constructor
+  // Initialize OCDB and containers used in the digitization
+
+  Init();
 }
 
 //____________________________________________________________________________ 
-AliADDigitizer::AliADDigitizer(AliDigitizationInput* digInput)
-  :AliDigitizer(digInput),
-   fNdigits(0),
-   fDigits(0)
+  AliADDigitizer::AliADDigitizer(AliAD *AD, DigiTask_t task)
+                    :AliDigitizer(),
+                    fCalibData(GetCalibData()),
+                    fNdigits(0),
+                     fDigits(0),
+                    fSignalShape(NULL),
+                     fPMResponse(NULL),
+                     fSinglePhESpectrum(NULL),
+                    fEvenOrOdd(kFALSE),
+                    fTask(task),
+                    fAD(AD)
+{
+  // constructor
+  // Initialize OCDB and containers used in the digitization
+
+  Init();
+}
+           
+//____________________________________________________________________________ 
+  AliADDigitizer::AliADDigitizer(AliDigitizationInput* digInput)
+                    :AliDigitizer(digInput),
+                    fCalibData(GetCalibData()),
+                    fNdigits(0),
+                     fDigits(0),
+                    fSignalShape(NULL),
+                     fPMResponse(NULL),
+                     fSinglePhESpectrum(NULL),
+                    fEvenOrOdd(kFALSE),
+                    fTask(kHits2Digits),
+                    fAD(NULL)
 {
   // constructor
+  // Initialize OCDB and containers used in the digitization
+
+  Init();
 }
            
 //____________________________________________________________________________ 
-AliADDigitizer::~AliADDigitizer()
+  AliADDigitizer::~AliADDigitizer()
 {
   // destructor
   
@@ -83,145 +123,524 @@ AliADDigitizer::~AliADDigitizer()
     fDigits=0; 
   }
 
+  if (fSignalShape) {
+    delete fSignalShape;
+    fSignalShape = NULL;
+  }
+  if (fPMResponse) {
+    delete fPMResponse;
+    fPMResponse = NULL;
+  }
+  if (fSinglePhESpectrum) {
+    delete fSinglePhESpectrum;
+    fSinglePhESpectrum = NULL;
+  }
+
+  for(Int_t i = 0 ; i < 16; ++i) {
+    if (fTime[i]) delete [] fTime[i];
+  }
 }
 
 //____________________________________________________________________________ 
 Bool_t AliADDigitizer::Init()
 {
-       fDigits = new TClonesArray ("AliADdigit",1000);
-       return kTRUE;
+  // Initialises the digitizer
+  // Initialize OCDB and containers used in the digitization
+
+  // check if the digitizer was already initialized
+  if (fSignalShape) return kTRUE;
+
+  fSignalShape = new TF1("ADSignalShape",this,&AliADDigitizer::SignalShape,0,200,6,"AliADDigitizer","SignalShape");
+  //  fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,2.9,6.40982,3.69339e-01);
+  //  fSignalShape->SetParameters(1.34330e+00,1.13007e+02,-4.95705e-01,
+  //                         3.68911e+00,1.01040e+00, 3.94675e-01);
+  fSignalShape->SetParameters(-1.07335e+00,2.16002e+01,-1.26133e-01,
+                             1.41619e+00,5.50334e-01,3.86111e-01);
+  fPMResponse = new TF1("ADPMResponse",this,&AliADDigitizer::PMResponse,-kPMRespTime,2.*kPMRespTime,0,"AliADDigitizer","PMResponse");
+  fSinglePhESpectrum = new TF1("ADSinglePhESpectrum",this,&AliADDigitizer::SinglePhESpectrum,0,20,0,"AliADDigitizer","SinglePhESpectrum");
+  
+  // Now get the CTP L0->L1 delay
+  AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/CTP/CTPtiming");
+  if (!entry) AliFatal("CTP timing parameters are not found in OCDB !");
+  AliCTPTimeParams *ctpParams = (AliCTPTimeParams*)entry->GetObject();
+  Float_t l1Delay = (Float_t)ctpParams->GetDelayL1L0()*25.0;
+
+  AliCDBEntry *entry1 = AliCDBManager::Instance()->Get("GRP/CTP/TimeAlign");
+  if (!entry1) AliFatal("CTP time-alignment is not found in OCDB !");
+  AliCTPTimeParams *ctpTimeAlign = (AliCTPTimeParams*)entry1->GetObject();
+  l1Delay += ((Float_t)ctpTimeAlign->GetDelayL1L0()*25.0);
+
+  AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("VZERO/Calib/TimeDelays");
+  if (!entry2) AliFatal("AD time delays are not found in OCDB !");
+  TH1F *delays = (TH1F*)entry2->GetObject();
+
+  AliCDBEntry *entry3 = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase");
+  if (!entry3) AliFatal("LHC clock-phase shift is not found in OCDB !");
+  AliLHCClockPhase *phase = (AliLHCClockPhase*)entry3->GetObject();
+
+  for(Int_t i = 0 ; i < 16; ++i) {
+
+    for(Int_t j = 0; j < kNClocks; ++j) fAdc[i][j] = 0;
+    fLeadingTime[i] = fTimeWidth[i] = 0;
+
+    fPmGain[i] = fCalibData->GetGain(i);
+
+    fAdcPedestal[i][0] = fCalibData->GetPedestal(i);
+    fAdcSigma[i][0]    = fCalibData->GetSigma(i); 
+    fAdcPedestal[i][1] = fCalibData->GetPedestal(i+16);
+    fAdcSigma[i][1]    = fCalibData->GetSigma(i+16); 
+
+    Int_t board = AliADCalibData::GetBoardNumber(i);
+    fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+
+                            (Float_t)kMaxTDCWidth*fCalibData->GetWidthResolution(board))/
+                           fCalibData->GetTimeResolution(board));
+    fNBinsLT[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0)/
+                             fCalibData->GetTimeResolution(board));
+    fBinSize[i] = fCalibData->GetTimeResolution(board);
+    fHptdcOffset[i] = (((Float_t)fCalibData->GetRollOver(board)-
+                       (Float_t)fCalibData->GetTriggerCountOffset(board))*25.0+
+                      fCalibData->GetTimeOffset(i)-
+                      l1Delay-
+                      phase->GetMeanPhase()+
+                      delays->GetBinContent(i+1)+
+                      kV0Offset);
+    fClockOffset[i] = (((Float_t)fCalibData->GetRollOver(board)-
+                       (Float_t)fCalibData->GetTriggerCountOffset(board))*25.0+
+                      fCalibData->GetTimeOffset(i)-
+                      l1Delay+
+                      kV0Offset);
+
+    fTime[i] = new Float_t[fNBins[i]];
+    memset(fTime[i],0,fNBins[i]*sizeof(Float_t));
+    
+   // AliWarning(Form("PMT %d,PM gain %f, fNBins %d, TimeBinSize %f,",i, fPmGain[i], fNBins[i],fBinSize[i]));
+    
+  }
+
+  return kTRUE;
+
 }
-//____________________________________________________________________________ 
 
+//____________________________________________________________________________ 
 void AliADDigitizer::Digitize(Option_t* /*option*/) 
 {   
    // Creates digits from hits
-       // tem. variables
-  Int_t modules[16]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
-  Int_t moduls[16]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
-  Int_t mods;
+  fNdigits = 0;  
+
+  if (fAD && !fDigInput) {
+    AliLoader *loader = fAD->GetLoader();
+    if (!loader) {
+      AliError("Can not get AD Loader via AliAD object!");
+      return;
+    }
+    AliRunLoader* runLoader = AliRunLoader::Instance();
+    for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); ++iEvent) {
+      runLoader->GetEvent(iEvent);
+      if (fTask == kHits2Digits) {
+       DigitizeHits();
+       DigitizeSDigits();
+       WriteDigits(loader);
+      }
+      else {
+       DigitizeHits();
+       WriteSDigits(loader);
+      }
+    }
+  }
+  else if (fDigInput) {
+      ReadSDigits();
+      DigitizeSDigits();
+      AliRunLoader *currentLoader = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
+      AliLoader *loader = currentLoader->GetLoader("ADLoader");
+      if (!loader) { 
+       AliError("Cannot get AD Loader via RunDigitizer!");
+       return;
+      }
+      WriteDigits(loader);
+  }
+  else {
+    AliFatal("Invalid digitization task! Exiting!");
+  }
+}
+
+//____________________________________________________________________________ 
+void AliADDigitizer::DigitizeHits()
+{
+  // Digitize the hits to the level of
+  // SDigits (fTime arrays)
+
+  for(Int_t i = 0 ; i < 16; ++i) {
+    memset(fTime[i],0,fNBins[i]*sizeof(Float_t));
+    fLabels[i][0] = fLabels[i][1] = fLabels[i][2] = -1;
+  }
+  Float_t integral = fPMResponse->Integral(-kPMRespTime,2.*kPMRespTime);
+  Float_t meansPhE = fSinglePhESpectrum->Mean(0,20);
+  
+     AliLoader* loader = fAD->GetLoader();
+     if (!loader) {
+       AliError("Can not get AD Loader!");
+       return;
+     }
+     loader->LoadHits();
+     TTree* treeH = loader->TreeH();
+     if (!treeH) {
+       AliError("Cannot get TreeH!");
+       return;
+     }
+     TClonesArray* hits = fAD->Hits();
+
+//  Now makes Digits from hits
+     Int_t nTracks = (Int_t) treeH->GetEntries();
+     for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+         fAD->ResetHits();
+         treeH->GetEvent(iTrack);
+         Int_t nHits = hits->GetEntriesFast();
+         for (Int_t iHit = 0; iHit < nHits; iHit++) {
+          AliADhit* hit = (AliADhit *)hits->UncheckedAt(iHit);
+          Int_t nPhot = hit->GetNphot();
+          Int_t pmt  = hit->GetCell();//One PM per cell in AD                          
+          if (pmt < 0) continue;
+          Int_t trackLabel = hit->GetTrack();
+          for(Int_t l = 0; l < 3; ++l) {
+            if (fLabels[pmt][l] < 0) {
+              fLabels[pmt][l] = trackLabel;
+              break;
+            }
+          }
+          Float_t dt_scintillator = gRandom->Gaus(0,kIntTimeRes);
+          Float_t t = dt_scintillator + hit->GetTof();
+          //if (pmt < 32) t += kV0CDelayCables;
+          t += fHptdcOffset[pmt];
+          Int_t nPhE;
+          Float_t prob = fCalibData->GetLightYields(pmt)*kPhotoCathodeEfficiency; // Optical losses included!
+          if (nPhot > 100)
+            nPhE = (Int_t)gRandom->Gaus(prob*Float_t(nPhot)+0.5,
+                                        sqrt(Float_t(nPhot)*prob*(1.-prob)));
+          else
+            nPhE = gRandom->Binomial(nPhot,prob);
+          Float_t charge = TMath::Qe()*fPmGain[pmt]*fBinSize[pmt]/integral;
+          
+          
+          for (Int_t iPhE = 0; iPhE < nPhE; ++iPhE) {       
+            Float_t tPhE = t + fSignalShape->GetRandom(0,fBinSize[pmt]*Float_t(fNBins[pmt]));
+            Float_t gainVar = fSinglePhESpectrum->GetRandom(0,20)/meansPhE;
+            Int_t firstBin = TMath::Max(0,(Int_t)((tPhE-kPMRespTime)/fBinSize[pmt]));
+            Int_t lastBin = TMath::Min(fNBins[pmt]-1,(Int_t)((tPhE+2.*kPMRespTime)/fBinSize[pmt]));
+            for(Int_t iBin = firstBin; iBin <= lastBin; ++iBin) {
+              Float_t tempT = fBinSize[pmt]*(0.5+iBin)-tPhE;
+              fTime[pmt][iBin] += gainVar*charge*fPMResponse->Eval(tempT);
+            }
+          }         // ph.e. loop
+         }           // hit loop
+     }               // track loop
+     loader->UnloadHits();
+}
 
-       // loaders
+//____________________________________________________________________________ 
+void AliADDigitizer::DigitizeSDigits()
+{
+  // Digitize the fTime arrays (SDigits) to the level of
+  // Digits (fAdc arrays)
+  for(Int_t i = 0 ; i < 16; ++i) {
+    for(Int_t j = 0; j < kNClocks; ++j) fAdc[i][j] = 0;
+    fLeadingTime[i] = fTimeWidth[i] = 0;
+  }
 
-       AliRunLoader* outRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
-       if (!outRunLoader)
-       {
-               Error("Exec","Can not get output Run Loader");
-               return;
+  Float_t maximum = 0.9*fSignalShape->GetMaximum(0,200); // Not exact, one needs to do this on the convoluted
+  Float_t integral2 = fSignalShape->Integral(0,200); // function. Anyway the effect is small <10% on the 2.5 ADC thr
+  for (Int_t ipmt = 0; ipmt < 16; ++ipmt) {
+    Float_t thr = fCalibData->GetCalibDiscriThr(ipmt,kFALSE)*kChargePerADC*maximum*fBinSize[ipmt]/integral2;
+    //Float_t thr = 1e-25;
+       
+    Bool_t ltFound = kFALSE, ttFound = kFALSE;
+    for (Int_t iBin = 0; iBin < fNBins[ipmt]; ++iBin) {
+      Float_t t = fBinSize[ipmt]*Float_t(iBin);
+      
+      if (fTime[ipmt][iBin] > thr) {
+       if (!ltFound && (iBin < fNBinsLT[ipmt])) {
+         ltFound = kTRUE;
+         fLeadingTime[ipmt] = t;
        }
-       AliLoader* outLoader = outRunLoader->GetLoader("ADLoader");
-       if (!outLoader)
-       {
-               Error("Exec","Can not get output AD Loader");
-               return;
+      }
+      else {
+       if (ltFound) {
+         if (!ttFound) {
+           ttFound = kTRUE;
+           fTimeWidth[ipmt] = t - fLeadingTime[ipmt];
+         }
        }
+      }
+      Float_t tadc = t - fClockOffset[ipmt];
+      Int_t clock = kNClocks/2 - Int_t(tadc/25.0);
+      if (clock >= 0 && clock < kNClocks)
+       fAdc[ipmt][clock] += fTime[ipmt][iBin]/kChargePerADC;
+    }
+    AliDebug(1,Form("Channel %d Offset %f Time %f",ipmt,fClockOffset[ipmt],fLeadingTime[ipmt]));
+    Int_t board = AliADCalibData::GetBoardNumber(ipmt);
+    if (ltFound && ttFound) {
+      fTimeWidth[ipmt] = fCalibData->GetWidthResolution(board)*
+       Float_t(Int_t(fTimeWidth[ipmt]/fCalibData->GetWidthResolution(board)));
+      if (fTimeWidth[ipmt] < Float_t(kMinTDCWidth)*fCalibData->GetWidthResolution(board))
+       fTimeWidth[ipmt] = Float_t(kMinTDCWidth)*fCalibData->GetWidthResolution(board);
+      if (fTimeWidth[ipmt] > Float_t(kMaxTDCWidth)*fCalibData->GetWidthResolution(board))
+       fTimeWidth[ipmt] = Float_t(kMaxTDCWidth)*fCalibData->GetWidthResolution(board);
+    }
+  }
+
+  fEvenOrOdd = gRandom->Integer(2);
+  for (Int_t j=0; j<16; ++j){
+    for (Int_t iClock = 0; iClock < kNClocks; ++iClock) {
+      Int_t integrator = (iClock + fEvenOrOdd) % 2;
+      AliDebug(1,Form("ADC %d %d %f",j,iClock,fAdc[j][iClock]));
+      fAdc[j][iClock]  += gRandom->Gaus(fAdcPedestal[j][integrator], fAdcSigma[j][integrator]);
+    }
+  }
+       
+}
+
+//____________________________________________________________________________ 
+void AliADDigitizer::ReadSDigits()
+{
+  // Read SDigits which are then to precessed
+  // in the following method
+  for(Int_t i = 0 ; i < 16; ++i) {
+    memset(fTime[i],0,fNBins[i]*sizeof(Float_t));
+    fLabels[i][0] = fLabels[i][1] = fLabels[i][2] = -1;
+  }
+
+  // Loop over input files
+  Int_t nFiles= fDigInput->GetNinputs();
+  for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
+    // Get the current loader 
+    AliRunLoader* currentLoader = 
+      AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(inputFile));
+
+    AliLoader *loader = currentLoader->GetLoader("ADLoader");
+    loader->LoadSDigits("READ");
   
-       outLoader->LoadDigits("update");
-       if (!outLoader->TreeD()) outLoader->MakeTree("D");
-       outLoader->MakeDigitsContainer();
-       TTree* treeD = outLoader->TreeD();
-       Int_t bufsize = 16000;
-       treeD->Branch("ADdigit",&fDigits, bufsize);
-
-       const Float_t eMin = 0;//1.52; //! MeVs, minimum energy
-       // loop over inputs
-
-       for (Int_t iInput=0; iInput < fDigInput->GetNinputs();iInput++)
-       {
-               AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(iInput));
-               AliLoader* loader = runLoader->GetLoader("ADLoader");
-               if (!loader)
-               {
-                       Error("Exec","Can not get AD Loader for input %d",iInput);
-                       continue;
-               }
-               if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
-               AliAD* ad = (AliAD*) runLoader->GetAliRun()->GetDetector("AD");
-               if (!ad)
-               {
-                       Error("Exec","No AD detector for input %d",iInput);
-                       continue;
-               }
-               loader->LoadHits();
-               TTree* treeH = loader->TreeH();
-               if (!treeH)
-               {
-                       Error("Exec","Cannot get TreeH for input %d",iInput);
-                       continue;
-               }
-               TClonesArray* hits = ad->Hits();
-
-               // here I loop over tracks
-               Int_t nTracks = (Int_t) treeH->GetEntries();
-               for (Int_t iTrack=0; iTrack < nTracks; iTrack++)
-               {
-                       ad->ResetHits();
-                       treeH->GetEvent(iTrack);
-                       Int_t nHits = hits->GetEntriesFast();
-                       // here comes the loop over AD hits
-                       for (Int_t iHit=0; iHit < nHits; iHit++)
-                       {
-                               AliADhit* hit = (AliADhit *)hits->UncheckedAt(iHit);
-                               Float_t eloss_mev = hit->GetEloss()*1000.0;
-                               Int_t module = hit->GetModule();
-                               //cout << "Module AD!!! " << module << endl;
-                               // at some point we shoukd have some calib objects to set real theresholds 
-                               // simple checking on hit, minimum energy at scintillator pad should be > 1.52 MeV's
-                               if (eloss_mev > eMin)
-                               {
-                                 modules[module] = 1;//cout << "energy: " << eloss_mev << endl;
-                               }else modules[module] = 0;
-                       }// end loop over hits
-               } // endo loop over tracks
-               for (Int_t i=0; i<16; i++)
-               {
-                       moduls[i] = modules[i];
-                       //cout << "iModule: " << i <<  " AD hits: " << moduls[i] << endl;
-               }
-
-               loader->UnloadHits();
-
-       } // end loop over inputs
-
-       // here I add the hits to the TreeD
-
-       Int_t tracks[3] = {-1,-1,-1};
-       for (Int_t i=0; i<16; i++)
-       {
-               if (moduls[i]==1)
-               {
-                       mods = i;
-                       AddDigit(tracks,mods,mods);
-               }
+    // Get the tree of summable digits
+    TTree* sdigitsTree = loader->TreeS();
+    if (!sdigitsTree)  {
+      AliError("No sdigit tree from digInput");
+      continue;
+    }
+
+    // Get the branch 
+    TBranch* sdigitsBranch = sdigitsTree->GetBranch("ADSDigit");
+    if (!sdigitsBranch) {
+      AliError("Failed to get sdigit branch");
+      return;
+    }
+
+    // Set the branch address
+    TClonesArray *sdigitsArray = NULL;
+    sdigitsBranch->SetAddress(&sdigitsArray);
+
+    // Sum contributions from the sdigits
+    // Get number of entries in the tree 
+    Int_t nentries  = Int_t(sdigitsBranch->GetEntries());
+    for (Int_t entry = 0; entry < nentries; ++entry)  {
+      sdigitsBranch->GetEntry(entry);
+      // Get the number of sdigits 
+      Int_t nsdigits = sdigitsArray->GetEntries();
+      for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) {
+       AliADSDigit* sDigit = static_cast<AliADSDigit*>(sdigitsArray->UncheckedAt(sdigit));
+       Int_t pmNumber = sDigit->PMNumber();
+       Int_t nbins = sDigit->GetNBins();
+       if (nbins != fNBins[pmNumber]) {
+         AliError(Form("Incompatible number of bins between digitizer (%d) and sdigit (%d) for PM %d! Skipping sdigit!",
+                       fNBins[pmNumber],nbins,pmNumber));
+         continue;
        }
-       treeD->Fill();
-       outLoader->WriteDigits("OVERWRITE");
-       outLoader->UnloadDigits();
-       ResetDigit();
+       // Sum the charges
+       Float_t *charges = sDigit->GetCharges();
+       for(Int_t iBin = 0; iBin < nbins; ++iBin) fTime[pmNumber][iBin] += charges[iBin];
+       //for(Int_t iBin = 0; iBin < nbins; ++iBin) AliWarning(Form("Charge %e ",fTime[pmNumber][iBin]));
+       // and the labels
+       Int_t *labels = sDigit->GetTracks();
+       Int_t j = 0;
+       for(Int_t i = 0; i < 3; ++i) {
+         if (fLabels[pmNumber][i] < 0) {
+           if (labels[j] < 0) break;
+           fLabels[pmNumber][i] = labels[j];
+           j++;
+         }
+       }
+      }
+    }
+    loader->UnloadSDigits();
+  }
+}
+
+
+//____________________________________________________________________________
+void AliADDigitizer::WriteDigits(AliLoader *loader)
+{
+  // Take fAdc arrays filled by the previous
+  // method and produce and add digits to the digit Tree
+
+  loader->LoadDigits("UPDATE");
+
+  if (!loader->TreeD()) loader->MakeTree("D");
+  loader->MakeDigitsContainer();
+  TTree* treeD  = loader->TreeD();
+  DigitsArray();
+  treeD->Branch("ADDigit", &fDigits); 
+  
+  Short_t *chargeADC = new Short_t[kNClocks];
+  for (Int_t i=0; i<16; i++) {      
+    for (Int_t j = 0; j < kNClocks; ++j) {
+      Int_t tempadc = Int_t(fAdc[i][j]);
+      if (tempadc > 1023) tempadc = 1023;
+      chargeADC[j] = tempadc;
+    }
+    AddDigit(i, fLeadingTime[i], fTimeWidth[i], Bool_t((10+fEvenOrOdd)%2), chargeADC, fLabels[i]);
+  }
+  delete [] chargeADC;
+
+  treeD->Fill();
+  loader->WriteDigits("OVERWRITE");  
+  loader->UnloadDigits();     
+  ResetDigits();
+}
+
+//____________________________________________________________________________
+void AliADDigitizer::WriteSDigits(AliLoader *loader)
+{
+  // Take fTime arrays filled by the previous
+  // method and produce and add sdigits to the sdigit Tree
+
+  loader->LoadSDigits("UPDATE");
+
+  if (!loader->TreeS()) loader->MakeTree("S");
+  loader->MakeSDigitsContainer();
+  TTree* treeS  = loader->TreeS();
+  SDigitsArray();
+  treeS->Branch("ADSDigit", &fDigits); 
+  
+  for (Int_t ipmt = 0; ipmt < 16; ++ipmt) {
+    AddSDigit(ipmt,fNBins[ipmt],fTime[ipmt],fLabels[ipmt]);
+  }
+
+  treeS->Fill();
+  loader->WriteSDigits("OVERWRITE");  
+  loader->UnloadSDigits();     
+  ResetDigits();
+}
+
+
+
+//____________________________________________________________________________
+void AliADDigitizer::AddDigit(Int_t pmnumber, Float_t time, Float_t width, Bool_t integrator, Short_t *chargeADC, Int_t *labels) 
+ { 
+// Adds Digit 
+  TClonesArray &ldigits = *fDigits;  
+        
+  new(ldigits[fNdigits++]) AliADdigit(pmnumber,time,width,integrator,chargeADC,labels);
+        
+}
+//____________________________________________________________________________
+void AliADDigitizer::AddSDigit(Int_t pmnumber, Int_t nbins, Float_t *charges, Int_t *labels) 
+ { 
+// Adds SDigit 
+  TClonesArray &ldigits = *fDigits;  
+        
+  new(ldigits[fNdigits++]) AliADSDigit(pmnumber,nbins,charges,labels);
+        
+}
+//____________________________________________________________________________
+void AliADDigitizer::ResetDigits()
+{
+
+// Clears Digits
+
+  fNdigits = 0;
+  if (fDigits) fDigits->Clear();
+}
+
+//____________________________________________________________________________
+AliADCalibData* AliADDigitizer::GetCalibData() const
+
+{
+/*/
+  AliCDBManager *man = AliCDBManager::Instance();
+
+  AliCDBEntry *entry=0;
+
+  entry = man->Get("AD/Calib/Data");
+
+  AliADCalibData *calibdata = 0;
+
+  if (entry) calibdata = (AliADCalibData*) entry->GetObject();
+  if (!calibdata)  AliFatal("No calibration data from calibration database !");
+/*/
+  AliADCalibData *calibdata = new AliADCalibData();
+  
+  return calibdata;
+
 }
 
 //____________________________________________________________________________
-void AliADDigitizer::AddDigit(Int_t* track, Int_t module, Float_t cell) 
-{  
-   // Adds Digit 
-   TClonesArray &ldigits = *fDigits;  
-   new(ldigits[fNdigits++]) AliADdigit(track,module,cell);
+double AliADDigitizer::SignalShape(double *x, double *par)
+{
+  // this function simulates the time
+  // of arrival of the photons at the
+  // photocathode
+  Double_t xx = x[0];
+  if (xx <= par[0]) return 0;
+  Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
+  if (xx <= par[3]) return a;
+  Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
+  Double_t f = a*b/(a+b);
+  AliDebug(100,Form("x=%f func=%f",xx,f));
+  return f;
 }
+
 //____________________________________________________________________________
-void AliADDigitizer::AddDigit(Int_t* modul,Float_t cell) 
-{  
-   // Adds Digit 
-   TClonesArray &ldigits = *fDigits;  
-   new(ldigits[fNdigits++]) AliADdigit(modul,cell);
+double AliADDigitizer::PMResponse(double *x, double * /* par */)
+{
+  // this function describes the
+  // PM time response to a single
+  // photoelectron
+  Double_t xx = x[0]+kPMRespTime;
+  return xx*xx*TMath::Exp(-xx*xx/(kPMRespTime*kPMRespTime));
 }
 
 //____________________________________________________________________________
-void AliADDigitizer::ResetDigit()
+double AliADDigitizer::SinglePhESpectrum(double *x, double * /* par */)
 {
-   // Clears Digits
-   fNdigits = 0;
-   if (fDigits) fDigits->Delete();
+  // this function describes the
+  // PM amplitude response to a single
+  // photoelectron
+  Double_t xx = x[0];
+  if (xx < 0) return 0;
+  return (TMath::Poisson(xx,kPMNbOfSecElec)+kPMTransparency*TMath::Poisson(xx,1.0));
+}
+//____________________________________________________________________
+TClonesArray* AliADDigitizer::DigitsArray() 
+{
+  // Initialize digit array if not already and
+  // return pointer to it. 
+  if (!fDigits) { 
+    fDigits = new TClonesArray("AliADdigit", 16);
+    fNdigits = 0;
+  }
+  return fDigits;
+}
+
+//____________________________________________________________________
+TClonesArray* AliADDigitizer::SDigitsArray() 
+{
+  // Initialize sdigit array if not already and
+  // return pointer to it. 
+  if (!fDigits) { 
+    fDigits = new TClonesArray("AliADSDigit", 16);
+    fNdigits = 0;
+  }
+  return fDigits;
 }
 
index 78e2dd0..803aafc 100644 (file)
 // --- AliRoot header files ---
 
 #include "AliDigitizer.h"
-
-// #include "AliADConst.h"
+#include "AliADConst.h"
 
 class TClonesArray;
+class TF1;
 class AliDigitizationInput;
 class AliCDBManager;
 class AliCDBStorage;
+class AliADCalibData;
+class AliAD;
 
 class AliADDigitizer: public AliDigitizer {
 
 public:
-                    AliADDigitizer() ;                              // default constructor
-                    AliADDigitizer(AliDigitizationInput* digInput); // constructor
-  virtual          ~AliADDigitizer() ;              // destructor
-
-  virtual Bool_t    Init();
-  virtual void      Digitize(Option_t* option=0);
-  void AddDigit(Int_t* track, Int_t module, Float_t cell);
-  void AddDigit(Int_t* modul, Float_t cell);   
+   enum DigiTask_t { 
+     kHits2Digits, 
+     kHits2SDigits
+   };
+  
+   AliADDigitizer() ;                       // default constructor
+   AliADDigitizer(AliAD *AD, DigiTask_t task);         // constructor
+   AliADDigitizer(AliDigitizationInput* digInput); // constructor
+   virtual ~AliADDigitizer() ;              // destructor
+
+   virtual Bool_t Init();
+   virtual void   Digitize(Option_t* option=0);
+
+   void DigitizeHits();
+   void DigitizeSDigits();
+   void WriteDigits(AliLoader *loader);
+   void WriteSDigits(AliLoader *loader);
+   void ReadSDigits();
+
+   void AddDigit(Int_t pmnumber, Float_t time, Float_t width, Bool_t integrator, Short_t *chargeADC, Int_t *labels);
+   void AddSDigit(Int_t pmnumber, Int_t nbins, Float_t *charges, Int_t *labels);
+   TClonesArray* DigitsArray(); 
+   TClonesArray* SDigitsArray(); 
+   void ResetDigits();
+                                               
+   AliADCalibData *GetCalibData() const;
 
-            void    ResetDigit();
-            
+   TF1*   GetSignalShape() const { return fSignalShape; }
+   TF1*   GetPMResponse() const { return fPMResponse; }
+   TF1*   GetSinglePhESpectrum() const { return fSinglePhESpectrum; }
+   double SignalShape(double *x, double *par);
+   double PMResponse(double *x, double *par);
+   double SinglePhESpectrum(double *x, double *par);
+
+ protected:
  
-private:
+   AliADCalibData *fCalibData;  //! calibration data
  
-                    AliADDigitizer(const AliADDigitizer& /*digitizer*/); 
-                    AliADDigitizer& operator = (const AliADDigitizer& /*digitizer*/); 
-  
+ private:
+   AliADDigitizer(const AliADDigitizer& /*digitizer*/); 
+      
+   AliADDigitizer& operator = (const AliADDigitizer& /*digitizer*/); 
+
+   Int_t    fNdigits;                //! Number of digits
+   TClonesArray *fDigits;            //! List of digits
+   
+   TF1*     fSignalShape;            // function which describes the PMT signal shape
+   TF1*     fPMResponse;             // function which describes the PM time response
+   TF1*     fSinglePhESpectrum;      // function which describes the single ph.e. PM response
+
+   Float_t  fAdc[16][kNClocks];      //! Container for ADC samples
+   Float_t  fLeadingTime[16];        //! Leading time container
+   Float_t  fTimeWidth[16];          //! Time width container
+   Float_t  fAdcPedestal[16][2];     //! Pedestals, one per integrator
+   Float_t  fAdcSigma[16][2];        //! Sigma of pedestals
+   Float_t  fPmGain[16];             //! PMT gains
+   Int_t    fNBins[16];              //! Number of bins in fTime container
+   Int_t    fNBinsLT[16];            //! Number of bins in fTime container (match window only)
+   Float_t  fBinSize[16];            //! Bin size in fTime container
+   Float_t  fHptdcOffset[16];        //! HPTDC time offsets channel by channel
+   Float_t  fClockOffset[16];        //! Clock offsets channel by channel
 
-           Int_t    fNdigits;           //! Number of digits
-    TClonesArray*   fDigits;            //! List of digits
+   Float_t *fTime[16];               //! Main container used in digitization
+   Int_t    fLabels[16][3];          //! Container for MC labels
+   Bool_t   fEvenOrOdd;              //! Choise of integrator in central ADC sample
 
+   DigiTask_t fTask;                 //! The task (to be) executed by the digitizer
+   AliAD  *fAD;                //! Pointer to AliDetector object
    ClassDef(AliADDigitizer,1)     // digitizer for AD
 
 };
index a87acfd..1ceac83 100644 (file)
@@ -98,7 +98,7 @@ void AliADReconstructor::FillESD(TTree* digitsTree, TTree* /*clustersTree*/,AliE
     
   for (Int_t d=0; d<nDigits; d++) {    
     AliADdigit* digit = (AliADdigit*) fDigitsArray->At(d);
-    Int_t module = digit->GetCell();
+    Int_t module = digit->PMNumber();
  //   printf("AD Module: %d\n",module);
     ADHits[module] = kTRUE;
   }  
index 7307fc5..a83f1ac 100644 (file)
@@ -26,7 +26,7 @@ class AliADSDigit: public AliDigit  {
                         AliADSDigit(const AliADSDigit& /*sdigit*/); 
                         AliADSDigit& operator = (const AliADSDigit& /*sdigit*/); 
 
-    Int_t               fPMNumber;      // PhotoMultiplier number (0 to 63)
+    Int_t               fPMNumber;      // PhotoMultiplier number (0 to 16)
     Int_t               fNBins;         // Number of charge bins
     Float_t*            fCharges;       //[fNBins] Array with charges
 
index 8057d0e..4d03126 100644 (file)
  **************************************************************************/
 
 #include "AliADdigit.h"
+#include "AliADConst.h"
 
 ClassImp(AliADdigit)
 
 //__________________________________________________________________________
 AliADdigit::AliADdigit()
    :AliDigit(),
-    fModule(0),
-    fCell(0)
+    fPMNumber(0),
+    fTime(0.),
+    fWidth(0.),
+    fIntegrator(0)
 
 {
   // Standard default
   // constructor 
+  for(Int_t iClock = 0; iClock < kNClocks; ++iClock) fChargeADC[iClock] = 0;
 }
 
 //__________________________________________________________________________
-AliADdigit::AliADdigit(Int_t module, Float_t cellPad)
-       : AliDigit(),
-       fModule(module),
-       fCell(cellPad)
-{
-}
+AliADdigit::AliADdigit(Int_t   PMnumber, Float_t time, 
+                             Float_t width,
+                            Bool_t integrator,
+                            Short_t *chargeADC,
+                            Int_t *labels)
+:AliDigit(),
+fPMNumber(PMnumber),
+fTime(time),
+fWidth(width),
+fIntegrator(integrator)
+{  
+  // Constructor
+  // Used in the digitizer
+  if (chargeADC) {
+    for(Int_t iClock = 0; iClock < kNClocks; ++iClock)
+      fChargeADC[iClock] = chargeADC[iClock];
+  }
+  else {
+    for(Int_t iClock = 0; iClock < kNClocks; ++iClock)
+      fChargeADC[iClock] = 0;
+  }
 
-//__________________________________________________________________________
-AliADdigit::AliADdigit(Int_t* tracks, Int_t  module, Float_t cellPad)
-       :AliDigit(tracks),
-       fModule(module),
-       fCell(cellPad)
-{
-}
-//__________________________________________________________________________
-AliADdigit::AliADdigit(Int_t* module, Float_t cellPad)
-       : AliDigit(module),
-         fModule(0),
-         fCell(cellPad)
-{
-}
-//__________________________________________________________________________
-AliADdigit::~AliADdigit()
-{
-  //
-  //
-  //
+  if (labels)
+    for(Int_t iTrack = 0; iTrack < 3; ++iTrack) fTracks[iTrack] = labels[iTrack];
 }
 
 
index 48d35d5..ab56788 100644 (file)
@@ -6,25 +6,34 @@
 /* $Id: AliADdigit.h  $ */
 
 #include "AliDigit.h"
+#include "AliADConst.h"
 
 //_____________________________________________________________________________
 class AliADdigit: public AliDigit  {
 
 public:
-                   AliADdigit();
-                  AliADdigit(Int_t* tracks, Int_t module, Float_t cell);
-                  AliADdigit(Int_t* module, Float_t cell);
-                  AliADdigit(Int_t module, Float_t cell);
-  virtual          ~AliADdigit();
-  virtual void      Print(const Option_t* option="") const;
-
-  Int_t GetModule() const {return fModule;}
-  Float_t GetCell() const {return fCell;}
+    AliADdigit();
+    AliADdigit(Int_t   PMnumber, Float_t time, 
+                  Float_t TimeWidth,
+                 Bool_t  Integrator,
+                 Short_t *chargeADC = 0,
+                 Int_t *labels = 0);
+    virtual ~AliADdigit() {};
+    virtual void Print(const Option_t* option="") const;
 
+    Int_t   PMNumber()   const {return fPMNumber;}    
+    Short_t ADC()        const {return fChargeADC[kNClocks/2];}
+    Float_t Time()       const {return fTime;}
+    Float_t Width()      const {return fWidth;} 
+    Bool_t  Integrator() const {return fIntegrator;}
+    Short_t ChargeADC(Int_t clock) const {return (clock >= 0 && clock < kNClocks) ? fChargeADC[clock] : 0;}
     
-private:
-  Int_t fModule; //! module producing the digit 
-  Float_t           fCell;                  // Time of Flight
+  protected:
+    Int_t   fPMNumber;      // PhotoMultiplier number (0 to 16)
+    Float_t fTime;          // Time of Flight
+    Float_t fWidth;         // Width of the time distribution
+    Bool_t  fIntegrator;    // Integrator used
+    Short_t fChargeADC[kNClocks]; // ADC samples as present in raw data
 
   ClassDef(AliADdigit,1)  // AD Digit class
 };
index 74de1bb..d91dc72 100644 (file)
@@ -426,16 +426,16 @@ void AliADv1::StepManager()
       
       
    // Get sensitive volumes id (scintillator pads)
-   static Int_t idADA = TVirtualMC::GetMC()->VolId( "ADApad" );
-   static Int_t idADC = TVirtualMC::GetMC()->VolId( "ADCpad" );
+   static Int_t idADA = gMC->VolId( "ADApad" );
+   static Int_t idADC = gMC->VolId( "ADCpad" );
    
    // We keep only charged tracks : 
-   // if ( !TVirtualMC::GetMC()->TrackCharge() || !TVirtualMC::GetMC()->IsTrackAlive() ) return;   
+   // if ( !gMC->TrackCharge() || !gMC->IsTrackAlive() ) return;   
    // We keep charged and non-charged tracks : 
-   if ( !TVirtualMC::GetMC()->IsTrackAlive() ) return;   
+   if ( !gMC->IsTrackAlive() ) return;   
    
    Int_t copy;
-   Int_t current_volid = TVirtualMC::GetMC()->CurrentVolID( copy );
+   Int_t current_volid = gMC->CurrentVolID( copy );
 
    // check is the track is in a sensitive volume
    if( current_volid != idADA && current_volid != idADC ) {
@@ -445,7 +445,7 @@ void AliADv1::StepManager()
    // First read the position, otherwise weird reults! //ecv
    Double_t s[3];
    Float_t  x[3];
-   TVirtualMC::GetMC()->TrackPosition( s[0], s[1], s[2] );
+   gMC->TrackPosition( s[0], s[1], s[2] );
    for ( Int_t j=0; j<3; j++ ) x[j] = s[j];
    
    // Set detectro type: ADA or ADC
@@ -453,11 +453,11 @@ void AliADv1::StepManager()
    
    // Get sector copy (1,2,3,4) ( 1 level up from pad )
    Int_t sect;
-   TVirtualMC::GetMC()->CurrentVolOffID( 1, sect );
+   gMC->CurrentVolOffID( 1, sect );
 
    // Get Detector copy (1,2) ( 2 levels up from pad )
    Int_t detc;
-   TVirtualMC::GetMC()->CurrentVolOffID( 2, detc );
+   gMC->CurrentVolOffID( 2, detc );
    
    // Sector number 
    // ADA1 = 10-14
@@ -477,8 +477,8 @@ void AliADv1::StepManager()
       photoCathodeEfficiency = fADAPhotoCathodeEfficiency;
    }
       
-   Float_t destep_ad = TVirtualMC::GetMC()->Edep();
-   Float_t step_ad   = TVirtualMC::GetMC()->TrackStep();
+   Float_t destep_ad = gMC->Edep();
+   Float_t step_ad   = gMC->TrackStep();
    Int_t  nPhotonsInStep_ad = Int_t( destep_ad / (lightYield_ad * 1e-9) ); 
    nPhotonsInStep_ad = gRandom->Poisson( nPhotonsInStep_ad );
    
@@ -491,10 +491,10 @@ void AliADv1::StepManager()
    eloss_ad   += destep_ad;
    tlength_ad += step_ad;  
  
-   if ( TVirtualMC::GetMC()->IsTrackEntering() ) { 
+   if ( gMC->IsTrackEntering() ) { 
       nPhotons_ad = nPhotonsInStep_ad;
       Double_t p[4];
-      TVirtualMC::GetMC()->TrackMomentum( p[0], p[1], p[2], p[3] );
+      gMC->TrackMomentum( p[0], p[1], p[2], p[3] );
       Float_t pt  = TMath::Sqrt( p[0]*p[0] + p[1]*p[1] + p[2]*p[2] ); 
       TParticle *par = gAlice->GetMCApp()->Particle(gAlice->GetMCApp()->GetCurrentTrackNumber());
       Int_t imo = par->GetFirstMother();
@@ -517,7 +517,7 @@ void AliADv1::StepManager()
       hits_ad[5]  = p[0];     // Px
       hits_ad[6]  = p[1];     // Py
       hits_ad[7]  = p[2];     // Pz
-      hits_ad[8]  = 1.0e09*TVirtualMC::GetMC()->TrackTime(); // in ns!
+      hits_ad[8]  = 1.0e09*gMC->TrackTime(); // in ns!
   
       tlength_ad = 0.0;
       eloss_ad   = 0.0; 
@@ -527,7 +527,7 @@ void AliADv1::StepManager()
    
    nPhotons_ad += nPhotonsInStep_ad;
 
-   if( TVirtualMC::GetMC()->IsTrackExiting() || TVirtualMC::GetMC()->IsTrackStop() || TVirtualMC::GetMC()->IsTrackDisappeared() ) {
+   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared() ) {
 
       // Set integer values
       vol_ad[3]  = nPhotons_ad;
@@ -568,7 +568,7 @@ void AliADv1::StepManager()
    }
        
    //   Do we need track reference ????
-   // if( TVirtualMC::GetMC()->IsTrackEntering() || TVirtualMC::GetMC()->IsTrackExiting() ) {
+   // if( gMC->IsTrackEntering() || gMC->IsTrackExiting() ) {
    //    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), 49);
    // }
 }
@@ -579,12 +579,6 @@ void AliADv1::AddHit(Int_t track, Int_t *vol, Float_t *hits)
        new(lhits[fNhits++]) AliADhit(fIshunt,track,vol,hits);
 }
 //_________________________________________________________
-void AliADv1::AddDigits(Int_t* track, Int_t module, Float_t time)
-{
-       TClonesArray &ldigits = *fDigits;
-       new(ldigits[fNdigits++]) AliADdigit(track,module,time);
-}
-//_________________________________________________________
 void AliADv1::MakeBranch(Option_t *option)
 {
 
index e9dd293..509d616 100644 (file)
@@ -33,7 +33,6 @@ public:
   virtual      TString Version() { return TString("v1"); }
   virtual       Int_t  IsVersion() const { return 1; }
   virtual      void   AddHit(Int_t track, Int_t *vol, Float_t *hits);
-  virtual      void   AddDigits(Int_t* track, Int_t module, Float_t time);
   virtual      void   MakeBranch(Option_t *option);
   virtual       void   CreateGeometry();
   virtual      void   Init();