]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALSDigitizer.cxx
Added quotas for TTL. Merging ignores extra files.
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
index b9d5cc2ce58df0b33d3df26e011e1e02b73e5af4..fd7087db36bca0d235bbfc6194ac17b64d9646ef 100644 (file)
 // 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".
 //
 //*-- 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 <TH1.h>
 #include <TBrowser.h>
-#include <Riostream.h>
+//#include <Riostream.h>
 #include <TMath.h>
+#include <TROOT.h>
 
 // --- Standard library ---
 #include "stdlib.h"
@@ -67,7 +71,7 @@
 #include "AliEMCALHit.h"
 #include "AliEMCALSDigitizer.h"
 #include "AliEMCALGeometry.h"
-#include "AliEMCALHistoUtilities.h"
+#include "AliEMCALSimParam.h"
 
 ClassImp(AliEMCALSDigitizer)
            
@@ -81,9 +85,7 @@ AliEMCALSDigitizer::AliEMCALSDigitizer()
     fSDigitsInRun(0),
     fFirstEvent(0),
     fLastEvent(0),
-    fSampling(0.),
-    fControlHists(0),
-    fHists(0)
+    fSampling(0.)
 {
   // ctor
   InitParameters();
@@ -100,14 +102,12 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName,
     fSDigitsInRun(0),
     fFirstEvent(0),
     fLastEvent(0),
-    fSampling(0.),
-    fControlHists(1),
-    fHists(0)
+    fSampling(0.)
 {
   // ctor
   Init();
   InitParameters() ; 
-  if(fControlHists) BookControlHists(1);
+
 }
 
 
@@ -123,9 +123,7 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd)
     fSDigitsInRun(sd.fSDigitsInRun),
     fFirstEvent(sd.fFirstEvent),
     fLastEvent(sd.fLastEvent),
-    fSampling(sd.fSampling),
-    fControlHists(sd.fControlHists),
-    fHists(sd.fHists)
+    fSampling(sd.fSampling)
 {
   //cpy ctor 
 }
@@ -135,7 +133,7 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd)
 AliEMCALSDigitizer::~AliEMCALSDigitizer() {
   //dtor
   AliLoader *emcalLoader=0;
-  if ((emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")))
+  if ((emcalLoader = AliRunLoader::Instance()->GetDetectorLoader("EMCAL")))
       emcalLoader->CleanSDigitizer();
 }
 
@@ -150,7 +148,7 @@ void AliEMCALSDigitizer::Init(){
 
   fInit = kTRUE; 
    
-  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
 
   if ( emcalLoader == 0 ) {
     Fatal("Init", "Could not obtain the AliEMCALLoader");
@@ -171,14 +169,46 @@ void AliEMCALSDigitizer::InitParameters()
   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         = 0;
-  fB         = 1.e+7;
+  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  = 0.; // 24-nov-04 - was 1.e-6;
-  Print("");
+  // 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"));
+
 }
 
 //____________________________________________________________________________
@@ -187,14 +217,26 @@ 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"));
+
     return ; 
   }
-  
+
   if(strstr(option,"tim"))
     gBenchmark->Start("EMCALSDigitizer");
 
-  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  AliRunLoader *rl = AliRunLoader::Instance();
   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 
   //switch off reloading of this task while getting event
@@ -228,7 +270,8 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
 
     Int_t nSdigits = 0 ;
     Int_t i;
-    AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); 
+
+       AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); 
     for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
       AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
       AliEMCALDigit * curSDigit = 0 ;
@@ -242,12 +285,14 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
          // Assign primary number only if deposited energy is significant
          curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
                                          hit->GetIparent(), hit->GetId(), 
-                                         Digitize(energy), hit->GetTime() ) ;
+                                         Digitize(energy), hit->GetTime(), 
+                                         -1, energy ) ;
          else
            curSDigit =  new AliEMCALDigit( -1               , 
                                            -1               ,
                                            hit->GetId(), 
-                                           Digitize(energy), hit->GetTime() ) ;
+                                           Digitize(energy), hit->GetTime(), 
+                                           -1, energy ) ;
       } else {
        Warning("Exec"," abs id %i is bad \n", hit->GetId());
        newsdigit = kFALSE;
@@ -274,20 +319,11 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
     
     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    
     
@@ -321,22 +357,47 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
 }
 
 //__________________________________________________________________
-
 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);
-    }
-    return (Int_t ) aSignal;
+  //
+  //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 (Int_t ) aSignal;
+}
 
 //__________________________________________________________________
+Float_t AliEMCALSDigitizer::Calibrate(Int_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(); 
@@ -375,10 +436,8 @@ Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
 void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
 {
   //Prints list of digits produced at the current pass of AliEMCALDigitizer
-  
-  
-  // AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
-  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+    
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
   const TClonesArray * sdigits = rl->SDigits() ; 
   
   printf("\n") ;  
@@ -395,17 +454,17 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
       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);
+      printf("%s",tempo);
       isum += digit->GetAmp();
       
       Int_t iprimary;
       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
-       sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
-       printf(tempo); 
+                 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
+       printf("%s",tempo); 
       }         
     }
     delete tempo ;
-    printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6);
+    printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, Calibrate(isum));
   } else printf("\n");
 }
 
@@ -413,40 +472,13 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
 void AliEMCALSDigitizer::Unload() const
 {
   // Unload Hits and SDigits from the folder
-  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+  AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->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){
-    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); 
-}