]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALDigitizer.cxx
new classes for track segments
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
index b2b4dc4ab969176e0bc2068a21b00f99d9623cf1..9c62d73a618a9990c22ea5012e5e7b292de273eb 100644 (file)
@@ -68,7 +68,6 @@
 #include "TBenchmark.h"
 
 // --- Standard library ---
-#include <Riostream.h>
 
 // --- AliRoot header files ---
 #include "AliRun.h"
@@ -90,7 +89,7 @@ ClassImp(AliEMCALDigitizer)
   AliEMCALDigitizer::AliEMCALDigitizer()
 {
   // ctor
-  InitParameters() ; 
+  //InitParameters() ; 
   fDefaultInit = kTRUE ; 
   fManager = 0 ;                     // We work in the standalong mode
  
@@ -105,9 +104,9 @@ AliEMCALDigitizer::AliEMCALDigitizer(const char *headerFile, const char *name, c
   SetName(name) ;
   fManager = 0 ;                     // We work in the standalong mode
   fSplitFile= 0 ; 
-  InitParameters() ; 
   fToSplit = toSplit ;
   Init() ;
+  InitParameters() ; 
   fDefaultInit = kFALSE ; 
 }
 
@@ -155,16 +154,19 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
 
   const AliEMCALGeometry *geom = gime->EMCALGeometry() ; 
   //Making digits with noise, first EMC
-  Int_t nEMC = 2*geom->GetNPhi()*geom->GetNZ();
-
-  Int_t absID ;
+  Int_t nEMC = 0 ; 
+  if (geom->GetNHCLayers() > 0 )
+    nEMC = 3*geom->GetNPhi()*geom->GetNZ(); //max number of digits possible (Preshower, ECAL, HCAL)
+  else 
+    nEMC = 2*geom->GetNPhi()*geom->GetNZ(); //max number of digits possible (Preshower, ECAL)
+ Int_t absID ;
   TString name      =  geom->GetName() ;
 
  // get first the sdigitizer from the tasks list (must have same name as the digitizer)
   const AliEMCALSDigitizer * sDigitizer = gime->SDigitizer(GetName()); 
   if ( !sDigitizer) {
-    cerr << "ERROR: AliEMCALDigitizer::Digitize -> SDigitizer with name " << GetName() << " not found " << endl ; 
-    abort() ; 
+    Fatal("Digitize", "SDigitizer with name %s not found", GetName() ); 
   }
 
 // loop through the sdigits posted to the White Board and add them to the noise
@@ -178,15 +180,15 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
     if ( (sdigits = (TClonesArray*)folder->FindObject(GetName()) ) ) {
       TString fileName(folder->GetName()) ;
       fileName.ReplaceAll("_","/") ;
-      //      cout << "INFO: AliEMCALDigitizer::Digitize -> Adding SDigits " 
-      //          << GetName() << " from " << fileName << endl ; 
+      if (gDebug)
+       Info("Digitize", "Adding SDigits %s from %s", GetName(), fileName.Data()) ; 
       sdigArray->AddAt(sdigits, input) ;
       input++ ;
     }
   }
 
   //Find the first tower with signal
-  Int_t nextSig = 200000 ; 
+  Int_t nextSig = nEMC + 1 ; 
   Int_t i;
   for(i=0; i<input; i++){
     sdigits = (TClonesArray *)sdigArray->At(i) ;
@@ -207,11 +209,13 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
 
   //Put Noise contribution
   for(absID = 1; absID <= nEMC; absID++){
-    Float_t noise = gRandom->Gaus(0., fPinNoise); 
-    new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
+    Float_t amp = 0 ; 
+    //    Float_t noise = TMath::Abs(gRandom->Gaus(0., fPinNoise)); 
+    //new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID, sDigitizer->Digitize(noise), TimeOfNoise() ) ;
+    // amplitude set to zero, noise will be added later
+    new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ) ;
     //look if we have to add signal?
     digit = (AliEMCALDigit *) digits->At(absID-1) ;
     if(absID==nextSig){
       //Add SDigits from all inputs    
       ticks->Clear() ;
@@ -253,13 +257,19 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
            curSDigit = 0 ;
        }
       }
-
+      // add fluctuations for photo-electron creation
+      amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
+      if (gDebug==1) 
+       Info("Digitize", "id = %d BEFORE amp = %f\n", digit->GetId(),amp) ; 
+      amp *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
+      if (gDebug==1) 
+       Info("Digitize", "After fluctuation amp = %f\n", amp) ; 
       //calculate and set time
       Float_t time = FrontEdgeTime(ticks) ;
       digit->SetTime(time) ;
 
       //Find next signal module
-      nextSig = 200000 ;
+      nextSig = nEMC + 1 ;
       for(i=0; i<input; i++){
        sdigits = ((TClonesArray *)sdigArray->At(i)) ;
        Int_t curNext = nextSig ;
@@ -270,6 +280,17 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
        if(curNext < nextSig) nextSig = curNext ;
       }
     }
+    // add the noise now
+    
+    if (geom->IsInECAL(digit->GetId())) 
+      amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
+    else if (geom->IsInPRE(digit->GetId())) 
+      amp += TMath::Abs(gRandom->Gaus(0., fPinNoise/100.)) ; // arbitrarely divide by 100, assuming that the gain of APD will be higher
+    else if (geom->IsInHCAL(digit->GetId())) 
+      amp += TMath::Abs(gRandom->Gaus(0., fPinNoise/10.)) ;  // arbitrarely divide by 10, assuming that the gain of APD will be higher
+   if (gDebug==1) 
+      Info("Digitize", "After noise amp = %f \n", amp) ; 
+    digit->SetAmp(sDigitizer->Digitize(amp)) ;  
   }
   
   ticks->Delete() ;
@@ -278,27 +299,30 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
   delete sdigArray ; //We should not delete its contents
 
   //remove digits below thresholds
-  for(absID = 0; absID < nEMC/2 ; absID++){
+  for(absID = 0; absID < nEMC; absID++){
     digit = dynamic_cast<AliEMCALDigit*>( digits->At(absID) ) ;
-    if(sDigitizer->Calibrate( digit->GetAmp() ) < fTowerDigitThreshold)
+    Float_t threshold = 0 ; 
+
+    if (geom->IsInECAL(digit->GetId())) 
+      threshold = fDigitThreshold ; 
+    else  if (geom->IsInPRE(digit->GetId()))
+      threshold = fDigitThreshold / 100. ; // arbitrary see before when noise is added
+    else  if (geom->IsInHCAL(digit->GetId()))
+      threshold = fDigitThreshold / 10. ; // arbitrary see before when noise is added    
+
+    if(sDigitizer->Calibrate( digit->GetAmp() ) <= threshold)
       digits->RemoveAt(absID) ;
-    else
+    else {
+      if (gDebug == 1)
+       Info("Digitize", "id = %d, amp = %f, noise = %f", absID, sDigitizer->Calibrate(digit->GetAmp()), fPinNoise) ; 
       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
-  }
-  
-    
-    for(absID = nEMC/2; absID < nEMC ; absID++){
-      digit = dynamic_cast<AliEMCALDigit*>( digits->At(absID) ) ;
-      if(sDigitizer->Calibrate( digit->GetAmp() ) < fPreShowerDigitThreshold)
-       digits->RemoveAt(absID) ;
-      else
-       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
     }
+  }
   
   digits->Compress() ;  
   
   Int_t ndigits = digits->GetEntriesFast() ; 
-  digits->Expand(ndigits) ;
+  //  digits->Expand(ndigits) ;
   
   //Set indexes in list of digits
   for (i = 0 ; i < ndigits ; i++) { 
@@ -314,17 +338,22 @@ void AliEMCALDigitizer::Digitize(const Int_t event)
 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t absId)
 { 
   Int_t channel = -999;
-  Int_t nphi = AliEMCALGetter::GetInstance()->EMCALGeometry()->GetNPhi() ; 
-  Int_t nz   = AliEMCALGetter::GetInstance()->EMCALGeometry()->GetNZ() ;
+  AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ; 
   
-  if(absId <= nphi*nz){  //digitize as tower
-    channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalTower)/fADCchannelTower ))  ;
-  if(channel > fNADCTower ) 
-    channel =  fNADCTower ;
-  } else {
-    channel =  static_cast<Int_t>(TMath::Ceil( (energy + fADCpedestalPreSho)/fADCchannelPreSho ))  ;
-  if(channel > fNADCPreSho ) 
-    channel =  fNADCPreSho ;
+  if(geom->IsInPRE(absId)){        //digitize as PRE section
+    channel =  static_cast<Int_t>(TMath::Ceil( (energy + fADCpedestalPRE)/fADCchannelPRE ))  ;
+    if(channel > fNADCPRE ) 
+      channel =  fNADCPRE ;
+  }
+  else if(geom->IsInECAL(absId)){  //digitize as ECAL section
+    channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
+    if(channel > fNADCEC ) 
+      channel =  fNADCEC ;
+  } 
+  else if(geom->IsInHCAL(absId)){  //digitize as HCAL section
+    channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalHC)/fADCchannelHC ))  ;
+    if(channel > fNADCHC ) 
+      channel =  fNADCHC ;
   }
   
   return channel ;
@@ -373,14 +402,12 @@ void AliEMCALDigitizer::Exec(Option_t *option)
       }
       
       if ( emcalfound ) {
-       cerr << "WARNING: AliEMCALDigitizer -> Digits branch with name " << GetName() 
-            << " already exits" << endl ;
+       Error( "Exec", "Digits branch with name %s already exits", GetName() ) ;
        return ; 
       }   
       if ( digitizerfound ) {
-       cerr << "WARNING: AliEMCALDigitizer -> Digitizer branch with name " << GetName() 
-            << " already exits" << endl ;
-       return ; 
+       Error( "Exec", "Digitizer branch with name %s already exit", GetName() ) ;
+              return ; 
       }
     }   
   }
@@ -400,7 +427,7 @@ void AliEMCALDigitizer::Exec(Option_t *option)
       for(input = 0 ; input < fManager->GetNinputs(); input ++){
        TTree * treeS = fManager->GetInputTreeS(input) ;
        if(!treeS){
-         cerr << "AliEMCALDigitizer -> No Input " << endl ;
+         Error( "Exec", "No Input") ;
          return ;
        }
        gime->ReadTreeS(treeS,input) ;
@@ -423,10 +450,8 @@ void AliEMCALDigitizer::Exec(Option_t *option)
   
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALDigitizer");
-    cout << "AliEMCALDigitizer:" << endl ;
-    cout << "  took " << gBenchmark->GetCpuTime("EMCALDigitizer") << " seconds for Digitizing " 
-        <<  gBenchmark->GetCpuTime("EMCALDigitizer")/nevents << " seconds per event " << endl ;
-    cout << endl ;
+    Info("Exec", "took %f seconds for Digitizing %f seconds per event", 
+        gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nevents ) ;
   }
   
 }
@@ -461,7 +486,7 @@ Bool_t AliEMCALDigitizer::Init()
   
   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), fToSplit) ; 
   if ( gime == 0 ) {
-    cerr << "ERROR: AliEMCALDigitizer::Init -> Could not obtain the Getter object !" << endl ; 
+    Error("Init", "Could not obtain the Getter object !" ) ; 
     return kFALSE;
   } 
   
@@ -509,20 +534,26 @@ Bool_t AliEMCALDigitizer::Init()
 //____________________________________________________________________________ 
 void AliEMCALDigitizer::InitParameters()
 {
-  fPinNoise = 0.00001 ;
-  fTowerDigitThreshold = 0.001 ;
+  fMeanPhotonElectron = 1250 ; // electrons per GeV
+  
+  fPinNoise           = 0.001 ; // noise equivalent GeV (random choice)
+  fDigitThreshold     = fPinNoise * 3; //2 sigma
   fTimeResolution     = 0.5e-9 ;
   fTimeSignalLength   = 1.0e-9 ;
-  fPreShowerDigitThreshold = fTowerDigitThreshold/25. ;
-  fADCchannelTower = 0.000220;       // width of one ADC channel in GeV
-  fADCpedestalTower = 0.005 ;      // GeV
-  fNADCTower = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
 
-  fADCchannelPreSho = 0.0000300;          // width of one ADC channel in Pre Shower
-  fADCpedestalPreSho = 0.005 ;         // 
-  fNADCPreSho = (Int_t) TMath::Power(2,12);      // number of channels in Pre ShowerADC
+  fADCchannelEC    = 0.000220;                     // width of one ADC channel in GeV
+  fADCpedestalEC   = 0.005 ;                       // GeV
+  fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
+
+  fADCchannelHC    = 0.000220;                     // width of one ADC channel in GeV
+  fADCpedestalHC   = 0.005 ;                       // GeV
+  fNADCHC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
 
-  fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
+  fADCchannelPRE   = 0.0000300;                    // width of one ADC channel in Pre Shower
+  fADCpedestalPRE  = 0.005 ;                       // GeV 
+  fNADCPRE         = (Int_t) TMath::Power(2,12);   // number of channels in Pre ShowerADC
+
+  fTimeThreshold      = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
  
 }
 
@@ -545,7 +576,7 @@ void AliEMCALDigitizer::MixWith(char* headerFile)
     Init() ;
   
   if(fManager){
-    cout << "Can not use this method under AliRunDigitizer " << endl ;
+    Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
     return ;
   } 
   
@@ -560,7 +591,7 @@ void AliEMCALDigitizer::MixWith(char* headerFile)
   path += "/" ; 
   path += GetName() ;
   if ( gROOT->FindObjectAny(path.Data()) ) {
-    cerr << "WARNING: AliEMCALDigitizer::MixWith -> Entry already exists, do not add" << endl ;
+    Error("MixWith", "Entry already exists, do not add" ) ;
     return;
   }
 
@@ -571,7 +602,7 @@ void AliEMCALDigitizer::MixWith(char* headerFile)
   if ( !file ) { 
     file = new TFile(headerFile, "READ") ; 
     if (!file) { 
-      cerr << "ERROR: AliEMCALDigitizer::MixWith -> File " << headerFile << " does not exist!" << endl ; 
+      Error("MixWith", "File %s does not exist!", headerFile) ; 
       return ; 
     }
   }
@@ -580,42 +611,40 @@ void AliEMCALDigitizer::MixWith(char* headerFile)
 
 //__________________________________________________________________
 void AliEMCALDigitizer::Print(Option_t* option)const {
+
+  TString message("\n") ; 
+
   if( strcmp(GetName(), "") != 0) {
-    
-    cout << "------------------- "<< GetName() << " -------------" << endl ;
+    message += "------------------- " ; 
+    message += GetName() ; 
+    message += " -------------" ;
     const Int_t nStreams = GetNInputStreams() ; 
     if (nStreams) {
       Int_t index = 0 ;  
-      for (index = 0 ; index < nStreams ; index++)  
-       cout << "Adding SDigits " << GetName() << " from " <<  fManager->GetInputFileName(index, 0) << endl ; 
-      
-      cout << endl ;
-      cout << "Writing digits to " <<   fManager->GetInputFileName(0, 0) << endl ;   
+      for (index = 0 ; index < nStreams ; index++) {  
+       message += "\nAdding SDigits " ; 
+       message += GetName() ; 
+       message += " from " ; 
+       message += fManager->GetInputFileName(index, 0) ;
+      } 
+     
+      message += "\nWriting digits to " ; 
+      message += fManager->GetInputFileName(0, 0) ;   
     } else { 
-//       AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;  
-//       gime->Folder("sdigits")  ;
-//       cout << "Digitizing sDigits from file(s): " <<endl ;
-//       TCollection * folderslist = gime->Folder("sdigits")->GetListOfFolders() ; 
-//       TIter next(folderslist) ; 
-//       TFolder * folder = 0 ; 
-      
-//       while ( (folder = (TFolder*)next()) ) {
-//     if ( folder->FindObject(GetName())  ) 
-//      cout << "Adding SDigits " << GetName() << " from " << GetSDigitsFileName() << endl ; 
-//      }
-      cout << endl ;
-      cout << "Writing digits to " << GetTitle() << endl ;
+      message += "\nWriting digits to " ;
+      message += GetTitle() ;
     }       
-    cout << endl ;
-    cout << "With following parameters: " << endl ;
-    cout << "     Electronics noise in EMC (fPinNoise) = " << fPinNoise << endl ;
-    cout << "  Threshold  in EMC  (fTowerDigitThreshold) = " << fTowerDigitThreshold  << endl;
-    cout << "  Threshold  in PreShower  (fPreShowerDigitThreshold) = " << fPreShowerDigitThreshold  << endl ; ;
-    cout << "---------------------------------------------------" << endl ;
+    message += "\nWith following parameters: " ;
+    message += "\n     Electronics noise in EMC (fPinNoise) = " ; 
+    message += fPinNoise ;
+    message += "\n  Threshold  in EMC  (fDigitThreshold) = " ; 
+    message += fDigitThreshold  ;
+    message += "\n---------------------------------------------------"  ;
   }
   else
-    cout << "AliEMCALDigitizer not initialized " << endl ;
+    message += "\nAliEMCALDigitizer not initialized " ;
+
+  Info("Print", message.Data() ) ; 
 }
 
 //__________________________________________________________________
@@ -624,28 +653,31 @@ void AliEMCALDigitizer::PrintDigits(Option_t * option){
   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
   TClonesArray * fDigits = gime->Digits() ;
 
-  cout << "AliEMCALDigitiser:"<< endl ;
-  cout << "       Number of entries in Digits list " << fDigits->GetEntriesFast() << endl ;
-  cout << endl ;
+  TString message("\n") ; 
+  message += "       Number of entries in Digits list " ; 
+  message += fDigits->GetEntriesFast() ;
+
   if(strstr(option,"all")){
     
     //loop over digits
     AliEMCALDigit * digit;
-    cout << "Digit Id " << " Amplitude " <<  " Index "  <<  " Nprim " << " Primaries list " <<  endl;      
+    message += "\n   Id  Amplitude    Time          Index Nprim: Primaries list \n" ;    
     Int_t index ;
+    char * tempo = new char[8192]; 
     for (index = 0 ; index < fDigits->GetEntries() ; index++) {
       digit = (AliEMCALDigit * )  fDigits->At(index) ;
-      cout << setw(8)  <<  digit->GetId() << " "  <<   setw(3)  <<  digit->GetAmp() <<   "  "  
-          << setw(6)  <<  digit->GetIndexInList() << "  "   
-          << setw(5)  <<  digit->GetNprimary() <<"  ";
-      
+      sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d : ",
+             digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
+      message += tempo ; 
       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) ) ; 
+       message += tempo ; 
+      }
+    }   
+    delete tempo ; 
   }
+  Info("PrintDigits", message.Data() ) ; 
 }
 
 //__________________________________________________________________