]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALDigitizer.cxx
follow modification in AliEMCALCalibData regarding the time calibration and bunch...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
index 5f1c50da82232606113106a396cfcac8a9b6027d..9e61fcf170192e80928c40456bea5c87d1560adf 100644 (file)
 //_________________________________________________________________________
 // 
 //////////////////////////////////////////////////////////////////////////////
-// Class performs digitization of Summable digits 
+// Class performs digitization of Summable digits from simulated data
 //  
-// In addition it performs mixing of summable digits from different events.
+// In addition it performs mixing/embedding of summable digits from different events.
 //
-// For each event two branches are created in TreeD:
+// For each event 3 branches are created in TreeD:
 //   "EMCAL" - list of digits
+//   "EMCALTRG" - list of trigger digits
 //   "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
 //
-// Note, that one cset title for new digits branch, and repeat digitization with
-// another set of parameters.
 //
-// Examples of use:
-// root[0] AliEMCALDigitizer * d = new AliEMCALDigitizer() ;
-// root[1] d->ExecuteTask()             
-// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-//                       //Digitizes SDigitis in all events found in file galice.root 
-//
-// root[2] AliEMCALDigitizer * d1 = new AliEMCALDigitizer("galice1.root") ;  
-//                       // Will read sdigits from galice1.root
-// root[3] d1->MixWith("galice2.root")       
-// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-//                       // Reads another portion of sdigits from galice2.root
-// root[3] d1->MixWith("galice3.root")       
-//                       // Reads another portion of sdigits from galice3.root
-// root[4] d->ExecuteTask("deb timing")    
-//                       // Reads SDigits from files galice1.root, galice2.root ....
-//                       // mixes them and stores produced Digits in file galice1.root          
-//                       // deb - prints number of produced digits
-//                       // deb all - prints list of produced digits
-//                       // timing  - prints time used for digitization
 ////////////////////////////////////////////////////////////////////////////////////
 //
 //*-- Author: Sahal Yacoob (LBL)
 // based on : AliEMCALDigitizer
 // Modif: 
 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
-//                           of new  IO (à la PHOS)
-///_________________________________________________________________________________
+//                           of new IO (a la PHOS)
+//  November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry 
+//  July 2011 GCB: Digitizer modified to accomodate embedding. 
+//                 Time calibration added. Decalibration possibility of energy and time added
+//_________________________________________________________________________________
 
 // --- ROOT system ---
-#include "TFile.h"
-#include "TTree.h"
-#include "TSystem.h"
-#include "TROOT.h"
-#include "TFolder.h"
-#include "TObjString.h"
-#include "TGeometry.h"
-#include "TBenchmark.h"
-
-// --- Standard library ---
+#include <TROOT.h>
+#include <TTree.h>
+#include <TSystem.h>
+#include <TBenchmark.h>
+#include <TBrowser.h>
+#include <TObjectTable.h>
+#include <TRandom.h>
+#include <TF1.h>
+#include <cassert>
 
 // --- AliRoot header files ---
+#include "AliLog.h"
 #include "AliRun.h"
-#include "AliHeader.h"
-#include "AliStream.h"
-#include "AliRunDigitizer.h"
+#include "AliDigitizationInput.h"
+#include "AliRunLoader.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCAL.h"
-#include "AliEMCALGetter.h"
+#include "AliEMCALLoader.h"
 #include "AliEMCALDigitizer.h"
 #include "AliEMCALSDigitizer.h"
 #include "AliEMCALGeometry.h"
-#include "AliEMCALTick.h"
+#include "AliEMCALCalibData.h"
+#include "AliEMCALSimParam.h"
+#include "AliEMCALRawDigit.h"
+
+  namespace
+  {
+    Double_t HeavisideTheta(Double_t x)
+    {
+      Double_t signal = 0.;
+      
+      if (x > 0.) signal = 1.;  
+      
+      return signal;  
+    }
+    
+    Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
+    {
+      Double_t v0 = par[0];
+      Double_t t0 = par[1];
+      Double_t tr = par[2];
+      
+      Double_t R1 = 1000.;
+      Double_t C1 = 33e-12;
+      Double_t R2 = 1800;
+      Double_t C2 = 22e-12;
+      
+      Double_t t  =   x[0];
+      
+      return (((0.8*(-((TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-t + t0)/(C1*R1))*
+                       TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2)) + 
+                    C1*C2*R1*R2*(1 - (C2*TMath::Power(TMath::E(),(-t + t0)/(C2*R2))*R2)/(-(C1*R1) + C2*R2)))*v0*
+               HeavisideTheta(t - t0))/tr 
+              - (0.8*(C1*C2*R1*R2 - 
+                      (TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C1*R1))*
+                       TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2) + 
+                      (C1*TMath::Power(C2,2)*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C2*R2))*
+                       R1*TMath::Power(R2,2))/(C1*R1 - C2*R2))*v0*
+                 HeavisideTheta(t - t0 - 1.25*tr))/tr)/(C2*R1));
+    }
+  }
 
 ClassImp(AliEMCALDigitizer)
-
-
+  
+  
 //____________________________________________________________________________ 
-  AliEMCALDigitizer::AliEMCALDigitizer()
+AliEMCALDigitizer::AliEMCALDigitizer()
+  : AliDigitizer("",""),
+    fDefaultInit(kTRUE),
+    fDigitsInRun(0),
+    fInit(0),
+    fInput(0),
+    fInputFileNames(0x0),
+    fEventNames(0x0),
+    fDigitThreshold(0),
+    fMeanPhotonElectron(0),
+    fGainFluctuations(0),
+    fPinNoise(0),
+    fTimeNoise(0),
+    fTimeDelay(0),
+    fTimeResolutionPar0(0),
+    fTimeResolutionPar1(0),
+    fADCchannelEC(0),
+    fADCpedestalEC(0),
+    fADCchannelECDecal(0),
+    fTimeChannel(0),
+    fTimeChannelDecal(0),
+    fNADCEC(0),
+    fEventFolderName(""),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fCalibData(0x0),
+    fSDigitizer(0x0)
 {
   // ctor
-  //InitParameters() ; 
-  fDefaultInit = kTRUE ; 
-  fManager = 0 ;                     // We work in the standalong mode
+  InitParameters() ; 
+  fDigInput = 0 ;                     // We work in the standalone mode
 }
 
 //____________________________________________________________________________ 
-AliEMCALDigitizer::AliEMCALDigitizer(const char *headerFile, const char *name, const Bool_t toSplit)
+AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
+  : AliDigitizer("EMCALDigitizer", alirunFileName),
+    fDefaultInit(kFALSE),
+    fDigitsInRun(0),
+    fInit(0),
+    fInput(0),
+    fInputFileNames(0), 
+    fEventNames(0), 
+    fDigitThreshold(0),
+    fMeanPhotonElectron(0),
+    fGainFluctuations(0),
+    fPinNoise(0),
+    fTimeNoise(0),
+    fTimeDelay(0),
+    fTimeResolutionPar0(0),
+    fTimeResolutionPar1(0),
+    fADCchannelEC(0),
+    fADCpedestalEC(0),
+    fADCchannelECDecal(0),
+    fTimeChannel(0),
+    fTimeChannelDecal(0),
+    fNADCEC(0),
+    fEventFolderName(eventFolderName),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fCalibData(0x0),
+    fSDigitizer(0x0)
 {
   // ctor
-
-  SetTitle(headerFile) ;
-  SetName(name) ;
-  fManager = 0 ;                     // We work in the standalong mode
-  fSplitFile= 0 ; 
-  fToSplit = toSplit ;
-  Init() ;
   InitParameters() ; 
-  fDefaultInit = kFALSE ; 
+  Init() ;
+  fDigInput = 0 ;                     // We work in the standalone mode
 }
 
 //____________________________________________________________________________ 
-AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * ard):AliDigitizer(ard)
+AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d) 
+  : AliDigitizer(d.GetName(),d.GetTitle()),
+    fDefaultInit(d.fDefaultInit),
+    fDigitsInRun(d.fDigitsInRun),
+    fInit(d.fInit),
+    fInput(d.fInput),
+    fInputFileNames(d.fInputFileNames),
+    fEventNames(d.fEventNames),
+    fDigitThreshold(d.fDigitThreshold),
+    fMeanPhotonElectron(d.fMeanPhotonElectron),
+    fGainFluctuations(d.fGainFluctuations),
+    fPinNoise(d.fPinNoise),
+    fTimeNoise(d.fTimeNoise),
+    fTimeDelay(d.fTimeDelay),
+    fTimeResolutionPar0(d.fTimeResolutionPar0),
+    fTimeResolutionPar1(d.fTimeResolutionPar1),
+    fADCchannelEC(d.fADCchannelEC),
+    fADCpedestalEC(d.fADCpedestalEC),
+    fADCchannelECDecal(d.fADCchannelECDecal),
+    fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
+    fNADCEC(d.fNADCEC),
+    fEventFolderName(d.fEventFolderName),
+    fFirstEvent(d.fFirstEvent),
+    fLastEvent(d.fLastEvent),
+    fCalibData(d.fCalibData),
+    fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
 {
-  // ctor
-  SetTitle(ard->GetInputFileName(0,0)) ;
+  // copyy ctor 
+}
+
+//____________________________________________________________________________ 
+AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
+  : AliDigitizer(rd,"EMCALDigitizer"),
+    fDefaultInit(kFALSE),
+    fDigitsInRun(0),
+    fInit(0),
+    fInput(0),
+    fInputFileNames(0),
+    fEventNames(0),
+    fDigitThreshold(0),
+    fMeanPhotonElectron(0),
+    fGainFluctuations(0),
+    fPinNoise(0.),
+    fTimeNoise(0.),
+    fTimeDelay(0.),
+    fTimeResolutionPar0(0.),
+    fTimeResolutionPar1(0.),
+    fADCchannelEC(0),
+    fADCpedestalEC(0),
+    fADCchannelECDecal(0),
+    fTimeChannel(0),
+    fTimeChannelDecal(0),    
+    fNADCEC(0),
+    fEventFolderName(0),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fCalibData(0x0),
+    fSDigitizer(0x0)
+{
+  // ctor Init() is called by RunDigitizer
+  fDigInput = rd ; 
+  fEventFolderName = fDigInput->GetInputFolderName(0) ;
+  SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
   InitParameters() ; 
-  fDefaultInit = kFALSE ; 
-  fSplitFile   = 0 ; 
-  if (ard->GetOutputFile()) {
-    SetName(ard->GetOutputFile().Data());
-    fToSplit = kTRUE ;
-  } else {
-    SetName("Default") ;
-    fToSplit = kFALSE ;
-  }
 }
 
 //____________________________________________________________________________ 
   AliEMCALDigitizer::~AliEMCALDigitizer()
 {
-  // dtor
-
-    fSplitFile = 0 ; 
+  //dtor
+  delete [] fInputFileNames ; 
+  delete [] fEventNames ; 
+  if (fSDigitizer) delete fSDigitizer;
 }
 
 //____________________________________________________________________________
-void AliEMCALDigitizer::Digitize(const Int_t event) 
+void AliEMCALDigitizer::Digitize(Int_t event) 
 { 
-
+  
   // Makes the digitization of the collected summable digits
   // for this it first creates the array of all EMCAL modules
-  // filled with noise (different for EMC, CPV and PPSD) and
-  // after that adds contributions from SDigits. This design 
-  // helps to avoid scanning over the list of digits to add 
-  // contribution of any new SDigit.
-
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  TClonesArray * digits = gime->Digits(GetName()) ; 
+  // filled with noise and after that adds contributions from 
+  // SDigits. This design helps to avoid scanning over the 
+  // list of digits to add  contribution of any new SDigit.
+  //
+  // JLK 26-Jun-2008
+  // Note that SDigit energy info is stored as an amplitude, so we
+  // must call the Calibrate() method of the SDigitizer to convert it
+  // back to an energy in GeV before adding it to the Digit
+  //
+  static int nEMC=0; //max number of digits possible
+  AliRunLoader *rl = AliRunLoader::Instance();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
   
-  digits->Clear() ;
-
-  const AliEMCALGeometry *geom = gime->EMCALGeometry() ; 
-  //Making digits with noise, first EMC
-  Int_t nEMC = 0 ; 
-  if (geom->GetNHCLayers() > 0 )
-    nEMC = 3*geom->GetNPhi()*geom->GetNZ(); //max number of digits possible (Preshower, ECAL, HCAL)
-  else 
-    nEMC = 2*geom->GetNPhi()*geom->GetNZ(); //max number of digits possible (Preshower, ECAL)
- Int_t absID ;
-  TString name      =  geom->GetName() ;
-
- // get first the sdigitizer from the tasks list (must have same name as the digitizer)
-  const AliEMCALSDigitizer * sDigitizer = gime->SDigitizer(GetName()); 
-  if ( !sDigitizer) {
-    Fatal("Digitize", "SDigitizer with name %s not found", GetName() ); 
-  }
-
-// loop through the sdigits posted to the White Board and add them to the noise
-  TCollection * folderslist = gime->SDigitsFolder()->GetListOfFolders() ; 
-  TIter next(folderslist) ; 
-  TFolder * folder = 0 ; 
-  TClonesArray * sdigits = 0 ;
-  Int_t input = 0 ;
-  TObjArray * sdigArray = new TObjArray(2) ;
-  while ( (folder = (TFolder*)next()) ) {
-    if ( (sdigits = (TClonesArray*)folder->FindObject(GetName()) ) ) {
-      TString fileName(folder->GetName()) ;
-      fileName.ReplaceAll("_","/") ;
-      if (gDebug)
-       Info("Digitize", "Adding SDigits %s from %s", GetName(), fileName.Data()) ; 
-      sdigArray->AddAt(sdigits, input) ;
-      input++ ;
+  if(emcalLoader){
+    Int_t readEvent = event ; 
+    // fDigInput is data member from AliDigitizer
+    if (fDigInput) 
+      readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ; 
+    AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
+                    readEvent, GetTitle(), fEventFolderName.Data())) ; 
+    rl->GetEvent(readEvent);
+    
+    TClonesArray * digits = emcalLoader->Digits() ; 
+    digits->Delete() ;  //correct way to clear array when memory is
+    //allocated by objects stored in array
+    
+    // Load Geometry
+    AliEMCALGeometry *geom = 0;
+    if (!rl->GetAliRun()) {
+      AliFatal("Could not get AliRun from runLoader");
     }
-  }
-
-  //Find the first tower with signal
-  Int_t nextSig = nEMC + 1 ; 
-  Int_t i;
-  for(i=0; i<input; i++){
-    sdigits = (TClonesArray *)sdigArray->At(i) ;
-    if ( !sdigits->GetEntriesFast() )
-      continue ; 
-    Int_t curNext = ((AliEMCALDigit *)sdigits->At(0))->GetId() ;
-     if(curNext < nextSig) 
-       nextSig = curNext ;
-  }
-
-  TArrayI index(input) ;
-  index.Reset() ;  //Set all indexes to zero
-
-  AliEMCALDigit * digit ;
-  AliEMCALDigit * curSDigit ;
-
-  TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
-
-  //Put Noise contribution
-  for(absID = 1; absID <= nEMC; absID++){
-    Float_t amp = 0 ; 
-    //    Float_t noise = TMath::Abs(gRandom->Gaus(0., fPinNoise)); 
-    //new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID, sDigitizer->Digitize(noise), TimeOfNoise() ) ;
-    // amplitude set to zero, noise will be added later
-    new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ) ;
-    //look if we have to add signal?
-    digit = (AliEMCALDigit *) digits->At(absID-1) ;
-    if(absID==nextSig){
-      //Add SDigits from all inputs    
-      ticks->Clear() ;
-      Int_t contrib = 0 ;
-      Float_t a = digit->GetAmp() ;
-      Float_t b = TMath::Abs( a /fTimeSignalLength) ;
-      //Mark the beginnign of the signal
-      new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);  
-      //Mark the end of the ignal     
-      new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
-      
- // loop over input
-      for(i = 0; i< input ; i++){  //loop over (possible) merge sources
-       if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
-         curSDigit = (AliEMCALDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;        
-       else
-         curSDigit = 0 ;
-       //May be several digits will contribute from the same input
-       while(curSDigit && curSDigit->GetId() == absID){           
-         //Shift primary to separate primaries belonging different inputs
-         Int_t primaryoffset ;
-         if(fManager)
-           primaryoffset = fManager->GetMask(i) ; 
-         else
-           primaryoffset = i ;
-         curSDigit->ShiftPrimary(primaryoffset) ;
-         
-         a = curSDigit->GetAmp() ;
-         b = a /fTimeSignalLength ;
-         new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);  
-         new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
-
-         *digit = *digit + *curSDigit ;  //add energies
-
-         index[i]++ ;
-         if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
-           curSDigit = (AliEMCALDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;      
-         else
-           curSDigit = 0 ;
-       }
-      }
-      // add fluctuations for photo-electron creation
-      amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
-      if (gDebug==1) 
-       Info("Digitize", "id = %d BEFORE amp = %f\n", digit->GetId(),amp) ; 
-      amp *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
-      if (gDebug==1) 
-       Info("Digitize", "After fluctuation amp = %f\n", amp) ; 
-      //calculate and set time
-      Float_t time = FrontEdgeTime(ticks) ;
-      digit->SetTime(time) ;
-
-      //Find next signal module
-      nextSig = nEMC + 1 ;
-      for(i=0; i<input; i++){
-       sdigits = ((TClonesArray *)sdigArray->At(i)) ;
-       Int_t curNext = nextSig ;
-       if(sdigits->GetEntriesFast() > index[i] ){
-         curNext = ((AliEMCALDigit *) sdigits->At(index[i]))->GetId() ;
-         
+    else{
+      AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
+      geom = emcal->GetGeometry();
+      nEMC = geom->GetNCells();
+      AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
+    }
+    
+    Int_t absID = -1 ;
+    
+    digits->Expand(nEMC) ;
+    
+    // RS create a digitizer on the fly
+    if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
+    fSDigitizer->SetEventRange(0, -1) ;
+    //
+    //take all the inputs to add together and load the SDigits
+    TObjArray * sdigArray = new TObjArray(fInput) ;
+    sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
+    Int_t i ;
+    Int_t embed = kFALSE;
+    for(i = 1 ; i < fInput ; i++){
+      TString tempo(fEventNames[i]) ; 
+      tempo += i ;
+        
+      AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ; 
+        
+      if (rl2==0) 
+       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ; 
+        
+      if (fDigInput) 
+       readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ; 
+      Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ; 
+      rl2->LoadSDigits();
+      //     rl2->LoadDigits();
+      rl2->GetEvent(readEvent);
+      if(rl2->GetDetectorLoader("EMCAL"))
+        {
+          AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
+          if(emcalLoader2) {
+            //Merge 2 simulated sdigits
+            if(emcalLoader2->SDigits()){
+              TClonesArray* sdigits2 = emcalLoader2->SDigits();
+              sdigArray->AddAt(sdigits2, i) ;
+              // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
+              // do not smear energy of embedded, do not add noise to any sdigits
+              if(sdigits2->GetEntriesFast()>0){
+               //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
+               AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
+               if(digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded){
+                 embed = kTRUE;
+               }
+              }
+            }
+            
+          }//loader2
+          else  AliFatal("EMCAL Loader of second event not available!");
+          
+        }
+      else Fatal("Digitize", "Loader of second input not found");
+    }// input loop
+      
+    //Find the first tower with signal
+    Int_t nextSig = nEMC + 1 ; 
+    TClonesArray * sdigits ;  
+    for(i = 0 ; i < fInput ; i++){
+      if(i > 0 && embed) continue;
+      sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
+      if(sdigits){
+       if (sdigits->GetEntriesFast() ){
+         AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
+         if(sd){
+           Int_t curNext = sd->GetId() ;
+           if(curNext < nextSig) 
+             nextSig = curNext ;
+           AliDebug(1,Form("input %i : #sdigits %i \n",
+                           i, sdigits->GetEntriesFast()));
+              
+         }//First entry is not NULL
+         else {
+           Info("Digitize", "NULL sdigit pointer");
+           continue;
+         }
+       }//There is at least one entry
+       else {
+         Info("Digitize", "NULL sdigits array");
+         continue;
        }
-       if(curNext < nextSig) nextSig = curNext ;
+      }// SDigits array exists
+      else {
+       Info("Digitizer","Null sdigit pointer");
+       continue;
       }
-    }
-    // add the noise now
+    }// input loop
+    AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
+      
+    TArrayI index(fInput) ;
+    index.Reset() ;  //Set all indexes to zero
+      
+    AliEMCALDigit * digit ;
+    AliEMCALDigit * curSDigit ;
+      
+    //Put Noise contribution, smear time and energy     
+    Float_t timeResolution = 0;
+    for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
+      Float_t energy = 0 ;    
+        
+      // amplitude set to zero, noise will be added later
+      Float_t noiseTime = 0.;
+      if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
+      new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
+      //look if we have to add signal?
+      digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
+        
+      if (digit) {
+          
+       if(absID==nextSig){
+
+         // Calculate time as time of the largest digit
+         Float_t time = digit->GetTime() ;
+         Float_t aTime= digit->GetAmplitude() ;
+            
+         // loop over input
+         Int_t nInputs = fInput;
+         if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
+         for(i = 0; i< nInputs ; i++){  //loop over (possible) merge sources
+           TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
+           if(sdtclarr) {
+             Int_t sDigitEntries = sdtclarr->GetEntriesFast();
+                
+             if(sDigitEntries > index[i] )
+               curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;      
+             else
+               curSDigit = 0 ;
+                
+             //May be several digits will contribute from the same input
+             while(curSDigit && (curSDigit->GetId() == absID)){           
+               //Shift primary to separate primaries belonging different inputs
+               Int_t primaryoffset ;
+               if(fDigInput)
+                 primaryoffset = fDigInput->GetMask(i) ; 
+               else
+                 primaryoffset = i ;
+               curSDigit->ShiftPrimary(primaryoffset) ;
+                  
+               if(curSDigit->GetAmplitude()>aTime) {
+                 aTime = curSDigit->GetAmplitude();
+                 time  = curSDigit->GetTime();
+               }
+                  
+               *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
+                  
+               index[i]++ ;
+                  
+               if( sDigitEntries > index[i] )
+                 curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
+               else
+                 curSDigit = 0 ;
+             }// while
+           }// source exists
+         }// loop over merging sources
+            
+            
+         //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
+         energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
+            
+         // add fluctuations for photo-electron creation
+    // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
+    Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
+    energy       *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;    
     
-    if (geom->IsInECAL(digit->GetId())) 
-      amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
-    else if (geom->IsInPRE(digit->GetId())) 
-      amp += TMath::Abs(gRandom->Gaus(0., fPinNoise/100.)) ; // arbitrarely divide by 100, assuming that the gain of APD will be higher
-    else if (geom->IsInHCAL(digit->GetId())) 
-      amp += TMath::Abs(gRandom->Gaus(0., fPinNoise/10.)) ;  // arbitrarely divide by 10, assuming that the gain of APD will be higher
-   if (gDebug==1) 
-      Info("Digitize", "After noise amp = %f \n", amp) ; 
-    digit->SetAmp(sDigitizer->Digitize(amp)) ;  
+         //calculate and set time
+         digit->SetTime(time) ;
+            
+         //Find next signal module
+         nextSig = nEMC + 1 ;
+         for(i = 0 ; i < nInputs ; i++){
+           sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
+              
+           if(sdigits){
+             Int_t curNext = nextSig ;
+             if(sdigits->GetEntriesFast() > index[i])
+                {
+                  AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
+                  if(tmpdigit)
+                   {
+                     curNext = tmpdigit->GetId() ;       
+                   }
+                }
+             if(curNext < nextSig) nextSig = curNext ;
+           }// sdigits exist
+         } // input loop       
+            
+       }//absID==nextSig
+          
+       // add the noise now, no need for embedded events with real data
+       if(!embed)
+         energy += gRandom->Gaus(0., fPinNoise) ;
+          
+
+       // JLK 26-June-2008
+       //Now digitize the energy according to the fSDigitizer method,
+       //which merely converts the energy to an integer.  Later we will
+       //check that the stored value matches our allowed dynamic ranges
+       digit->SetAmplitude(fSDigitizer->Digitize(energy)) ; 
+          
+       //Set time resolution with final energy
+       timeResolution = GetTimeResolution(energy);
+       digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
+       AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
+                        absID, energy, nextSig));
+       //Add delay to time
+       digit->SetTime(digit->GetTime()+fTimeDelay) ;
+          
+      }// Digit pointer exists
+      else{
+       Info("Digitizer","Digit pointer is null");
+      }
+    } // for(absID = 0; absID < nEMC; absID++)
+      
+      //ticks->Delete() ;
+      //delete ticks ;
+
+      //Embed simulated amplitude (and time?) to real data digits
+    if(embed){
+      for(Int_t input = 1; input<fInput; input++){
+       TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
+       if(!realDigits){
+         AliDebug(1,"No sdigits to merge\n");
+         continue;
+       }
+       for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++){
+         AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
+         if(digit2){
+           digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
+           if(digit){
+                
+             // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
+             // and multiply to get a big int.
+             Float_t time2 = digit2->GetTime(); 
+             Float_t e2    = digit2->GetAmplitude();
+             CalibrateADCTime(e2,time2,digit2->GetId());
+             Float_t amp2 = fSDigitizer->Digitize(e2);
+                
+             Float_t e0    = digit ->GetAmplitude();
+             if(e0>0.01)
+               AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
+                               digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
+                               digit2->GetId(),amp2,                   digit2->GetType(), time2           ));
+                
+             // Sum amplitudes, change digit amplitude
+             digit->SetAmplitude( digit->GetAmplitude() + amp2);
+             // Rough assumption, assign time of the largest of the 2 digits, 
+             // data or signal to the final digit.
+             if(amp2 > digit->GetAmplitude())  digit->SetTime(time2);
+             //Mark the new digit as embedded
+             digit->SetType(AliEMCALDigit::kEmbedded);
+                
+             if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
+               AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
+                               digit->GetId(), digit->GetAmplitude(), digit->GetType()));
+           }//digit2 
+         }//digit2
+       }//loop digit 2
+      }//input loop
+    }//embed
+      
+    //JLK is it better to call Clear() here?
+    delete sdigArray ; //We should not delete its contents
+      
+    //remove digits below thresholds
+    // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
+    // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3,
+    // before merge in the same loop real data digits if available
+    Float_t energy = 0;
+    Float_t time   = 0;
+    for(i = 0 ; i < nEMC ; i++){
+      digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
+      if(digit){
+          
+       //Then get the energy in GeV units.
+       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
+       //Then digitize using the calibration constants of the ocdb
+       Float_t ampADC = energy;
+       DigitizeEnergyTime(ampADC, time, digit->GetId())  ;       
+       if(ampADC < fDigitThreshold)
+         digits->RemoveAt(i) ;
+          
+      }// digit exists
+    } // digit loop
+      
+    digits->Compress() ;  
+      
+    Int_t ndigits = digits->GetEntriesFast() ; 
+
+    //JLK 26-June-2008
+    //After we have done the summing and digitizing to create the
+    //digits, now we want to calibrate the resulting amplitude to match
+    //the dynamic range of our real data.  
+    for (i = 0 ; i < ndigits ; i++) { 
+      digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
+      if(digit){
+       digit->SetIndexInList(i) ; 
+       energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
+       time   = digit->GetTime();
+       Float_t ampADC = energy;
+       DigitizeEnergyTime(ampADC, time, digit->GetId());
+       digit->SetAmplitude(ampADC) ;
+       digit->SetTime(time);
+       // printf("digit amplitude set at end: i %d, amp %f\n",i,digit->GetAmplitude());
+      }// digit exists
+    }//Digit loop
+      
   }
+  else AliFatal("EMCALLoader is NULL!") ;
   
-  ticks->Delete() ;
-  delete ticks ;
-
-  delete sdigArray ; //We should not delete its contents
-
-  //remove digits below thresholds
-  for(absID = 0; absID < nEMC; absID++){
-    digit = dynamic_cast<AliEMCALDigit*>( digits->At(absID) ) ;
-    Float_t threshold = 0 ; 
-
-    if (geom->IsInECAL(digit->GetId())) 
-      threshold = fDigitThreshold ; 
-    else  if (geom->IsInPRE(digit->GetId()))
-      threshold = fDigitThreshold / 100. ; // arbitrary see before when noise is added
-    else  if (geom->IsInHCAL(digit->GetId()))
-      threshold = fDigitThreshold / 10. ; // arbitrary see before when noise is added    
-
-    if(sDigitizer->Calibrate( digit->GetAmp() ) <= threshold)
-      digits->RemoveAt(absID) ;
-    else {
-      if (gDebug == 1)
-       Info("Digitize", "id = %d, amp = %f, noise = %f", absID, sDigitizer->Calibrate(digit->GetAmp()), fPinNoise) ; 
-      digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
-    }
-  }
+}
+
+//_____________________________________________________________________
+void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, const Int_t absId)
+{ 
+  // JLK 26-June-2008
+  // Returns digitized value of the energy in a cell absId
+  // using the calibration constants stored in the OCDB
+  // or default values if no CalibData object is found.
+  // This effectively converts everything to match the dynamic range
+  // of the real data we will collect
+  //
+  // Load Geometry
+  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
   
-  digits->Compress() ;  
+  if (geom==0){
+    AliFatal("Did not get geometry from EMCALLoader");
+    return;
+  }
   
-  Int_t ndigits = digits->GetEntriesFast() ; 
-  //  digits->Expand(ndigits) ;
+  Int_t iSupMod = -1;
+  Int_t nModule = -1;
+  Int_t nIphi   = -1;
+  Int_t nIeta   = -1;
+  Int_t iphi    = -1;
+  Int_t ieta    = -1;
+  Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+  if(!bCell)
+    Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
+  geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
   
-  //Set indexes in list of digits
-  for (i = 0 ; i < ndigits ; i++) { 
-    digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
-    digit->SetIndexInList(i) ; 
-    Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ;
-    digit->SetAmp(DigitizeEnergy(energy,digit->GetId()) ) ;
+  if(fCalibData) {
+    fADCpedestalEC     = fCalibData->GetADCpedestal     (iSupMod,ieta,iphi);
+    fADCchannelEC      = fCalibData->GetADCchannel      (iSupMod,ieta,iphi);
+    fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
+    fTimeChannel       = fCalibData->GetTimeChannel     (iSupMod,ieta,iphi,0); // Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
+    fTimeChannelDecal  = fCalibData->GetTimeChannelDecal(iSupMod,ieta,iphi);      
   }
+  
+  //Apply calibration to get ADC counts and partial decalibration as especified in OCDB
+  energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal   ;
+  time  += fTimeChannel-fTimeChannelDecal;
+  
+  if(energy > fNADCEC ) 
+    energy =  fNADCEC ; 
 }
 
-//____________________________________________________________________________
+//_____________________________________________________________________
+void AliEMCALDigitizer::DecalibrateEnergy(Double_t& energy, const Int_t absId)
+{ 
+       // Load Geometry
+       const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+       
+       if (geom==0){
+               AliFatal("Did not get geometry from EMCALLoader");
+               return;
+       }
+       
+       Int_t iSupMod = -1;
+       Int_t nModule = -1;
+       Int_t nIphi   = -1;
+       Int_t nIeta   = -1;
+       Int_t iphi    = -1;
+       Int_t ieta    = -1;
+       Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+       if(!bCell)
+               Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
+       geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+       
+       if(fCalibData) {
+               fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
+       }
+
+       energy /= fADCchannelEC/0.0162;
+}
 
-Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t absId)
+//_____________________________________________________________________
+void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
 { 
-  Int_t channel = -999;
-  AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ; 
+  // Returns the energy in a cell absId with a given adc value
+  // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
+  // Used in case of embedding, transform ADC counts from real event into energy
+  // so that we can add the energy of the simulated sdigits which are in energy
+  // units.
+  // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
+  // or with time out of window
   
-  if(geom->IsInPRE(absId)){        //digitize as PRE section
-    channel =  static_cast<Int_t>(TMath::Ceil( (energy + fADCpedestalPRE)/fADCchannelPRE ))  ;
-    if(channel > fNADCPRE ) 
-      channel =  fNADCPRE ;
+  // Load Geometry
+  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+  
+  if (!geom){
+    AliFatal("Did not get geometry from EMCALLoader");
+    return;
   }
-  else if(geom->IsInECAL(absId)){  //digitize as ECAL section
-    channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
-    if(channel > fNADCEC ) 
-      channel =  fNADCEC ;
-  } 
-  else if(geom->IsInHCAL(absId)){  //digitize as HCAL section
-    channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalHC)/fADCchannelHC ))  ;
-    if(channel > fNADCHC ) 
-      channel =  fNADCHC ;
+  
+  Int_t iSupMod = -1;
+  Int_t nModule = -1;
+  Int_t nIphi   = -1;
+  Int_t nIeta   = -1;
+  Int_t iphi    = -1;
+  Int_t ieta    = -1;
+  Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+  if(!bCell)
+    Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
+  geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+  
+  if(fCalibData) {
+    fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
+    fADCchannelEC  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
+    fTimeChannel   = fCalibData->GetTimeChannel(iSupMod,ieta,iphi,0);// Assign bunch crossing number equal to 0 (has simulation different bunch crossings?)
   }
   
-  return channel ;
+  adc   = adc * fADCchannelEC - fADCpedestalEC;
+  time -= fTimeChannel;
+  
+  
 }
 
+
 //____________________________________________________________________________
-void AliEMCALDigitizer::Exec(Option_t *option) 
+void AliEMCALDigitizer::Digitize(Option_t *option) 
 { 
-  // Managing method
+  // Steering method to process digitization for events
+  // in the range from fFirstEvent to fLastEvent.
+  // This range is optionally set by SetEventRange().
+  // if fLastEvent=-1, then process events until the end.
+  // by default fLastEvent = fFirstEvent (process only one event)
+
+  if (!fInit) { // to prevent overwrite existing file
+    Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
+    return ;
+  } 
 
-  if(strcmp(GetName(), "") == 0 )   
-    Init() ;
   if (strstr(option,"print")) {
-    Print("");
+
+    Print();
     return ; 
   }
   
   if(strstr(option,"tim"))
     gBenchmark->Start("EMCALDigitizer");
 
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
-  
-  Int_t nevents ;
-  
-  TTree * treeD ;
-  
-  if(fManager){
-    treeD = fManager->GetTreeD() ;
-    nevents = 1 ;    // Will process only one event
-    
-    //Check, if this branch already exits
-    if (treeD) { 
-      TObjArray * lob = (TObjArray*)treeD->GetListOfBranches() ;
-      TIter next(lob) ; 
-      TBranch * branch = 0 ;  
-      Bool_t emcalfound = kFALSE, digitizerfound = kFALSE ; 
-      
-      while ( (branch = (TBranch*)next()) && (!emcalfound || !digitizerfound) ) {
-       if ( (strcmp(branch->GetName(), "EMCAL")==0) && 
-            (strcmp(branch->GetTitle(), GetName())==0) ) 
-         emcalfound = kTRUE ;
-       
-       else if ( (strcmp(branch->GetName(), "AliEMCALDigitizer")==0) && 
-                 (strcmp(branch->GetTitle(), GetName())==0) ) 
-         digitizerfound = kTRUE ; 
-      }
-      
-      if ( emcalfound ) {
-       Error( "Exec", "Digits branch with name %s already exits", GetName() ) ;
-       return ; 
-      }   
-      if ( digitizerfound ) {
-       Error( "Exec", "Digitizer branch with name %s already exit", GetName() ) ;
-              return ; 
-      }
-    }   
+  AliRunLoader *rl = AliRunLoader::Instance();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
+  Int_t nEvents = 0;
+  if(!emcalLoader){
+    AliFatal("Did not get the  Loader");
   }
-  else { //EMCAL standalone
-    if(gime->BranchExists("Digits") ) 
-      return ;
-    nevents=gime->MaxEvent() ;
-  }
-  
-  Int_t ievent ;
-
-  for(ievent = 0; ievent < nevents; ievent++){
+  else{
     
-    if(fManager){
-
-      Int_t input ;
-      for(input = 0 ; input < fManager->GetNinputs(); input ++){
-       TTree * treeS = fManager->GetInputTreeS(input) ;
-       if(!treeS){
-         Error( "Exec", "No Input") ;
-         return ;
-       }
-       gime->ReadTreeS(treeS,input) ;
-      }
-
+    if (fLastEvent == -1)  {
+      fLastEvent = rl->GetNumberOfEvents() - 1 ;
     }
-    else 
-      gime->Event(ievent,"S") ; 
+    else if (fDigInput) 
+      fLastEvent = fFirstEvent ; // what is this ??
     
-    Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
-    
-    WriteDigits(ievent) ;
-    
-    if(strstr(option,"deb"))
-      PrintDigits(option);
-    
-    //increment the total number of Digits per run 
-    fDigitsInRun += gime->Digits()->GetEntriesFast() ;  
-  }
+    nEvents = fLastEvent - fFirstEvent + 1;
+    Int_t ievent  = -1;
+
+    TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit",    32*96);
+    TClonesArray* digitsTRG = new TClonesArray("AliEMCALRawDigit", 32*96);
+
+    rl->LoadSDigits("EMCAL");
+    for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+      
+      rl->GetEvent(ievent);
+      
+      Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
+      
+      WriteDigits() ;
+      
+      //Trigger Digits
+      //-------------------------------------
+     
+
+      Digits2FastOR(digitsTMP, digitsTRG);  
+      
+      WriteDigits(digitsTRG);
+      
+      (emcalLoader->TreeD())->Fill();
+      
+      emcalLoader->WriteDigits(   "OVERWRITE");
+      
+      Unload();
+      
+      digitsTRG  ->Delete();
+      digitsTMP  ->Delete();
+      
+      //-------------------------------------
+      
+      if(strstr(option,"deb"))
+        PrintDigits(option);
+      if(strstr(option,"table")) gObjectTable->Print();
+      
+      //increment the total number of Digits per run 
+      fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;  
+    }//loop
+       
+  }//loader exists
   
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALDigitizer");
-    Info("Exec", "took %f seconds for Digitizing %f seconds per event", 
-        gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nevents ) ;
-  }
-  
+    Float_t cputime   = gBenchmark->GetCpuTime("EMCALDigitizer");
+    Float_t avcputime = cputime;
+    if(nEvents==0) avcputime = 0 ;
+    AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
+  } 
 }
 
-//____________________________________________________________________________ 
-Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks) 
-{ // 
-  ticks->Sort() ; //Sort in accordance with times of ticks
-  TIter it(ticks) ;
-  AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
-  Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
-  
-  AliEMCALTick * t ;  
-  while((t=(AliEMCALTick*) it.Next())){
-    if(t->GetTime() < time)  //This tick starts before crossing
-      *ctick+=*t ;
-    else
-      return time ;
-    
-    time = ctick->CrossingTime(fTimeThreshold) ;    
+//__________________________________________________________________
+Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
+{  
+  // From F. Blanco
+  Float_t res = -1;
+  if (energy > 0) {
+    res = TMath::Sqrt(fTimeResolutionPar0 + 
+                     fTimeResolutionPar1/(energy*energy) );
   }
-  return time ;
+  // parametrization above is for ns. Convert to seconds:
+  return res*1e-9;
 }
 
 //____________________________________________________________________________ 
-Bool_t AliEMCALDigitizer::Init()
+void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
 {
-  // Makes all memory allocations
-
-  if( strcmp(GetTitle(), "") == 0 )
-    SetTitle("galice.root") ;
+  // FEE digits afterburner to produce TRG digits 
+  // we are only interested in the FEE digit deposited energy
+  // to be converted later into a voltage value
   
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), fToSplit) ; 
-  if ( gime == 0 ) {
-    Error("Init", "Could not obtain the Getter object !" ) ; 
-    return kFALSE;
-  } 
+  // push the FEE digit to its associated FastOR (numbered from 0:95)
+  // TRU is in charge of summing module digits
   
-  //const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
-
-  //fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
+  AliRunLoader *runLoader = AliRunLoader::Instance();
   
-  // Post Digits to the white board
-  gime->PostDigits(GetName() ) ;   
+  AliRun* run = runLoader->GetAliRun();
   
-  // Post Digitizer to the white board
-  gime->PostDigitizer(this) ;
-  fSplitFile = 0 ;
-  if(fToSplit){
-    // construct the name of the file as /path/EMCAL.SDigits.root
-    //First - extract full path if necessary
-    TString digitsFileName(GetTitle()) ;
-    Ssiz_t islash = digitsFileName.Last('/') ;
-    if(islash<digitsFileName.Length())
-      digitsFileName.Remove(islash+1,digitsFileName.Length()) ;
-    else
-      digitsFileName="" ;
-    // Next - append the file name 
-    digitsFileName+="EMCAL.Digits." ;
-    if((strcmp(GetName(),"Default")!=0)&&(strcmp(GetName(),"")!=0)){
-      digitsFileName+=GetName() ;
-      digitsFileName+="." ;
-    }
-    digitsFileName+="root" ;
-    // Finally - check if the file already opened or open the file
-    fSplitFile = static_cast<TFile*>(gROOT->GetFile(digitsFileName.Data()));   
-    if(!fSplitFile)
-      fSplitFile =  TFile::Open(digitsFileName.Data(),"update") ;
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
+  if(!emcalLoader){
+    AliFatal("Did not get the  Loader");
   }
+  else {
+    AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(run->GetDetector("EMCAL"));
+    if(emcal){
+      AliEMCALGeometry* geom = emcal->GetGeometry();
+      
+      // build FOR from simulated digits
+      // and xfer to the corresponding TRU input (mapping)
+      
+      TClonesArray* digits = emcalLoader->SDigits();
+      
+         AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
+               
+      TIter NextDigit(digits);
+      while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
+      {
+        Int_t id = digit->GetId();
+        
+        Int_t trgid;
+        if (geom && geom->GetFastORIndexFromCellIndex(id , trgid)) 
+        {
+          AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
+          
+          AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
+          
+          if (!d)
+          {
+            new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
+            d = (AliEMCALDigit*)digitsTMP->At(trgid);
+            d->SetId(trgid);
+          }    
+          else
+          {
+            *d = *d + *digit;
+          }
+        }
+      }
+      
+      if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
+      
+         TF1 g4noiseF("g4noise","[0]*exp(-0.5*((x-[1])/[2])**4)",-10,10);
+         g4noiseF.SetParameter( 0,   1 );  //
+         g4noiseF.SetParameter( 1,   0 );  // 
+         g4noiseF.SetParameter( 2,  .2 );  // 
+               
+      Int_t    nSamples = 32;
+      Int_t *timeSamples = new Int_t[nSamples];
+      
+      NextDigit = TIter(digitsTMP);
+      while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
+      {
+        if (digit)
+        {
+          Int_t     id = digit->GetId();
+          Float_t time = 50.e-9;
+          
+          Double_t depositedEnergy = 0.;
+          for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
+          
+          if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
+          
+          // FIXME: Check digit time!
+          if (depositedEnergy)
+          {
+//                     depositedEnergy += depositedEnergy * g4noiseF.GetRandom();
+                         
+                   DecalibrateEnergy(depositedEnergy, id);  
+                         
+            DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
+            
+            for (Int_t j=0;j<nSamples;j++) 
+            {
+              timeSamples[j] = ((j << 12) & 0xFF000) | (timeSamples[j] & 0xFFF);
+            }
+            
+            new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALRawDigit(id, timeSamples, nSamples);
+                         
+                       if (AliDebugLevel()) ((AliEMCALRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
+          }
+        }
+      }
+      
+      delete [] timeSamples;
+    }// AliEMCAL exists
+    else AliFatal("Could not get AliEMCAL");
+  }// loader exists
   
-  //Mark that we will use current header file
-  if(!fManager){
-    gime->PostSDigits(GetName(),GetTitle()) ;
-    gime->PostSDigitizer(GetName(),GetTitle()) ;
-  }
-  return kTRUE ;    
 }
 
 //____________________________________________________________________________ 
-void AliEMCALDigitizer::InitParameters()
+void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
 {
-  fMeanPhotonElectron = 1250 ; // electrons per GeV
-  
-  Warning("InitParameters", "No noise added\n") ; 
-  fPinNoise           = 0. ; // 0.001 ; // noise equivalent GeV (random choice)
-  fDigitThreshold     = fPinNoise * 3; //2 sigma
-  fTimeResolution     = 0.5e-9 ;
-  fTimeSignalLength   = 1.0e-9 ;
-
-  fADCchannelEC    = 0.000220;                     // width of one ADC channel in GeV
-  fADCpedestalEC   = 0.005 ;                       // GeV
-  fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
-
-  fADCchannelHC    = 0.000220;                     // width of one ADC channel in GeV
-  fADCpedestalHC   = 0.005 ;                       // GeV
-  fNADCHC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
-
-  fADCchannelPRE   = 0.0000300;                    // width of one ADC channel in Pre Shower
-  fADCpedestalPRE  = 0.005 ;                       // GeV 
-  fNADCPRE         = (Int_t) TMath::Power(2,12);   // number of channels in Pre ShowerADC
-
-  fTimeThreshold      = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
+       // parameters:  
+       // id: 0..95
+       const Int_t    reso = 11;      // 11-bit resolution ADC
+       const Double_t vFSR = 1;       // Full scale input voltage range
+//     const Double_t dNe  = 125;     // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
+       const Double_t dNe  = 125/1.3; // F-ALTRO max V. FEE: factor 4 
+       const Double_t vA   = .136e-6; // CSP output range: 0.136uV/e-
+       const Double_t rise = 50e-9;   // rise time (10-90%) of the FastOR signal before shaping
+       
+       const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
+       
+       Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
+               
+       TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
+       signalF.SetParameter( 0,   vV ); 
+       signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
+       signalF.SetParameter( 2, rise );
+       
+       for (Int_t iTime=0; iTime<nSamples; iTime++) 
+       {
+               // FIXME: add noise (probably not simply Gaussian) according to DA measurements
+               // probably plan an access to OCDB
+               
+               timeSamples[iTime] = int((TMath::Power(2, reso) / vFSR) * signalF.Eval(iTime * kTimeBinWidth) + 0.5);
+       }
 }
 
-//__________________________________________________________________
-void AliEMCALDigitizer::MixWith(char* headerFile)
+//____________________________________________________________________________ 
+Bool_t AliEMCALDigitizer::Init()
 {
-  // Allows to produce digits by superimposing background and signal event.
-  // It is assumed, that headers file with SIGNAL events is opened in 
-  // the constructor. 
-  // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed 
-  // Thus we avoid writing (changing) huge and expensive 
-  // backgound files: all output will be writen into SIGNAL, i.e. 
-  // opened in constructor file. 
-  //
-  // One can open as many files to mix with as one needs.
-  // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
-  // can be mixed.
-
-  if( strcmp(GetName(), "") == 0 )
-    Init() ;
+  // Makes all memory allocations
+  fInit = kTRUE ; 
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
   
-  if(fManager){
-    Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
-    return ;
+  if ( emcalLoader == 0 ) {
+    Fatal("Init", "Could not obtain the AliEMCALLoader");  
+    return kFALSE;
   } 
   
-  // check if the specified SDigits do not already exist on the White Board:
-  // //Folders/RunMC/Event/Data/EMCAL/SDigits/headerFile/sdigitsname
-
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  TString path = gime->SDigitsFolder()->GetName() ; 
-
-  // before it was ???? "Folders/RunMC/Event/Data/EMCAL/SDigits" ; 
-  path += headerFile ; 
-  path += "/" ; 
-  path += GetName() ;
-  if ( gROOT->FindObjectAny(path.Data()) ) {
-    Error("MixWith", "Entry already exists, do not add" ) ;
-    return;
+  fFirstEvent = 0 ; 
+  fLastEvent = fFirstEvent ; 
+  
+  if(fDigInput)
+    fInput = fDigInput->GetNinputs() ; 
+  else 
+    fInput           = 1 ;
+
+  fInputFileNames  = new TString[fInput] ; 
+  fEventNames      = new TString[fInput] ; 
+  fInputFileNames[0] = GetTitle() ; 
+  fEventNames[0]     = fEventFolderName.Data() ; 
+  Int_t index ; 
+  for (index = 1 ; index < fInput ; index++) {
+    fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0); 
+    TString tempo = fDigInput->GetInputFolderName(index) ;
+    fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput 
   }
+    
+  //Calibration instance
+  fCalibData = emcalLoader->CalibData();
+  return fInit ;    
+}
 
-  gime->PostSDigits(GetName(),headerFile) ;
+//____________________________________________________________________________ 
+void AliEMCALDigitizer::InitParameters()
+{ 
+  // Parameter initialization for digitizer
   
-  // check if the requested file is already open or exist and if SDigits Branch exist
-  TFile * file = (TFile*)gROOT->FindObject(headerFile); 
-  if ( !file ) { 
-    file = new TFile(headerFile, "READ") ; 
-    if (!file) { 
-      Error("MixWith", "File %s does not exist!", headerFile) ; 
-      return ; 
-    }
+  // Get the parameters from the OCDB via the loader
+  AliRunLoader *rl = AliRunLoader::Instance();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
+  AliEMCALSimParam * simParam = 0x0;
+  if(emcalLoader) simParam = emcalLoader->SimulationParameters();
+       
+  if(!simParam){
+    simParam = AliEMCALSimParam::GetInstance();
+    AliWarning("Simulation Parameters not available in OCDB?");
   }
   
+  fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400;  // electrons per GeV 
+  fGainFluctuations   = simParam->GetGainFluctuations()   ; // 15.0; 
+  
+  fPinNoise           = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data 
+  if (fPinNoise < 0.0001 ) 
+    Warning("InitParameters", "No noise added\n") ; 
+  fTimeNoise          = simParam->GetTimeNoise(); // 1.28E-5 s
+  fDigitThreshold     = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
+  fTimeResolutionPar0 = simParam->GetTimeResolutionPar0(); 
+  fTimeResolutionPar1 = simParam->GetTimeResolutionPar1(); 
+  fTimeDelay          = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
+  
+  // These defaults are normally not used. 
+  // Values are read from calibration database instead
+  fADCchannelEC       = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
+  fADCpedestalEC      = 0.0 ;   // GeV
+  fADCchannelECDecal  = 1.0;    // No decalibration by default
+  fTimeChannel        = 0.0;    // No time calibration by default
+  fTimeChannelDecal   = 0.0;    // No time decalibration by default
+
+  fNADCEC             = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
+  
+  AliDebug(2,Form("Mean Photon Electron %d, Gain Fluct. %2.1f; Noise: APD %f, Time %f; Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
+                  fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
+  
 }
 
 //__________________________________________________________________
-void AliEMCALDigitizer::Print(Option_t* option)const {
-
-  TString message("\n") ; 
-
-  if( strcmp(GetName(), "") != 0) {
-    message += "------------------- " ; 
-    message += GetName() ; 
-    message += " -------------" ;
-    const Int_t nStreams = GetNInputStreams() ; 
-    if (nStreams) {
-      Int_t index = 0 ;  
-      for (index = 0 ; index < nStreams ; index++) {  
-       message += "\nAdding SDigits " ; 
-       message += GetName() ; 
-       message += " from " ; 
-       message += fManager->GetInputFileName(index, 0) ;
-      } 
-     
-      message += "\nWriting digits to " ; 
-      message += fManager->GetInputFileName(0, 0) ;   
-    } else { 
-      message += "\nWriting digits to " ;
-      message += GetTitle() ;
-    }       
-    message += "\nWith following parameters: " ;
-    message += "\n     Electronics noise in EMC (fPinNoise) = " ; 
-    message += fPinNoise ;
-    message += "\n  Threshold  in EMC  (fDigitThreshold) = " ; 
-    message += fDigitThreshold  ;
-    message += "\n---------------------------------------------------"  ;
+void AliEMCALDigitizer::Print1(Option_t * option)
+{ // 19-nov-04 - just for convenience
+  Print(); 
+  PrintDigits(option);
+}
+
+//__________________________________________________________________
+void AliEMCALDigitizer::Print (Option_t * ) const
+{
+  // Print Digitizer's parameters
+  printf("Print: \n------------------- %s -------------", GetName() ) ; 
+  if( strcmp(fEventFolderName.Data(), "") != 0 ){
+    printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
+    
+    Int_t nStreams ; 
+    if (fDigInput) 
+      nStreams =  GetNInputStreams() ;
+    else 
+      nStreams = fInput ; 
+    
+    AliRunLoader *rl=0;
+    
+    Int_t index = 0 ;  
+    for (index = 0 ; index < nStreams ; index++) {  
+      TString tempo(fEventNames[index]) ; 
+      tempo += index ;
+      if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
+        rl = AliRunLoader::Open(fInputFileNames[index], tempo) ; 
+      AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
+      if(emcalLoader){
+        TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
+        if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
+          fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
+        printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ; 
+      }//loader
+    }
+    
+    AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+    
+    if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
+    else printf("\nNULL LOADER");
+    
+    printf("\nWith following parameters:\n") ;
+    printf("    Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
+    printf("    Threshold  in Tower  (fDigitThreshold) = %d\n", fDigitThreshold)  ;
+    printf("---------------------------------------------------\n")  ;
   }
   else
-    message += "\nAliEMCALDigitizer not initialized " ;
-
-  Info("Print", message.Data() ) ; 
+    printf("Print: AliEMCALDigitizer not initialized") ; 
 }
 
 //__________________________________________________________________
-void AliEMCALDigitizer::PrintDigits(Option_t * option){
+void AliEMCALDigitizer::PrintDigits(Option_t * option)
+{
+  //utility method for printing digit information
+  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  if(emcalLoader){
+    TClonesArray * digits  = emcalLoader->Digits() ;
+    TClonesArray * sdigits = emcalLoader->SDigits() ;
     
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  TClonesArray * fDigits = gime->Digits() ;
-
-  TString message("\n") ; 
-  message += "       Number of entries in Digits list " ; 
-  message += fDigits->GetEntriesFast() ;
-
-  if(strstr(option,"all")){
+    printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ; 
+    printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
     
-    //loop over digits
-    AliEMCALDigit * digit;
-    message += "\n   Id  Amplitude    Time          Index Nprim: Primaries list \n" ;    
-    Int_t index ;
-    char * tempo = new char[8192]; 
-    for (index = 0 ; index < fDigits->GetEntries() ; index++) {
-      digit = (AliEMCALDigit * )  fDigits->At(index) ;
-      sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d : ",
-             digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
-      message += tempo ; 
-      Int_t iprimary;
-      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
-       sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
-       message += tempo ; 
-      }
-    }   
-    delete tempo ; 
-  }
-  Info("PrintDigits", message.Data() ) ; 
+    if(strstr(option,"all")){  
+      //loop over digits
+      AliEMCALDigit * digit;
+      printf("\nEMC digits (with primaries):\n")  ;
+      printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
+      Int_t index ;
+      for (index = 0 ; index < digits->GetEntries() ; index++) {
+        digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
+        if(digit){
+          printf("\n%6d  %8f    %6.5e %4d      %2d : ",
+                 digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
+          Int_t iprimary;
+          for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
+            printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
+          }
+        }// digit exists
+      }// loop
+    }
+    printf("\n");
+  }// loader exists
+  else printf("NULL LOADER, cannot print\n");
 }
 
 //__________________________________________________________________
 Float_t AliEMCALDigitizer::TimeOfNoise(void)
-{  // Calculates the time signal generated by noise
-  //to be rewritten, now returns just big number
-  return 1. ;
+{  
+  // Calculates the time signal generated by noise
+  //printf("Time noise %e\n",fTimeNoise);
+  return gRandom->Rndm() * fTimeNoise;
+}
 
+//__________________________________________________________________
+void AliEMCALDigitizer::Unload() 
+{  
+  // Unloads the SDigits and Digits
+  AliRunLoader *rl=0;
+    
+  Int_t i ; 
+  for(i = 1 ; i < fInput ; i++){
+    TString tempo(fEventNames[i]) ; 
+    tempo += i;
+    if ((rl = AliRunLoader::GetRunLoader(tempo))) 
+      rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ; 
+  }
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  if(emcalLoader)emcalLoader->UnloadDigits() ; 
+  else AliFatal("Did not get the loader");
 }
 
 //_________________________________________________________________________________________
-void AliEMCALDigitizer::WriteDigits(Int_t event)
-{
-
-  // Makes TreeD in the output file. 
+void AliEMCALDigitizer::WriteDigits()
+{ // Makes TreeD in the output file. 
   // Check if branch already exists: 
   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
   //      already existing branches. 
   //   else creates branch with Digits, named "EMCAL", title "...",
   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
   //      and names of files, from which digits are made.
+  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  
+  if(emcalLoader){
+    const TClonesArray * digits = emcalLoader->Digits() ; 
+    TTree * treeD = emcalLoader->TreeD(); 
+    if ( !treeD ) {
+      emcalLoader->MakeDigitsContainer();
+      treeD = emcalLoader->TreeD(); 
+    }
+    
+    // -- create Digits branch
+    Int_t bufferSize = 32000 ;    
+    TBranch * digitsBranch = 0;
+    if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
+      digitsBranch->SetAddress(&digits);
+      AliWarning("Digits Branch already exists. Not all digits will be visible");
+    }
+    else
+      treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
+    //digitsBranch->SetTitle(fEventFolderName);
+    
+    // treeD->Fill() ;
+    /*  
+     emcalLoader->WriteDigits("OVERWRITE");
+     emcalLoader->WriteDigitizer("OVERWRITE");
+     
+     Unload() ; 
+     */
+    
+  }// loader exists
+  else AliFatal("Loader not available");
+}
 
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  const TClonesArray * digits = gime->Digits(GetName()) ; 
- TTree * treeD ;
-
- if(fManager)
-   treeD = fManager->GetTreeD() ;
- else {
-   if ((gAlice->TreeD() == 0) || (fSplitFile)) // we should not create TreeD if it is already here 
-     gAlice->MakeTree("D", fSplitFile); // We overwrite TreeD in split file in the case of second reconstruction
-   if(fSplitFile)
-     fSplitFile->cd() ;
-   treeD = gAlice->TreeD();
- }
- // -- create Digits branch
- Int_t bufferSize = 32000 ;    
- TBranch * digitsBranch = treeD->Branch("EMCAL",&digits,bufferSize);
- digitsBranch->SetTitle(GetName());
- // -- Create Digitizer branch
- Int_t splitlevel = 0 ;
- const  AliEMCALDigitizer * d = gime->Digitizer(GetName()) ;
- TBranch * digitizerBranch = treeD->Branch("AliEMCALDigitizer", "AliEMCALDigitizer", &d,bufferSize,splitlevel); 
- digitizerBranch->SetTitle(GetName());
- digitsBranch->Fill() ;
- digitizerBranch->Fill() ; 
- treeD->AutoSave() ;  
+//__________________________________________________________________
+void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
+{ // overloaded method
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  if(emcalLoader){
+    
+    TTree* treeD = emcalLoader->TreeD(); 
+    if (!treeD) 
+      {
+       emcalLoader->MakeDigitsContainer();
+       treeD = emcalLoader->TreeD(); 
+      }
+    
+    // -- create Digits branch
+    Int_t bufferSize = 32000;
+    
+    if (TBranch* triggerBranch = treeD->GetBranch(branchName)) 
+      {
+       triggerBranch->SetAddress(&digits);
+      }
+    else
+      {
+       treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
+      }
+    
+    // treeD->Fill();
+  }// loader exists
+  else AliFatal("Loader not available");
 }