]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALSDigitizer.cxx
Adding the full covariance matrix for the ITS space-points
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
index 2c212a9d7cc6cea47b7d4c87cc2ef3f7420e08a8..676cd717f92fd712192451e1bc86a93800a989c8 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 "TBenchmark.h"
+#include "TH1.h"
+#include "TBrowser.h"
+#include "Riostream.h"
+//#include "TObjectTable.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 "AliEMCALHistoUtilities.h"
 
 ClassImp(AliEMCALSDigitizer)
-
            
 //____________________________________________________________________________ 
-  AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","") 
+  AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","") 
 {
   // ctor
-  fA = 0;
-  fB = 10000000. ;
-  fPrimThreshold = 0.01 ;
-  fNevents = 0 ;     
-  fSDigits = 0 ;
-  fHits = 0 ;
-  fIsInitialized = kFALSE ;
-
+  fFirstEvent = fLastEvent  = fControlHists = 0;  
+  fDefaultInit = kTRUE ; 
+  fHists = 0;
 }
 
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliEMCALSDigitizer","")
+AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, 
+                                      const char * eventFolderName):
+  TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
+  fEventFolderName(eventFolderName)
 {
   // ctor
-  fA = 0;
-  fB = 10000000.;
-  fPrimThreshold = 0.01 ;
-  fNevents = 0 ;      
-  fSDigitsTitle = sDigitsTitle ;
-  fHeadersFile = headerFile ;
-  fSDigits = new TClonesArray("AliEMCALDigit",30000);
-  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 ;
+  fFirstEvent = fLastEvent  = fControlHists = 0 ; // runs one event by defaut  
+  Init();
+  InitParameters() ; 
+  fDefaultInit = kFALSE ; 
+  fHists = 0;
+  fControlHists = 1;
+  if(fControlHists) BookControlHists(1);
 }
 
+
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::~AliEMCALSDigitizer()
-{
-  // dtor
-  if(fSDigits)
-    delete fSDigits ;
-  if(fHits)
-    delete fHits ;
+AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) {
+  //cpy ctor 
+
+  fFirstEvent    = sd.fFirstEvent ; 
+  fLastEvent     = sd.fLastEvent ;
+  fA             = sd.fA ;
+  fB             = sd.fB ;
+  fECPrimThreshold  = sd.fECPrimThreshold ;
+  fSDigitsInRun  = sd.fSDigitsInRun ;
+  SetName(sd.GetName()) ; 
+  SetTitle(sd.GetTitle()) ; 
+  fEventFolderName = sd.fEventFolderName;
 }
+
+
+//____________________________________________________________________________ 
+AliEMCALSDigitizer::~AliEMCALSDigitizer() {
+  AliLoader *emcalLoader=0;
+  if ((emcalLoader = AliRunLoader::GetRunLoader()->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::GetRunLoader()->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",30000);
-    
-    // add Task to //root/Tasks folder
-    TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
-    roottasks->Add(this) ; 
-    
-    fIsInitialized = kTRUE ;
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::InitParameters()
+{
+  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+  if (geom->GetSampling() == 0.) {
+    Fatal("InitParameters", "Sampling factor not set !") ; 
   }
+  // Initializes parameters
+  fA         = 0;
+  fB         = 1.e+7;
+  fSampling  = geom->GetSampling();
+
+ // threshold for deposit energy of hit
+  fECPrimThreshold  = 0.; // 24-nov-04 - was 1.e-6;
+  Print("");
 }
+
 //____________________________________________________________________________
-void AliEMCALSDigitizer::Exec(Option_t *option) { 
-  //Collects all hits in the same active volume into digit
+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() ; 
+    return ; 
+  }
   
-  if(!fIsInitialized)
-    Init() ;
-
   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 ;
+
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  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 ;
     }
-    
-    //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 ;
+
+  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();
     }
-    
-    fSDigits->Clear();
+    TClonesArray * hits = emcalLoader->Hits() ; 
+    TClonesArray * sdigits = emcalLoader->SDigits() ;
+    sdigits->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);
-      AliEMCALv0 * EMCAL = (AliEMCALv0 *) gAlice->GetDetector("EMCAL") ;
-      AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance( EMCAL->GetGeometry()->GetName(), EMCAL->GetGeometry()->GetTitle() );
+    //Now make SDigits from hits, for EMCAL it is the same, so just copy    
+    Int_t nPrim =  static_cast<Int_t>((emcalLoader->TreeH())->GetEntries()) ; 
+    // This is not true: there is only one list of hits (MvL jan 2006)
+    // Attention nPrim is the number of primaries tracked by Geant 
+    // and this number could be different to the number of Primaries in TreeK; 
+    Int_t iprim ;
+    for ( iprim = 0 ; iprim < nPrim ; iprim++ ) { 
+      //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
+      rl->Stack()->Particle(iprim) ;
       Int_t i;
-      for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
-       AliEMCALHit * hit = (AliEMCALHit*)fHits->At(i) ;
-        AliEMCALDigit  * curSDigit = 0  ;
-        AliEMCALDigit * sdigit ;
-        Bool_t newsdigit = kTRUE; 
-// Assign primary number only if contribution is significant
-        if( hit->GetEnergy() > fPrimThreshold)
-   curSDigit =  new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1 ) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
-        else
-   curSDigit =  new AliEMCALDigit( -1               , -1               ,(((hit->GetId()/geom->GetNPhi())%geom->GetNZ() + 1 ) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
-     //cout << "Hit ID = " <<hit->GetId() << endl ; 
-     //cout << "ID for detector = " << curSDigit->GetId() << endl ;  
-         //  cout << hit->GetEnergy() << " - hit energy   -   Digit Energy - " << curSDigit->GetAmp() << endl;
-      for(Int_t check= 0; check < nSdigits ; check++) {
-          sdigit = (AliEMCALDigit *)fSDigits->At(check);
-          if( sdigit->GetId() == curSDigit->GetId())   
-            {// cout << "SDigit - Get Amp  " << sdigit->GetAmp() << endl ; 
-             *sdigit = *sdigit + *curSDigit ;
-              newsdigit = kFALSE;
-          //  cout << " and after addition " << sdigit->GetAmp() << endl ; 
-           }
-          }
-     if (newsdigit) 
-         { new((*fSDigits)[nSdigits])  AliEMCALDigit(*curSDigit);
-         nSdigits++ ;  
-       //cout << "Detector nsdigits = " << nSdigits << endl ;
-             }
-        newsdigit = kTRUE;  
-         
-
-        if( hit->GetEnergy() > fPrimThreshold)
-      curSDigit =  new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), ((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ() + 1) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
-        else
-    curSDigit =  new AliEMCALDigit( -1               , -1               ,((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
-      if((hit->GetId()/geom->GetNPhi()) < (2*geom->GetNZ())) 
-       { //cout << "ID for Preshower = " << curSDigit->GetId()  << endl ;
-        for(Int_t check= 0; check < nSdigits; check++) {
-          sdigit = (AliEMCALDigit *)fSDigits->At(check);
-          if( sdigit->GetId() == curSDigit->GetId())   
-           { 
-             *sdigit = *sdigit + *curSDigit ;
-             newsdigit = kFALSE ;
-            }
-           }
-     if (newsdigit) 
-         { new((*fSDigits)[nSdigits])  AliEMCALDigit(*curSDigit);
+      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; 
+
+        // 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() ) ;
+         else
+           curSDigit =  new AliEMCALDigit( -1               , 
+                                           -1               ,
+                                           hit->GetId(), 
+                                           Digitize(energy), hit->GetTime() ) ;
+       } else {
+          Warning("Exec"," abs id %i is bad \n", hit->GetId());
+          newsdigit = kFALSE;
+         curSDigit = 0;
+        }
+
+       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 (newsdigit) {
+         new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
          nSdigits++ ;  
-       //cout << "Preshower nsdigits = " << nSdigits << endl ;
-           }
-        newsdigit=kTRUE;       
-      } 
-     } 
-    } // loop over tracks
-    fSDigits->Sort() ;
+       }
+       delete curSDigit ; 
+      }  // loop over all hits (hit = deposited energy/entering particle)
+    } // loop over iprim
+    sdigits->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) ;     
-    }
+    nSdigits = sdigits->GetEntriesFast() ;
+    fSDigitsInRun += nSdigits ;  
 
-    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++){
-      
-      if(emcalNotFound){
-       sdigitsBranch=(TBranch *) branches->At(ibranch) ;
-       if( (strcmp("EMCAL",sdigitsBranch->GetName())==0 ) &&
-           (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) )
-         emcalNotFound = kFALSE ;
-      }
-      if(sdigitizerNotFound){
-       sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
-       if( (strcmp(sdigitizerBranch->GetName(),"AliEMCALSDigitizer") == 0)&&
-           (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) )
-         sdigitizerNotFound = kFALSE ;
-      }
-    }
+    Int_t nPrimarymax = -1 ; 
+    Int_t i ;
+    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) ;
 
-    if(!(sdigitizerNotFound && emcalNotFound)){
-      cout << "AliEMCALSdigitizer error:" << endl ;
-      cout << "Can not overwrite existing branches: do not write" << endl ;
-      return ;
+      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);
     
-    //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()) ;
+    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { // for what 
+      if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > nPrimarymax)
+       nPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
     }
+
+    // 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();
     
-    TDirectory *cwd = gDirectory;
+    emcalLoader->WriteSDigits("OVERWRITE");
     
-    //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) ;
+    //NEXT - SDigitizer
+    emcalLoader->WriteSDigitizer("OVERWRITE");  // why in event cycle ?
     
     if(strstr(option,"deb"))
-      PrintSDigits(option) ;
-    
+      PrintSDigits(option) ;  
   }
+
+  Unload();
+  
+  emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
   
   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 ;
+    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 ;
+
+
+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  " << fSDigitsTitle.Data() << endl ;
-  cout << "   with digitization parameters  A = " << fA << endl ;
-  cout << "                                 B = " << fB << endl ;
-  cout << "   Threshold for Primary assignment= " << fPrimThreshold << 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{
-  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
   
-  cout << "AliEMCALSDigitizer: " << endl ;
-  cout << "       Number of entries in SDigits list  " << fSDigits->GetEntriesFast() << endl ;
-  cout << endl ;
   
-  if(strstr(option,"all")){// print all digits
+  // 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 " <<  " 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, 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  %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 << setw(5)  <<  digit->GetPrimary(iprimary+1) << "  ";
-      cout << endl;     
+      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
+       sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
+       printf(tempo); 
+      }         
     }
-    
+    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)
+{ // 22-nov-04
+  gROOT->cd();
+  const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance() ;
+  if(var>=1){
+    printf(" AliEMCALSDigitizer::BookControlHists() in action \n");
+    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;
+}
+
+void AliEMCALSDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
+{
+  AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt); 
 }