]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALRecPoint.cxx
Additional protection in Digitize, which was moved to the implementation file
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecPoint.cxx
index 763294077f75d1e3433ef456cb742e649252b6f1..e1940d602892bc9bffc291997ba6e98302ba0310 100644 (file)
 
 // --- AliRoot header files ---
 #include "AliGenerator.h"
+#include "AliRunLoader.h"
+#include "AliRun.h"
+#include "AliEMCAL.h"
+#include "AliEMCALLoader.h"
 #include "AliEMCALGeometry.h"
+#include "AliEMCALHit.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALRecPoint.h"
 
@@ -57,8 +62,11 @@ AliEMCALRecPoint::AliEMCALRecPoint()
   fTime = 0. ;
   //  fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
   fCoreRadius = 10;        //HG Check this
-  fGeom = AliEMCALGeometry::GetInstance();
-  fGeom->GetTransformationForSM(); // Global <-> Local
+
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  //fGeomPtr = AliEMCALGeometry::GetInstance();
+  fGeomPtr->GetTransformationForSM(); // Global <-> Local
 }
 
 //____________________________________________________________________________
@@ -79,8 +87,10 @@ AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) : AliRecPoint(opt)
   fTime = -1. ;
   //fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
   fCoreRadius = 10;        //HG Check this
-  fGeom = AliEMCALGeometry::GetInstance();
-  fGeom->GetTransformationForSM(); // Global <-> Local
+  //fGeomPtr = AliEMCALGeometry::GetInstance();
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  fGeomPtr->GetTransformationForSM(); // Global <-> Local
 }
 //____________________________________________________________________________
 AliEMCALRecPoint::~AliEMCALRecPoint()
@@ -108,7 +118,7 @@ void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
     fTimeList =  new Float_t[fMaxDigit]; 
   if(fAbsIdList == 0) {
     fAbsIdList =  new Int_t[fMaxDigit];
-    fSuperModuleNumber = fGeom->GetSuperModuleNumber(digit.GetId());
+    fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId());
   }
 
   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
@@ -173,11 +183,11 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d
 
   areNeighbours = kFALSE ;
 
-  fGeom->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta);
-  fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, relid1[0],relid1[1]);
+  fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta);
+  fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, relid1[0],relid1[1]);
 
-  fGeom->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
-  fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, relid2[0],relid2[1]);
+  fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
+  fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, relid2[0],relid2[1]);
   
   rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
   coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
@@ -205,7 +215,7 @@ Int_t AliEMCALRecPoint::Compare(const TObject * obj) const
   TVector3 locpos2;  
   clu->GetLocalPosition(locpos2);  
 
-  Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
+  Int_t rowdif = (Int_t)(TMath::Ceil(locpos1.X()/delta)-TMath::Ceil(locpos2.X()/delta)) ;
   if (rowdif> 0) 
     rv = 1 ;
   else if(rowdif < 0) 
@@ -337,7 +347,7 @@ void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
   //  printf("eval axis done\n");
   EvalDispersion(logWeight, digits) ;
   //  printf("eval dispersion done\n");
// EvalCoreEnergy(logWeight, digits);
 //EvalCoreEnergy(logWeight, digits);
  // printf("eval energy done\n");
   EvalTime(digits) ;
   //  printf("eval time done\n");
@@ -365,7 +375,7 @@ void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
 
-    fGeom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
+    fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
     w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 
     if(w>0.0) {
@@ -397,7 +407,7 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit
   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
 
-    fGeom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
+    fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
     // printf(" Id %i : Local x,y,z %f %f %f \n", digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
 
     if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
@@ -448,7 +458,7 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit
 void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
 {
   // This function calculates energy in the core, 
-  // i.e. within a radius rad = 3cm around the center. Beyond this radius
+  // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius
   // in accordance with shower profile the energy deposition 
   // should be less than 2%
 
@@ -463,12 +473,14 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
   
   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
-    Float_t etai = 0. ;
-    Float_t phii = 0. ;
-    fGeom->PosInAlice(digit->GetId(), etai, phii);
-    phii = phii * kDeg2Rad;
+    
+    Float_t ieta = 0.0;
+    Float_t iphi = 0.0;
+    //fGeomPtr->PosInAlice(digit->GetId(), ieta, iphi);
+    fGeomPtr->EtaPhiFromIndex(digit->GetId(),ieta, iphi) ;
+    iphi = iphi * kDeg2Rad;
   
-    Float_t distance = TMath::Sqrt((etai-fLocPos.X())*(etai-fLocPos.X())+(phii-fLocPos.Y())*(phii-fLocPos.Y())) ;
+    Float_t distance = TMath::Sqrt((ieta-fLocPos.X())*(ieta-fLocPos.X())+(iphi-fLocPos.Y())*(iphi-fLocPos.Y())) ;
     if(distance < fCoreRadius)
       fCoreEnergy += fEnergyList[iDigit] ;
   }
@@ -489,7 +501,7 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
   const Float_t kDeg2Rad = TMath::DegToRad();
   AliEMCALDigit * digit ;
 
-  TString gn(fGeom->GetName());
+  TString gn(fGeomPtr->GetName());
 
   Int_t iDigit;
   
@@ -502,12 +514,12 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
     // for this geometry does not exist
       int nSupMod=0, nTower=0, nIphi=0, nIeta=0;
       int iphi=0, ieta=0;
-      fGeom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
-      fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
+      fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
+      fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
       etai=(Float_t)ieta;
       phii=(Float_t)iphi;
     } else {
-      fGeom->EtaPhiFromIndex(digit->GetId(), etai, phii);
+      fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii);
       phii = phii * kDeg2Rad;
     }
 
@@ -673,7 +685,7 @@ void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
   // returns the position of the cluster in the global reference system of ALICE
   // These are now the Cartesian X, Y and Z
   //  cout<<" geom "<<geom<<endl;
-  fGeom->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
+  fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
 }
 
 //____________________________________________________________________________
@@ -759,6 +771,66 @@ Int_t  AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit **  maxAt, Float_t *
   }
   return iDigitN ;
 }
+
+//____________________________________________________________________________
+Int_t AliEMCALRecPoint::GetPrimaryIndex() const  
+{
+  // Get the primary track index in TreeK which deposits the most energy 
+  // in Digits which forms RecPoint. Kinematics, Hits and Digits must be 
+  // loaded before the call of the method.
+
+  AliRunLoader *rl = AliRunLoader::GetRunLoader(); 
+  if (!rl) 
+    AliError(Form(" No Runloader ")) ; 
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
+    (rl->GetDetectorLoader("EMCAL"));
+
+  // Get the list of digits forming this RecPoint
+  Int_t  nDigits   = fMulDigit   ;
+  Int_t *digitList = fDigitsList ;
+  
+  // Find the digit with maximum amplitude
+  AliEMCALDigit *digit = 0;
+  TClonesArray *digits = emcalLoader->Digits();
+  Int_t maxAmp = 0;
+  Int_t bestDigitIndex = -1;
+  for (Int_t iDigit=0; iDigit<nDigits; iDigit++) {
+    digit = static_cast<AliEMCALDigit *>(digits->At(digitList[iDigit]));
+    if (digit->GetAmp() > maxAmp) {
+      maxAmp = digit->GetAmp();
+      bestDigitIndex = iDigit;
+    }
+  }
+
+  digit = static_cast<AliEMCALDigit *>(digits->At(digitList[bestDigitIndex]));  
+
+  // Get the list of hits producing this digit,
+  // find which hit has deposited more energy 
+  // and find the primary track.
+
+  AliEMCALHit *hit = 0;
+  TClonesArray *hits = emcalLoader->Hits();
+
+  Double_t maxedep  =  0;
+  Int_t    maxtrack = -1;
+  Int_t    nHits    = hits ->GetEntries();
+  Int_t    id       = digit->GetId();
+  for (Int_t iHit=0; iHit<nHits; iHit++) {
+    hit = static_cast<AliEMCALHit*> (hits->At(iHit)) ;
+    if(hit->GetId() == id){
+      Double_t edep  = hit->GetEnergy();
+      Int_t    track = hit->GetIparent();//Primary();
+      if(edep > maxedep){
+       maxedep  = edep;
+       maxtrack = track;
+      }
+    }
+  }
+  if (maxtrack != -1) return maxtrack; 
+  return -12345;                       // no track found :(
+}
+
 //____________________________________________________________________________
 void AliEMCALRecPoint::EvalTime(TClonesArray * digits){
   // time is set to the time of the digit with the maximum energy