Improved the mechanism to divert SDigits into a separate file.
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
index 7420072..c29b59a 100644 (file)
@@ -61,6 +61,7 @@
 #include "TROOT.h"
 #include "TFolder.h"
 #include "TObjString.h"
+#include "TGeometry.h"
 #include "TBenchmark.h"
 // --- Standard library ---
 #include <iomanip.h>
@@ -68,6 +69,7 @@
 // --- AliRoot header files ---
 
 #include "AliRun.h"
+#include "AliHeader.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALHit.h"
 #include "AliEMCALTick.h"
@@ -86,21 +88,22 @@ ClassImp(AliEMCALDigitizer)
   // ctor
 
   fSDigitizer = 0 ;
-  fNinputs = 1 ;
-  fPinNoise = 0.01 ;
-  fEMCDigitThreshold = 0.01 ;
-  fTimeResolution     = 0.5e-9 ;
-  fTimeSignalLength   = 1.0e-9 ;
-  fPreShowerDigitThreshold = fEMCDigitThreshold/100. ;
-  fADCchannelTower = 0.0030;        // width of one ADC channel in GeV
-  fADCpedestalTower = 0.005 ;      // pedestal of ADC
-  fNADCTower = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
-
-  fADCchannelPreSho = 0.0030 ;          // width of one ADC channel in Pre Shower
-  fADCpedestalPreSho = 0.005 ;         // pedestal of ADC
-  fNADCPreSho = (Int_t) TMath::Power(2,16);      // number of channels in Pre Shower ADC
-  fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
+  fNinputs = 0 ;
+  fPinNoise = 0.0 ;
+  fTowerDigitThreshold = 0.0 ;
+  fTimeResolution     = 0. ;
+  fTimeSignalLength   = 0. ;
+  fPreShowerDigitThreshold = 0. ;
+  fADCchannelTower = 0.0;      // width of one ADC channel in GeV
+  fADCpedestalTower = 0. ;      // pedestal of ADC
+  fNADCTower = 0;  // number of channels in Tower ADC
+
+  fADCchannelPreSho  = 0.0;          // width of one ADC channel in Pre Shower
+  fADCpedestalPreSho = 0.0 ;         // pedestal of ADC
+  fNADCPreSho = 0;      // number of channels in Pre Shower ADC
+  fTimeThreshold = 0.0; //Means 1 MeV in terms of SDigits amplitude
   fManager = 0 ;
+  fSplitFile= 0 ; 
 
 
 
@@ -112,29 +115,29 @@ Bool_t AliEMCALDigitizer::Init()
 
   fSDigitizer = 0 ;
   fNinputs = 1 ;
-  fPinNoise = 0.01 ;
-  fEMCDigitThreshold = 0.01 ;
+  fPinNoise = 0.00001 ;
+  fTowerDigitThreshold = 0.001 ;
   fTimeResolution     = 0.5e-9 ;
   fTimeSignalLength   = 1.0e-9 ;
-  fPreShowerDigitThreshold = fEMCDigitThreshold/100. ;
+  fPreShowerDigitThreshold = fTowerDigitThreshold/25. ;
   fInitialized = kFALSE ;
-  fADCchannelTower = 0.0030;        // width of one ADC channel in GeV
-  fADCpedestalTower = 0.005 ;      //
+  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.0030 ;          // width of one ADC channel in Pre Shower
+  fADCchannelPreSho = 0.0000300;          // width of one ADC channel in Pre Shower
   fADCpedestalPreSho = 0.005 ;         // 
-  fNADCPreSho = (Int_t) TMath::Power(2,16);      // number of channels in Pre ShowerADC
+  fNADCPreSho = (Int_t) TMath::Power(2,12);      // number of channels in Pre ShowerADC
 
   fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
  
 
 
-if(fManager)
+  if(fManager)
     SetTitle("aliroot") ;
   else if (strcmp(GetTitle(),"")==0) 
-   SetTitle("galice.root") ;
-
+    SetTitle("galice.root") ;
+  
   if( strcmp(GetName(), "") == 0 )
     SetName("Default") ;
   
@@ -167,9 +170,10 @@ if(fManager)
 //____________________________________________________________________________ 
 AliEMCALDigitizer::AliEMCALDigitizer(const char *headerFile,const char *name)
 {
-   SetName(name) ;
+  SetName(name) ;
   SetTitle(headerFile) ;
   fManager = 0 ;                     // We work in the standalong mode
+  fSplitFile= 0 ; 
   Init() ;
   
 
@@ -188,7 +192,12 @@ AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * ard):AliDigitizer(ard)
 {
   // dtor
 
+  if (fSplitFile)       
+    if ( fSplitFile->IsOpen() )
+      fSplitFile->Close() ;
+
 }
+
 //____________________________________________________________________________
 void AliEMCALDigitizer::Reset() { 
   //sets current event number to the first simulated event
@@ -230,7 +239,6 @@ void AliEMCALDigitizer::Digitize(const Int_t event) {
     cerr << "ERROR: AliEMCALDigitizer::Digitize -> SDigitizer with name " << GetName() << " not found " << endl ; 
     abort() ; 
   }
-
 // loop through the sdigits posted to the White Board and add them to the noise
   TCollection * folderslist = gime->SDigitsFolder()->GetListOfFolders() ; 
   TIter next(folderslist) ; 
@@ -240,13 +248,16 @@ void AliEMCALDigitizer::Digitize(const Int_t event) {
   TObjArray * sdigArray = new TObjArray(2) ;
   while ( (folder = (TFolder*)next()) ) 
     if ( (sdigits = (TClonesArray*)folder->FindObject(GetName()) ) ) {
+      TString fileName(folder->GetName()) ;
+      fileName.ReplaceAll("_","/") ;
       cout << "INFO: AliEMCALDigitizer::Digitize -> Adding SDigits " 
-          << GetName() << " from " << folder->GetName() << endl ; 
+          << GetName() << " from " << fileName << endl ; 
       sdigArray->AddAt(sdigits, input) ;
       input++ ;
     }
 
-  //Find the first crystall with signal
+
+  //Find the first tower with signal
   Int_t nextSig = 200000 ; 
   Int_t i;
   for(i=0; i<input; i++){
@@ -268,13 +279,12 @@ void AliEMCALDigitizer::Digitize(const Int_t event) {
 
   //Put Noise contribution
   for(absID = 1; absID <= nEMC; absID++){
-    Float_t noise = gRandom->Gaus(0., fPinNoise) ; 
+    Float_t noise = gRandom->Gaus(0., fPinNoise); 
     new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
     //look if we have to add signal?
+    digit = (AliEMCALDigit *) digits->At(absID-1) ;
     if(absID==nextSig){
-      //Add SDigits from all inputs 
-      digit = (AliEMCALDigit *) digits->At(absID-1) ;
-      
+      //Add SDigits from all inputs    
       ticks->Clear() ;
       Int_t contrib = 0 ;
       Float_t a = digit->GetAmp() ;
@@ -342,13 +352,15 @@ void AliEMCALDigitizer::Digitize(const Int_t event) {
 
   //remove digits below thresholds
   for(absID = 0; absID < nEMC/2 ; absID++){
-    if(sDigitizer->Calibrate(((AliEMCALDigit*)digits->At(absID))->GetAmp()) < fEMCDigitThreshold)
+
+    if(sDigitizer->Calibrate(((AliEMCALDigit*)digits->At(absID))->GetAmp()) < fTowerDigitThreshold)
       digits->RemoveAt(absID) ;
     else
       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
   }
   
   for(absID = nEMC/2; absID < nEMC ; absID++){
+
     if(sDigitizer->Calibrate(((AliEMCALDigit*)digits->At(absID))->GetAmp()) < fPreShowerDigitThreshold)
       digits->RemoveAt(absID) ;
     else
@@ -364,7 +376,7 @@ void AliEMCALDigitizer::Digitize(const Int_t event) {
   
   //Set indexes in list of digits
   //Int_t i ;
-  for (i = 0 ; i < ndigits ; i++) { 
+ for (i = 0 ; i < ndigits ; i++) { 
     AliEMCALDigit * digit = (AliEMCALDigit *) digits->At(i) ; 
     digit->SetIndexInList(i) ; 
     Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ;
@@ -380,12 +392,15 @@ Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t absId)
   Int_t nphi = AliEMCALGetter::GetInstance()->EMCALGeometry()->GetNPhi() ; 
   Int_t nz   = AliEMCALGetter::GetInstance()->EMCALGeometry()->GetNZ() ;
   
-  if(absId <= nphi*nz) //digitize as tower
-    channel = (Int_t) TMath::Ceil( (energy + fADCpedestalTower)/fADCchannelTower )  ;
-  if(channel > fNADCTower ) channel =  fNADCTower ;
-  else 
-    channel =  (Int_t) TMath::Ceil( (energy + fADCpedestalPreSho)/fADCchannelPreSho )  ;
-  if(channel > fNADCPreSho ) channel =  fNADCPreSho ;
+  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 ;
+  }
   
   return channel ;
 }
@@ -419,38 +434,36 @@ if(strcmp(GetName(), "") == 0 )
     nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
     treeD=gAlice->TreeD() ;
   }
-  if(treeD == 0 ){
-    cerr << " AliEMCALDigitizer :: Can not find TreeD " << endl ;
-    return ;
-  }
 
   //Check, if this branch already exits
-  TObjArray * lob = (TObjArray*)treeD->GetListOfBranches() ;
-  TIter next(lob) ; 
-  TBranch * branch = 0 ;  
-  Bool_t emcalfound = kFALSE, digitizerfound = kFALSE ; 
-  
-  while ( (branch = (TBranch*)next()) && (!emcalfound || !digitizerfound) ) {
-    if ( (strcmp(branch->GetName(), "EMCAL")==0) && 
-        (strcmp(branch->GetTitle(), GetName())==0) ) 
-      emcalfound = kTRUE ;
+  if (treeD) {
+    TObjArray * lob = (TObjArray*)treeD->GetListOfBranches() ;
+    TIter next(lob) ; 
+    TBranch * branch = 0 ;  
+    Bool_t emcalfound = kFALSE, digitizerfound = kFALSE ; 
     
-    else if ( (strcmp(branch->GetName(), "AliEMCALDigitizer")==0) && 
-             (strcmp(branch->GetTitle(), GetName())==0) ) 
-      digitizerfound = kTRUE ; 
+    while ( (branch = (TBranch*)next()) && (!emcalfound || !digitizerfound) ) {
+      if ( (strcmp(branch->GetName(), "EMCAL")==0) && 
+          (strcmp(branch->GetTitle(), GetName())==0) ) 
+       emcalfound = kTRUE ;
+      
+      else if ( (strcmp(branch->GetName(), "AliEMCALDigitizer")==0) && 
+               (strcmp(branch->GetTitle(), GetName())==0) ) 
+       digitizerfound = kTRUE ; 
+    }
+    
+    if ( emcalfound ) {
+      cerr << "WARNING: AliEMCALDigitizer -> Digits branch with name " << GetName() 
+          << " already exits" << endl ;
+      return ; 
+    }   
+    if ( digitizerfound ) {
+      cerr << "WARNING: AliEMCALDigitizer -> Digitizer branch with name " << GetName() 
+          << " already exits" << endl ;
+      return ; 
+    }   
   }
-
-  if ( emcalfound ) {
-    cerr << "WARNING: AliEMCALDigitizer -> Digits branch with name " << GetName() 
-        << " already exits" << endl ;
-    return ; 
-  }   
-  if ( digitizerfound ) {
-    cerr << "WARNING: AliEMCALDigitizer -> Digitizer branch with name " << GetName() 
-        << " already exits" << endl ;
-    return ; 
-  }   
-
   Int_t ievent ;
 
   for(ievent = 0; ievent < nevents; ievent++){
@@ -466,8 +479,8 @@ if(strcmp(GetName(), "") == 0 )
        gime->ReadTreeS(treeS,input) ;
       }
     }
-    else
-      gime->Event(ievent,"S") ;
+    else 
+      gime->Event(ievent,"S") ; 
     
     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
     
@@ -480,6 +493,10 @@ if(strcmp(GetName(), "") == 0 )
     fDigitsInRun += gime->Digits()->GetEntriesFast() ;  
   }
   
+  if (fSplitFile) 
+    if ( fSplitFile->IsOpen() ) 
+      fSplitFile->Close() ; 
+
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALDigitizer");
     cout << "AliEMCALDigitizer:" << endl ;
@@ -556,7 +573,7 @@ void AliEMCALDigitizer::Print(Option_t* option)const {
       cout << endl ;
       cout << "With following parameters: " << endl ;
       cout << "     Electronics noise in EMC (fPinNoise) = " << fPinNoise << endl ;
-      cout << "  Threshold  in EMC  (fEMCDigitThreshold) = " << fEMCDigitThreshold  << endl;
+      cout << "  Threshold  in EMC  (fTowerDigitThreshold) = " << fTowerDigitThreshold  << endl;
       cout << "  Threshold  in PreShower  (fPreShowerDigitThreshold) = " << fPreShowerDigitThreshold  << endl ; ;
       cout << "---------------------------------------------------" << endl ;
     }
@@ -610,53 +627,28 @@ void AliEMCALDigitizer::WriteDigits(Int_t event)
   const TClonesArray * digits = gime->Digits(GetName()) ; 
  TTree * treeD ;
 
-  if(fManager)
-    treeD = fManager->GetTreeD() ;
- else
-    treeD = gAlice->TreeD();
-  
-  // create new branches
-  // -- generate file name if necessary
-  char * file =0;
-  if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
-    file = new char[strlen(gAlice->GetBaseFile())+20] ;
-    sprintf(file,"%s/EMCAL.Digits.root",gAlice->GetBaseFile()) ;
-  }
-
-  TDirectory *cwd = gDirectory;
-  // -- create Digits branch
-  Int_t bufferSize = 32000 ;    
-  TBranch * digitsBranch = treeD->Branch("EMCAL",&digits,bufferSize);
-  digitsBranch->SetTitle(GetName());
-  if (file) {
-    digitsBranch->SetFile(file);
-    TIter next( digitsBranch->GetListOfBranches());
-    TBranch * sbr ;
-    while ((sbr=(TBranch*)next())) {
-      sbr->SetFile(file);
-    }   
-    cwd->cd();
-  } 
-    
-  // -- Create Digitizer branch
-  Int_t splitlevel = 0 ;
-  AliEMCALDigitizer * d = gime->Digitizer(GetName()) ;
-  TBranch * digitizerBranch = treeD->Branch("AliEMCALDigitizer", "AliEMCALDigitizer", &d,bufferSize,splitlevel); 
-  digitizerBranch->SetTitle(GetName());
-  if (file) {
-    digitizerBranch->SetFile(file);
-    TIter next( digitizerBranch->GetListOfBranches());
-    TBranch * sbr;
-    while ((sbr=(TBranch*)next())) {
-      sbr->SetFile(file);
-    }   
-    cwd->cd();
-  }
-  
-  digitsBranch->Fill() ;      
-  digitizerBranch->Fill() ;
-
-  treeD->Write(0,kOverwrite) ;  
+ if(fManager)
+   treeD = fManager->GetTreeD() ;
+ else {
+   if (!gAlice->TreeD() ) 
+     gAlice->MakeTree("D",fSplitFile);
+   treeD = gAlice->TreeD();
+ }
+ // -- create Digits branch
+ Int_t bufferSize = 32000 ;    
+ TBranch * digitsBranch = treeD->Branch("EMCAL",&digits,bufferSize);
+ digitsBranch->SetTitle(GetName());
+ // -- Create Digitizer branch
+ Int_t splitlevel = 0 ;
+ AliEMCALDigitizer * d = gime->Digitizer(GetName()) ;
+ TBranch * digitizerBranch = treeD->Branch("AliEMCALDigitizer", "AliEMCALDigitizer", &d,bufferSize,splitlevel); 
+ digitizerBranch->SetTitle(GetName());
+ digitsBranch->Fill() ;
+ digitizerBranch->Fill() ; 
+ treeD->AutoSave() ;  
  
 }
 //____________________________________________________________________________ 
@@ -694,3 +686,55 @@ void AliEMCALDigitizer::SetSDigitsBranch(const char* title)
 
   AliEMCALGetter::GetInstance()->SDigits()->SetName(title) ; 
 }
+
+//__________________________________________________________________
+void AliEMCALDigitizer::SetSplitFile(const TString splitFileName)
+{
+  // Diverts the Digits in a file separate from the hits file
+  
+  // I guess it is not going to work if we do merging
+  if (fManager) {
+    cerr << "ERROR: AliEMCALDigitizer::SetSplitFile -> Not yet available in case of merging activated " << endl ;  
+    return ; 
+  }
+
+  TDirectory * cwd = gDirectory ;
+  if ( !(gAlice->GetTreeDFileName() == splitFileName) ) {
+    if (gAlice->GetTreeDFile() )  
+      gAlice->GetTreeDFile()->Close() ; 
+  }
+
+  fSplitFile = gAlice->InitTreeFile("D",splitFileName.Data());
+  fSplitFile->cd() ; 
+  if ( !fSplitFile->Get("gAlice") ) 
+    gAlice->Write();
+  
+  TTree *treeE  = gAlice->TreeE();
+  if (!treeE) {
+    cerr << "ERROR: AliEMCALDigitizer::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: AliEMCALDigitizer::SetSplitFile -> AliceGeom was not found in the input file "<<endl;
+      abort() ;
+    }
+    AliceGeom->Write();
+  }
+  
+  gAlice->MakeTree("D",fSplitFile);
+  cwd->cd() ; 
+  cout << "INFO: AliEMCALDigitizer::SetSPlitMode -> Digits will be stored in " << splitFileName.Data() << endl ; 
+}