]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALSDigitizer.cxx
Possibility to
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
index 0961a29c87ff624328fa5442739b5987dcfba4d5..09b9d35974810ab81529b803d2c57629accf959a 100644 (file)
@@ -1,4 +1,4 @@
-/**************************************************************************
+/*************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
 //_________________________________________________________________________
 // This is a TTask that makes SDigits out of Hits
 // A Summable Digits is the sum of all hits originating 
-// from one primary in one active cell
+// from one in one tower of the EMCAL 
 // A threshold for assignment of the primary to SDigit is applied 
+//
+// JLK 26-Jun-2008 Added explanation:
+// SDigits need to hold the energy sum of the hits, but AliEMCALDigit
+// can (should) only store amplitude.  Therefore, the SDigit energy is
+// "digitized" before being stored and must be "calibrated" back to an
+// energy before SDigits are summed to form true Digits
+//
 // SDigits are written to TreeS, branch "EMCAL"
 // AliEMCALSDigitizer with all current parameters is written 
 // to TreeS branch "AliEMCALSDigitizer".
 //             deb all  - print # and list of produced SDigits
 //             tim - print benchmarking information
 //
-//*-- Author :  Dmitri Peressounko (SUBATECH & KI) 
+//*-- Author : Sahal Yacoob (LBL)
+// based on  : AliPHOSSDigitzer 
 //////////////////////////////////////////////////////////////////////////////
 
-
 // --- ROOT system ---
-#include "TFile.h"
-#include "TTask.h"
-#include "TTree.h"
-#include "TSystem.h"
-#include "TROOT.h"
-#include "TFolder.h"
-#include "TBenchmark.h"
+#include <TBenchmark.h>
+#include <TBrowser.h>
+//#include <Riostream.h>
+#include <TMath.h>
+#include <TROOT.h>
+
 // --- Standard library ---
-#include <iomanip.h>
+#include "stdlib.h"
 
 // --- AliRoot header files ---
-#include "AliRun.h"
+#include "AliLog.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"
 #include "AliEMCALDigit.h"
+#include "AliEMCALLoader.h"
 #include "AliEMCALHit.h"
 #include "AliEMCALSDigitizer.h"
 #include "AliEMCALGeometry.h"
-#include "AliEMCALv1.h"
+#include "AliEMCALSimParam.h"
 
 ClassImp(AliEMCALSDigitizer)
-
            
 //____________________________________________________________________________ 
-  AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","") 
+AliEMCALSDigitizer::AliEMCALSDigitizer()
+  : TTask("",""),
+    fA(0.),fB(0.),fECPrimThreshold(0.),
+    fDefaultInit(kTRUE),
+    fEventFolderName(0),
+    fInit(0),
+    fSDigitsInRun(0),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fSampling(0.)
 {
   // ctor
-  fA = 0;
-  fB = 10000000. ;
-  fPrimThreshold = 0.001 ;
-  fNevents = 0 ;     
-  fSDigits = 0 ;
-  fHits = 0 ;
-  fIsInitialized = kFALSE ;
-
+  InitParameters();
 }
 
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliEMCALSDigitizer","")
+AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, 
+                                      const char * eventFolderName)
+  : TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
+    fA(0.),fB(0.),fECPrimThreshold(0.),
+    fDefaultInit(kFALSE),
+    fEventFolderName(eventFolderName),
+    fInit(0),
+    fSDigitsInRun(0),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fSampling(0.)
 {
   // ctor
-  fA = 0;
-  fB = 10000000.;
-  fPrimThreshold = 0.001 ;
-  fNevents = 0 ;      
-  fSDigitsTitle = sDigitsTitle ;
-  fHeadersFile = headerFile ;
-  fSDigits = new TClonesArray("AliEMCALDigit",1000);
-  fHits    = new TClonesArray("AliEMCALHit",1000);
-
-  TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
-  
-  //File was not opened yet
-  if(file == 0){
-    if(fHeadersFile.Contains("rfio"))
-      file =   TFile::Open(fHeadersFile,"update") ;
-    else
-      file = new TFile(fHeadersFile.Data(),"update") ;
-    gAlice = (AliRun *) file->Get("gAlice") ;
-  }
-  
-  //add Task to //root/Tasks folder
-  TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
-  roottasks->Add(this) ; 
-  
-  fIsInitialized = kTRUE ;
+  Init();
+  InitParameters() ; 
+
 }
 
+
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::~AliEMCALSDigitizer()
+AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) 
+  : TTask(sd.GetName(),sd.GetTitle()),
+    fA(sd.fA),
+    fB(sd.fB),
+    fECPrimThreshold(sd.fECPrimThreshold),
+    fDefaultInit(sd.fDefaultInit),
+    fEventFolderName(sd.fEventFolderName),
+    fInit(sd.fInit),
+    fSDigitsInRun(sd.fSDigitsInRun),
+    fFirstEvent(sd.fFirstEvent),
+    fLastEvent(sd.fLastEvent),
+    fSampling(sd.fSampling)
 {
-  // dtor
-  if(fSDigits)
-    delete fSDigits ;
-  if(fHits)
-    delete fHits ;
+  //cpy ctor 
 }
+
+
+//____________________________________________________________________________ 
+AliEMCALSDigitizer::~AliEMCALSDigitizer() {
+  //dtor
+  AliLoader *emcalLoader=0;
+  if ((emcalLoader = AliRunLoader::Instance()->GetDetectorLoader("EMCAL")))
+      emcalLoader->CleanSDigitizer();
+}
+
 //____________________________________________________________________________ 
 void AliEMCALSDigitizer::Init(){
-  //Initialization can not be done in the default constructor
+  // Initialization: open root-file, allocate arrays for hits and sdigits,
+  // attach task SDigitizer to the list of EMCAL tasks
+  // 
+  // Initialization can not be done in the default constructor
+  //============================================================= YS
+  //  The initialisation is now done by the getter
 
-  if(!fIsInitialized){
+  fInit = kTRUE; 
+   
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 
-    if(fHeadersFile.IsNull())
-      fHeadersFile="galice.root" ;
+  if ( emcalLoader == 0 ) {
+    Fatal("Init", "Could not obtain the AliEMCALLoader");
+    return ;
+  } 
+  
+  emcalLoader->PostSDigitizer(this);
+  emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
+  
+}
 
-    TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
-    
-    //if file was not opened yet, read gAlice
-    if(file == 0){
-      file = new TFile(fHeadersFile.Data(),"update") ;
-      gAlice = (AliRun *) file->Get("gAlice") ;
-    }
-    
-    fHits    = new TClonesArray("AliEMCALHit",1000);
-    fSDigits = new TClonesArray("AliEMCALDigit",1000);
-    
-    // add Task to //root/Tasks folder
-    TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
-    roottasks->Add(this) ; 
-    
-    fIsInitialized = kTRUE ;
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::InitParameters()
+{
+  //initialize parameters for sdigitization
+
+  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+  if (geom->GetSampling() == 0.) {
+    Fatal("InitParameters", "Sampling factor not set !") ; 
   }
+
+  // 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?");
+  }
+       
+  //
+  //JLK 26-Jun-2008 THIS SHOULD HAVE BEEN EXPLAINED AGES AGO:
+  //
+  //In order to be able to store SDigit Energy info into
+  //AliEMCALDigit, we need to convert it temporarily to an ADC amplitude
+  //and later when summing SDigits to form digits, convert it back to
+  //energy.  These fA and fB parameters accomplish this through the
+  //Digitize() and Calibrate() methods
+  //
+  // Initializes parameters
+  fA         = simParam->GetA(); //0;
+  fB         = simParam->GetB(); //1.e+6;  // Changed 24 Apr 2007. Dynamic range now 2 TeV
+  fSampling  = geom->GetSampling();
+
+  // threshold for deposit energy of hit
+  fECPrimThreshold  = simParam->GetECPrimaryThreshold();//0.05;// GeV // 22-may-07 was 0// 24-nov-04 - was 1.e-6;
+  
+  AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
+  AliDebug(2,Form("   fInit                                 %i\n", int(fInit)));
+  AliDebug(2,Form("   fFirstEvent                           %i\n", fFirstEvent));
+  AliDebug(2,Form("   fLastEvent                            %i\n", fLastEvent));
+  AliDebug(2,Form("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()));
+  AliDebug(2,Form("   with digitization parameters       A = %f\n", fA));
+  AliDebug(2,Form("                                      B = %f\n", fB));
+  AliDebug(2,Form("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold));
+  AliDebug(2,Form("   Sampling                             = %f\n", fSampling));
+  AliDebug(2,Form("---------------------------------------------------\n"));
+
 }
+
 //____________________________________________________________________________
-void AliEMCALSDigitizer::Exec(Option_t *option) { 
-  //Collects all hits in the same active volume into digit
-  
-  if(!fIsInitialized)
-    Init() ;
+void AliEMCALSDigitizer::Exec(Option_t *option) 
+{ 
+  // Collects all hit of the same tower into digits
+  TString o(option); o.ToUpper();
+  if (strstr(option, "print") ) {
+
+    AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
+    AliDebug(2,Form("   fInit                                 %i\n", int(fInit)));
+    AliDebug(2,Form("   fFirstEvent                           %i\n", fFirstEvent));
+    AliDebug(2,Form("   fLastEvent                            %i\n", fLastEvent));
+    AliDebug(2,Form("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()));
+    AliDebug(2,Form("   with digitization parameters       A = %f\n", fA));
+    AliDebug(2,Form("                                      B = %f\n", fB));
+    AliDebug(2,Form("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold));
+    AliDebug(2,Form("   Sampling                             = %f\n", fSampling));
+    AliDebug(2,Form("---------------------------------------------------\n"));
+
+    return ; 
+  }
 
   if(strstr(option,"tim"))
     gBenchmark->Start("EMCALSDigitizer");
-  
-  fNevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
-  
-  Int_t ievent ;
-  for(ievent = 0; ievent < fNevents; ievent++){
-    gAlice->GetEvent(ievent) ;
-    gAlice->SetEvent(ievent) ;
-    
-    if(gAlice->TreeH()==0){
-      cout << "AliEMCALSDigitizer: There is no Hit Tree" << endl;
-      return ;
-    }
-    
-    //set address of the hits 
-    TBranch * branch = gAlice->TreeH()->GetBranch("EMCAL");
-    if (branch) 
-      branch->SetAddress(&fHits);
-    else{
-      cout << "ERROR in AliEMCALSDigitizer: "<< endl ;
-      cout << "      no branch EMCAL in TreeH"<< endl ;
-      cout << "      do nothing " << endl ;
-      return ;
+
+  AliRunLoader *rl = AliRunLoader::Instance();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
+
+  //switch off reloading of this task while getting event
+  if (!fInit) { // to prevent overwrite existing file
+    AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
+    return ;
     }
-    
-    fSDigits->Clear();
-    Int_t nSdigits = 0 ;
-    
-    
-    //Now made SDigits from hits, for EMCAL it is the same, so just copy    
-    Int_t itrack ;
-    for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
-      
-      //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
-      branch->GetEntry(itrack,0);
-      AliEMCAL * EMCAL = (AliEMCAL *) gAlice->GetDetector("EMCAL") ;
-      AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance( EMCAL->GetGeometry()->GetName(), EMCAL->GetGeometry()->GetTitle() );
-      Int_t i;
-      for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
-       AliEMCALHit * hit = (AliEMCALHit*)fHits->At(i) ;
-        AliEMCALDigit curSDigit = AliEMCALDigit(1,1,1,1);
-        AliEMCALDigit *sdigit ;
-       // Assign primary number only if contribution is significant
-     
-        if( hit->GetEnergy() > fPrimThreshold)
-          curSDigit =  AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
-        else
-          curSDigit =  AliEMCALDigit( -1               , -1               ,(((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
-       
-        for(Int_t check= 0; check < nSdigits; check++) {
-          sdigit = (AliEMCALDigit *)fSDigits->At(check);
-          if( (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()) == ((sdigit->GetId()/geom->GetNPhi())%geom->GetNZ()))  && ((hit->GetId()%geom->GetNPhi()) == (sdigit->GetId()%geom->GetNPhi()))) 
-           { 
-             *sdigit = *sdigit + curSDigit ;
-           }
-     else 
-         { new((*fSDigits)[nSdigits])  AliEMCALDigit(curSDigit);
-         nSdigits++ ; } 
-       }
-  
 
-        if( hit->GetEnergy() > fPrimThreshold)
-          curSDigit =  AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), ((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
-        else
-          curSDigit =  AliEMCALDigit( -1               , -1               ,((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
-      if((hit->GetId()/geom->GetNPhi()) < (2*geom->GetNZ())) 
-       {
-        for(Int_t check= 0; check < nSdigits; check++) {
-          sdigit = (AliEMCALDigit *)fSDigits->At(check);
-          if( (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()) == ((sdigit->GetId()/geom->GetNPhi())%geom->GetNZ()))  && ((hit->GetId()%geom->GetNPhi()) == (sdigit->GetId()%geom->GetNPhi()))) 
-           { 
-             *sdigit = *sdigit + curSDigit ;
-           }
-     else 
-         { new((*fSDigits)[nSdigits])  AliEMCALDigit(curSDigit);
-         nSdigits++ ; } 
-       }
-      } 
-     } 
-    } // loop over tracks
-    
-    fSDigits->Sort() ;
-    
-    nSdigits = fSDigits->GetEntriesFast() ;
-    fSDigits->Expand(nSdigits) ;
-    
-    Int_t i ;
-    for (i = 0 ; i < nSdigits ; i++) { 
-      AliEMCALDigit * digit = (AliEMCALDigit *) fSDigits->At(i) ; 
-      digit->SetIndexInList(i) ;     
+  if (fLastEvent == -1) 
+    fLastEvent = rl->GetNumberOfEvents() - 1 ;
+  else {
+    fLastEvent = TMath::Min(fLastEvent, rl->GetNumberOfEvents()-1);
+  }
+  Int_t nEvents   = fLastEvent - fFirstEvent + 1;
+
+  Int_t ievent;
+  Float_t energy=0.; // de * fSampling - 23-nov-04
+  rl->LoadKinematics();
+  rl->LoadHits("EMCAL");
+
+  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+    rl->GetEvent(ievent);
+    TTree * treeS = emcalLoader->TreeS();
+    if ( !treeS ) { 
+      emcalLoader->MakeSDigitsContainer();
+      treeS = emcalLoader->TreeS();
     }
+    TClonesArray * hits = emcalLoader->Hits() ; 
+    TClonesArray * sdigits = emcalLoader->SDigits() ;
+    sdigits->Clear();
 
-    if(gAlice->TreeS() == 0)
-      gAlice->MakeTree("S") ;
-    
-    //check, if this branch already exits?
-    TBranch * sdigitsBranch = 0;
-    TBranch * sdigitizerBranch = 0;
-    
-    TObjArray * branches = gAlice->TreeS()->GetListOfBranches() ;
-    Int_t ibranch;
-    Bool_t emcalNotFound = kTRUE ;
-    Bool_t sdigitizerNotFound = kTRUE ;
-    
-    for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
+    Int_t nSdigits = 0 ;
+    Int_t i;
+
+       AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); 
+    for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
+      AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
+      AliEMCALDigit * curSDigit = 0 ;
+      AliEMCALDigit * sdigit = 0 ;
+      Bool_t newsdigit = kTRUE; 
       
-      if(emcalNotFound){
-       sdigitsBranch=(TBranch *) branches->At(ibranch) ;
-       if( (strcmp("EMCAL",sdigitsBranch->GetName())==0 ) &&
-           (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) )
-         emcalNotFound = kFALSE ;
+      // hit->GetId() - Absolute Id number EMCAL segment
+      if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId())
+       energy = hit->GetEnergy() * fSampling; // 23-nov-04
+       if(energy >  fECPrimThreshold )
+         // Assign primary number only if deposited energy is significant
+         curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
+                                         hit->GetIparent(), hit->GetId(), 
+                                         Digitize(energy), hit->GetTime(),kFALSE,
+                                         -1, energy ) ;
+         else
+           curSDigit =  new AliEMCALDigit( -1, 
+                                           -1,
+                                           hit->GetId(), 
+                                           Digitize(energy), hit->GetTime(),kFALSE,
+                                           -1, energy ) ;
+      } else {
+       Warning("Exec"," abs id %i is bad \n", hit->GetId());
+       newsdigit = kFALSE;
+       curSDigit = 0;
       }
-      if(sdigitizerNotFound){
-       sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
-       if( (strcmp(sdigitizerBranch->GetName(),"AliEMCALSDigitizer") == 0)&&
-           (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) )
-         sdigitizerNotFound = kFALSE ;
+      
+      if(curSDigit != 0){
+       for(Int_t check= 0; check < nSdigits ; check++) {
+         sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
+         
+         if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?              
+           *sdigit = *sdigit + *curSDigit;
+           newsdigit = kFALSE;
+         }
+       }
       }
-    }
-
-    if(!(sdigitizerNotFound && emcalNotFound)){
-      cout << "AliEMCALSdigitizer error:" << endl ;
-      cout << "Can not overwrite existing branches: do not write" << endl ;
-      return ;
-    }
+      if (newsdigit) {
+       new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
+       nSdigits++ ;  
+      }
+      delete curSDigit ; 
+    }  // loop over all hits (hit = deposited energy/entering particle)
+    sdigits->Sort() ;
     
-    //Make (if necessary) branches    
-    char * file =0;
-    if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
-      file = new char[strlen(gAlice->GetBaseFile())+20] ;
-      sprintf(file,"%s/EMCAL.SDigits.root",gAlice->GetBaseFile()) ;
+    nSdigits = sdigits->GetEntriesFast() ;
+    fSDigitsInRun += nSdigits ;  
+
+    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
+      AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
+      sdigit->SetIndexInList(i) ;
     }
     
-    TDirectory *cwd = gDirectory;
+    // Now write SDigits    
     
-    //First list of sdigits
     Int_t bufferSize = 32000 ;    
-    sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&fSDigits,bufferSize);
-    sdigitsBranch->SetTitle(fSDigitsTitle.Data());
-    if (file) {
-      sdigitsBranch->SetFile(file);
-      TIter next( sdigitsBranch->GetListOfBranches());
-      TBranch * subbr;
-      while ((subbr=(TBranch*)next())) {
-       subbr->SetFile(file);
-      }   
-      cwd->cd();
-    } 
-      
-    //second - SDigitizer
-    Int_t splitlevel = 0 ;
-    AliEMCALSDigitizer * sd = this ;
-    sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
-                                              &sd,bufferSize,splitlevel); 
-    sdigitizerBranch->SetTitle(fSDigitsTitle.Data());
-    if (file) {
-      sdigitizerBranch->SetFile(file);
-      TIter next( sdigitizerBranch->GetListOfBranches());
-      TBranch * subbr ;
-      while ((subbr=(TBranch*)next())) {
-       subbr->SetFile(file);
-      }   
-      cwd->cd();
-      delete file;
-    }
-
-    sdigitsBranch->Fill() ;
-    sdigitizerBranch->Fill() ;
-    gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
+    TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
+    if (sdigitsBranch)
+      sdigitsBranch->SetAddress(&sdigits);
+    else
+      treeS->Branch("EMCAL",&sdigits,bufferSize);
     
-    if(strstr(option,"deb"))
-      PrintSDigits(option) ;
+    treeS->Fill();
     
+    emcalLoader->WriteSDigits("OVERWRITE");
+    
+    //NEXT - SDigitizer
+    emcalLoader->WriteSDigitizer("OVERWRITE");  // why in event cycle ?
+    
+    if(strstr(option,"deb"))
+      PrintSDigits(option) ;  
   }
   
-  if(strstr(option,"tim")){
-    gBenchmark->Stop("EMCALSDigitizer");
-    cout << "AliEMCALSDigitizer:" << endl ;
-    cout << "   took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing " 
-        <<  gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
-    cout << endl ;
-  }
+  Unload();
   
+  emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
   
+  if(strstr(option,"tim")){
+    gBenchmark->Stop("EMCALSDigitizer"); 
+    printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n", 
+          gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; 
+  }
 }
+
 //__________________________________________________________________
-void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
-  //Seting title to branch SDigits 
-  if(!fSDigitsTitle.IsNull())
-    cout << "AliEMCALSdigitizer: changing SDigits file from " <<fSDigitsTitle.Data() << " to " << title << endl ;
-  fSDigitsTitle=title ;
+Float_t AliEMCALSDigitizer::Digitize(Float_t energy)const {
+  // Digitize the energy
+  //
+  //JLK 26-Jun-2008 EXPLANATION LONG OVERDUE:
+  //
+  //We have to digitize the SDigit energy so that it can be stored in
+  //AliEMCALDigit, which has only an ADC amplitude member and
+  //(rightly) no energy member.  This method converts the energy to an
+  //integer which can be re-converted back to an energy with the
+  //Calibrate(energy) method when it is time to create Digits from SDigits 
+  //
+  Double_t aSignal = fA + energy*fB;
+  if (TMath::Abs(aSignal)>2147483647.0) { 
+    //PH 2147483647 is the max. integer
+    //PH This apparently is a problem which needs investigation
+    AliWarning(Form("Too big or too small energy %f",aSignal));
+    aSignal = TMath::Sign((Double_t)2147483647,aSignal);
+  }
+
+  return (Float_t ) aSignal;
 }
+
 //__________________________________________________________________
-void AliEMCALSDigitizer::Print(Option_t* option)const{
-  cout << "------------------- "<< GetName() << " -------------" << endl ;
-  cout << "   Writing SDigitis to branch with title  " << fSDigitsTitle.Data() << endl ;
-  cout << "   with digitization parameters  A = " << fA << endl ;
-  cout << "                                 B = " << fB << endl ;
-  cout << "   Threshold for Primary assignment= " << fPrimThreshold << endl ; 
-  cout << "---------------------------------------------------"<<endl ;
-  
+Float_t AliEMCALSDigitizer::Calibrate(Float_t amp)const {
+  //
+  // Convert the amplitude back to energy in GeV
+  //
+  //JLK 26-Jun-2008 EXPLANATION LONG OVERDUE:
+  //
+  //We have to digitize the SDigit energy with the method Digitize() 
+  //so that it can be stored in AliEMCALDigit, which has only an ADC 
+  //amplitude member and (rightly) no energy member.  This method is
+  //just the reverse of Digitize(): it converts the stored amplitude 
+  //back to an energy value in GeV so that the SDigit energies can be 
+  //summed before adding noise and creating digits out of them
+  //
+  return (Float_t)(amp - fA)/fB;
+
 }
+
+//__________________________________________________________________
+void AliEMCALSDigitizer::Print1(Option_t * option)
+{
+  Print(); 
+  PrintSDigits(option);
+}
+
+//__________________________________________________________________
+void AliEMCALSDigitizer::Print(Option_t *option) const
+{ 
+  // Prints parameters of SDigitizer
+  printf("Print: \n------------------- %s ------------- option %s\n", GetName() , option) ; 
+  printf("   fInit                                 %i\n", int(fInit));
+  printf("   fFirstEvent                           %i\n", fFirstEvent);
+  printf("   fLastEvent                            %i\n", fLastEvent);
+  printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
+  printf("   with digitization parameters       A = %f\n", fA) ; 
+  printf("                                      B = %f\n", fB) ;
+  printf("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold)  ;
+  printf("   Sampling                             = %f\n", fSampling);
+  printf("---------------------------------------------------\n") ;
+}
+
 //__________________________________________________________________
-Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
-  if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
+Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
+{
+  // Equal operator.
+  // SDititizers are equal if their pedestal, slope and threshold are equal
+  if( (fA==sd.fA)&&(fB==sd.fB)&&
+      (fECPrimThreshold==sd.fECPrimThreshold))
     return kTRUE ;
   else
     return kFALSE ;
 }
+
 //__________________________________________________________________
-void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
+void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
+{
   //Prints list of digits produced at the current pass of AliEMCALDigitizer
+    
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  const TClonesArray * sdigits = rl->SDigits() ; 
   
-  cout << "AliEMCALSDigitizer: " << endl ;
-  cout << "       Number of entries in SDigits list  " << fSDigits->GetEntriesFast() << endl ;
-  cout << endl ;
-  
-  if(strstr(option,"all")){// print all digits
+  printf("\n") ;  
+  printf("event %i", rl->GetRunLoader()->GetEventNumber());
+  printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast()); 
+  if(strstr(option,"all")||strstr(option,"EMC")){
     
     //loop over digits
     AliEMCALDigit * digit;
-    cout << "SDigit Id " << " Amplitude " <<  " Index "  <<  " Nprim " << " Primaries list " <<  endl;    
-    Int_t index ;
-    for (index = 0 ; index < fSDigits->GetEntries() ; index++) {
-      digit = (AliEMCALDigit * )  fSDigits->At(index) ;
-      cout << setw(8)  <<  digit->GetId() << " "  <<   setw(3)  <<  digit->GetAmp() <<   "  "  
-          << setw(6)  <<  digit->GetIndexInList() << "  "   
-          << setw(5)  <<  digit->GetNprimary() <<"  ";
+    printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
+       Int_t   index = 0;
+    Float_t isum  = 0.;
+    char * tempo = new char[8192]; 
+    for (index = 0 ; index < sdigits->GetEntries() ; index++) {
+      digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
+      sprintf(tempo, "\n%6d  %8f    %6.5e %4d      %2d :",
+             digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
+      printf("%s",tempo);
+      isum += digit->GetAmplitude();
       
       Int_t iprimary;
-      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
-       cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "  ";
-      cout << endl;     
+      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
+                 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
+       printf("%s",tempo); 
+      }         
     }
-    
-  }
+    delete tempo ;
+    printf("\n** Sum %2.3f : %10.3f GeV/c **\n ", isum, Calibrate(isum));
+  } else printf("\n");
+}
+
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::Unload() const
+{
+  // Unload Hits and SDigits from the folder
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  rl->UnloadHits() ; 
+  rl->UnloadSDigits() ; 
+}
+
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::Browse(TBrowser* b)
+{
+  TTask::Browse(b);
 }