// --- 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"
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
}
//____________________________________________________________________________
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()
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
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] ) ;
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)
// 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");
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) {
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 ));
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%
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] ;
}
const Float_t kDeg2Rad = TMath::DegToRad();
AliEMCALDigit * digit ;
- TString gn(fGeom->GetName());
+ TString gn(fGeomPtr->GetName());
Int_t iDigit;
// 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;
}
// 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);
}
//____________________________________________________________________________
}
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