]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDdigitizer.cxx
Add online gain tables (Annika)
[u/mrichter/AliRoot.git] / TRD / AliTRDdigitizer.cxx
index 0df0e23016dd2cc7856c71ae2cf34d609a62fab4..8ca1249e15e6c485f87d5cf688e935680bb1f27e 100644 (file)
 
 /* $Id$ */
 
-///////////////////////////////////////////////////////////////////////////////
-//                                                                           //
-//  Creates and handles digits from TRD hits                                 //
-//  Author: C. Blume (C.Blume@gsi.de)                                        //
-//                                                                           //
-//  The following effects are included:                                      //
-//      - Diffusion                                                          //
-//      - ExB effects                                                        //
-//      - Gas gain including fluctuations                                    //
-//      - Pad-response (simple Gaussian approximation)                       //
-//      - Time-response                                                      //
-//      - Electronics noise                                                  //
-//      - Electronics gain                                                   //
-//      - Digitization                                                       //
-//      - ADC threshold                                                      //
-//  The corresponding parameter can be adjusted via the various              //
-//  Set-functions. If these parameters are not explicitly set, default       //
-//  values are used (see Init-function).                                     //
-//  As an example on how to use this class to produce digits from hits       //
-//  have a look at the macro hits2digits.C                                   //
-//  The production of summable digits is demonstrated in hits2sdigits.C      //
-//  and the subsequent conversion of the s-digits into normal digits is      //
-//  explained in sdigits2digits.C.                                           //
-//                                                                           //
-///////////////////////////////////////////////////////////////////////////////
-
-#include <stdlib.h>
-
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  Creates and handles digits from TRD hits                              //
+//                                                                        //
+//  Authors: C. Blume (blume@ikf.uni-frankfurt.de)                        //
+//           C. Lippmann                                                  //
+//           B. Vulpescu                                                  //
+//                                                                        //
+//  The following effects are included:                                   //
+//      - Diffusion                                                       //
+//      - ExB effects                                                     //
+//      - Gas gain including fluctuations                                 //
+//      - Pad-response (simple Gaussian approximation)                    //
+//      - Time-response                                                   //
+//      - Electronics noise                                               //
+//      - Electronics gain                                                //
+//      - Digitization                                                    //
+//      - Zero suppression                                                //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include <TGeoManager.h>
+#include <TList.h>
 #include <TMath.h>
-#include <TVector.h>
 #include <TRandom.h>
-#include <TROOT.h>
 #include <TTree.h>
-#include <TFile.h>
-#include <TF1.h>
-#include <TList.h>
-#include <TTask.h>
 
 #include "AliRun.h"
+#include "AliMC.h"
 #include "AliRunLoader.h"
 #include "AliLoader.h"
 #include "AliConfig.h"
-#include "AliMagF.h"
 #include "AliRunDigitizer.h"
 #include "AliRunLoader.h"
 #include "AliLoader.h"
+#include "AliLog.h"
 
 #include "AliTRD.h"
 #include "AliTRDhit.h"
 #include "AliTRDdigitizer.h"
-#include "AliTRDdataArrayI.h"
-#include "AliTRDdataArrayF.h"
-#include "AliTRDsegmentArray.h"
+#include "AliTRDarrayDictionary.h"
+#include "AliTRDarrayADC.h"
+#include "AliTRDarraySignal.h"
 #include "AliTRDdigitsManager.h"
 #include "AliTRDgeometry.h"
-#include "AliTRDparameter.h"
+#include "AliTRDpadPlane.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDSimParam.h"
+#include "AliTRDCommonParam.h"
+#include "AliTRDfeeParam.h"
+#include "AliTRDmcmSim.h"
+#include "AliTRDdigitsParam.h"
+
+#include "Cal/AliTRDCalROC.h"
+#include "Cal/AliTRDCalDet.h"
+#include "Cal/AliTRDCalOnlineGainTableROC.h"
 
 ClassImp(AliTRDdigitizer)
 
 //_____________________________________________________________________________
 AliTRDdigitizer::AliTRDdigitizer()
+  :AliDigitizer()
+  ,fRunLoader(0)
+  ,fDigitsManager(0)
+  ,fSDigitsManager(0)
+  ,fSDigitsManagerList(0)
+  ,fTRD(0)
+  ,fGeo(0)
+  ,fMcmSim(new AliTRDmcmSim)
+  ,fEvent(0)
+  ,fMasks(0)
+  ,fCompress(kTRUE)
+  ,fSDigits(kFALSE)
+  ,fMergeSignalOnly(kFALSE)
 {
   //
   // AliTRDdigitizer default constructor
   //
-
-  fDigitsManager      = 0;
-  fSDigitsManagerList = 0;
-  fSDigitsManager     = 0;
-  fTRD                = 0;
-  fGeo                = 0;
-  fPar                = 0;
-  fMasks              = 0;
-  fEvent              = 0;
-  fCompress           = kTRUE;
-  fDebug              = 0;
-  fSDigits            = kFALSE;
-  fSDigitsScale       = 0.0;
-  fMergeSignalOnly    = kFALSE;
-  fSimpleSim          = kFALSE;
-  fSimpleDet          = 0;
+  
 }
 
 //_____________________________________________________________________________
-AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
-                :AliDigitizer(name,title)
+AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)              
+  :AliDigitizer(name,title)
+  ,fRunLoader(0)
+  ,fDigitsManager(0)
+  ,fSDigitsManager(0)
+  ,fSDigitsManagerList(0)
+  ,fTRD(0)
+  ,fGeo(0)
+  ,fMcmSim(new AliTRDmcmSim)
+  ,fEvent(0)
+  ,fMasks(0)
+  ,fCompress(kTRUE)
+  ,fSDigits(kFALSE)
+  ,fMergeSignalOnly(kFALSE)
 {
   //
   // AliTRDdigitizer constructor
   //
 
-  fDigitsManager      = NULL;
-  fSDigitsManager     = NULL;
-  fSDigitsManagerList = NULL;
-  fTRD                = NULL;
-  fGeo                = NULL;
-
-  //NewIO: These data members probably are not needed anymore
-  fDigitsManager      = 0;
-  fSDigitsManagerList = 0;
-  fSDigitsManager     = 0;
-  fTRD                = 0;
-  fGeo                = 0;
-  fPar                = 0;
-  //End NewIO comment
-  fMasks              = 0;
-  fEvent              = 0;
-  fCompress           = kTRUE;
-  fDebug              = 0;
-  fSDigits            = kFALSE;
-  fMergeSignalOnly    = kFALSE;
-  fSimpleSim          = kFALSE;
-  fSimpleDet          = 0;
-  // For the summable digits
-  fSDigitsScale       = 100.;
-
 }
 
 //_____________________________________________________________________________
 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
-                                , const Text_t *name, const Text_t *title)
-                :AliDigitizer(manager,name,title)
+                               , const Text_t *name, const Text_t *title)
+  :AliDigitizer(manager,name,title)
+  ,fRunLoader(0)
+  ,fDigitsManager(0)
+  ,fSDigitsManager(0)
+  ,fSDigitsManagerList(0)
+  ,fTRD(0)
+  ,fGeo(0)
+  ,fMcmSim(new AliTRDmcmSim)
+  ,fEvent(0)
+  ,fMasks(0)
+  ,fCompress(kTRUE)
+  ,fSDigits(kFALSE)
+  ,fMergeSignalOnly(kFALSE)
 {
   //
   // AliTRDdigitizer constructor
   //
 
-  fDigitsManager      = 0;
-  fSDigitsManagerList = 0;
-  fSDigitsManager     = 0;
-  fTRD                = 0;
-  fGeo                = 0;
-  fPar                = 0;
-  fMasks              = 0;
-  fEvent              = 0;
-  fCompress           = kTRUE;
-  fDebug              = 0;
-  fSDigits            = kFALSE;
-  fMergeSignalOnly    = kFALSE;
-  fSimpleSim          = kFALSE;
-  fSimpleDet          = 0;
-  // For the summable digits
-  fSDigitsScale       = 100.;
-
 }
 
 //_____________________________________________________________________________
 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
-                :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
+  :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
+  ,fRunLoader(0)
+  ,fDigitsManager(0)
+  ,fSDigitsManager(0)
+  ,fSDigitsManagerList(0)
+  ,fTRD(0)
+  ,fGeo(0)
+  ,fMcmSim(new AliTRDmcmSim)
+  ,fEvent(0)
+  ,fMasks(0)
+  ,fCompress(kTRUE)
+  ,fSDigits(kFALSE)
+  ,fMergeSignalOnly(kFALSE)
 {
   //
   // AliTRDdigitizer constructor
   //
 
-
-
-  fDigitsManager      = 0;
-  fSDigitsManagerList = 0;
-  fSDigitsManager     = 0;
-  fTRD                = 0;
-  fGeo                = 0;
-  fPar                = 0;
-
-  fMasks              = 0;
-  fEvent              = 0;
-  fCompress           = kTRUE;
-  fDebug              = 0;
-  fSDigits            = kFALSE;
-  fMergeSignalOnly    = kFALSE;
-  fSimpleSim          = kFALSE;
-  fSimpleDet          = 0;
-
-  // For the summable digits
-  fSDigitsScale       = 100.;
-
 }
 
 //_____________________________________________________________________________
-AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d):AliDigitizer(d)
+AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
+  :AliDigitizer(d)
+  ,fRunLoader(0)
+  ,fDigitsManager(0)
+  ,fSDigitsManager(0)
+  ,fSDigitsManagerList(0)
+  ,fTRD(0)
+  ,fGeo(0)
+  ,fMcmSim(new AliTRDmcmSim)
+  ,fEvent(0)
+  ,fMasks(0)
+  ,fCompress(d.fCompress)
+  ,fSDigits(d.fSDigits)
+  ,fMergeSignalOnly(d.fMergeSignalOnly)
 {
   //
   // AliTRDdigitizer copy constructor
   //
 
-  ((AliTRDdigitizer &) d).Copy(*this);
-
 }
 
 //_____________________________________________________________________________
@@ -213,25 +192,35 @@ AliTRDdigitizer::~AliTRDdigitizer()
   // AliTRDdigitizer destructor
   //
 
-
   if (fDigitsManager) {
     delete fDigitsManager;
-    fDigitsManager = 0;
+    fDigitsManager      = 0;
   }
 
   if (fSDigitsManager) {
-    delete fSDigitsManager;
-    fSDigitsManager = 0;
+    // s-digitsmanager will be deleted via list
+    fSDigitsManager     = 0;
   }
 
   if (fSDigitsManagerList) {
+    fSDigitsManagerList->Delete();
     delete fSDigitsManagerList;
     fSDigitsManagerList = 0;
   }
 
   if (fMasks) {
     delete [] fMasks;
-    fMasks = 0;
+    fMasks       = 0;
+  }
+
+  if (fMcmSim) {
+    delete fMcmSim;
+    fMcmSim = 0;
+  }
+
+  if (fGeo) {
+    delete fGeo;
+    fGeo = 0;
   }
 
 }
@@ -243,38 +232,37 @@ AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
   // Assignment operator
   //
 
-  if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
+  if (this != &d) {
+    ((AliTRDdigitizer &) d).Copy(*this);
+  }
+
   return *this;
 
 }
 
 //_____________________________________________________________________________
-void AliTRDdigitizer::Copy(TObject &d)
+void AliTRDdigitizer::Copy(TObject &d) const
 {
   //
   // Copy function
   //
 
-  ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
-  ((AliTRDdigitizer &) d).fSDigitsManager     = 0;
+  ((AliTRDdigitizer &) d).fRunLoader          = 0;
   ((AliTRDdigitizer &) d).fDigitsManager      = 0;
+  ((AliTRDdigitizer &) d).fSDigitsManager     = 0;
+  ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
   ((AliTRDdigitizer &) d).fTRD                = 0;
   ((AliTRDdigitizer &) d).fGeo                = 0;
-  ((AliTRDdigitizer &) d).fMasks              = 0;
   ((AliTRDdigitizer &) d).fEvent              = 0;
-  ((AliTRDdigitizer &) d).fPar                = 0;
+  ((AliTRDdigitizer &) d).fMasks              = 0;
   ((AliTRDdigitizer &) d).fCompress           = fCompress;
-  ((AliTRDdigitizer &) d).fDebug              = fDebug  ;
   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
-  ((AliTRDdigitizer &) d).fSDigitsScale       = fSDigitsScale;
   ((AliTRDdigitizer &) d).fMergeSignalOnly    = fMergeSignalOnly;
-  ((AliTRDdigitizer &) d).fSimpleSim          = fSimpleSim;
-  ((AliTRDdigitizer &) d).fSimpleDet          = fSimpleDet;
-                                       
+
 }
 
 //_____________________________________________________________________________
-void AliTRDdigitizer::Exec(Option_t* option)
+void AliTRDdigitizer::Exec(const Option_t * const option)
 {
   //
   // Executes the merging
@@ -286,115 +274,120 @@ void AliTRDdigitizer::Exec(Option_t* option)
 
   TString optionString = option;
   if (optionString.Contains("deb")) {
-    fDebug = 1;
-    if (optionString.Contains("2")) {
-      fDebug = 2;
-    }
-    printf("<AliTRDdigitizer::Exec> ");
-    printf("Called with debug option %d\n",fDebug);
+    AliLog::SetClassDebugLevel("AliTRDdigitizer",1);
+    AliInfo("Called with debug option");
   }
 
   // The AliRoot file is already connected by the manager
-  AliRunLoader* inrl;
+  AliRunLoader *inrl = 0x0;
   
-  if (gAlice) 
-   {
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::Exec> ");
-      printf("AliRun object found on file.\n");
-    }
-   }
+  if (gAlice) {
+    AliDebug(1,"AliRun object found on file.");
+  }
   else {
     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
     inrl->LoadgAlice();
     gAlice = inrl->GetAliRun();
-    if (!gAlice)
-     {
-       printf("<AliTRDdigitizer::Exec> ");
-       printf("Could not find AliRun object.\n");
-       return;
-     }
+    if (!gAlice) {
+      AliError("Could not find AliRun object.");
+      return;
+    }
   }
                                                                            
   Int_t nInput = fManager->GetNinputs();
-  fMasks = new Int_t[nInput];
+  fMasks       = new Int_t[nInput];
   for (iInput = 0; iInput < nInput; iInput++) {
     fMasks[iInput] = fManager->GetMask(iInput);
   }
 
+  //
   // Initialization
-  InitDetector();
+  //
 
-  for (iInput = 0; iInput < nInput; iInput++) {
+  AliRunLoader *orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+
+  if (InitDetector()) {
+
+    AliLoader *ogime = orl->GetLoader("TRDLoader");
 
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::Exec> ");
-      printf("Add input stream %d\n",iInput);
+    TTree *tree = 0;
+    if (fSDigits) { 
+      // If we produce SDigits
+      tree = ogime->TreeS();
+      if (!tree) {
+       ogime->MakeTree("S");
+       tree = ogime->TreeS();
+      }
+    }
+    else {
+      // If we produce Digits
+      tree = ogime->TreeD();
+      if (!tree) {
+       ogime->MakeTree("D");
+        tree = ogime->TreeD();
+      }
     }
 
-    // check if the input tree exists
+    MakeBranch(tree);
+
+  }
+  for (iInput = 0; iInput < nInput; iInput++) {
+
+    AliDebug(1,Form("Add input stream %d",iInput));
+
+    // Check if the input tree exists
     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
-    AliLoader* gime = inrl->GetLoader("TRDLoader");
-
-    TTree * treees =  gime->TreeS();
-    if (treees == 0x0) 
-     {
-      if (gime->LoadSDigits())
-       {
-         Error("Exec","Error Occured while loading S. Digits for input %d.",iInput);
-         return;
-       }
-      treees =  gime->TreeS();
-     }
+    AliLoader *gime = inrl->GetLoader("TRDLoader");
+
+    TTree *treees = gime->TreeS();
+    if (treees == 0x0) {
+      if (gime->LoadSDigits()) {
+        AliError(Form("Error Occured while loading S. Digits for input %d.",iInput));
+        return;
+      }
+      treees = gime->TreeS();
+    }
     
     if (treees == 0x0) {
-      printf("<AliTRDdigitizer::Exec> ");
-      printf("Input stream %d does not exist\n",iInput);
+      AliError(Form("Input stream %d does not exist",iInput));
       return;
     } 
 
     // Read the s-digits via digits manager
     sdigitsManager = new AliTRDdigitsManager();
-    sdigitsManager->SetDebug(fDebug);
     sdigitsManager->SetSDigits(kTRUE);
     
-    AliRunLoader* rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
-    AliLoader* gimme = rl->GetLoader("TRDLoader");
-    if (!gimme->TreeS()) gimme->LoadSDigits();
-    sdigitsManager->ReadDigits(gimme->TreeS());
+    AliRunLoader *rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
+    AliLoader *gimme = rl->GetLoader("TRDLoader");
+    if (!gimme->TreeS()) 
+      {
+       gimme->LoadSDigits();
+      }
 
+    sdigitsManager->ReadDigits(gimme->TreeS());
+   
     // Add the s-digits to the input list 
     AddSDigitsManager(sdigitsManager);
 
   }
 
   // Convert the s-digits to normal digits
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::Exec> ");
-    printf("Do the conversion\n");
-  }
+  AliDebug(1,"Do the conversion");
   SDigits2Digits();
 
   // Store the digits
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::Exec> ");
-    printf("Write the digits\n");
-  }
-  
-  AliRunLoader* orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
-  AliLoader* ogime = orl->GetLoader("TRDLoader");
-
-  fDigitsManager->MakeBranch(ogime->TreeD());
-
+  AliDebug(1,"Write the digits");
   fDigitsManager->WriteDigits();
 
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::Exec> ");
-    printf("Done\n");
-  }
+  // Write parameters
+  orl->CdGAFile();
 
+  // Clean up
   DeleteSDigitsManager();
 
+  AliDebug(1,"Done");
+
 }
 
 //_____________________________________________________________________________
@@ -403,65 +396,121 @@ Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
   //
   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
   //
-
   // Connect the AliRoot file containing Geometry, Kine, and Hits
+  //  
+
+  TString evfoldname = AliConfig::GetDefaultEventFolderName();
+
+  fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
+  if (!fRunLoader) {
+    fRunLoader = AliRunLoader::Open(file,evfoldname,"UPDATE");
+  }  
+  if (!fRunLoader) {
+    AliError(Form("Can not open session for file %s.",file));
+    return kFALSE;
+  }
+   
+  if (!fRunLoader->GetAliRun()) {
+    fRunLoader->LoadgAlice();
+  }
+  gAlice = fRunLoader->GetAliRun();
   
-  fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,"UPDATE");
+  if (gAlice) {
+    AliDebug(1,"AliRun object found on file.");
+  }
+  else {
+    AliError("Could not find AliRun object.");
+    return kFALSE;
+  }
+
+  fEvent = nEvent;
+
+  AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
+  if (!loader) {
+    AliError("Can not get TRD loader from Run Loader");
+    return kFALSE;
+  }
   
-  if (!fRunLoader)
-   {
-     Error("Open","Can not open session for file %s.",file);
-     return kFALSE;
-   }
+  if (InitDetector()) {
+    TTree *tree = 0;
+    if (fSDigits) { 
+      // If we produce SDigits
+      tree = loader->TreeS();
+      if (!tree) {
+        loader->MakeTree("S");
+        tree = loader->TreeS();
+      }
+    }
+    else {
+      // If we produce Digits
+      tree = loader->TreeD();
+      if (!tree) {
+        loader->MakeTree("D");
+        tree = loader->TreeD();
+      }
+    }
+    return MakeBranch(tree);
+  }
+  else {
+    return kFALSE;
+  }
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::Open(AliRunLoader * const runLoader, Int_t nEvent)
+{
+  //
+  // Opens a ROOT-file with TRD-hits and reads in the hit-tree
+  //
+  // Connect the AliRoot file containing Geometry, Kine, and Hits
+  //  
+
+  fRunLoader = runLoader;
+  if (!fRunLoader) {
+    AliError("RunLoader does not exist");
+    return kFALSE;
+  }
    
-  fRunLoader->LoadgAlice();
+  if (!fRunLoader->GetAliRun()) {
+    fRunLoader->LoadgAlice();
+  }
   gAlice = fRunLoader->GetAliRun();
   
   if (gAlice) {
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::Open> ");
-      printf("AliRun object found on file.\n");
-    }
+    AliDebug(1,"AliRun object found on file.");
   }
   else {
-    printf("<AliTRDdigitizer::Open> ");
-    printf("Could not find AliRun object.\n");
+    AliError("Could not find AliRun object.");
     return kFALSE;
   }
 
   fEvent = nEvent;
 
-  // Import the Trees for the event nEvent in the file
-  fRunLoader->GetEvent(fEvent);
-  
-  AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
-  if (!loader)
-   {
-     Error("Open","Can not get TRD loader from Run Loader");
-     return kFALSE;
-   }
+  AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
+  if (!loader) {
+    AliError("Can not get TRD loader from Run Loader");
+    return kFALSE;
+  }
   
   if (InitDetector()) {
-    TTree* tree = 0;
-    if (fSDigits)
-     { 
-     //if we produce SDigits
-       tree = loader->TreeS();
-       if (!tree)
-        {
-         loader->MakeTree("S");
-         tree = loader->TreeS();
-        }
-     }
-    else
-     {//if we produce Digits
-       tree = loader->TreeD();
-       if (!tree)
-        {
-         loader->MakeTree("D");
-         tree = loader->TreeD();
-        }
-     }
+    TTree *tree = 0;
+    if (fSDigits) { 
+      // If we produce SDigits
+      tree = loader->TreeS();
+      if (!tree) {
+        loader->MakeTree("S");
+        tree = loader->TreeS();
+      }
+    }
+    else {
+      // If we produce Digits
+      tree = loader->TreeD();
+      if (!tree) {
+        loader->MakeTree("D");
+        tree = loader->TreeD();
+      }
+    }
     return MakeBranch(tree);
   }
   else {
@@ -480,39 +529,40 @@ Bool_t AliTRDdigitizer::InitDetector()
   // Get the pointer to the detector class and check for version 1
   fTRD = (AliTRD *) gAlice->GetDetector("TRD");
   if (!fTRD) {
-    printf("<AliTRDdigitizer::InitDetector> ");
-    printf("No TRD module found\n");
+    AliFatal("No TRD module found");
     exit(1);
   }
   if (fTRD->IsVersion() != 1) {
-    printf("<AliTRDdigitizer::InitDetector> ");
-    printf("TRD must be version 1 (slow simulator).\n");
+    AliFatal("TRD must be version 1 (slow simulator)");
     exit(1);
   }
 
   // Get the geometry
-  fGeo = fTRD->GetGeometry();
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::InitDetector> ");
-    printf("Geometry version %d\n",fGeo->IsVersion());
-  }
+  fGeo = new AliTRDgeometry();
 
   // Create a digits manager
+  if (fDigitsManager) {
+    delete fDigitsManager;
+  }
   fDigitsManager = new AliTRDdigitsManager();
   fDigitsManager->SetSDigits(fSDigits);
   fDigitsManager->CreateArrays();
   fDigitsManager->SetEvent(fEvent);
-  fDigitsManager->SetDebug(fDebug);
 
   // The list for the input s-digits manager to be merged
-  fSDigitsManagerList = new TList();
+  if (fSDigitsManagerList) {
+    fSDigitsManagerList->Delete();
+  } 
+  else {
+    fSDigitsManagerList = new TList();
+  }
 
   return kTRUE;
 
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDdigitizer::MakeBranch(TTreetree) const
+Bool_t AliTRDdigitizer::MakeBranch(TTree *tree) const
 {
   // 
   // Create the branches for the digits array
@@ -522,6 +572,28 @@ Bool_t AliTRDdigitizer::MakeBranch(TTree* tree) const
 
 }
 
+//_____________________________________________________________________________
+void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
+{
+  //
+  // Add a digits manager for s-digits to the input list.
+  //
+
+  fSDigitsManagerList->Add(man);
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdigitizer::DeleteSDigitsManager()
+{
+  //
+  // Removes digits manager from the input list.
+  //
+
+  fSDigitsManagerList->Delete();
+
+}
+
 //_____________________________________________________________________________
 Bool_t AliTRDdigitizer::MakeDigits()
 {
@@ -529,755 +601,948 @@ Bool_t AliTRDdigitizer::MakeDigits()
   // Creates digits.
   //
 
-  ///////////////////////////////////////////////////////////////
-  // Parameter 
-  ///////////////////////////////////////////////////////////////
+  AliDebug(1,"Start creating digits");
 
-  // Converts number of electrons to fC
-  const Double_t kEl2fC  = 1.602E-19 * 1.0E15; 
+  if (!fGeo) {
+    AliError("No geometry defined");
+    return kFALSE;
+  }
 
-  ///////////////////////////////////////////////////////////////
+  AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliFatal("Could not get calibration object");
+    return kFALSE;
+  }
 
-  // Number of pads included in the pad response
-  const Int_t kNpad  = 3;
+  const Int_t kNdet  = AliTRDgeometry::Ndet();
 
-  // Number of track dictionary arrays
-  const Int_t kNDict = AliTRDdigitsManager::kNDict;
+  Float_t **hits = new Float_t*[kNdet];
+  Int_t    *nhit = new Int_t[kNdet];
+
+  AliTRDarraySignal *signals = 0x0;
 
-  // Half the width of the amplification region
-  const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
-
-  Int_t   iRow, iCol, iTime, iPad;
-  Int_t   iDict  = 0;
-  Int_t   nBytes = 0;
-
-  Int_t   totalSizeDigits = 0;
-  Int_t   totalSizeDict0  = 0;
-  Int_t   totalSizeDict1  = 0;
-  Int_t   totalSizeDict2  = 0;
-
-  Int_t   timeTRDbeg = 0;
-  Int_t   timeTRDend = 1;
-
-  Float_t pos[3];
-  Float_t rot[3];
-  Float_t xyz[3];
-  Float_t padSignal[kNpad];
-  Float_t signalOld[kNpad];
-
-  AliTRDdataArrayF *signals = 0;
-  AliTRDdataArrayI *digits  = 0;
-  AliTRDdataArrayI *dictionary[kNDict];
-
-  // Create a default parameter class if none is defined
-  if (!fPar) {
-    fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::MakeDigits> ");
-      printf("Create the default parameter object\n");
+  // Check the number of time bins from simParam against OCDB,
+  // if OCDB value is not supposed to be used.
+  // As default, the value from OCDB is taken
+  if (AliTRDSimParam::Instance()->GetNTBoverwriteOCDB()) {
+    if (calibration->GetNumberOfTimeBinsDCS() != AliTRDSimParam::Instance()->GetNTimeBins()) {
+      AliWarning(Form("Number of time bins is different to OCDB value [SIM=%d, OCDB=%d]"
+                     ,AliTRDSimParam::Instance()->GetNTimeBins()
+                     ,calibration->GetNumberOfTimeBinsDCS()));
     }
+    // Save the values for the raw data headers
+    fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(AliTRDSimParam::Instance()->GetNTimeBins());
+  }
+  else {
+    // Save the values for the raw data headers
+    fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(calibration->GetNumberOfTimeBinsDCS());
+  }
+
+  // Save the values for the raw data headers
+  fDigitsManager->GetDigitsParam()->SetADCbaselineAll(AliTRDSimParam::Instance()->GetADCbaseline());
+  // Sort all hits according to detector number
+  if (!SortHits(hits,nhit)) {
+    AliError("Sorting hits failed");
+    delete [] hits;
+    delete [] nhit;
+    return kFALSE;
   }
 
-  // Create a container for the amplitudes
-  AliTRDsegmentArray *signalsArray 
-                     = new AliTRDsegmentArray("AliTRDdataArrayF"
-                                             ,AliTRDgeometry::Ndet());
+  // Loop through all detectors
+  for (Int_t det = 0; det < kNdet; det++) {
+
+    // Detectors that are switched off, not installed, etc.
+    if (( calibration->IsChamberInstalled(det)) &&
+       (!calibration->IsChamberMasked(det))    &&
+        ( fGeo->ChamberInGeometry(det))         &&
+        (nhit[det] > 0)) {
+
+      signals = new AliTRDarraySignal();
+         
+      // Convert the hits of the current detector to detector signals
+      if (!ConvertHits(det,hits[det],nhit[det],signals)) {
+       AliError(Form("Conversion of hits failed for detector=%d",det));
+        delete [] hits;
+        delete [] nhit;
+        delete signals;
+        signals = 0x0;
+        return kFALSE;
+      }
+
+      // Convert the detector signals to digits or s-digits
+      if (!ConvertSignals(det,signals)) {
+       AliError(Form("Conversion of signals failed for detector=%d",det));
+        delete [] hits;
+        delete [] nhit;
+        delete signals;
+        signals = 0x0;
+       return kFALSE;
+      }
+
+      // Delete the signals array
+      delete signals;
+      signals = 0x0;
 
-  if (fPar->TRFOn()) {
-    timeTRDbeg = ((Int_t) (-fPar->GetTRFlo() / fPar->GetTimeBinSize())) - 1;
-    timeTRDend = ((Int_t) ( fPar->GetTRFhi() / fPar->GetTimeBinSize())) - 1;
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::MakeDigits> ");
-      printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
+    } // if: detector status
+
+    delete [] hits[det];
+
+  } // for: detector
+
+  if (!fSDigits) {
+    if (AliDataLoader *trklLoader 
+          = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
+      if (trklLoader->Tree())
+        trklLoader->WriteData("OVERWRITE");
     }
   }
+    
+  delete [] hits;
+  delete [] nhit;
 
-  Float_t elAttachProp = fPar->GetElAttachProp() / 100.; 
+  return kTRUE;
 
-  if (!fGeo) {
-    printf("<AliTRDdigitizer::MakeDigits> ");
-    printf("No geometry defined\n");
-    return kFALSE;
-  }
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::SortHits(Float_t **hits, Int_t *nhit)
+{
+  //
+  // Read all the hits and sorts them according to detector number
+  // in the output array <hits>.
+  //
+
+  AliDebug(1,"Start sorting hits");
+
+  const Int_t kNdet = AliTRDgeometry::Ndet();
+  // Size of the hit vector
+  const Int_t kNhit = 6;
 
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::MakeDigits> ");
-    printf("Start creating digits.\n");
+  Float_t *xyz      = 0;
+  Int_t    nhitTrk  = 0;
+
+  Int_t   *lhit     = new Int_t[kNdet];
+
+  for (Int_t det = 0; det < kNdet; det++) {
+    lhit[det] = 0;
+    nhit[det] = 0;
+    hits[det] = 0;
   }
 
-  AliLoader* gimme = fRunLoader->GetLoader("TRDLoader");
-  if (!gimme->TreeH()) gimme->LoadHits();
-  TTree* hitTree = gimme->TreeH();
-  if (hitTree == 0x0)
-    {
-      Error("MakeDigits","Can not get TreeH");
-      return kFALSE;
-    }
+  AliLoader *gimme = fRunLoader->GetLoader("TRDLoader");
+  if (!gimme->TreeH()) {
+    gimme->LoadHits();
+  }
+  TTree *hitTree = gimme->TreeH();
+  if (hitTree == 0x0) {
+    AliError("Can not get TreeH");
+    delete [] lhit;
+    return kFALSE;
+  }
   fTRD->SetTreeAddress();
-  
+
   // Get the number of entries in the hit tree
   // (Number of primary particles creating a hit somewhere)
-  Int_t nTrack = 1;
-  if (!fSimpleSim) {
-    nTrack = (Int_t) hitTree->GetEntries();
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::MakeDigits> ");
-      printf("Found %d primary particles\n",nTrack);
-    } 
-  }
+  Int_t nTrk = (Int_t) hitTree->GetEntries();
+  AliDebug(1,Form("Found %d tracks",nTrk));
 
-  Int_t detectorOld = -1;
-  Int_t countHits   =  0; 
+  // Loop through all the tracks in the tree
+  for (Int_t iTrk = 0; iTrk < nTrk; iTrk++) {
 
-  // Loop through all entries in the tree
-  for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+    gAlice->GetMCApp()->ResetHits();
+    hitTree->GetEvent(iTrk);
 
-    if (!fSimpleSim) {   
-      gAlice->ResetHits();
-      nBytes += hitTree->GetEvent(iTrack);
+    if (!fTRD->Hits()) {
+      AliError(Form("No hits array for track = %d",iTrk));
+      continue;
     }
 
+    // Number of hits for this track
+    nhitTrk = fTRD->Hits()->GetEntriesFast();
+
+    Int_t hitCnt = 0;
     // Loop through the TRD hits
-    Int_t iHit = 0;
     AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
     while (hit) {
-      countHits++;
-      iHit++;
-
-              pos[0]      = hit->X();
-              pos[1]      = hit->Y();
-              pos[2]      = hit->Z();
-      Float_t q           = hit->GetCharge();
-      Int_t   track       = hit->Track();
-      Int_t   detector    = hit->GetDetector();
-      Int_t   plane       = fGeo->GetPlane(detector);
-      Int_t   sector      = fGeo->GetSector(detector);
-      Int_t   chamber     = fGeo->GetChamber(detector);
-      Int_t   nRowMax     = fPar->GetRowMax(plane,chamber,sector);
-      Int_t   nColMax     = fPar->GetColMax(plane);
-      Int_t   nTimeMax    = fPar->GetTimeMax();
-      Int_t   nTimeBefore = fPar->GetTimeBefore();
-      Int_t   nTimeAfter  = fPar->GetTimeAfter();
-      Int_t   nTimeTotal  = fPar->GetTimeTotal();
-      Float_t row0        = fPar->GetRow0(plane,chamber,sector);
-      Float_t col0        = fPar->GetCol0(plane);
-      Float_t time0       = fPar->GetTime0(plane);
-      Float_t rowPadSize  = fPar->GetRowPadSize(plane,chamber,sector);
-      Float_t colPadSize  = fPar->GetColPadSize(plane);
-      Float_t timeBinSize = fPar->GetTimeBinSize();
-      Float_t divideRow   = 1.0 / rowPadSize;
-      Float_t divideCol   = 1.0 / colPadSize;
-      Float_t divideTime  = 1.0 / timeBinSize;
-
-      if (fDebug > 1) {
-        printf("Analyze hit no. %d ",iHit);
-        printf("-----------------------------------------------------------\n");
-        hit->Dump();
-        printf("plane = %d, sector = %d, chamber = %d\n"
-              ,plane,sector,chamber);
-        printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n" 
-              ,nRowMax,nColMax,nTimeMax);
-        printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
-             ,nTimeBefore,nTimeAfter,nTimeTotal);
-        printf("row0 = %f, col0 = %f, time0 = %f\n"
-              ,row0,col0,time0);
-        printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
-              ,rowPadSize,colPadSize,timeBinSize); 
-      }
-       
-      // Don't analyze test hits and switched off detectors
-      if ((CheckDetector(plane,chamber,sector)) &&
-          (((Int_t) q) != 0)) {
-
-        if (detector != detectorOld) {
 
-          if (fDebug > 1) {
-            printf("<AliTRDdigitizer::MakeDigits> ");
-            printf("Get new container. New det = %d, Old det = %d\n"
-                  ,detector,detectorOld);
-         }
-          // Compress the old one if enabled
-          if ((fCompress) && (detectorOld > -1)) {
-            if (fDebug > 1) {
-              printf("<AliTRDdigitizer::MakeDigits> ");
-              printf("Compress the old container ...");
-           }
-            signals->Compress(1,0);
-            for (iDict = 0; iDict < kNDict; iDict++) {
-              dictionary[iDict]->Compress(1,0);
-           }
-            if (fDebug > 1) printf("done\n");
-         }
-         // Get the new container
-          signals = (AliTRDdataArrayF *) signalsArray->At(detector);
-          if (signals->GetNtime() == 0) {
-            // Allocate a new one if not yet existing
-            if (fDebug > 1) {
-              printf("<AliTRDdigitizer::MakeDigits> ");
-              printf("Allocate a new container ... ");
-           }
-            signals->Allocate(nRowMax,nColMax,nTimeTotal);
-         }
-          else if (fSimpleSim) {
-            // Clear an old one for the simple simulation
-            if (fDebug > 1) {
-              printf("<AliTRDdigitizer::MakeDigits> ");
-              printf("Clear a old container ... ");
-            }
-            signals->Clear();
-          }
-          else {
-           // Expand an existing one
-            if (fCompress) {
-              if (fDebug > 1) {
-                printf("<AliTRDdigitizer::MakeDigits> ");
-                printf("Expand an existing container ... ");
-             }
-              signals->Expand();
-           }
+      hitCnt++;
+
+      // Don't analyze test hits
+      if (((Int_t) hit->GetCharge()) != 0) {
+
+        Int_t   trk  = hit->Track();
+        Int_t   det  = hit->GetDetector();
+        Int_t   q    = hit->GetCharge();
+        Float_t x    = hit->X();
+        Float_t y    = hit->Y();
+        Float_t z    = hit->Z();
+        Float_t time = hit->GetTime();
+
+        if (nhit[det] == lhit[det]) {
+          // Inititialization of new detector
+          xyz        = new Float_t[kNhit*(nhitTrk+lhit[det])];
+          if (hits[det]) {
+            memcpy(xyz,hits[det],sizeof(Float_t)*kNhit*lhit[det]);
+            delete [] hits[det];
          }
-         // The same for the dictionary
-          if (!fSimpleSim) {       
-            for (iDict = 0; iDict < kNDict; iDict++) {       
-              dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
-              if (dictionary[iDict]->GetNtime() == 0) {
-                dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
-             }
-              else {
-                if (fCompress) dictionary[iDict]->Expand();
-             }
-           }
-          }      
-          if (fDebug > 1) printf("done\n");
-          detectorOld = detector;
-        }
-
-        // Rotate the sectors on top of each other
-        if (fSimpleSim) {
-          rot[0] = pos[0];
-          rot[1] = pos[1];
-          rot[2] = pos[2];
-        }
+          lhit[det] += nhitTrk;
+          hits[det]  = xyz;
+       }
         else {
-          fGeo->Rotate(detector,pos,rot);
+          xyz        = hits[det];
        }
+        xyz[nhit[det]*kNhit+0] = x;
+        xyz[nhit[det]*kNhit+1] = y;
+        xyz[nhit[det]*kNhit+2] = z;
+        xyz[nhit[det]*kNhit+3] = q;
+        xyz[nhit[det]*kNhit+4] = trk;
+        xyz[nhit[det]*kNhit+5] = time;
+        nhit[det]++;
 
-        // The driftlength. It is negative if the hit is in the 
-        // amplification region.
-        Float_t driftlength = time0 - rot[0];
-
-        // Take also the drift in the amplification region into account
-        // The drift length is at the moment still the same, regardless of
-        // the position relativ to the wire. This non-isochronity needs still
-        // to be implemented.
-        Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
-        if (fPar->ExBOn()) driftlengthL /= TMath::Sqrt(fPar->GetLorentzFactor());
-
-        // Loop over all electrons of this hit
-        // TR photons produce hits with negative charge
-        Int_t nEl = ((Int_t) TMath::Abs(q));
-        for (Int_t iEl = 0; iEl < nEl; iEl++) {
-
-          xyz[0] = rot[0];
-          xyz[1] = rot[1];
-          xyz[2] = rot[2];
-
-          // Electron attachment
-          if (fPar->ElAttachOn()) {
-            if (gRandom->Rndm() < (driftlengthL * elAttachProp)) 
-              continue;
-          }
+      } // if: charge != 0
 
-          // Apply the diffusion smearing
-          if (fPar->DiffusionOn()) {
-            if (!(fPar->Diffusion(driftlengthL,xyz))) continue;
-         }
+      hit = (AliTRDhit *) fTRD->NextHit();   
 
-          // Apply E x B effects (depends on drift direction)
-          if (fPar->ExBOn()) { 
-            if (!(fPar->ExB(driftlength+kAmWidth,xyz))) continue;   
-         }
+    } // for: hits of one track
 
-          // The electron position after diffusion and ExB in pad coordinates 
-          // The pad row (z-direction)
-          Float_t rowDist   = xyz[2] - row0;
-          Int_t   rowE      = ((Int_t) (rowDist * divideRow));
-          if ((rowE < 0) || (rowE >= nRowMax)) continue;   
-          Float_t rowOffset = ((((Float_t) rowE) + 0.5) * rowPadSize) - rowDist;
-
-          // The pad column (rphi-direction)
-          Float_t col0tilt  = fPar->Col0Tilted(col0,rowOffset,plane);
-          Float_t colDist   = xyz[1] - col0tilt;
-          Int_t   colE      = ((Int_t) (colDist * divideCol));
-          if ((colE < 0) || (colE >= nColMax)) continue;   
-          Float_t colOffset = ((((Float_t) colE) + 0.5) * colPadSize) - colDist;    
-
-          // The time bin (negative for hits in the amplification region)
-         // In the amplification region the electrons drift from both sides
-         // to the middle (anode wire plane)
-          Float_t timeDist   = time0 - xyz[0];
-          Float_t timeOffset = 0;
-          Int_t   timeE      = 0;
-          if (timeDist > 0) {
-           // The time bin
-            timeE      = ((Int_t) (timeDist * divideTime));
-            // The distance of the position to the middle of the timebin
-            timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
-         }
-          else {
-           // Difference between half of the amplification gap width and
-           // the distance to the anode wire
-            Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
-            // The time bin
-            timeE      = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
-            // The distance of the position to the middle of the timebin
-            timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
-         }
-          // Apply the gas gain including fluctuations
-          Float_t ggRndm = 0.0;
-          do {
-            ggRndm = gRandom->Rndm();
-         } while (ggRndm <= 0);
-          Int_t signal = (Int_t) (-fPar->GetGasGain() * TMath::Log(ggRndm));
-
-          // Apply the pad response 
-          if (fPar->PRFOn()) {
-           // The distance of the electron to the center of the pad 
-           // in units of pad width
-            Float_t dist = - colOffset * divideCol;
-            if (!(fPar->PadResponse(signal,dist,plane,padSignal))) continue;
-         }
-         else {
-            padSignal[0] = 0.0;
-            padSignal[1] = signal;
-            padSignal[2] = 0.0;
-         }
+  } // for: tracks
 
-         // Sample the time response inside the drift region
-         // + additional time bins before and after.
-          // The sampling is done always in the middle of the time bin
-          for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg,        -nTimeBefore) 
-                   ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter ) 
-                   ;iTimeBin++) {
-
-           // Apply the time response
-            Float_t timeResponse = 1.0;
-            Float_t crossTalk    = 0.0;
-            Float_t time         = (iTimeBin - timeE) * timeBinSize + timeOffset;
-            if (fPar->TRFOn()) {
-              timeResponse = fPar->TimeResponse(time);
-           }
-            if (fPar->CTOn()) {
-              crossTalk    = fPar->CrossTalk(time);
-            }
+  delete [] lhit;
 
-            signalOld[0] = 0.0;
-            signalOld[1] = 0.0;
-            signalOld[2] = 0.0;
-
-            for (iPad = 0; iPad < kNpad; iPad++) {
-
-              Int_t colPos = colE + iPad - 1;
-              if (colPos <        0) continue;
-              if (colPos >= nColMax) break;
-
-              // Add the signals
-              // Note: The time bin number is shifted by nTimeBefore to avoid negative
-              // time bins. This has to be subtracted later.
-              Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
-              signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
-              if( colPos != colE ) {
-                signalOld[iPad] += padSignal[iPad] * (timeResponse + crossTalk);
-              } 
-              else {
-                signalOld[iPad] += padSignal[iPad] * timeResponse;
-              }
-              signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
-
-              // Store the track index in the dictionary
-              // Note: We store index+1 in order to allow the array to be compressed
-              if ((signalOld[iPad] > 0) && (!fSimpleSim)) { 
-                for (iDict = 0; iDict < kNDict; iDict++) {
-                  Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
-                                                                      ,colPos
-                                                                      ,iCurrentTimeBin);
-                  if (oldTrack == track+1) break;
-                  if (oldTrack ==       0) {
-                    dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
-                    break;
-                  }
-                }
-              }
-
-           } // Loop: pads
-
-         } // Loop: time bins
-
-        } // Loop: electrons of a single hit
-
-      } // If: detector and test hit
+  return kTRUE;
 
-      hit = (AliTRDhit *) fTRD->NextHit();   
+}
 
-    } // Loop: hits of one primary track
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::ConvertHits(Int_t det
+                                  , const Float_t * const hits
+                                  , Int_t nhit
+                                  , AliTRDarraySignal *signals)
+{
+  //
+  // Converts the detectorwise sorted hits to detector signals
+  //
 
-  } // Loop: primary tracks
+  AliDebug(1,Form("Start converting hits for detector=%d (nhits=%d)",det,nhit));
 
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::MakeDigits> ");
-    printf("Finished analyzing %d hits\n",countHits);
-  }
+  // Number of pads included in the pad response
+  const Int_t kNpad      = 3;
+  // Number of track dictionary arrays
+  const Int_t kNdict     = AliTRDdigitsManager::kNDict;
+  // Size of the hit vector
+  const Int_t kNhit      = 6;
+
+  // Width of the amplification region
+  const Float_t kAmWidth = AliTRDgeometry::AmThick();
+  // Width of the drift region
+  const Float_t kDrWidth = AliTRDgeometry::DrThick();
+  // Drift + amplification region 
+  const Float_t kDrMin   =          - 0.5 * kAmWidth;
+  const Float_t kDrMax   = kDrWidth + 0.5 * kAmWidth;
+  
+  Int_t    iPad          = 0;
+  Int_t    dict          = 0;
+  Int_t    timeBinTRFend = 1;
 
-  // The coupling factor
-  Float_t coupling = fPar->GetPadCoupling() 
-                   * fPar->GetTimeCoupling();
+  Double_t pos[3];
+  Double_t loc[3];
+  Double_t padSignal[kNpad];
+  Double_t signalOld[kNpad];
 
-  // The conversion factor
-  Float_t convert  = kEl2fC
-                   * fPar->GetChipGain();
+  AliTRDarrayDictionary *dictionary[kNdict];
 
-  // Loop through all chambers to finalize the digits
-  Int_t iDetBeg = 0;
-  Int_t iDetEnd = AliTRDgeometry::Ndet();
-  if (fSimpleSim) {
-    iDetBeg = fSimpleDet;
-    iDetEnd = iDetBeg + 1;
-  }
-  for (Int_t iDet = iDetBeg; iDet < iDetEnd; iDet++) {
+  AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
+  AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
+  AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
 
-    Int_t plane       = fGeo->GetPlane(iDet);
-    Int_t sector      = fGeo->GetSector(iDet);
-    Int_t chamber     = fGeo->GetChamber(iDet);
-    Int_t nRowMax     = fPar->GetRowMax(plane,chamber,sector);
-    Int_t nColMax     = fPar->GetColMax(plane);
-    Int_t nTimeMax    = fPar->GetTimeMax();
-    Int_t nTimeTotal  = fPar->GetTimeTotal();
+  if (!commonParam) {
+    AliFatal("Could not get common parameterss");
+    return kFALSE;
+  }
+  if (!simParam) {
+    AliFatal("Could not get simulation parameters");
+    return kFALSE;
+  }  
+  if (!calibration) {
+    AliFatal("Could not get calibration object");  
+    return kFALSE;
+  }
 
-    Double_t *inADC  = new Double_t[nTimeTotal];
-    Double_t *outADC = new Double_t[nTimeTotal];
+  // Get the detector wise calibration objects
+  AliTRDCalROC       *calVdriftROC      = 0;
+  Float_t             calVdriftDetValue = 0.0;
+  const AliTRDCalDet *calVdriftDet      = calibration->GetVdriftDet();  
+  AliTRDCalROC       *calT0ROC          = 0;
+  Float_t             calT0DetValue     = 0.0;
+  const AliTRDCalDet *calT0Det          = calibration->GetT0Det();  
+
+  if (simParam->TRFOn()) {
+    timeBinTRFend = ((Int_t) (simParam->GetTRFhi() 
+                  * commonParam->GetSamplingFrequency())) - 1;
+  }
 
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::MakeDigits> ");
-      printf("Digitization for chamber %d\n",iDet);
+  Int_t   nTimeTotal   = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
+  Float_t samplingRate = commonParam->GetSamplingFrequency();
+  Float_t elAttachProp = simParam->GetElAttachProp() / 100.0; 
+
+  AliTRDpadPlane *padPlane = fGeo->GetPadPlane(det);
+  Int_t   layer   = fGeo->GetLayer(det);         //update
+  Float_t row0    = padPlane->GetRow0ROC();
+  Int_t   nRowMax = padPlane->GetNrows();
+  Int_t   nColMax = padPlane->GetNcols();
+
+  // Create a new array for the signals
+  signals->Allocate(nRowMax,nColMax,nTimeTotal);
+
+  // Create a new array for the dictionary
+  for (dict = 0; dict < kNdict; dict++) {       
+    dictionary[dict] = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
+    dictionary[dict]->Allocate(nRowMax,nColMax,nTimeTotal);
+  }      
+
+  // Loop through the hits in this detector
+  for (Int_t hit = 0; hit < nhit; hit++) {
+
+    pos[0]          = hits[hit*kNhit+0];
+    pos[1]          = hits[hit*kNhit+1];
+    pos[2]          = hits[hit*kNhit+2];
+    Float_t q       = hits[hit*kNhit+3];
+    Float_t hittime = hits[hit*kNhit+5];
+    Int_t   track   = ((Int_t) hits[hit*kNhit+4]);
+
+    Int_t   inDrift = 1;
+
+    // Find the current volume with the geo manager
+    gGeoManager->SetCurrentPoint(pos);
+    gGeoManager->FindNode();      
+    if (strstr(gGeoManager->GetPath(),"/UK")) {
+      inDrift = 0;
     }
 
-    // Add a container for the digits of this detector
-    digits = fDigitsManager->GetDigits(iDet);        
-    // Allocate memory space for the digits buffer
-    if (digits->GetNtime() == 0) {
-      digits->Allocate(nRowMax,nColMax,nTimeTotal);
+    // Get the calibration objects
+    calVdriftROC      = calibration->GetVdriftROC(det);
+    calVdriftDetValue = calVdriftDet->GetValue(det);
+    calT0ROC          = calibration->GetT0ROC(det);
+    calT0DetValue     = calT0Det->GetValue(det);
+
+    // Go to the local coordinate system:
+    // loc[0] - col  direction in amplification or driftvolume
+    // loc[1] - row  direction in amplification or driftvolume
+    // loc[2] - time direction in amplification or driftvolume
+    gGeoManager->MasterToLocal(pos,loc);
+    if (inDrift) {
+      // Relative to middle of amplification region
+      loc[2] = loc[2] - kDrWidth/2.0 - kAmWidth/2.0;
+    } 
+
+    // The driftlength [cm] (w/o diffusion yet !).
+    // It is negative if the hit is between pad plane and anode wires.
+    Double_t driftlength = -1.0 * loc[2];
+
+    // Stupid patch to take care of TR photons that are absorbed
+    // outside the chamber volume. A real fix would actually need
+    // a more clever implementation of the TR hit generation
+    if (q < 0.0) {
+      if ((loc[1] < padPlane->GetRowEndROC()) ||
+          (loc[1] > padPlane->GetRow0ROC())) {
+        continue;
+      }
+      if ((driftlength < kDrMin) ||
+          (driftlength > kDrMax)) {
+        continue;
+      }
     }
-    else if (fSimpleSim) {
-      digits->Clear();
+
+    // Get row and col of unsmeared electron to retrieve drift velocity
+    // The pad row (z-direction)
+    Int_t    rowE         = padPlane->GetPadRowNumberROC(loc[1]);
+    if (rowE < 0) {
+      continue;
     }
-    // Get the signal container
-    signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
-    if (signals->GetNtime() == 0) {
-      // Create missing containers
-      signals->Allocate(nRowMax,nColMax,nTimeTotal);      
+    Double_t rowOffset    = padPlane->GetPadRowOffsetROC(rowE,loc[1]);
+    // The pad column (rphi-direction)
+    Double_t offsetTilt   = padPlane->GetTiltOffset(rowOffset);
+    Int_t    colE         = padPlane->GetPadColNumber(loc[0]+offsetTilt);
+    if (colE < 0) {
+      continue;          
     }
-    else {
-      // Expand the container if neccessary
-      if (fCompress) signals->Expand();
+    Double_t colOffset    = 0.0;
+
+    // Normalized drift length
+    Float_t  driftvelocity  = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
+    Double_t absdriftlength = TMath::Abs(driftlength);
+    if (commonParam->ExBOn()) {
+      absdriftlength /= TMath::Sqrt(GetLorentzFactor(driftvelocity));
     }
-    // Create the missing dictionary containers
-    if (!fSimpleSim) {    
-      for (iDict = 0; iDict < kNDict; iDict++) {       
-        dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
-        if (dictionary[iDict]->GetNtime() == 0) {
-          dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
+
+    // Loop over all electrons of this hit
+    // TR photons produce hits with negative charge
+    Int_t nEl = ((Int_t) TMath::Abs(q));
+    for (Int_t iEl = 0; iEl < nEl; iEl++) {
+
+      // Now the real local coordinate system of the ROC
+      // column direction: locC
+      // row direction:    locR 
+      // time direction:   locT
+      // locR and locC are identical to the coordinates of the corresponding
+      // volumina of the drift or amplification region.
+      // locT is defined relative to the wire plane (i.e. middle of amplification
+      // region), meaning locT = 0, and is negative for hits coming from the
+      // drift region. 
+      Double_t locC = loc[0];
+      Double_t locR = loc[1];
+      Double_t locT = loc[2];
+
+      // Electron attachment
+      if (simParam->ElAttachOn()) {
+        if (gRandom->Rndm() < (absdriftlength * elAttachProp)) {
+          continue;
         }
+      }
+          
+      // Apply the diffusion smearing
+      if (simParam->DiffusionOn()) {
+        if (!(Diffusion(driftvelocity,absdriftlength,locR,locC,locT))) {
+          continue;
+       }
+      }
+
+      // Apply E x B effects (depends on drift direction)
+      if (commonParam->ExBOn()) { 
+        if (!(ExB(driftvelocity,driftlength,locC))) {
+          continue;
+       }
+      }
+
+      // The electron position after diffusion and ExB in pad coordinates.
+      // The pad row (z-direction)
+      rowE       = padPlane->GetPadRowNumberROC(locR);
+      if (rowE < 0) continue;
+      rowOffset  = padPlane->GetPadRowOffsetROC(rowE,locR);
+
+      // The pad column (rphi-direction)
+      offsetTilt = padPlane->GetTiltOffset(rowOffset);
+      colE       = padPlane->GetPadColNumber(locC+offsetTilt);
+      if (colE < 0) continue;         
+      colOffset  = padPlane->GetPadColOffset(colE,locC+offsetTilt);
+         
+      // Also re-retrieve drift velocity because col and row may have changed
+      driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
+      Float_t t0    = calT0DetValue     + calT0ROC->GetValue(colE,rowE);
+
+      // Convert the position to drift time [mus], using either constant drift velocity or
+      // time structure of drift cells (non-isochronity, GARFIELD calculation).
+      // Also add absolute time of hits to take pile-up events into account properly
+      Double_t drifttime;
+      if (simParam->TimeStructOn()) {
+       // Get z-position with respect to anode wire
+        Double_t zz  =  row0 - locR + padPlane->GetAnodeWireOffset();
+       zz -= ((Int_t)(2 * zz)) / 2.0;
+        if (zz > 0.25) {
+          zz  = 0.5 - zz;
+        }
+        // Use drift time map (GARFIELD)
+        drifttime = commonParam->TimeStruct(driftvelocity,0.5*kAmWidth-1.0*locT,zz)
+                  + hittime;
       } 
-    }
+      else {
+       // Use constant drift velocity
+        drifttime = TMath::Abs(locT) / driftvelocity
+                  + hittime;
+      }
 
-    Int_t nDigits = 0;
+      // Apply the gas gain including fluctuations
+      Double_t ggRndm = 0.0;
+      do {
+        ggRndm = gRandom->Rndm();
+      } while (ggRndm <= 0);
+      Double_t signal = -(simParam->GetGasGain()) * TMath::Log(ggRndm);
+
+      // Apply the pad response 
+      if (simParam->PRFOn()) {
+        // The distance of the electron to the center of the pad 
+        // in units of pad width
+        Double_t dist = (colOffset - 0.5*padPlane->GetColSize(colE))
+                      / padPlane->GetColSize(colE);
+        // This is a fixed parametrization, i.e. not dependent on
+        // calibration values !
+        if (!(calibration->PadResponse(signal,dist,layer,padSignal))) continue;
+      }
+      else {
+        padSignal[0] = 0.0;
+        padSignal[1] = signal;
+        padSignal[2] = 0.0;
+      }
 
-    // Don't create noise in detectors that are switched off
-    if (CheckDetector(plane,chamber,sector)) {
+      // The time bin (always positive), with t0 distortion
+      Double_t timeBinIdeal = drifttime * samplingRate + t0;
+      // Protection 
+      if (TMath::Abs(timeBinIdeal) > 2*nTimeTotal) {
+        timeBinIdeal = 2 * nTimeTotal;
+      }
+      Int_t    timeBinTruncated = ((Int_t) timeBinIdeal);
+      // The distance of the position to the middle of the timebin
+      Double_t timeOffset       = ((Float_t) timeBinTruncated 
+                                + 0.5 - timeBinIdeal) / samplingRate;
+          
+      // Sample the time response inside the drift region
+      // + additional time bins before and after.
+      // The sampling is done always in the middle of the time bin
+      for (Int_t iTimeBin = TMath::Max(timeBinTruncated,0)
+          ;iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal)
+         ;iTimeBin++) {
+
+        // Apply the time response
+        Double_t timeResponse = 1.0;
+        Double_t crossTalk    = 0.0;
+        Double_t time         = (iTimeBin - timeBinTruncated) / samplingRate + timeOffset;
+
+        if (simParam->TRFOn()) {
+          timeResponse = simParam->TimeResponse(time);
+        }
+        if (simParam->CTOn()) {
+          crossTalk    = simParam->CrossTalk(time);
+        }
 
-      // Create the digits for this chamber
-      for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
-        for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
+        signalOld[0] = 0.0;
+        signalOld[1] = 0.0;
+        signalOld[2] = 0.0;
 
-         // Create summable digits
-          if (fSDigits) {
+        for (iPad = 0; iPad < kNpad; iPad++) {
 
-            for (iTime = 0; iTime < nTimeTotal; iTime++) {         
-              Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
-              signalAmp *= fSDigitsScale;
-              signalAmp  = TMath::Min(signalAmp,(Float_t) 1.0e9);
-              Int_t adc  = (Int_t) signalAmp;
-              if (adc > 0) nDigits++;
-              digits->SetDataUnchecked(iRow,iCol,iTime,adc);
-           }
+          Int_t colPos = colE + iPad - 1;
+          if (colPos <        0) continue;
+          if (colPos >= nColMax) break;
 
-         }
-         // Create normal digits
+          // Add the signals
+          signalOld[iPad] = signals->GetData(rowE,colPos,iTimeBin);
+
+          if (colPos != colE) {
+           // Cross talk added to non-central pads
+            signalOld[iPad] += padSignal[iPad] 
+                            * (timeResponse + crossTalk);
+          } 
           else {
+           // W/o cross talk at central pad
+            signalOld[iPad] += padSignal[iPad] 
+                            * timeResponse;
+          }
 
-            for (iTime = 0; iTime < nTimeTotal; iTime++) {         
-              Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
-              // Pad and time coupling
-              signalAmp *= coupling;
-              // Add the noise, starting from minus ADC baseline in electrons
-              Double_t baselineEl = fPar->GetADCbaseline() * (fPar->GetADCinRange()
-                                                           / fPar->GetADCoutRange()) 
-                                                           / convert;
-              signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fPar->GetNoise())
-                                     ,-baselineEl);
-              // Convert to mV
-              signalAmp *= convert;
-              // Add ADC baseline in mV
-              signalAmp += fPar->GetADCbaseline() * (fPar->GetADCinRange()
-                                                   / fPar->GetADCoutRange());
-             // Convert to ADC counts. Set the overflow-bit fADCoutRange if the 
-             // signal is larger than fADCinRange
-              Int_t adc  = 0;
-              if (signalAmp >= fPar->GetADCinRange()) {
-                adc = ((Int_t) fPar->GetADCoutRange());
-             }
-              else {
-                adc = ((Int_t) (signalAmp * (fPar->GetADCoutRange() 
-                                           / fPar->GetADCinRange())));
+          signals->SetData(rowE,colPos,iTimeBin,signalOld[iPad]);
+
+          // Store the track index in the dictionary
+          // Note: We store index+1 in order to allow the array to be compressed
+         // Note2: Taking out the +1 in track
+          if (signalOld[iPad] > 0.0) { 
+            for (dict = 0; dict < kNdict; dict++) {
+             Int_t oldTrack = dictionary[dict]->GetData(rowE,colPos,iTimeBin); 
+             if (oldTrack == track) break;
+             if (oldTrack ==  -1 ) {
+               dictionary[dict]->SetData(rowE,colPos,iTimeBin,track);
+               break;
              }
-              inADC[iTime]  = adc;
-              outADC[iTime] = adc;
-           }
-
-           // Apply the tail cancelation via the digital filter
-            if (fPar->TCOn()) {
-              DeConvExp(inADC,outADC,nTimeTotal,fPar->GetTCnexp());
-           }
-
-            for (iTime = 0; iTime < nTimeTotal; iTime++) {   
-              // Store the amplitude of the digit if above threshold
-              if (outADC[iTime] > fPar->GetADCthreshold()) {
-                if (fDebug > 2) {
-                  printf("  iRow = %d, iCol = %d, iTime = %d, adc = %f\n"
-                        ,iRow,iCol,iTime,outADC[iTime]);
-               }
-                nDigits++;
-                digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
-             }
-           }
+            }
+          }
 
-         }
+       } // Loop: pads
 
-        }
-      }
+      } // Loop: time bins
 
-    }
+    } // Loop: electrons of a single hit
 
-    // Compress the arrays
-    if (!fSimpleSim) {  
-      digits->Compress(1,0);
-      for (iDict = 0; iDict < kNDict; iDict++) {
-        dictionary[iDict]->Compress(1,0);
-      }
+  } // Loop: hits
 
-      totalSizeDigits += digits->GetSize();
-      totalSizeDict0  += dictionary[0]->GetSize();
-      totalSizeDict1  += dictionary[1]->GetSize();
-      totalSizeDict2  += dictionary[2]->GetSize();
-
-      Float_t nPixel = nRowMax * nColMax * nTimeMax;
-      if (fDebug > 0) {
-        printf("<AliTRDdigitizer::MakeDigits> ");
-        printf("Found %d digits in detector %d (%3.0f).\n"
-              ,nDigits,iDet
-              ,100.0 * ((Float_t) nDigits) / nPixel);
-      } 
+  AliDebug(2,Form("Finished analyzing %d hits",nhit));
 
-      if (fCompress) signals->Compress(1,0);
+  return kTRUE;
 
-    }
+}
 
-    delete [] inADC;
-    delete [] outADC;
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::ConvertSignals(Int_t det, AliTRDarraySignal *signals)
+{
+  //
+  // Convert signals to digits
+  //
 
-  }
+  AliDebug(1,Form("Start converting the signals for detector %d",det));
 
-  if (signalsArray) {
-    delete signalsArray;
-    signalsArray = 0;
+  if (fSDigits) {
+    // Convert the signal array to s-digits
+    if (!Signal2SDigits(det,signals)) {
+      return kFALSE;
+    }
   }
-
-  if (fDebug > 0) {
-    printf("<AliTRDdigitizer::MakeDigits> ");
-    printf("Total number of analyzed hits = %d\n",countHits);
-    if (!fSimpleSim) {    
-      printf("<AliTRDdigitizer::MakeDigits> ");
-      printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
-                                                        ,totalSizeDict0
-                                                        ,totalSizeDict1
-                                                        ,totalSizeDict2);        
+  else {
+    // Convert the signal array to digits
+    if (!Signal2ADC(det,signals)) {
+      return kFALSE;
     }
+    // Run digital processing for digits
+    RunDigitalProcessing(det);
   }
 
+  // Compress the arrays
+  CompressOutputArrays(det);   
+
   return kTRUE;
 
 }
 
 //_____________________________________________________________________________
-void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
+Bool_t AliTRDdigitizer::Signal2ADC(Int_t det, AliTRDarraySignal *signals)
 {
   //
-  // Add a digits manager for s-digits to the input list.
+  // Converts the sampled electron signals to ADC values for a given chamber
   //
 
-  fSDigitsManagerList->Add(man);
+  AliDebug(1,Form("Start converting signals to ADC values for detector=%d",det));
+
+  AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliFatal("Could not get calibration object");
+    return kFALSE;
+  }
+
+  AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
+  if (!simParam) {
+    AliFatal("Could not get simulation parameters");
+    return kFALSE;
+  }  
+
+  // Converts number of electrons to fC
+  const Double_t kEl2fC = 1.602e-19 * 1.0e15; 
+
+  // Coupling factor
+  Double_t coupling     = simParam->GetPadCoupling() 
+                        * simParam->GetTimeCoupling();
+  // Electronics conversion factor
+  Double_t convert      = kEl2fC 
+                        * simParam->GetChipGain();
+  // ADC conversion factor
+  Double_t adcConvert   = simParam->GetADCoutRange()
+                        / simParam->GetADCinRange();
+  // The electronics baseline in mV
+  Double_t baseline     = simParam->GetADCbaseline() 
+                        / adcConvert;
+  // The electronics baseline in electrons
+  Double_t baselineEl   = baseline
+                        / convert;
+
+  Int_t row  = 0;
+  Int_t col  = 0;
+  Int_t time = 0;
+
+  Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
+  Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
+  Int_t nTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
+  if (fSDigitsManager->GetDigitsParam()->GetNTimeBins(det)) {
+    nTimeTotal = fSDigitsManager->GetDigitsParam()->GetNTimeBins(det);
+  }
+  else {
+    AliFatal("Could not get number of time bins");
+    return kFALSE;
+  }
+
+  // The gain factor calibration objects
+  const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();  
+  AliTRDCalROC       *calGainFactorROC      = 0x0;
+  Float_t             calGainFactorDetValue = 0.0;
+
+  AliTRDarrayADC     *digits = 0x0;
+
+  if (!signals) {
+    AliError(Form("Signals array for detector %d does not exist\n",det));
+    return kFALSE;
+  }
+  if (signals->HasData()) {
+    // Expand the container if neccessary
+    signals->Expand();   
+  }
+  else {
+    // Create missing containers
+    signals->Allocate(nRowMax,nColMax,nTimeTotal);      
+  }
+
+  // Get the container for the digits of this detector
+  if (fDigitsManager->HasSDigits()) {
+    AliError("Digits manager has s-digits");
+    return kFALSE;
+  }
+
+  digits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det);
+  // Allocate memory space for the digits buffer
+  if (!digits->HasData()) {
+    digits->Allocate(nRowMax,nColMax,nTimeTotal);
+  }
+
+  // Get the calibration objects
+  calGainFactorROC      = calibration->GetGainFactorROC(det);
+  calGainFactorDetValue = calGainFactorDet->GetValue(det);
+
+  // Get the online gain factors
+  AliTRDCalOnlineGainTableROC *onlineGainFactorROC 
+    = calibration->GetOnlineGainTableROC(det);
+
+  // Create the digits for this chamber
+  for (row  = 0; row  <  nRowMax; row++ ) {
+    for (col  = 0; col  <  nColMax; col++ ) {
+
+      // Check whether pad is masked
+      // Bridged pads are not considered yet!!!
+      if (calibration->IsPadMasked(det,col,row) || 
+          calibration->IsPadNotConnected(det,col,row)) {
+        continue;
+      }
+
+      // The gain factors
+      Float_t padgain = calGainFactorDetValue 
+                      * calGainFactorROC->GetValue(col,row);
+      if (padgain <= 0) {
+        AliError(Form("Not a valid gain %f, %d %d %d",padgain,det,col,row));
+      }
+
+      for (time = 0; time < nTimeTotal; time++) {
+
+       // Get the signal amplitude
+        Float_t signalAmp = signals->GetData(row,col,time);
+        // Pad and time coupling
+        signalAmp *= coupling;
+       // Gain factors
+       signalAmp *= padgain;
+
+        // Add the noise, starting from minus ADC baseline in electrons
+        signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,simParam->GetNoise())
+                               ,-baselineEl);
+
+        // Convert to mV
+        signalAmp *= convert;
+        // Add ADC baseline in mV
+        signalAmp += baseline;
+
+       // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
+       // signal is larger than fADCinRange
+        Short_t adc  = 0;
+        if (signalAmp >= simParam->GetADCinRange()) {
+          adc = ((Short_t) simParam->GetADCoutRange());
+       }
+        else {
+          adc = TMath::Nint(signalAmp * adcConvert);
+       }
+
+        // Saving all digits
+       digits->SetData(row,col,time,adc);
+
+      } // for: time
+
+    } // for: col
+  } // for: row
+
+  return kTRUE;
 
 }
 
 //_____________________________________________________________________________
-void AliTRDdigitizer::DeleteSDigitsManager()
+Bool_t AliTRDdigitizer::Signal2SDigits(Int_t det, AliTRDarraySignal *signals)
 {
   //
-  // Removes digits manager from the input list.
+  // Converts the sampled electron signals to s-digits
   //
 
-  fSDigitsManagerList->Delete();
+  AliDebug(1,Form("Start converting signals to s-digits for detector=%d",det));
+
+  AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliFatal("Could not get calibration object");
+    return kFALSE;
+  }
+
+  Int_t row  = 0;
+  Int_t col  = 0;
+  Int_t time = 0;
+
+  Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
+  Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
+  Int_t nTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
+
+  // Get the container for the digits of this detector
+  if (!fDigitsManager->HasSDigits()) {
+    AliError("Digits manager has no s-digits");
+    return kFALSE;
+  }
+
+  AliTRDarraySignal *digits = (AliTRDarraySignal *) fDigitsManager->GetSDigits(det);
+  // Allocate memory space for the digits buffer
+  if (!digits->HasData()) {
+    digits->Allocate(nRowMax,nColMax,nTimeTotal);
+  }
+
+  // Create the sdigits for this chamber
+  for (row  = 0; row  <  nRowMax; row++ ) {
+    for (col  = 0; col  <  nColMax; col++ ) {
+      for (time = 0; time < nTimeTotal; time++) {         
+        digits->SetData(row,col,time,signals->GetData(row,col,time));
+      } // for: time
+    } // for: col
+  } // for: row
+  
+  return kTRUE;
 
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDdigitizer::ConvertSDigits()
+Bool_t AliTRDdigitizer::Digits2SDigits(AliTRDdigitsManager * const manDig
+                                     , AliTRDdigitsManager * const manSDig)
 {
   //
-  // Converts s-digits to normal digits
+  // Converts digits into s-digits. Needed for embedding into real data.
   //
 
-  // Number of track dictionary arrays
-  const Int_t    kNDict = AliTRDdigitsManager::kNDict;
-
-  // Converts number of electrons to fC
-  const Double_t kEl2fC = 1.602E-19 * 1.0E15; 
+  AliDebug(1,"Start converting digits to s-digits");
 
-  Int_t iDict = 0;
-  Int_t iRow;
-  Int_t iCol;
-  Int_t iTime;
-
-  if (!fPar) {    
-    fPar = new AliTRDparameter("TRDparameter","Standard parameter");
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::ConvertSDigits> ");
-      printf("Create the default parameter object\n");
-    }
+  if (!fGeo) {
+    fGeo = new AliTRDgeometry();
   }
 
-  Double_t sDigitsScale = 1.0 / GetSDigitsScale();
-  Double_t noise        = fPar->GetNoise();
-  Double_t padCoupling  = fPar->GetPadCoupling();
-  Double_t timeCoupling = fPar->GetTimeCoupling();
-  Double_t chipGain     = fPar->GetChipGain();
-  Double_t coupling     = padCoupling * timeCoupling;
-  Double_t convert      = kEl2fC * chipGain;
-  Double_t adcInRange   = fPar->GetADCinRange();
-  Double_t adcOutRange  = fPar->GetADCoutRange();
-  Int_t    adcThreshold = fPar->GetADCthreshold();
-  Int_t    adcBaseline  = fPar->GetADCbaseline();   
-
-  AliTRDdataArrayI *digitsIn;
-  AliTRDdataArrayI *digitsOut;
-  AliTRDdataArrayI *dictionaryIn[kNDict];
-  AliTRDdataArrayI *dictionaryOut[kNDict];
-
-  // Loop through the detectors
-  for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
+  AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliFatal("Could not get calibration object");
+    return kFALSE;
+  }
 
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::ConvertSDigits> ");
-      printf("Convert detector %d to digits.\n",iDet);
-    }
+  AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
+  if (!simParam) {
+    AliFatal("Could not get simulation parameters");
+    return kFALSE;
+  }  
 
-    Int_t plane      = fGeo->GetPlane(iDet);
-    Int_t sector     = fGeo->GetSector(iDet);
-    Int_t chamber    = fGeo->GetChamber(iDet);
-    Int_t nRowMax    = fPar->GetRowMax(plane,chamber,sector);
-    Int_t nColMax    = fPar->GetColMax(plane);
-    Int_t nTimeTotal = fPar->GetTimeTotal();
-
-    Double_t *inADC  = new Double_t[nTimeTotal];
-    Double_t *outADC = new Double_t[nTimeTotal];
-
-    digitsIn  = fSDigitsManager->GetDigits(iDet);
-    digitsIn->Expand();
-    digitsOut = fDigitsManager->GetDigits(iDet);
-    digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
-    for (iDict = 0; iDict < kNDict; iDict++) {
-      dictionaryIn[iDict]  = fSDigitsManager->GetDictionary(iDet,iDict);
-      dictionaryIn[iDict]->Expand();
-      dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
-      dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
+  // Converts number of electrons to fC
+  const Double_t kEl2fC = 1.602e-19 * 1.0e15; 
+
+  // Coupling factor
+  Double_t coupling     = simParam->GetPadCoupling() 
+                        * simParam->GetTimeCoupling();
+  // Electronics conversion factor
+  Double_t convert      = kEl2fC 
+                        * simParam->GetChipGain();
+  // ADC conversion factor
+  Double_t adcConvert   = simParam->GetADCoutRange()
+                        / simParam->GetADCinRange();
+  // The electronics baseline in mV
+  Double_t baseline     = simParam->GetADCbaseline() 
+                        / adcConvert;
+  // The electronics baseline in electrons
+  //Double_t baselineEl   = baseline
+  //                      / convert;
+
+  // The gainfactor calibration objects
+  //const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();  
+  //AliTRDCalROC       *calGainFactorROC      = 0;
+  //Float_t             calGainFactorDetValue = 0.0;
+
+  Int_t row  = 0;
+  Int_t col  = 0;
+  Int_t time = 0;
+
+  for (Int_t det = 0; det < AliTRDgeometry::Ndet(); det++) {
+
+    Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
+    Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
+    Int_t nTimeTotal = manDig->GetDigitsParam()->GetNTimeBins(det);
+
+    // Get the calibration objects
+    //calGainFactorROC      = calibration->GetGainFactorROC(det);
+    //calGainFactorDetValue = calGainFactorDet->GetValue(det);
+
+    // Get the digits
+    AliTRDarrayADC *digits = (AliTRDarrayADC *) manDig->GetDigits(det);
+
+    if (!manSDig->HasSDigits()) {
+      AliError("SDigits manager has no s-digits");
+      return kFALSE;
     }
+    // Get the s-digits
+    AliTRDarraySignal     *sdigits = (AliTRDarraySignal *)     manSDig->GetSDigits(det);
+    AliTRDarrayDictionary *tracks0 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,0);
+    AliTRDarrayDictionary *tracks1 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,1);
+    AliTRDarrayDictionary *tracks2 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,2);
+    // Allocate memory space for the digits buffer
+    sdigits->Allocate(nRowMax,nColMax,nTimeTotal);
+    tracks0->Allocate(nRowMax,nColMax,nTimeTotal);
+    tracks1->Allocate(nRowMax,nColMax,nTimeTotal);
+    tracks2->Allocate(nRowMax,nColMax,nTimeTotal);
+
+    // Keep the digits param
+    manSDig->GetDigitsParam()->SetNTimeBinsAll(manDig->GetDigitsParam()->GetNTimeBins(0));
+    manSDig->GetDigitsParam()->SetADCbaselineAll(manDig->GetDigitsParam()->GetADCbaseline(0));
+
+    if (digits->HasData()) {
+
+      digits->Expand();
+
+      // Create the sdigits for this chamber
+      for (row  = 0; row  <  nRowMax; row++ ) {
+        for (col  = 0; col  <  nColMax; col++ ) {
+
+          // The gain factors
+          //Float_t padgain = calGainFactorDetValue 
+          //                * calGainFactorROC->GetValue(col,row);
+
+          for (time = 0; time < nTimeTotal; time++) {
+
+            Short_t  adcVal = digits->GetData(row,col,time);
+            Double_t signal = (Double_t) adcVal;
+            // ADC -> signal in mV
+            signal /= adcConvert;
+            // Subtract baseline in mV
+            signal -= baseline;
+            // Signal in mV -> signal in #electrons
+            signal /= convert;
+            // Gain factor
+            //signal /= padgain; // Not needed for real data
+            // Pad and time coupling
+            signal /= coupling;
+
+            sdigits->SetData(row,col,time,signal);
+            tracks0->SetData(row,col,time,0);
+            tracks1->SetData(row,col,time,0);
+            tracks2->SetData(row,col,time,0);
+
+          } // for: time
+
+        } // for: col
+      } // for: row
+  
+    } // if: has data
 
-    for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
-      for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
-
-        for (iTime = 0; iTime < nTimeTotal; iTime++) {         
-          Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
-          signal *= sDigitsScale;
-          // Pad and time coupling
-          signal *= coupling;
-          // Add the noise, starting from minus ADC baseline in electrons
-          Double_t baselineEl = adcBaseline * (adcInRange / adcOutRange) / convert;
-          signal  = TMath::Max((Double_t) gRandom->Gaus(signal,noise),-baselineEl);
-          // Convert to mV
-          signal *= convert;
-          // add ADC baseline in mV
-          signal += adcBaseline * (adcInRange / adcOutRange);
-         // Convert to ADC counts. Set the overflow-bit adcOutRange if the 
-         // signal is larger than adcInRange
-          Int_t adc  = 0;
-          if (signal >= adcInRange) {
-            adc = ((Int_t) adcOutRange);
-         }
-          else {
-            adc = ((Int_t) (signal * (adcOutRange / adcInRange)));
-         }
-          inADC[iTime]  = adc;
-          outADC[iTime] = adc;
-       }
+    sdigits->Compress(0);
+    tracks0->Compress();
+    tracks1->Compress();
+    tracks2->Compress();
 
-       // Apply the tail cancelation via the digital filter
-        if (fPar->TCOn()) {
-          DeConvExp(inADC,outADC,nTimeTotal,fPar->GetTCnexp());
-       }
+    // No compress just remove
+    manDig->RemoveDigits(det);
+    manDig->RemoveDictionaries(det);      
 
-        for (iTime = 0; iTime < nTimeTotal; iTime++) {   
-          // Store the amplitude of the digit if above threshold
-          if (outADC[iTime] > adcThreshold) {
-            digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
-           // Copy the dictionary
-            for (iDict = 0; iDict < kNDict; iDict++) { 
-              Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
-              dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
-           }
-         }
-       }
+  } // for: det
 
-      }
-    }
+  return kTRUE;
 
-    if (fCompress) {
-      digitsIn->Compress(1,0);
-      digitsOut->Compress(1,0);
-      for (iDict = 0; iDict < kNDict; iDict++) {
-        dictionaryIn[iDict]->Compress(1,0);
-        dictionaryOut[iDict]->Compress(1,0);
-      }
-    }
+}
 
-    delete [] inADC;
-    delete [] outADC;
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::SDigits2Digits()
+{
+  //
+  // Merges the input s-digits and converts them to normal digits
+  //
 
-  }    
+  if (!MergeSDigits()) {
+    return kFALSE;
+  }
 
-  return kTRUE;
+  return ConvertSDigits();
 
 }
 
@@ -1294,43 +1559,46 @@ Bool_t AliTRDdigitizer::MergeSDigits()
   // Number of track dictionary arrays
   const Int_t kNDict = AliTRDdigitsManager::kNDict;
 
-  if (!fPar) {
-    fPar = new AliTRDparameter("TRDparameter","Standard parameter");
-    if (fDebug > 0) {
-      printf("<AliTRDdigitizer::MergeSDigits> ");
-      printf("Create the default parameter object\n");
-    }
+  AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
+  if (!simParam) {
+    AliFatal("Could not get simulation parameters");
+    return kFALSE;
   }
-
+  
+  AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliFatal("Could not get calibration object");
+    return kFALSE;
+  }
+  
   Int_t iDict = 0;
   Int_t jDict = 0;
 
-  AliTRDdataArrayI *digitsA;
-  AliTRDdataArrayI *digitsB;
-  AliTRDdataArrayI *dictionaryA[kNDict];
-  AliTRDdataArrayI *dictionaryB[kNDict];
+  AliTRDarraySignal     *digitsA;
+  AliTRDarraySignal     *digitsB;
+  AliTRDarrayDictionary *dictionaryA[kNDict];
+  AliTRDarrayDictionary *dictionaryB[kNDict];
 
+  AliTRDdigitsManager   *mergeSDigitsManager = 0x0;
   // Get the first s-digits
   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
-  if (!fSDigitsManager) return kFALSE;
+  if (!fSDigitsManager) { 
+    AliError("No SDigits manager");
+    return kFALSE;
+  }
 
   // Loop through the other sets of s-digits
-  AliTRDdigitsManager *mergeSDigitsManager;
-  mergeSDigitsManager = (AliTRDdigitsManager *) 
-                        fSDigitsManagerList->After(fSDigitsManager);
-
-  if (fDebug > 0) {
-    if (mergeSDigitsManager) {
-      printf("<AliTRDdigitizer::MergeSDigits> ");
-      printf("Merge %d input files.\n",fSDigitsManagerList->GetSize());
-    }
-    else {
-      printf("<AliTRDdigitizer::MergeSDigits> ");
-      printf("Only one input file.\n");
-    }
-  }
+  mergeSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->After(fSDigitsManager);
 
+  if (mergeSDigitsManager) {
+    AliDebug(1,Form("Merge %d input files.",fSDigitsManagerList->GetSize()));
+  }
+  else {
+    AliDebug(1,"Only one input file.");
+  }
+  
   Int_t iMerge = 0;
+
   while (mergeSDigitsManager) {
 
     iMerge++;
@@ -1338,131 +1606,241 @@ Bool_t AliTRDdigitizer::MergeSDigits()
     // Loop through the detectors
     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
 
-      Int_t plane      = fGeo->GetPlane(iDet);
-      Int_t sector     = fGeo->GetSector(iDet);
-      Int_t chamber    = fGeo->GetChamber(iDet);
-      Int_t nRowMax    = fPar->GetRowMax(plane,chamber,sector);
-      Int_t nColMax    = fPar->GetColMax(plane);
-      Int_t nTimeTotal = fPar->GetTimeTotal();
+      Int_t nTimeTotal = fSDigitsManager->GetDigitsParam()->GetNTimeBins(iDet);
+      if (mergeSDigitsManager->GetDigitsParam()->GetNTimeBins(iDet) != nTimeTotal) {
+        AliError(Form("Mismatch in the number of time bins [%d,%d] in detector %d"
+                     ,nTimeTotal
+                    ,mergeSDigitsManager->GetDigitsParam()->GetNTimeBins(iDet)
+                    ,iDet));
+        return kFALSE;
+      }
 
+      Int_t nRowMax = fGeo->GetPadPlane(iDet)->GetNrows();
+      Int_t nColMax = fGeo->GetPadPlane(iDet)->GetNcols();
+         
       // Loop through the pixels of one detector and add the signals
-      digitsA = fSDigitsManager->GetDigits(iDet);
-      digitsB = mergeSDigitsManager->GetDigits(iDet);
-      digitsA->Expand();
-      digitsB->Expand();
+      digitsA = (AliTRDarraySignal *) fSDigitsManager->GetSDigits(iDet);    
+      digitsB = (AliTRDarraySignal *) mergeSDigitsManager->GetSDigits(iDet); 
+      digitsA->Expand();  
+      if (!digitsA->HasData()) continue;
+      digitsB->Expand();    
+      if (!digitsB->HasData()) continue;
+         
       for (iDict = 0; iDict < kNDict; iDict++) {
-        dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
-        dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
-        dictionaryA[iDict]->Expand();
+       dictionaryA[iDict] = (AliTRDarrayDictionary *) fSDigitsManager->GetDictionary(iDet,iDict);
+       dictionaryB[iDict] = (AliTRDarrayDictionary *) mergeSDigitsManager->GetDictionary(iDet,iDict);
+       dictionaryA[iDict]->Expand();  
         dictionaryB[iDict]->Expand();
       }
 
       // Merge only detectors that contain a signal
       Bool_t doMerge = kTRUE;
       if (fMergeSignalOnly) {
-        if (digitsA->GetOverThreshold(0) == 0) {
-          doMerge = kFALSE;
+        if (digitsA->GetOverThreshold(0) == 0) {                             
+         doMerge = kFALSE;
        }
       }
-
+         
       if (doMerge) {
-
-        if (fDebug > 0) {
-          printf("<AliTRDdigitizer::MergeSDigits> ");
-          printf("Merge detector %d of input no.%d\n",iDet,iMerge+1);
+             
+       AliDebug(1,Form("Merge detector %d of input no.%d",iDet,iMerge+1));
+             
+       for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
+         for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
+           for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
+               
+             // Add the amplitudes of the summable digits 
+             Float_t ampA = digitsA->GetData(iRow,iCol,iTime);
+             Float_t ampB = digitsB->GetData(iRow,iCol,iTime);
+             ampA += ampB;
+             digitsA->SetData(iRow,iCol,iTime,ampA);
+
+             // Add the mask to the track id if defined.
+             for (iDict = 0; iDict < kNDict; iDict++) {
+               Int_t trackB = dictionaryB[iDict]->GetData(iRow,iCol,iTime);
+               if ((fMasks) && (trackB > 0))  {
+                 for (jDict = 0; jDict < kNDict; jDict++) { 
+                   Int_t trackA = dictionaryA[iDict]->GetData(iRow,iCol,iTime); 
+                   if (trackA == 0) {
+                     trackA = trackB + fMasks[iMerge];
+                     dictionaryA[iDict]->SetData(iRow,iCol,iTime,trackA);  
+                   } // if:  track A == 0
+                 } // for: jDict
+               } // if:  fMasks and trackB > 0
+             } // for: iDict
+
+           } // for: iTime
+         } // for: iCol
+        } // for: iRow
+
+      } // if:  doMerge
+
+      mergeSDigitsManager->RemoveDigits(iDet);
+      mergeSDigitsManager->RemoveDictionaries(iDet);
+  
+      if (fCompress) {
+        digitsA->Compress(0); 
+        for (iDict = 0; iDict < kNDict; iDict++) {                                     
+         dictionaryA[iDict]->Compress();
         }
+      }
+      
+    } // for: detectors    
+      
+    // The next set of s-digits
+    mergeSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->After(mergeSDigitsManager);
+    
+  } // while: mergeDigitsManagers
+  
+  return kTRUE;
 
-        for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
-          for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
-            for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {         
+}
 
-             // Add the amplitudes of the summable digits 
-              Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
-              Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
-              ampA += ampB;
-              digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
-
-            // Add the mask to the track id if defined.
-              for (iDict = 0; iDict < kNDict; iDict++) {
-                Int_t trackB = dictionaryB[iDict]->GetDataUnchecked(iRow,iCol,iTime);
-                if ((fMasks) && (trackB > 0)) {
-                  for (jDict = 0; jDict < kNDict; jDict++) { 
-                    Int_t trackA = dictionaryA[iDict]->GetDataUnchecked(iRow,iCol,iTime);
-                    if (trackA == 0) {
-                      trackA = trackB + fMasks[iMerge];
-                      dictionaryA[iDict]->SetDataUnchecked(iRow,iCol,iTime,trackA);
-                   }
-                 }
-               }
-             }
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::ConvertSDigits()
+{
+  //
+  // Converts s-digits to normal digits
+  //
 
-           }
-         }
-        }
+  AliTRDarraySignal *digitsIn = 0x0;
 
-      }
+  if (!fSDigitsManager->HasSDigits()) {
+    AliError("No s-digits in digits manager");
+    return kFALSE;
+  }
 
-      if (fCompress) {
-        digitsA->Compress(1,0);
-        digitsB->Compress(1,0);
-        for (iDict = 0; iDict < kNDict; iDict++) {
-          dictionaryA[iDict]->Compress(1,0);
-          dictionaryB[iDict]->Compress(1,0);
-        }
-      }
+  // Loop through the detectors
+  for (Int_t det = 0; det < AliTRDgeometry::Ndet(); det++) {
 
-    }    
+    // Get the merged s-digits (signals)
+    digitsIn = (AliTRDarraySignal *) fSDigitsManager->GetSDigits(det);
+    if (!digitsIn->HasData()) {
+      AliDebug(2,Form("No digits for det=%d",det));
+      continue;
+    }
+    
+    // Convert the merged sdigits to digits
+    if (!Signal2ADC(det,digitsIn)) {
+      continue;
+    }
 
-    // The next set of s-digits
-    mergeSDigitsManager = (AliTRDdigitsManager *) 
-                          fSDigitsManagerList->After(mergeSDigitsManager);
+    // Copy the dictionary information to the output array
+    if (!CopyDictionary(det)) {
+      continue;
+    }
+
+    // Delete 
+    fSDigitsManager->RemoveDigits(det);
+    fSDigitsManager->RemoveDictionaries(det);
+
+    // Run digital processing
+    RunDigitalProcessing(det);
+
+    // Compress the arrays
+    CompressOutputArrays(det);
 
+  } // for: detector numbers
+
+  if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
+    if (trklLoader->Tree())
+      trklLoader->WriteData("OVERWRITE");
   }
 
+  // Save the values for the raw data headers
+  if (AliTRDSimParam::Instance()->GetNTBoverwriteOCDB()) {
+    fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(AliTRDSimParam::Instance()->GetNTimeBins());
+  }
+  else {
+    fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS());
+  }
+  fDigitsManager->GetDigitsParam()->SetADCbaselineAll(AliTRDSimParam::Instance()->GetADCbaseline());
+
   return kTRUE;
 
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDdigitizer::SDigits2Digits()
+Bool_t AliTRDdigitizer::CopyDictionary(Int_t det)
 {
   //
-  // Merges the input s-digits and converts them to normal digits
+  // Copies the dictionary information from the s-digits arrays
+  // to the output arrays
   //
 
-  if (!MergeSDigits()) return kFALSE;
+  AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliFatal("Could not get calibration object");
+    return kFALSE;
+  }
 
-  return ConvertSDigits();
+  AliDebug(1,Form("Start copying dictionaries for detector=%d",det));
+
+  const Int_t kNDict = AliTRDdigitsManager::kNDict;
+  AliTRDarrayDictionary *dictionaryIn[kNDict];
+  AliTRDarrayDictionary *dictionaryOut[kNDict];
+
+  Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
+  Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
+  Int_t nTimeTotal = fSDigitsManager->GetDigitsParam()->GetNTimeBins(det);
+
+  Int_t row  = 0;
+  Int_t col  = 0;
+  Int_t time = 0;
+  Int_t dict = 0;
+
+  for (dict = 0; dict < kNDict; dict++) {
+
+    dictionaryIn[dict]  = (AliTRDarrayDictionary *) fSDigitsManager->GetDictionary(det,dict);
+    dictionaryIn[dict]->Expand();
+    dictionaryOut[dict] = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
+    dictionaryOut[dict]->Allocate(nRowMax,nColMax,nTimeTotal);
+
+    for (row = 0; row < nRowMax; row++) {
+      for (col = 0; col < nColMax; col++) {
+        for (time = 0; time < nTimeTotal; time++) {
+         Int_t track = dictionaryIn[dict]->GetData(row,col,time);
+         dictionaryOut[dict]->SetData(row,col,time,track);
+       } // for: time
+      } // for: col
+    } // for: row
+    
+  } // for: dictionaries
+  
+  return kTRUE;
 
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
+void AliTRDdigitizer::CompressOutputArrays(Int_t det)
 {
   //
-  // Checks whether a detector is enabled
+  // Compress the output arrays
   //
 
-  if (fSimpleSim) return kTRUE; 
+  const Int_t kNDict = AliTRDdigitsManager::kNDict;
+  AliTRDarrayDictionary *dictionary = 0x0;
 
-  if ((fTRD->GetSensChamber() >=       0) &&
-      (fTRD->GetSensChamber() != chamber)) return kFALSE;
-  if ((fTRD->GetSensPlane()   >=       0) &&
-      (fTRD->GetSensPlane()   !=   plane)) return kFALSE;
-  if ( fTRD->GetSensSector()  >=       0) {
-    Int_t sens1 = fTRD->GetSensSector();
-    Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
-    sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
-           * AliTRDgeometry::Nsect();
-    if (sens1 < sens2) {
-      if ((sector < sens1) || (sector >= sens2)) return kFALSE;
+  if (fCompress) {
+
+    if (!fSDigits) {
+      AliTRDarrayADC *digits = 0x0;  
+      digits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det);
+      digits->Compress();
     }
-    else {
-      if ((sector < sens1) && (sector >= sens2)) return kFALSE;
+
+    if (fSDigits) {
+      AliTRDarraySignal *digits = 0x0; 
+      digits = (AliTRDarraySignal *) fDigitsManager->GetSDigits(det);
+      digits->Compress(0);
     }
-  }
 
-  return kTRUE;
+    for (Int_t dict = 0; dict < kNDict; dict++) {
+      dictionary = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
+      dictionary->Compress();
+    }
+
+  }
 
 }
 
@@ -1473,111 +1851,151 @@ Bool_t AliTRDdigitizer::WriteDigits() const
   // Writes out the TRD-digits and the dictionaries
   //
 
+  // Write parameters
+  fRunLoader->CdGAFile();
+
   // Store the digits and the dictionary in the tree
   return fDigitsManager->WriteDigits();
 
 }
 
 //_____________________________________________________________________________
-void AliTRDdigitizer::DeConvExp(Double_t *source, Double_t *target
-                              , Int_t n, Int_t nexp) 
+void AliTRDdigitizer::InitOutput(Int_t iEvent)
 {
   //
-  // Does the deconvolution by the digital filter.
-  //
-  // Author:        Marcus Gutfleisch, KIP Heidelberg
-  // Optimized for: New TRF from Venelin Angelov, simulated with CADENCE
-  //                Pad-ground capacitance = 25 pF
-  //                Pad-pad cross talk capacitance = 6 pF
-  //                For 10 MHz digitization, corresponding to 20 time bins
-  //                in the drift region
+  // Initializes the output branches
   //
 
-  Double_t rates[2];
-  Double_t coefficients[2];
+  fEvent = iEvent;
+   
+  if (!fRunLoader) {
+    AliError("Run Loader is NULL");
+    return;  
+  }
 
-  /* initialize (coefficient = alpha, rates = lambda) */
-  
-  if( nexp == 1 ) {
-    rates[0] = 0.466998;
-    /* no rescaling */
-    coefficients[0] = 1.0;
-  }
-  if( nexp == 2 ) {
-    rates[0] = 0.8988162;
-    coefficients[0] = 0.11392069;
-    rates[1] = 0.3745688;
-    coefficients[1] = 0.8860793;
-    /* no rescaling */
-    Float_t sumc = coefficients[0]+coefficients[1];
-    coefficients[0] /= sumc;
-    coefficients[1] /= sumc;
+  AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
+  if (!loader) {
+    AliError("Can not get TRD loader from Run Loader");
+    return;
   }
-      
-  Int_t i, k;
-  Double_t reminder[2];
-  Double_t correction, result;
 
-  /* attention: computation order is important */
-  correction=0.0;
-  for ( k = 0; k < nexp; k++ ) reminder[k]=0.0;
-    
-  for ( i = 0; i < n; i++ ) {
-    result = ( source[i] - correction );    /* no rescaling */
-    target[i] = result;
-    
-    for ( k = 0; k < nexp; k++ ) reminder[k] = rates[k] 
-                             * ( reminder[k] + coefficients[k] * result);
-      
-    correction=0.0;
-    for ( k = 0; k < nexp; k++ ) correction += reminder[k];
+  TTree *tree = 0;
+  
+  if (fSDigits) { 
+    // If we produce SDigits
+    tree = loader->TreeS();
+    if (!tree) {
+      loader->MakeTree("S");
+      tree = loader->TreeS();
+    }
   }
+  else {
+    // If we produce Digits
+    tree = loader->TreeD();
+    if (!tree) {
+      loader->MakeTree("D");
+      tree = loader->TreeD();
+    }
+  }
+  fDigitsManager->SetEvent(iEvent);
+  fDigitsManager->MakeBranch(tree);
+
+}
   
+//_____________________________________________________________________________
+Int_t AliTRDdigitizer::Diffusion(Float_t vdrift, Double_t absdriftlength
+                               , Double_t &lRow, Double_t &lCol, Double_t &lTime)
+{
+  //
+  // Applies the diffusion smearing to the position of a single electron.
+  // Depends on absolute drift length.
+  //
+  
+  Float_t diffL = 0.0;
+  Float_t diffT = 0.0;
+
+  if (AliTRDCommonParam::Instance()->GetDiffCoeff(diffL,diffT,vdrift)) {
+
+    Float_t driftSqrt = TMath::Sqrt(absdriftlength);
+    Float_t sigmaT    = driftSqrt * diffT;
+    Float_t sigmaL    = driftSqrt * diffL;
+    lRow  = gRandom->Gaus(lRow ,sigmaT);
+    lCol  = gRandom->Gaus(lCol ,sigmaT * GetLorentzFactor(vdrift));
+    lTime = gRandom->Gaus(lTime,sigmaL * GetLorentzFactor(vdrift));
+
+    return 1;
+
+  }
+  else {
+
+    return 0;
+
+  }
+
 }
 
 //_____________________________________________________________________________
-void AliTRDdigitizer::InitOutput(Int_t iEvent)
+Float_t AliTRDdigitizer::GetLorentzFactor(Float_t vd)
 {
   //
-  // Initializes the output branches
+  // Returns the Lorentz factor
   //
 
-  fEvent = iEvent;
-   
-  if (!fRunLoader)
-   {
-     Error("InitOutput","Run Loader is NULL");
-     return;  
-   }
-  AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
-  if (!loader)
-   {
-     Error("Open","Can not get TRD loader from Run Loader");
-     return;
-   }
-
-  TTree* tree = 0;
+  Double_t omegaTau      = AliTRDCommonParam::Instance()->GetOmegaTau(vd);
+  Double_t lorentzFactor = 1.0;
+  if (AliTRDCommonParam::Instance()->ExBOn()) {
+    lorentzFactor = 1.0 / (1.0 + omegaTau*omegaTau);
+  }
+
+  return lorentzFactor;
+
+}
   
-  if (fSDigits)
-   { 
-   //if we produce SDigits
-    tree = loader->TreeS();
-    if (!tree)
-     {
-      loader->MakeTree("S");
-      tree = loader->TreeS();
-     }
-   }
-  else
-   {//if we produce Digits
-     tree = loader->TreeD();
-     if (!tree)
-      {
-       loader->MakeTree("D");
-       tree = loader->TreeD();
+//_____________________________________________________________________________
+Int_t AliTRDdigitizer::ExB(Float_t vdrift, Double_t driftlength, Double_t &lCol)
+{
+  //
+  // Applies E x B effects to the position of a single electron.
+  // Depends on signed drift length.
+  //
+
+  lCol = lCol 
+       + AliTRDCommonParam::Instance()->GetOmegaTau(vdrift) 
+       * driftlength;
+
+  return 1;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDdigitizer::RunDigitalProcessing(Int_t det)
+{
+  //
+  // Run the digital processing in the TRAP
+  //
+
+  AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
+
+  AliTRDarrayADC *digits = fDigitsManager->GetDigits(det);
+  if (!digits)
+    return;
+
+  //Call the methods in the mcm class using the temporary array as input  
+  for(Int_t rob = 0; rob < digits->GetNrow() / 2; rob++)
+  {
+    for(Int_t mcm = 0; mcm < 16; mcm++)
+    {
+      fMcmSim->Init(det, rob, mcm); 
+      fMcmSim->SetDataByPad(digits, fDigitsManager);
+      fMcmSim->Filter();
+      if (feeParam->GetTracklet()) {
+        fMcmSim->Tracklet();
+        fMcmSim->StoreTracklets();
       }
-   }
-  fDigitsManager->SetEvent(iEvent);
-  fDigitsManager->MakeBranch(tree);
+      fMcmSim->ZSMapping();
+      fMcmSim->WriteData(digits);
+    }
+  }
 
 }
+