]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALSDigitizer.cxx
Coverity fixes
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
index 9793f86ee53df91009bb4e4828377a8f035e4a5e..a483a5bae41d6f643ee9eaf666f02e01cb607e0c 100644 (file)
 /* $Id$ */
 
 //_________________________________________________________________________
-// This is a TTask that makes SDigits out of Hits
+// This is a Class that makes SDigits out of Hits
 // A Summable Digits is the sum of all hits originating 
 // 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".
 // User case:
 // root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root")
 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-// root [1] s->ExecuteTask()
+// root [1] s->Digitize()
 //             // Makes SDigitis for all events stored in galice.root
 // root [2] s->SetPedestalParameter(0.001)
 //             // One can change parameters of digitization
 // root [3] s->SetSDigitsBranch("Redestal 0.001")
 //             // and write them into the new branch
-// root [4] s->ExeciteTask("deb all tim")
+// root [4] s->Digitize("deb all tim")
 //             // available parameters:
 //             deb - print # of produced SDigitis
 //             deb all  - print # and list of produced SDigits
 //
 //*-- 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 "TBenchmark.h"
+#include <TBenchmark.h>
+#include <TBrowser.h>
+//#include <Riostream.h>
+#include <TMath.h>
+#include <TROOT.h>
 
 // --- Standard library ---
 #include "stdlib.h"
 
 // --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"
 #include "AliEMCALDigit.h"
-#include "AliEMCALGetter.h"
+#include "AliEMCALLoader.h"
 #include "AliEMCALHit.h"
 #include "AliEMCALSDigitizer.h"
 #include "AliEMCALGeometry.h"
+#include "AliEMCALSimParam.h"
 
 ClassImp(AliEMCALSDigitizer)
            
 //____________________________________________________________________________ 
-  AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","") 
+AliEMCALSDigitizer::AliEMCALSDigitizer()
+  : TNamed("",""),
+    fA(0.),fB(0.),fECPrimThreshold(0.),
+    fDefaultInit(kTRUE),
+    fEventFolderName(0),
+    fInit(0),
+    fSDigitsInRun(0),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fSampling(0.),
+       fHits(0)
 {
   // ctor
-  InitParameters() ; 
-  fDefaultInit = kTRUE ; 
+  InitParameters();
 }
 
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, const char * eventFolderName):
-  TTask("EMCAL"+AliConfig::fgkSDigitizerTaskName, alirunFileName),
-  fEventFolderName(eventFolderName)
+AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, 
+                                      const char * eventFolderName)
+  : TNamed("EMCALSDigitizer", alirunFileName),
+    fA(0.),fB(0.),fECPrimThreshold(0.),
+    fDefaultInit(kFALSE),
+    fEventFolderName(eventFolderName),
+    fInit(0),
+    fSDigitsInRun(0),
+    fFirstEvent(0),
+    fLastEvent(0),
+    fSampling(0.),
+    fHits(0)
 {
   // ctor
   Init();
   InitParameters() ; 
-  fDefaultInit = kFALSE ; 
+
 }
 
 
 //____________________________________________________________________________ 
-AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) {
+AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) 
+  : TNamed(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),
+    fHits(0)
+{
   //cpy ctor 
+}
+
+//_____________________________________________________________________
+AliEMCALSDigitizer& AliEMCALSDigitizer::operator = (const AliEMCALSDigitizer& source)
+{ // assignment operator; use copy ctor
+  if (&source == this) return *this;
 
-  fA             = sd.fA ;
-  fB             = sd.fB ;
-  fPREPrimThreshold = sd.fPREPrimThreshold ;
-  fECPrimThreshold  = sd.fECPrimThreshold ;
-  fHCPrimThreshold  = sd.fHCPrimThreshold ;
-  fSDigitsInRun  = sd.fSDigitsInRun ;
-  SetName(sd.GetName()) ; 
-  SetTitle(sd.GetTitle()) ; 
-  fEventFolderName = sd.fEventFolderName;
+  new (this) AliEMCALSDigitizer(source);
+  return *this;
 }
 
+//____________________________________________________________________________ 
+AliEMCALSDigitizer::~AliEMCALSDigitizer() {
+  //dtor
+       
+  if(fHits){
+         fHits->Clear();
+         delete fHits;
+  }
+}
 
 //____________________________________________________________________________ 
 void AliEMCALSDigitizer::Init(){
-  // Initialization: open root-file, allocate arrays for hits and sdigits,
-  // attach task SDigitizer to the list of EMCAL tasks
+  // Initialization: open root-file, allocate arrays for hits and sdigits
   // 
   // Initialization can not be done in the default constructor
   //============================================================= YS
   //  The initialisation is now done by the getter
 
-  fInit = kTRUE 
+  fInit = kTRUE; 
    
-  AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());  
-  if ( gime == 0 ) {
-    Fatal("Init", "Could not obtain the Getter objectfor file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+
+  if ( emcalLoader == 0 ) {
+    Fatal("Init", "Could not obtain the AliEMCALLoader");
     return ;
   } 
   
-  TString opt("SDigits") ; 
-  if(gime->VersionExists(opt) ) { 
-    Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
-    fInit = kFALSE ; 
-  }
-  
-  gime->PostSDigitizer(this);
-  gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
-  
 }
 
 //____________________________________________________________________________ 
 void AliEMCALSDigitizer::InitParameters()
 {
-  // Initializes parameters
-  fA                      = 0;
-  fB                      = 10000000.;
+  //initialize parameters for sdigitization
 
-  AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
-  const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
+  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
   if (geom->GetSampling() == 0.) {
-    Error("InitParameters", "Sampling factor not set !") ; 
-    abort() ;
+    Fatal("InitParameters", "Sampling factor not set !") ; 
   }
-  else
-    Info("InitParameters", "Sampling factor set to %f\n", geom->GetSampling()) ; 
-  
-  // this threshold corresponds approximately to 100 MeV
-  fECPrimThreshold     = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNECLayers() ;
-  fPREPrimThreshold    = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNPRLayers() ; 
-  fHCPrimThreshold     = fECPrimThreshold/5. ; // 5 is totally arbitrary
-
-}
 
-//____________________________________________________________________________
-void AliEMCALSDigitizer::Exec(Option_t *option) 
-{ 
-  // Collects all hits in the section (PRE/ECAL/HCAL) of the same tower into digit
-  
-  if (strstr(option, "print") ) {
-    Print() ; 
-    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?");
   }
-  
-  if(strstr(option,"tim"))
-    gBenchmark->Start("EMCALSDigitizer");
+       
+  //
+  //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();
 
-  AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
-  //switch off reloading of this task while getting event
-  if (!fInit) { // to prevent overwrite existing file
-    Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
-    return ;
-  }
+  // threshold for deposit energy of hit
+  fECPrimThreshold  = simParam->GetECPrimaryThreshold();//0.05;// GeV // 22-may-07 was 0// 24-nov-04 - was 1.e-6;
   
-  Int_t nevents = gime->MaxEvent() ; 
-  Int_t ievent ;   
-  for(ievent = 0; ievent < nevents; ievent++){     
-  
-    gime->Event(ievent,"H") ;  
-
-    TTree * treeS = gime->TreeS(); 
-    TClonesArray * hits = gime->Hits() ; 
-    TClonesArray * sdigits = gime->SDigits() ;
-    sdigits->Clear();
-    Int_t nSdigits = 0 ;
-    //Now make SDigits from hits, for EMCAL it is the same, so just copy    
-    Int_t nPrim =  static_cast<Int_t>((gime->TreeH())->GetEntries()) ; 
-    // 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
-      gime->Track(iprim) ;
-      Int_t i;
-      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; 
+  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"));
 
-       // Assign primary number only if deposited energy is significant
-
-       AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
+}
 
-       if( geom->IsInPRE(hit->GetId()) )  
-         if( hit->GetEnergy() > fPREPrimThreshold )
-           curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
-                                           hit->GetIparent(), hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
-         else
-           curSDigit =  new AliEMCALDigit( -1               , 
-                                           -1               ,
-                                           hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
-       else if( geom->IsInECA(hit->GetId()) )
-         if( hit->GetEnergy() >  fECPrimThreshold )
-           curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
-                                           hit->GetIparent(), hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
-         else
-           curSDigit =  new AliEMCALDigit( -1               , 
-                                           -1               ,
-                                           hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
-       else if( geom->IsInHCA(hit->GetId()) )
-         if( hit->GetEnergy() >  fHCPrimThreshold )
-           
-           curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
-                                           hit->GetIparent(), hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
-         else
-           curSDigit =  new AliEMCALDigit( -1               , 
-                                           -1               ,
-                                           hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
+//____________________________________________________________________________
+void AliEMCALSDigitizer::Digitize(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 ; 
+       }
+       
        
-       Int_t check = 0 ;
-       for(check= 0; check < nSdigits ; check++) {
-         sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
-         if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL/HCAL/preshower tower ?              
-           *sdigit = *sdigit + *curSDigit;
-           newsdigit = kFALSE;
-         }
+       if(strstr(option,"tim"))
+               gBenchmark->Start("EMCALSDigitizer");
+       
+       AliRunLoader *rl = AliRunLoader::Instance();
+       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
+       
+       if (!fInit) { // to prevent overwrite existing file
+         AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
+         return ;
        }
-       if (newsdigit) { 
-         new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
-         nSdigits++ ;  
+       
+       if (fLastEvent == -1) 
+               fLastEvent = rl->GetNumberOfEvents() - 1 ;
+       else {
+               fLastEvent = TMath::Min(fLastEvent, rl->GetNumberOfEvents()-1);
        }
-       delete curSDigit ; 
-      }  // loop over all hits (hit = deposited energy/layer/entering particle)
-    } // loop over iprim
-    
-    sdigits->Sort() ;
-    
-    nSdigits = sdigits->GetEntriesFast() ;
-    fSDigitsInRun += nSdigits ;  
-    sdigits->Expand(nSdigits) ;
-
-    Int_t nPrimarymax = -1 ; 
-    Int_t i ;
-    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
-      AliEMCALDigit * 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()) ;
-    }
-    
-    // Now write SDigits    
+       Int_t nEvents   = fLastEvent - fFirstEvent + 1;
+       
+       Float_t energy=0.; // de * fSampling - 23-nov-04
+       rl->LoadKinematics();
+       rl->LoadHits("EMCAL");
+  
+  Int_t ievent;
+       for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+               rl->GetEvent(ievent);
+               TTree * treeS = emcalLoader->TreeS();
+               if ( !treeS ) { 
+                       emcalLoader->MakeSDigitsContainer();
+                       treeS = emcalLoader->TreeS();
+               }
     
-    //First list of sdigits
+               TClonesArray * sdigits = emcalLoader->SDigits() ;
+               sdigits->Clear("C");
+               
+               Int_t nSdigits = 0 ;
+               Int_t iHit, iTrack, iSDigit;
+               
+               AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); 
+               
+               TTree *treeH = emcalLoader->TreeH();    
+               if (treeH) {
+                       Int_t nTrack = treeH->GetEntries();  // TreeH has array of hits for every primary
+                       TBranch * branchH = treeH->GetBranch("EMCAL");
+                       //if(fHits)fHits->Clear();
+                       branchH->SetAddress(&fHits);
+                       for (iTrack = 0; iTrack < nTrack; iTrack++) {
+                               branchH->GetEntry(iTrack); 
+        
+        if(fHits){
+          
+          Int_t nHit = fHits->GetEntriesFast();
+          for(iHit = 0; iHit< nHit;iHit++){
+            
+            AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(fHits->At(iHit)) ;
+            AliEMCALDigit * curSDigit = 0 ;
+            AliEMCALDigit * sdigit = 0 ;
+            Bool_t newsdigit = kTRUE; 
+            
+            // hit->GetId() - Absolute Id number EMCAL segment
+            if(hit){
+              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, 0,0,energy ) ;
+                else
+                  curSDigit =  new AliEMCALDigit(-1, 
+                                                 -1,
+                                                 hit->GetId(), 
+                                                 Digitize(energy), hit->GetTime(),kFALSE,
+                                                 -1, 0,0,energy ) ;
+              } else {
+                Warning("Digitize"," 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){
+                    if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?              
+                      *sdigit = *sdigit + *curSDigit;
+                      newsdigit = kFALSE;
+                    }
+                  }// sdigit exists
+                  else {
+                    AliWarning("Sdigit do not exist");
+                    newsdigit = kFALSE;  
+                  }// sdigit does not exist
+                }//sdigit loop
+              }// currsdigit exists
+              
+              if (newsdigit) {
+                new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
+                nSdigits++ ;  
+              }
+              delete curSDigit ;
+              
+            }// hit exists
+            else AliFatal("Hit is NULL!");
+            
+          }  // loop over all hits (hit = deposited energy/entering particle)
+          
+        }//fHits is not NULL
+        else AliFatal("fHit is NULL!");
+        
+                               sdigits->Sort() ;
+                               
+                               nSdigits = sdigits->GetEntriesFast() ;
+                               fSDigitsInRun += nSdigits ;  
+                               
+                               for (iSDigit = 0 ; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) { 
+                                       AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(iSDigit)) ;
+                                       if(sdigit)sdigit->SetIndexInList(iSDigit) ;
+          else AliFatal("sdigit is NULL!");
+                               }       
+                               if(fHits)fHits->Clear();
+                       }//track loop
+               }// tree exists
+               
+               // 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) ;  
+       }
+       
+       Unload();
+               
+       if(strstr(option,"tim")){
+               gBenchmark->Stop("EMCALSDigitizer"); 
+               printf("\n Digitize: took %f seconds for SDigitizing %f seconds per event\n", 
+           gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; 
+       }
+}
 
-    Int_t bufferSize = 32000 ;    
-    TBranch * sdigitsBranch = treeS->Branch("EMCAL",&sdigits,bufferSize);
-    sdigitsBranch->Fill() ;
+//__________________________________________________________________
+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);
+  }
 
-    gime->WriteSDigits("OVERWRITE");
-    
-    //NEXT - SDigitizer
+  return (Float_t ) aSignal;
+}
 
-    gime->WriteSDigitizer("OVERWRITE");
-    
-    if(strstr(option,"deb"))
-      PrintSDigits(option) ;  
-  }
+//__________________________________________________________________
+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;
 
-  Unload();
-  
-  gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
-  
-  if(strstr(option,"tim")){
-    gBenchmark->Stop("EMCALSDigitizer"); 
-    Info("Exec", "took %f seconds for SDigitizing %f seconds per event", 
-        gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ; 
-  }
 }
 
+//__________________________________________________________________
+void AliEMCALSDigitizer::Print1(Option_t * option)
+{
+  Print(); 
+  PrintSDigits(option);
+}
 
 //__________________________________________________________________
-void AliEMCALSDigitizer::Print()const
+void AliEMCALSDigitizer::Print(Option_t *option) const
 { 
   // Prints parameters of SDigitizer
-  Info("Print", "\n------------------- %s -------------", GetName() ) ; 
+  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 PRE Primary assignment= %f\n", fPREPrimThreshold)  ; 
-  printf("   Threshold for EC Primary assignment= %f\n", fECPrimThreshold)  ; 
-  printf("   Threshold for HC Primary assignment= %f\n", fHCPrimThreshold)  ; 
+  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") ;
-
 }
 
 //__________________________________________________________________
@@ -319,60 +465,70 @@ 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) &&
-      (fHCPrimThreshold==sd.fHCPrimThreshold) &&
-      (fPREPrimThreshold==sd.fPREPrimThreshold))
+      (fECPrimThreshold==sd.fECPrimThreshold))
     return kTRUE ;
   else
     return kFALSE ;
 }
-
 //__________________________________________________________________
 void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
 {
   //Prints list of digits produced at the current pass of AliEMCALDigitizer
   
-  
-  AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
-  const TClonesArray * sdigits = gime->SDigits() ; 
-  
-  TString message("\n") ;  
-  message += "event " ; 
-  message += gAlice->GetEvNumber() ;
-  message += "\n      Number of entries in SDigits list " ;
-  message += sdigits->GetEntriesFast() ; 
-  if(strstr(option,"all")||strstr(option,"EMC")){
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  if(rl){
+    const TClonesArray * sdigits = rl->SDigits() ; 
     
-    //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 < 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()) ;  
-      message += tempo ; 
+    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")){
       
-      Int_t iprimary;
-      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
-       sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
-       message += tempo ; 
-      }         
-    }
-    delete tempo ;
+      //loop over digits
+      AliEMCALDigit * digit;
+      printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
+      Int_t   index = 0;
+      Float_t isum  = 0.;
+      const Int_t bufferSize = 8192;
+      char * tempo = new char[bufferSize]; 
+      for (index = 0 ; index < sdigits->GetEntries() ; index++) {
+        digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
+        if(digit){
+          snprintf(tempo, bufferSize,"\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++) {
+            snprintf(tempo,bufferSize, "%d ",digit->GetPrimary(iprimary+1) ) ; 
+            printf("%s",tempo); 
+          }  
+        } //sdigit exists
+        else AliFatal("SDigit is NULL!");
+      }//loop
+      delete [] tempo ;
+      printf("\n** Sum %2.3f : %10.3f GeV/c **\n ", isum, Calibrate(isum));
+    } else printf("\n");
   }
-  Info("PrintSDigits", message.Data() ) ; 
+  else AliFatal("EMCALLoader is NULL!");
 }
 
 //____________________________________________________________________________ 
 void AliEMCALSDigitizer::Unload() const
 {
   // Unload Hits and SDigits from the folder
-  AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
-  AliEMCALLoader * loader = gime->EmcalLoader() ; 
-  loader->UnloadHits() ; 
-  loader->UnloadSDigits() ; 
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
+  if(rl){
+    rl->UnloadHits() ; 
+    rl->UnloadSDigits() ;
+  }
+  else AliFatal("EMCALLoader is NULL!");
+}
+
+//____________________________________________________________________________ 
+void AliEMCALSDigitizer::Browse(TBrowser* b)
+{
+  TNamed::Browse(b);
 }