]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALSDigitizer.cxx
Load pythia libraries.
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
index 052545f52f189031687d59af9e55fc4dfb545fca..476610b447e30ede56ceb81cad2140837f74cc31 100644 (file)
 //
 //*-- Author : Sahal Yacoob (LBL)
 // based on  : AliPHOSSDigitzer 
+// Modif: 
+//  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
+//                           of new  IO (à la PHOS)
 //////////////////////////////////////////////////////////////////////////////
 
-
 // --- ROOT system ---
-#include "TFile.h"
-#include "TTask.h"
-#include "TTree.h"
-#include "TSystem.h"
-#include "TROOT.h"
-#include "TFolder.h"
-#include "TGeometry.h"
-#include "TBenchmark.h"
+#include <TBenchmark.h>
+#include <TH1.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 "AliHeader.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 "AliEMCALGetter.h"
+#include "AliEMCALHistoUtilities.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.),
+    fControlHists(0),
+    fHists(0)
 {
   // ctor
-  fA = fB =  fNevents = 0 ; 
-  fTowerPrimThreshold = fPreShowerPrimThreshold = fPhotonElectronFactor  = 0. ;
-  fHits = fSDigits = fSDigits = 0 ;
-  fSplitFile = 0 ; 
-  fIsInitialized = kFALSE ;
-
-
+  InitParameters();
 }
 
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask(sDigitsTitle, headerFile)
+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.),
+    fControlHists(1),
+    fHists(0)
 {
   // ctor
-  fA = 0;
-  fB = 10000000.;
-  fTowerPrimThreshold = 0.01 ;
-  fPreShowerPrimThreshold = 0.0001 ; 
-  fNevents = 0 ; 
-  fPhotonElectronFactor = 5000. ; // photoelectrons per GeV 
-  fSplitFile = 0 ; 
   Init();
+  InitParameters() ; 
+  if(fControlHists) BookControlHists(1);
 }
 
+
 //____________________________________________________________________________ 
-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),
+    fControlHists(sd.fControlHists),
+    fHists(sd.fHists)
 {
-  // dtor
+  //cpy ctor 
+}
+
 
-  if (fSplitFile) 
-    if ( fSplitFile->IsOpen() ) 
-      fSplitFile->Close() ; 
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  gime->CloseFile() ; 
+//____________________________________________________________________________ 
+AliEMCALSDigitizer::~AliEMCALSDigitizer() {
+  //dtor
+  AliLoader *emcalLoader=0;
+  if ((emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")))
+      emcalLoader->CleanSDigitizer();
 }
 
 //____________________________________________________________________________ 
 void AliEMCALSDigitizer::Init(){
-
   // Initialization: open root-file, allocate arrays for hits and sdigits,
   // attach task SDigitizer to the list of EMCAL tasks
   // 
@@ -120,397 +149,332 @@ void AliEMCALSDigitizer::Init(){
   //============================================================= YS
   //  The initialisation is now done by the getter
 
-  if( strcmp(GetTitle(), "") == 0 )
-    SetTitle("galice.root") ;
-    
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), "update") ;     
-  if ( gime == 0 ) {
-    cerr << "ERROR: AliEMCALSDigitizer::Init -> Could not obtain the Getter object !" << endl ; 
+  fInit = kTRUE; 
+   
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+
+  if ( emcalLoader == 0 ) {
+    Fatal("Init", "Could not obtain the AliEMCALLoader");
     return ;
   } 
   
-  gime->PostSDigits( GetName(), GetTitle() ) ; 
+  emcalLoader->PostSDigitizer(this);
+  emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
+  
+}
 
-  TString sdname(GetName() );
-  sdname.Append(":") ;
-  sdname.Append(GetTitle() ) ;
-  SetName(sdname.Data()) ;
-  gime->PostSDigitizer(this) ;
- }
-//____________________________________________________________________________
-void AliEMCALSDigitizer::Exec(Option_t *option) { 
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::InitParameters()
+{
+  //initialize parameters for sdigitization
 
-// Collects all hits in the same active volume into digit
-  if( strcmp(GetName(), "") == 0 )
-    Init() ;
+  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+  if (geom->GetSampling() == 0.) {
+    Fatal("InitParameters", "Sampling factor not set !") ; 
+  }
+  // Initializes parameters
+  fA         = 0;
+  fB         = 1.e+6;  // Changed 24 Apr 2007. Dynamic range now 2 TeV
+  fSampling  = geom->GetSampling();
+
+  // threshold for deposit energy of hit
+  fECPrimThreshold  = 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"));
+
+  //  Print("");
+
+}
+
+//____________________________________________________________________________
+void AliEMCALSDigitizer::Exec(Option_t *option) 
+{ 
+  // Collects all hit of the same tower into digits
+  TString o(option); o.ToUpper();
   if (strstr(option, "print") ) {
-    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"));
+
+    //    Print();
     return ; 
   }
   
   if(strstr(option,"tim"))
     gBenchmark->Start("EMCALSDigitizer");
 
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 
-  //Check, if this branch already exits
-  gAlice->GetEvent(0) ;
-  TString sdname(GetName()) ;
-  sdname.Remove(sdname.Index(GetTitle())-1) ;
-
-  if(gAlice->TreeS() ) {
-    TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
-    TIter next(lob) ; 
-    TBranch * branch = 0 ;  
-    Bool_t emcalfound = kFALSE, sdigitizerfound = kFALSE ; 
-    
-     while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
-      TString thisName( GetName() ) ; 
-      TString branchName( branch->GetTitle() ) ; 
-      branchName.Append(":") ; 
-      if ( (strcmp(branch->GetName(), "EMCAL")==0) && thisName.BeginsWith(branchName) )  
-       emcalfound = kTRUE ;
-      
-      else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && thisName.BeginsWith(branchName) )  
-       sdigitizerfound = kTRUE ; 
+  //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 ;
     }
-    
-    if ( emcalfound || sdigitizerfound ) {
-      cerr << "WARNING: AliEMCALSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName() 
-          << " already exits" << endl ;
-      return ; 
-    }   
-  }  
-
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;  
-  Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
-  
-  Int_t ievent ;
-  for(ievent = 0; ievent < nevents; ievent++){
-    gime->Event(ievent,"H") ;
-    const TClonesArray * fHits = gime->Hits() ;
-    TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
+
+  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();
-    Int_t nSdigits = 0 ;
 
-  //Collects all hits in the same active volume into digit
-  
-  
-    
-    
-    //Now made SDigits from hits, for EMCAL it is the same, so just copy    
-    //  Int_t itrack ;
-    // for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
-    //gime->Track(itrack);  
-      //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
-     
-      Int_t i;
-      for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { // loop over all hits (hit = deposited energy/layer/entering particle)
-       AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(fHits->At(i)) ;
-        AliEMCALDigit * curSDigit = 0 ;
-        AliEMCALDigit * sdigit = 0 ;
-        Bool_t newsdigit = kTRUE; 
-
-       
-
-       // Assign primary number only if deposited energy is significant
-         
-       if( (!hit->IsInPreShower() && hit->GetEnergy() > fTowerPrimThreshold) || 
-           (hit->IsInPreShower() && hit->GetEnergy() > fPreShowerPrimThreshold)) {
-         curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
-                                         hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()), 
-                                         Digitize(hit->GetEnergy()), 
-                                         hit->GetTime()) ;
-       } else {
-         curSDigit =  new AliEMCALDigit( -1               , 
-                                         -1               ,
-                                         Layer2TowerID(hit->GetId(),hit->IsInPreShower()), 
-                                         Digitize(hit->GetEnergy()), 
-                                         hit->GetTime() ) ;    
-       }
-       Int_t check = 0 ;
-       for(check= 0; check < nSdigits ; check++) {
-          sdigit = (AliEMCALDigit *)sdigits->At(check);
-          if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same tower or the same preshower ?              
-             *sdigit = *sdigit + *curSDigit;
-   
-             newsdigit = kFALSE;
-         }
-       }
-       if (newsdigit) { 
-         new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
-         nSdigits++ ;  
-       }
-      }  // loop over all hits (hit = deposited energy/layer/entering particle)
-      //    } // loop over tracks
+    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; 
       
-      sdigits->Sort() ;
+      // 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(), 
+                                         -1, energy ) ;
+         else
+           curSDigit =  new AliEMCALDigit( -1               , 
+                                           -1               ,
+                                           hit->GetId(), 
+                                           Digitize(energy), hit->GetTime(), 
+                                           -1, energy ) ;
+      } else {
+       Warning("Exec"," abs id %i is bad \n", hit->GetId());
+       newsdigit = kFALSE;
+       curSDigit = 0;
+      }
       
-      nSdigits = sdigits->GetEntriesFast() ;
-      if (nSdigits > 0) { 
-       sdigits->Expand(nSdigits) ;
-       
-       //  Int_t i ;
-       const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
-       
-       Int_t lastPreShowerIndex = nSdigits - 1 ;
-       if (!(dynamic_cast<AliEMCALDigit *>(sdigits->At(lastPreShowerIndex))->IsInPreShower()))
-         lastPreShowerIndex = -2; 
-       Int_t firstPreShowerIndex = 100000 ; 
-       Int_t index ; 
-       AliEMCALDigit * sdigit = 0 ;
-       for ( index = 0; index < nSdigits ; index++) {
-         sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ) ;
-         if (sdigit->IsInPreShower() ){ 
-           firstPreShowerIndex = index ;
-           break ;
-         }
-       }
-       
-       AliEMCALDigit * preshower ;
-       AliEMCALDigit * tower ;
-       Int_t lastIndex = lastPreShowerIndex +1 ; 
-       
-       
-       for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) {
-         preshower = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ); 
-         Bool_t towerFound = kFALSE ;
-         Int_t jndex ; 
-         for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
-           tower  = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) ); 
-           if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {          
-             Float_t towerEnergy  = static_cast<Float_t>(tower->GetAmp()) ; 
-             Float_t preshoEnergy = static_cast<Float_t>(preshower->GetAmp()) ; 
-             towerEnergy +=preshoEnergy ; 
-             *tower = *tower + *preshower    ; // and add preshower multiplied by layer ratio to tower
-             tower->SetAmp(static_cast<Int_t>(TMath::Ceil(towerEnergy))) ; 
-             towerFound = kTRUE ;
-           }
-         }
-         if ( !towerFound ) { 
-           
-           new((*sdigits)[lastIndex])  AliEMCALDigit(*preshower);
-           AliEMCALDigit * temp = dynamic_cast<AliEMCALDigit *>(sdigits->At(lastIndex)) ;
-           temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ;
-           lastIndex++ ; 
+      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;
          }
        }
-       
-       sdigits->Sort() ;
-       Int_t NPrimarymax = -1 ; 
-       for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
-         sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
-         sdigit->SetIndexInList(i) ;
-       }
-       
-       for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
-         if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
-           NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
-       }
       }
-       if(gAlice->TreeS() == 0)
-         gAlice->MakeTree("S",fSplitFile);
-
-      //First list of sdigits
-      Int_t bufferSize = 32000 ;    
-      TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
-      sdigitsBranch->SetTitle(sdname);
-      
-      //second - SDigitizer
-      Int_t splitlevel = 0 ;
-      AliEMCALSDigitizer * sd = this ;
-      TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
-                                                          &sd,bufferSize,splitlevel); 
-      sdigitizerBranch->SetTitle(sdname);
-      
-      sdigitsBranch->Fill() ; 
-      sdigitizerBranch->Fill() ; 
-      gAlice->TreeS()->AutoSave() ;
-           
-      if(strstr(option,"deb"))
-       PrintSDigits(option) ;
+      if (newsdigit) {
+       new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
+       nSdigits++ ;  
+      }
+      delete curSDigit ; 
+    }  // loop over all hits (hit = deposited energy/entering particle)
+    sdigits->Sort() ;
+    
+    nSdigits = sdigits->GetEntriesFast() ;
+    fSDigitsInRun += nSdigits ;  
+    
+    Double_t e=0.,esum=0.;
+    AliEMCALHistoUtilities::FillH1(fHists, 0, double(sdigits->GetEntriesFast()));
+    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
+      AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
+      sdigit->SetIndexInList(i) ;
       
+      AliEMCALHistoUtilities::FillH1(fHists, 2, double(sdigit->GetAmp()));
+      e = double(Calibrate(sdigit->GetAmp()));
+      esum += e;
+      AliEMCALHistoUtilities::FillH1(fHists, 3, e);
+      AliEMCALHistoUtilities::FillH1(fHists, 4, double(sdigit->GetId()));
+    }
+    if(esum>0.) AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
+    
+    // Now write SDigits    
+    
+    Int_t bufferSize = 32000 ;    
+    TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
+    if (sdigitsBranch)
+      sdigitsBranch->SetAddress(&sdigits);
+    else
+      treeS->Branch("EMCAL",&sdigits,bufferSize);
+    
+    treeS->Fill();
+    
+    emcalLoader->WriteSDigits("OVERWRITE");
+    
+    //NEXT - SDigitizer
+    emcalLoader->WriteSDigitizer("OVERWRITE");  // why in event cycle ?
+    
+    if(strstr(option,"deb"))
+      PrintSDigits(option) ;  
   }
   
-  if (fSplitFile) 
-    if ( fSplitFile->IsOpen() ) 
-      fSplitFile->Close() ; 
-  
-  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);
   
-}
-//__________________________________________________________________
-void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
-  // Setting title to branch SDigits 
-
-  TString stitle(title) ;
-
-  // check if branch with title already exists
-  TBranch * sdigitsBranch    = 
-    static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ; 
-  TBranch * sdigitizerBranch =  
-    static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
-  const char * sdigitsTitle    = sdigitsBranch ->GetTitle() ;  
-  const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
-  if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
-    cerr << "ERROR: AliEMCALSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
-    return ;
+  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 ) ; 
   }
-  
-  cout << "AliEMCALSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
-
-  SetName(title) ; 
-    
-  // Post to the WhiteBoard
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  gime->PostSDigits( title, GetTitle()) ; 
-
-
 }
 
 //__________________________________________________________________
-void AliEMCALSDigitizer::SetSplitFile(const TString splitFileName) 
-{
-  // Diverts the SDigits in a file separate from the hits file
-  
-  TDirectory * cwd = gDirectory ;
 
-  if ( !(gAlice->GetTreeSFileName() == splitFileName) ) {
-    if (gAlice->GetTreeSFile() )  
-      gAlice->GetTreeSFile()->Close() ; 
-  }
-
-  fSplitFile = gAlice->InitTreeFile("S",splitFileName.Data());
-  fSplitFile->cd() ; 
-  if ( !fSplitFile->Get("gAlice") ) 
-    gAlice->Write();
-  
-  TTree *treeE  = gAlice->TreeE();
-  if (!treeE) {
-    cerr << "ERROR: AliEMCALSDigitizer::SetSPlitFile -> No TreeE found "<<endl;
-   abort() ;
-  }      
-  
-  // copy TreeE
-  if ( !fSplitFile->Get("TreeE") ) {
-    AliHeader *header = new AliHeader();
-    treeE->SetBranchAddress("Header", &header);
-    treeE->SetBranchStatus("*",1);
-    TTree *treeENew =  treeE->CloneTree();
-    treeENew->Write();
-  }
-  
-  // copy AliceGeom
-  if ( !fSplitFile->Get("AliceGeom") ) {
-    TGeometry *AliceGeom = static_cast<TGeometry*>(cwd->Get("AliceGeom"));
-    if (!AliceGeom) {
-      cerr << "ERROR: AliEMCALSDigitizer::SetSPlitFile -> AliceGeom was not found in the input file "<<endl;
-      abort() ;
+Int_t AliEMCALSDigitizer::Digitize(Float_t energy)const {
+  // Digitize the energy
+    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);
     }
-    AliceGeom->Write();
+    return (Int_t ) aSignal;
   }
-  
-  gAlice->MakeTree("S",fSplitFile);
-  cwd->cd() ; 
-  cout << "INFO: AliEMCALSDigitizer::SetSPlitFile -> SDigits will be stored in " << splitFileName.Data() << endl ; 
+
+//__________________________________________________________________
+
+void AliEMCALSDigitizer::Print1(Option_t * option)
+{
+  Print(); 
+  PrintSDigits(option);
 }
 
 //__________________________________________________________________
-void AliEMCALSDigitizer::Print(Option_t* option)const{
-  cout << "------------------- "<< GetName() << " -------------" << endl ;
-  cout << "   Writing SDigitis to branch with title  " << GetName() << endl ;
-  cout << "   with digitization parameters  A               = " << fA << endl ;
-  cout << "                                 B               = " << fB << endl ;
-  cout << "   Threshold for Primary assignment in Tower     = " << fTowerPrimThreshold << endl ; 
-  cout << "   Threshold for Primary assignment in PreShower = " << fPreShowerPrimThreshold << endl ; 
-  cout << "---------------------------------------------------"<<endl ;
-  
+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{
+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)&&
-      (fTowerPrimThreshold==sd.fTowerPrimThreshold) &&
-      (fPreShowerPrimThreshold==sd.fPreShowerPrimThreshold))
+      (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
   
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  TString sdname(GetName()) ;
-  sdname.Remove(sdname.Index(GetTitle())-1) ;
-  const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ; 
   
-  cout << "AliEMCALSDigitiser: event " << gAlice->GetEvNumber() << endl ;
-  cout << "      Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
-  cout << endl ;
+  // AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+  const TClonesArray * sdigits = rl->SDigits() ; 
+  
+  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 " <<  "     Time " << "     Index "  <<  " Nprim " << " Primaries list " <<  endl;    
-    Int_t index ;
+    printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
+    Int_t index, isum=0;
+    char * tempo = new char[8192]; 
     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
-      digit = (AliEMCALDigit * )  sdigits->At(index) ;
-      cout << setw(6)  <<  digit->GetId() << "   "  <<         setw(10)  <<  digit->GetAmp() <<   "    "  << digit->GetTime()
-          << setw(6)  <<  digit->GetIndexInList() << "    "   
-          << setw(5)  <<  digit->GetNprimary() <<"    ";
+      digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
+      sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
+             digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
+      printf(tempo);
+      isum += digit->GetAmp();
       
       Int_t iprimary;
-      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
-       cout << " "  <<  digit->GetPrimary(iprimary+1) << "  ";
-      cout << endl;     
+      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
+       sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
+       printf(tempo); 
+      }         
     }
-    cout <<endl;
+    delete tempo ;
+    printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6);
+  } else printf("\n");
+}
+
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::Unload() const
+{
+  // Unload Hits and SDigits from the folder
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+  rl->UnloadHits() ; 
+  rl->UnloadSDigits() ; 
+}
+
+void AliEMCALSDigitizer::Browse(TBrowser* b)
+{
+  if(fHists) b->Add(fHists);
+  TTask::Browse(b);
+}
+
+TList *AliEMCALSDigitizer::BookControlHists(int var)
+{ 
+  //book histograms for monitoring sdigitization
+  // 22-nov-04
+  gROOT->cd();
+  const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance() ;
+  if(var>=1){
+    AliDebug(1, " BookControlHists() in action ");
+    new TH1F("HSDigiN",  "#EMCAL  sdigits ", 1001, -0.5, 1000.5);
+    new TH1F("HSDigiSumEnergy","Sum.EMCAL energy", 1000, 0.0, 100.);
+    new TH1F("HSDigiAmp",  "EMCAL sdigits amplitude", 1000, 0., 2.e+9);
+    new TH1F("HSDigiEnergy","EMCAL cell energy", 1000, 0.0, 100.);
+    new TH1F("HSDigiAbsId","EMCAL absID for sdigits",
+    geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
   }
+
+  fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalSDigiControlHists", kFALSE);
+  //  fHists = 0; ??
+
+  return fHists;
 }
-//________________________________________________________________________
-Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower){
-  // Method to Transform from Hit Id to Digit Id
-  // This function should be one to one
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
-  const AliEMCALGeometry * geom = gime->EMCALGeometry();
-  Int_t ieta  = ((ihit-1)/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index
-  Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index
-  Int_t it = -10;
-  Int_t ipre = 0;
-
-  if (preshower)ipre = 1;
-  if (iphi > 0 && ieta >= 0){
-    it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ();
-    return it;
-  }else{
-    cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Eta number = "
-        << ieta << "Phi number = " << iphi << endl ;
-    return it;
-  } // end if iphi>0 && ieta>0
+
+void AliEMCALSDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
+{
+  AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt); 
 }
-//_______________________________________________________________________________________
-// void AliEMCALSDigitizer::TestTowerID(void)
-// {
-//   Int_t j;
-
-//   Bool_t preshower = kFALSE;
-//   for (j = 0 ; j < 10 ; j++){  // loop over hit id
-//     Int_t i;
-//    for (i = 0 ; i <= 2 ; i++){  // loop over 
-//      Int_t k = i*96*144+j*144+1;
-//       cout << " Hit Index = " << k << "   " << j*10 << "   TOWERID = " <<  Layer2TowerID(k, preshower) << endl ;
-//     }
-//   }
-// }