// EMCAL includes
#include "AliEMCALRecoUtils.h"
#include "AliEMCALGeometry.h"
-#include "AliEMCALTrack.h"
+#include "AliTrackerBase.h"
#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
#include "AliEMCALPIDUtils.h"
+
ClassImp(AliEMCALRecoUtils)
//______________________________________________
fUseRunCorrectionFactors(kFALSE), fRunCorrectionFactorsSet(kFALSE),
fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
- fRejectExoticCluster(kFALSE), fPIDUtils(), fAODFilterMask(32),
+ fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
+ fExoticCellFraction(0.97), fExoticCellDiffTime(1e6), fExoticCellMinAmplitude(0.5),
+ fPIDUtils(), fAODFilterMask(32),
fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kTRUE), fCutEtaPhiSeparate(kFALSE),
- fCutR(0.1), fCutEta(0.025), fCutPhi(0.05),
- fMass(0.139), fStep(10),
+ fCutR(0.05), fCutEta(0.025), fCutPhi(0.05),
+ fClusterWindow(100), fMass(0.139),
+ fStepSurface(20.), fStepCluster(5.),
fTrackCutsType(kLooseCut), fCutMinTrackPt(0), fCutMinNClusterTPC(-1),
fCutMinNClusterITS(-1), fCutMaxChi2PerClusterTPC(1e10), fCutMaxChi2PerClusterITS(1e10),
fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
- fRejectExoticCluster(reco.fRejectExoticCluster), fPIDUtils(reco.fPIDUtils),
- fAODFilterMask(reco.fAODFilterMask),
+ fRejectExoticCluster(reco.fRejectExoticCluster), fRejectExoticCells(reco.fRejectExoticCells),
+ fExoticCellFraction(reco.fExoticCellFraction), fExoticCellDiffTime(reco.fExoticCellDiffTime),
+ fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
+ fPIDUtils(reco.fPIDUtils), fAODFilterMask(reco.fAODFilterMask),
fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
- fMass(reco.fMass), fStep(reco.fStep),
+ fClusterWindow(reco.fClusterWindow),
+ fMass(reco.fMass), fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
+
fRejectExoticCluster = reco.fRejectExoticCluster;
-
+ fRejectExoticCells = reco.fRejectExoticCells;
+ fExoticCellFraction = reco.fExoticCellFraction;
+ fExoticCellDiffTime = reco.fExoticCellDiffTime;
+ fExoticCellMinAmplitude = reco.fExoticCellMinAmplitude;
+
fPIDUtils = reco.fPIDUtils;
fAODFilterMask = reco.fAODFilterMask;
fCutR = reco.fCutR;
fCutEta = reco.fCutEta;
fCutPhi = reco.fCutPhi;
+ fClusterWindow = reco.fClusterWindow;
fMass = reco.fMass;
- fStep = reco.fStep;
- fRejectExoticCluster = reco.fRejectExoticCluster;
+ fStepSurface = reco.fStepSurface;
+ fStepCluster = reco.fStepCluster;
fTrackCutsType = reco.fTrackCutsType;
fCutMinTrackPt = reco.fCutMinTrackPt;
}
-//__________________________________________________
+//_____________________________________
AliEMCALRecoUtils::~AliEMCALRecoUtils()
{
//Destructor.
InitTrackCuts();
}
+//_______________________________________________________________________________
+Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
+ Float_t & amp, Double_t & time,
+ AliVCaloCells* cells)
+{
+ // Reject cell if criteria not passed and calibrate it
+
+ AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
+
+ if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE;
+
+ Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
+ geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
+ geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
+
+ // Do not include bad channels found in analysis,
+ if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) {
+ return kFALSE;
+ }
+
+ //Recalibrate energy
+ amp = cells->GetCellAmplitude(absID);
+ if(IsRecalibrationOn())
+ amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+
+
+ // Recalibrate time
+ time = cells->GetCellTime(absID);
+
+ RecalibrateCellTime(absID,bc,time);
+
+ return kTRUE;
+}
+
//_______________________________________________________________
Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
{
}
+//_______________________________________________________________________
+Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
+{
+ // Look to cell neighbourhood and reject if it seems exotic
+ // Do before recalibrating the cells
+ if(!fRejectExoticCells) return kFALSE;
+
+ AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+
+ Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
+ geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
+ geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
+
+ //Get close cells index, energy and time, not in corners
+ Int_t absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
+ Int_t absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
+ Int_t absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
+ Int_t absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
+
+ Float_t ecell = 0, ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
+ Double_t tcell = 0, tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
+ Bool_t accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
+
+ accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
+
+ if(!accept) return kTRUE; // reject this cell
+
+ if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
+
+ accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
+ accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
+ accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
+ accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
+
+ /*
+ printf("Cell absID %d \n",absID);
+ printf("\t accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
+ accept1,accept2,accept3,accept4);
+ printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
+ absID,absID1,absID2,absID3,absID4);
+ printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
+ ecell,ecell1,ecell2,ecell3,ecell4);
+ printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
+ tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
+ TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
+ */
+
+ if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
+ if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
+ if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
+ if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
+
+ Float_t eCross = ecell1+ecell2+ecell3+ecell4;
+
+ //printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell);
+
+ if(1-eCross/ecell > fExoticCellFraction) {
+ AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",absID,ecell,eCross,1-eCross/ecell));
+ return kTRUE;
+ }
+
+ return kFALSE;
+
+}
+
//_________________________________________________
-Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
- // Check if the cluster has high energy but small number of cells
- // The criteria comes from Gustavo's study
- //
+Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster, AliVCaloCells *cells, const Int_t bc) {
+ // Check if the cluster highest energy tower is exotic
if(!cluster){
AliInfo("Cluster pointer null!");
return kFALSE;
}
- Int_t nc = cluster->GetNCells() ;
+ if(!fRejectExoticCluster) return kFALSE;
- if ( nc > 8 ) return kFALSE ; // Good cluster, needed for 3x3 clusterizer
- else if ( nc < 1 + cluster->E()/3. ) return kTRUE ; // Bad cluster
- else return kFALSE ; // Good cluster
+ // Get highest energy tower
+ AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
+ Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
+ Bool_t shared = kFALSE;
+ GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
+
+ return IsExoticCell(absId,cells,bc);
}
}
Float_t energy = cluster->E();
-
+
switch (fNonLinearityFunction) {
case kPi0MC:
break;
}
-
+
return energy;
}
}
//________________________________________________________________
-void AliEMCALRecoUtils::RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells * cells, Int_t bc){
+void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc){
// Recalibrate the cells time and energy, considering the recalibration map and the energy
// of the cells that compose the cluster.
// bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
- if(!IsRecalibrationOn()) return;
+ if(!IsRecalibrationOn() && !IsTimeRecalibrationOn()) return;
if(!cells){
AliInfo("Cells pointer null!");
fCellsRecalibrated = kTRUE;
- Int_t absId =-1;
- Int_t icol =-1, irow =-1, imod = 1;
- Int_t iTower =-1, iIeta =-1, iIphi =-1;
-
- Int_t nEMcell = cells->GetNumberOfCells() ;
+ Int_t absId =-1;
+ Bool_t accept = kFALSE;
+ Float_t ecell = 0;
+ Double_t tcell = 0;
+ Int_t nEMcell = cells->GetNumberOfCells() ;
for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
- absId = cells->GetCellNumber(iCell);
-
- // Energy
- Float_t factor = 1;
- if(IsRecalibrationOn()){
- geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
- if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
- geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
- factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
- }
-
- Float_t cellE = cells->GetAmplitude(iCell) * factor ;
+ absId = cells->GetCellNumber(iCell);
- //Time
- Double_t celltime = cells->GetCellTime(absId);
- RecalibrateCellTime(absId, bc, celltime);
+ accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
+ if(!accept) {
+ ecell = 0;
+ tcell = 0;
+ }
//Set new values
- cells->SetCell(iCell,cells->GetCellNumber(iCell),cellE, celltime);
+ cells->SetCell(iCell,absId,ecell, tcell);
}
}
-//________________________________________________________________
-void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime){
+//_________________________________________________________________________________________________
+void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime)
+{
// Recalibrate time of cell with absID considering the recalibration map
// bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
+
+ if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0){
- if(!fCellsRecalibrated && IsTimeRecalibrationOn()){
-// printf("cell time org %g, ",celltime);
-
- Double_t timeBCoffset = 0.;
- if( bc%4 ==0 || bc%4==1) timeBCoffset = 100.*1.e-9; //in ns
-
- Double_t celloffset = GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;
-
-// printf("absId %d, time %f bc %d-%d: bc0 %f, bc1 %f, bc2 %f, bc3 %f \n", absId, celltime*1.e9,bc, bc%4,
-// GetEMCALChannelTimeRecalibrationFactor(0,absId),GetEMCALChannelTimeRecalibrationFactor(1,absId),
-// GetEMCALChannelTimeRecalibrationFactor(2,absId),GetEMCALChannelTimeRecalibrationFactor(3,absId));
+ celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; ;
- celltime -= timeBCoffset ;
- celltime -= celloffset ;
-// printf("new %g\n",celltime);
}
}
AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+ TObjArray *clusterArray = 0x0;
+ if(!clusterArr)
+ {
+ clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
+ for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
+ {
+ AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
+ if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
+ clusterArray->AddAt(cluster,icl);
+ }
+ }
+
Int_t matched=0;
Double_t cv[21];
for (Int_t i=0; i<21;i++) cv[i]=0;
AliExternalTrackParam *trackParam = 0;
//If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
+ AliESDtrack *esdTrack = 0;
+ AliAODTrack *aodTrack = 0;
if(esdevent)
{
- AliESDtrack *esdTrack = esdevent->GetTrack(itr);
- if(!esdTrack || !IsAccepted(esdTrack)) continue;
+ esdTrack = esdevent->GetTrack(itr);
+ if(!esdTrack) continue;
+ if(!IsAccepted(esdTrack)) continue;
if(esdTrack->Pt()<fCutMinTrackPt) continue;
+ Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
+ if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
}
//If the input event is AOD, the starting point for extrapolation is at vertex
- //AOD tracks are selected according to its bit.
+ //AOD tracks are selected according to its filterbit.
else if(aodevent)
{
- AliAODTrack *aodTrack = aodevent->GetTrack(itr);
+ aodTrack = aodevent->GetTrack(itr);
if(!aodTrack) continue;
if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
if(aodTrack->Pt()<fCutMinTrackPt) continue;
+ Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
+ if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
Double_t pos[3],mom[3];
aodTrack->GetXYZ(pos);
aodTrack->GetPxPyPz(mom);
else
{
printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
+ if(clusterArray)
+ {
+ clusterArray->Clear();
+ delete clusterArray;
+ }
return;
}
if(!trackParam) continue;
-
- Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
+
+ //Extrapolate the track to EMCal surface
+ AliExternalTrackParam emcalParam(*trackParam);
+ Float_t eta, phi;
+ if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi))
+ {
+ if(aodevent && trackParam) delete trackParam;
+ continue;
+ }
+
+ if(esdevent)
+ {
+ esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
+ }
+
+ if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
+ {
+ if(aodevent && trackParam) delete trackParam;
+ continue;
+ }
+
+
+ //Find matched clusters
Int_t index = -1;
- if(!clusterArr){// get clusters from event
- for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
+ Float_t dEta = -999, dPhi = -999;
+ if(!clusterArr)
{
- AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
- if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
- AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation
- Float_t tmpEta=-999, tmpPhi=-999;
- if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
- if(fCutEtaPhiSum)
- {
- Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
- if(tmpR<dRMax)
- {
- dRMax=tmpR;
- dEtaMax=tmpEta;
- dPhiMax=tmpPhi;
- index=icl;
- }
- }
- else if(fCutEtaPhiSeparate)
- {
- if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
- {
- dEtaMax = tmpEta;
- dPhiMax = tmpPhi;
- index=icl;
- }
- }
- else
- {
- printf("Error: please specify your cut criteria\n");
- printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
- printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
- if(aodevent && trackParam) delete trackParam;
- return;
- }
- }//cluster loop
- }
- else { // external cluster array, not from ESD event
- for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
+ index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
+ }
+ else
{
- AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
- if(!cluster){
- AliInfo("Cluster not found!!!");
- continue;
- }
- if(!cluster->IsEMCAL()) continue;
- AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
- Float_t tmpEta=-999, tmpPhi=-999;
- if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
- if(fCutEtaPhiSum)
- {
- Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
- if(tmpR<dRMax)
- {
- dRMax=tmpR;
- dEtaMax=tmpEta;
- dPhiMax=tmpPhi;
- index=icl;
- }
- }
- else if(fCutEtaPhiSeparate)
- {
- if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
- {
- dEtaMax = tmpEta;
- dPhiMax = tmpPhi;
- index=icl;
- }
- }
- else
- {
- printf("Error: please specify your cut criteria\n");
- printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
- printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
- if(aodevent && trackParam) delete trackParam;
- return;
- }
- }//cluster loop
- }// external list of clusters
+ index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
+ }
if(index>-1)
{
- fMatchedTrackIndex ->AddAt(itr,matched);
- fMatchedClusterIndex->AddAt(index,matched);
- fResidualEta ->AddAt(dEtaMax,matched);
- fResidualPhi ->AddAt(dPhiMax,matched);
+ fMatchedTrackIndex ->AddAt(itr,matched);
+ fMatchedClusterIndex ->AddAt(index,matched);
+ fResidualEta ->AddAt(dEta,matched);
+ fResidualPhi ->AddAt(dPhi,matched);
matched++;
}
if(aodevent && trackParam) delete trackParam;
}//track loop
+
+ if(clusterArray)
+ {
+ clusterArray->Clear();
+ delete clusterArray;
+ }
AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
- fMatchedTrackIndex ->Set(matched);
- fMatchedClusterIndex->Set(matched);
- fResidualPhi ->Set(matched);
- fResidualEta ->Set(matched);
+ fMatchedTrackIndex ->Set(matched);
+ fMatchedClusterIndex ->Set(matched);
+ fResidualPhi ->Set(matched);
+ fResidualEta ->Set(matched);
}
//________________________________________________________________________________
-Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
+Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi)
{
//
// This function returns the index of matched cluster to input track
// Returns -1 if no match is found
-
-
- Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
Int_t index = -1;
-
+ Double_t phiV = track->Phi()*TMath::RadToDeg();
+ if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
-
- if(!trackParam) return index;
+ if(!trackParam) return index;
+ AliExternalTrackParam emcalParam(*trackParam);
+ Float_t eta, phi;
+ if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) return index;
+ if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
+
+ TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
+
for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
{
AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
- if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
- AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
- Float_t tmpEta=-999, tmpPhi=-999;
- if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
- if(fCutEtaPhiSum)
+ if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
+ clusterArr->AddAt(cluster,icl);
+ }
+
+ index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
+ clusterArr->Clear();
+ delete clusterArr;
+
+ return index;
+}
+
+//________________________________________________________________________________
+Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, TObjArray * clusterArr, Float_t &dEta, Float_t &dPhi)
+{
+ dEta=-999, dPhi=-999;
+ Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
+ Int_t index = -1;
+ Float_t tmpEta=-999, tmpPhi=-999;
+
+ Double_t exPos[3] = {0.,0.,0.};
+ if(!emcalParam->GetXYZ(exPos)) return index;
+
+ Float_t clsPos[3] = {0.,0.,0.};
+ for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
{
- Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
- if(tmpR<dRMax)
+ AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
+ if(!cluster || !cluster->IsEMCAL()) continue;
+ cluster->GetPosition(clsPos);
+ Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
+ if(dR > fClusterWindow) continue;
+
+ AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
+ if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
+ if(fCutEtaPhiSum)
+ {
+ Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
+ if(tmpR<dRMax)
{
dRMax=tmpR;
dEtaMax=tmpEta;
dPhiMax=tmpPhi;
index=icl;
}
- }
- else if(fCutEtaPhiSeparate)
- {
- if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
+ }
+ else if(fCutEtaPhiSeparate)
+ {
+ if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
{
dEtaMax = tmpEta;
dPhiMax = tmpPhi;
index=icl;
}
+ }
+ else
+ {
+ printf("Error: please specify your cut criteria\n");
+ printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
+ printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
+ return index;
+ }
}
- else
- {
- printf("Error: please specify your cut criteria\n");
- printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
- printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
- return -1;
- }
- }//cluster loop
+
+ dEta=dEtaMax;
+ dPhi=dPhiMax;
+
return index;
}
-//________________________________________________________________________________
-Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, Double_t emcalR, Double_t mass, Double_t step, Float_t &eta, Float_t &phi)
+{
+ eta = -999, phi = -999;
+ if(!trkParam) return kFALSE;
+ if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
+ Double_t trkPos[3] = {0.,0.,0.};
+ if(!trkParam->GetXYZ(trkPos)) return kFALSE;
+ TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+ eta = trkPosVec.Eta();
+ phi = trkPosVec.Phi();
+ if(phi<0)
+ phi += 2*TMath::Pi();
+
+ return kTRUE;
+}
+
+
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, Float_t *clsPos, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
{
//
- //Return the residual by extrapolating a track to a cluster
+ //Return the residual by extrapolating a track param to a global position
//
- if(!trkParam || !cluster) return kFALSE;
- Float_t clsPos[3];
- Double_t trkPos[3];
- cluster->GetPosition(clsPos); //Has been recalculated
+ tmpEta = -999;
+ tmpPhi = -999;
+ if(!trkParam) return kFALSE;
+ Double_t trkPos[3] = {0.,0.,0.};
TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
- if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.8, -1)) return kFALSE;
- trkParam->GetXYZ(trkPos); //Get the extrapolated global position
+ if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
+ if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
// track cluster matching
- tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
- tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
+ tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
+ tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
return kTRUE;
}
+
+//
+//------------------------------------------------------------------------------
+Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
+{
+ //
+ //Return the residual by extrapolating a track param to a cluster
+ //
+ tmpEta = -999;
+ tmpPhi = -999;
+ if(!cluster || !trkParam) return kFALSE;
+
+ Float_t clsPos[3] = {0.,0.,0.};
+ cluster->GetPosition(clsPos);
+
+ return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
+}
+
+//
+//------------------------------------------------------------------------------
+Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
+{
+ //
+ //Return the residual by extrapolating a track param to a clusterfStepCluster
+ //
+
+ return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
+}
+
+
//________________________________________________________________________________
void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
{
return pos;
}
-//__________________________________________________________
-Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
+//___________________________________________________________________________________
+Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom,
+ AliVCaloCells* cells,const Int_t bc)
{
// check if the cluster survives some quality cut
//
//
Bool_t isGood=kTRUE;
- if(!cluster || !cluster->IsEMCAL()) return kFALSE;
+
+ if(!cluster || !cluster->IsEMCAL()) return kFALSE;
+
if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
+
if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
- if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
+
+ if(IsExoticCluster(cluster, cells,bc)) return kFALSE;
return isGood;
}
}
}
+
+//__________________________________________________________________
+void AliEMCALRecoUtils::SetClusterMatchedToTrack(AliESDEvent *event)
+{
+ // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
+
+ Int_t nTracks = event->GetNumberOfTracks();
+ for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
+ AliESDtrack* track = event->GetTrack(iTrack);
+ if (!track) {
+ AliWarning(Form("Could not receive track %d", iTrack));
+ continue;
+ }
+ Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
+ track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
+ if(matchClusIndex != -1)
+ track->SetStatus(AliESDtrack::kEMCALmatch);
+ else
+ track->ResetStatus(AliESDtrack::kEMCALmatch);
+ }
+ AliDebug(2,"Track matched to closest cluster");
+}
+
+//_________________________________________________________________
+void AliEMCALRecoUtils::SetTracksMatchedToCluster(AliESDEvent *event)
+{
+ // Checks if EMC clusters are matched to ESD track.
+ // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
+
+ for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) {
+ AliESDCaloCluster *cluster = event->GetCaloCluster(iClus);
+ if (!cluster->IsEMCAL())
+ continue;
+
+ Int_t nTracks = event->GetNumberOfTracks();
+ TArrayI arrayTrackMatched(nTracks);
+
+ // Get the closest track matched to the cluster
+ Int_t nMatched = 0;
+ Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
+ if (matchTrackIndex != -1) {
+ arrayTrackMatched[nMatched] = matchTrackIndex;
+ nMatched++;
+ }
+
+ // Get all other tracks matched to the cluster
+ for(Int_t iTrk=0; iTrk<nTracks; ++iTrk) {
+ AliESDtrack* track = event->GetTrack(iTrk);
+ if(iTrk == matchTrackIndex) continue;
+ if(track->GetEMCALcluster() == iClus){
+ arrayTrackMatched[nMatched] = iTrk;
+ ++nMatched;
+ }
+ }
+
+ //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
+
+ arrayTrackMatched.Set(nMatched);
+ cluster->AddTracksMatched(arrayTrackMatched);
+
+ Float_t eta= -999, phi = -999;
+ if (matchTrackIndex != -1)
+ GetMatchedResiduals(iClus, eta, phi);
+ cluster->SetTrackDistance(phi, eta);
+ }
+
+ AliDebug(2,"Cluster matched to tracks");
+}
+
+
//___________________________________________________
void AliEMCALRecoUtils::Print(const Option_t *) const
{
printf("Matching criteria: ");
if(fCutEtaPhiSum)
{
- printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
+ printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
}
else if(fCutEtaPhiSeparate)
{
- printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
+ printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
}
else
{
printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
}
- printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
+ printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
+ printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
printf("Track cuts: \n");
printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);