Correct the check of an already existing branch with the same name in Exec()
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
index 5008236c4d4b541949cde58bb683db9f6deafa93..2d528944961f510042665cbe733848e3557d5dbd 100644 (file)
@@ -73,15 +73,13 @@ ClassImp(AliEMCALSDigitizer)
   AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","") 
 {
   // ctor
-  fA = 0;
-  fB = 0. ;
-  fPrimThreshold = 0. ;
-  fNevents = 0 ;     
-  fSDigits = 0 ;
-  fHits = 0 ;
-  fLayerRatio = 0. ;
+  fA = fB =  fNevents = 0 ; 
+  fTowerPrimThreshold = fPreShowerPrimThreshold = fPhotonElectronFactor  = 0. ;
+  fHits = fSDigits = fSDigits = 0 ;
+
   fIsInitialized = kFALSE ;
 
+
 }
 
 //____________________________________________________________________________ 
@@ -90,9 +88,10 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigi
   // ctor
   fA = 0;
   fB = 10000000.;
-  fPrimThreshold = 0.01 ;
+  fTowerPrimThreshold = 0.01 ;
+  fPreShowerPrimThreshold = 0.0001 ; 
   fNevents = 0 ; 
-  fLayerRatio = 5./6. ;
+  fPhotonElectronFactor = 5000. ; // photoelectrons per GeV 
   Init();
 }
 
@@ -129,10 +128,9 @@ if( strcmp(GetTitle(), "") == 0 )
   TString sdname(GetName() );
   sdname.Append(":") ;
   sdname.Append(GetTitle() ) ;
-  SetName(sdname) ;
+  SetName(sdname.Data()) ;
   gime->PostSDigitizer(this) ;
  
-  // create a folder on the white board //YSAlice/WhiteBoard/SDigits/EMCAL/headerFile/sdigitsTitle
  
  }
 //____________________________________________________________________________
@@ -150,19 +148,25 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
   
   if(strstr(option,"tim"))
     gBenchmark->Start("EMCALSDigitizer");
+
+
   //Check, if this branch already exits
   gAlice->GetEvent(0) ;
+  TString sdname(GetName()) ;
+  sdname.Remove(sdname.Index(GetTitle())-1) ;
+
   if(gAlice->TreeS() ) {
     TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
     TIter next(lob) ; 
     TBranch * branch = 0 ;  
     Bool_t emcalfound = kFALSE, sdigitizerfound = kFALSE ; 
     
-    while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
+     while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
       if ( (strcmp(branch->GetName(), "EMCAL")==0) && (strcmp(branch->GetTitle(), GetName())==0) ) 
        emcalfound = kTRUE ;
       
-      else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) ) 
+      else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && (strcmp(branch->GetTitle(), sdname)==0) ) 
        sdigitizerfound = kTRUE ; 
     }
     
@@ -172,8 +176,7 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
       return ; 
     }   
   }  
-  TString sdname(GetName()) ;
-  sdname.Remove(sdname.Index(GetTitle())-1) ;
+
   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;  
   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
   
@@ -195,10 +198,7 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
     // for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
     //gime->Track(itrack);  
       //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
-      
-      
-
+     
       Int_t i;
       for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { // loop over all hits (hit = deposited energy/layer/entering particle)
        AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(fHits->At(i)) ;
@@ -206,25 +206,23 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
         AliEMCALDigit * sdigit = 0 ;
         Bool_t newsdigit = kTRUE; 
 
+       
+
        // Assign primary number only if deposited energy is significant
          
-       Float_t preShowerFactor ;
-       if (hit->IsInPreShower() ) 
-         preShowerFactor = fLayerRatio ; 
-       else 
-         preShowerFactor = 1 ; 
-       
-       if( hit->GetEnergy() > fPrimThreshold)
+       if( (!hit->IsInPreShower() && hit->GetEnergy() > fTowerPrimThreshold) || 
+           (hit->IsInPreShower() && hit->GetEnergy() > fPreShowerPrimThreshold)) {
          curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
                                          hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()), 
-                                         Digitize( hit->GetEnergy() * preShowerFactor ) 
+                                         Digitize(hit->GetEnergy())
                                          hit->GetTime()) ;
-        else
+       } else {
          curSDigit =  new AliEMCALDigit( -1               , 
                                          -1               ,
                                          Layer2TowerID(hit->GetId(),hit->IsInPreShower()), 
-                                         Digitize( hit->GetEnergy() * preShowerFactor ) 
+                                         Digitize(hit->GetEnergy())
                                          hit->GetTime() ) ;    
+       }
        Int_t check = 0 ;
        for(check= 0; check < nSdigits ; check++) {
           sdigit = (AliEMCALDigit *)sdigits->At(check);
@@ -249,8 +247,10 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
       //  Int_t i ;
      const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
      
-     Int_t lastPreShowerIndex = nSdigits - 1 ; 
-     Int_t firstPreShowerIndex = -1 ; 
+     Int_t lastPreShowerIndex = nSdigits - 1 ;
+     if (!(dynamic_cast<AliEMCALDigit *>(sdigits->At(lastPreShowerIndex))->IsInPreShower()))
+       lastPreShowerIndex = -2; 
+     Int_t firstPreShowerIndex = 100000 ; 
      Int_t index ; 
      AliEMCALDigit * sdigit = 0 ;
      for ( index = 0; index < nSdigits ; index++) {
@@ -272,8 +272,12 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
       Int_t jndex ; 
       for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
        tower  = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) ); 
-       if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {
-         *tower = *tower + *preshower ; // and add preshower to tower
+       if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {      
+         Float_t towerEnergy  = static_cast<Float_t>(tower->GetAmp()) ; 
+         Float_t preshoEnergy = static_cast<Float_t>(preshower->GetAmp()) ; 
+         towerEnergy +=preshoEnergy ; 
+         *tower = *tower + *preshower    ; // and add preshower multiplied by layer ratio to tower
+         tower->SetAmp(static_cast<Int_t>(TMath::Ceil(towerEnergy))) ; 
          towerFound = kTRUE ;
        }
       }
@@ -285,17 +289,20 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
        lastIndex++ ; 
       }
     }
-
+    
     sdigits->Sort() ;
-     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
-       sdigit->SetIndexInList(i) ;     
-      }
-      
-      if(gAlice->TreeS() == 0)
-       gAlice->MakeTree("S") ;
-      
-      
-      
+    Int_t NPrimarymax = -1 ; 
+    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
+      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()) ;
+       }
+       if(gAlice->TreeS() == 0)
+       gAlice->MakeTree("S") ;      
       
       //Make (if necessary) branches    
       char * file =0;
@@ -310,7 +317,6 @@ void AliEMCALSDigitizer::Exec(Option_t *option) {
       Int_t bufferSize = 32000 ;    
       TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
       sdigitsBranch->SetTitle(sdname);
-      cout << " AliEMCALSDigitizer::Exec sdname " << sdname << endl ; 
       
       if (file) {
        sdigitsBranch->SetFile(file);
@@ -391,15 +397,18 @@ void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
 void AliEMCALSDigitizer::Print(Option_t* option)const{
   cout << "------------------- "<< GetName() << " -------------" << endl ;
   cout << "   Writing SDigitis to branch with title  " << GetName() << endl ;
-  cout << "   with digitization parameters  A = " << fA << endl ;
-  cout << "                                 B = " << fB << endl ;
-  cout << "   Threshold for Primary assignment= " << fPrimThreshold << endl ; 
+  cout << "   with digitization parameters  A               = " << fA << endl ;
+  cout << "                                 B               = " << fB << endl ;
+  cout << "   Threshold for Primary assignment in Tower     = " << fTowerPrimThreshold << endl ; 
+  cout << "   Threshold for Primary assignment in PreShower = " << fPreShowerPrimThreshold << endl ; 
   cout << "---------------------------------------------------"<<endl ;
   
 }
 //__________________________________________________________________
 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
-  if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
+  if( (fA==sd.fA)&&(fB==sd.fB)&&
+      (fTowerPrimThreshold==sd.fTowerPrimThreshold) &&
+      (fPreShowerPrimThreshold==sd.fPreShowerPrimThreshold))
     return kTRUE ;
   else
     return kFALSE ;
@@ -437,21 +446,37 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
   }
 }
 //________________________________________________________________________
-Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t s, Bool_t preshower)
-{
+Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower){
+  // Method to Transform from Hit Id to Digit Id
+  // This function should be one to one
   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
   const AliEMCALGeometry * geom = gime->EMCALGeometry();
-  Int_t a  = (s/geom->GetNPhi())%geom->GetNZ()+1; // Phi Tower Index
-  Int_t b = (s-1)%(geom->GetNPhi())+1; //Eta Tower Index
-  Int_t t = -10;
+  Int_t ieta  = (ihit/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index
+  Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index
+  Int_t it = -10;
+  Int_t ipre = 0;
 
-if (a > 0 && b > 0)
-{
-  t = a*b + preshower*geom->GetNPhi()*geom->GetNZ();
-  return t;
-}
-else
-{cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Phi number = "
-      << b << "Eta number = " << a << endl ;
-return t;}
+  if (preshower)ipre = 1;
+  if (iphi > 0 && ieta >= 0){
+    it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ();
+    return it;
+  }else{
+    cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Eta number = "
+        << ieta << "Phi number = " << iphi << endl ;
+    return it;
+  } // end if iphi>0 && ieta>0
 }
+//_______________________________________________________________________________________
+// void AliEMCALSDigitizer::TestTowerID(void)
+// {
+//   Int_t j;
+
+//   Bool_t preshower = kFALSE;
+//   for (j = 0 ; j < 10 ; j++){  // loop over hit id
+//     Int_t i;
+//    for (i = 0 ; i <= 2 ; i++){  // loop over 
+//      Int_t k = i*96*144+j*144+1;
+//       cout << " Hit Index = " << k << "   " << j*10 << "   TOWERID = " <<  Layer2TowerID(k, preshower) << endl ;
+//     }
+//   }
+// }