1 //_________________________________________________________________________
2 // Utility Class for transverse energy studies
3 // Selection class for EMCAL
5 //*-- Authors: Oystein Djuvsland (Bergen)
6 //_________________________________________________________________________
9 #include "AliAnalysisEtSelectorEmcal.h"
10 #include "AliAnalysisEtCuts.h"
11 #include "AliEMCALTrack.h"
12 #include "TParticle.h"
16 AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal(AliAnalysisEtCuts* cuts):AliAnalysisEtSelector(cuts)
20 AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal():AliAnalysisEtSelector()
24 AliAnalysisEtSelectorEmcal::~AliAnalysisEtSelectorEmcal()
29 void AliAnalysisEtSelectorEmcal::Init()
31 AliAnalysisEtSelector::Init();
34 Int_t AliAnalysisEtSelectorEmcal::Init(const AliESDEvent* event)
37 AliAnalysisEtSelector::Init(event);
38 Printf("Initializing selector for run: %d", event->GetRunNumber());
43 TRefArray* AliAnalysisEtSelectorEmcal::GetClusters()
46 if(!fClusterArray) fClusterArray = new TRefArray;
50 fEvent->GetEMCALClusters(fClusterArray);
54 Printf("Could not initialize cluster array");
59 Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const AliESDCaloCluster& cl) const
64 return TMath::Sin(cp.Theta())*cl.E() > fCuts->GetReconstructedEmcalClusterEnergyCut();
67 Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const TParticle& p) const
69 return TMath::Sin(p.Theta())*p.Energy() > fCuts->GetReconstructedEmcalClusterEnergyCut();
71 Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(Double_t e) const
73 return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
76 Bool_t AliAnalysisEtSelectorEmcal::PassDistanceToBadChannelCut(const AliESDCaloCluster& ) const
81 Bool_t AliAnalysisEtSelectorEmcal::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
84 Int_t nTracksMatched = cluster.GetNTracksMatched();
85 if(nTracksMatched == 0){
89 Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
90 if(trackMatchedIndex < 0){
93 AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
94 if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random. It takes ~460 MeV to reach the EMCal
96 //Float_t recoE = cluster.E();
99 // //cluster.GetPosition(pos);
100 // Int_t trackMatchIdx = cluster.GetTrackMatchedIndex();
101 // //Double_t distance = 9999.0;
102 // if(trackMatchIdx>-1)
105 // //distance = CalcTrackClusterDistance(pos, &trackMatchIdx);
110 //return distance > fCuts->GetEmcalTrackDistanceCut();
113 Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const TParticle& part)
115 float myphi = part.Phi();
116 myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
117 return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
118 && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
119 && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
122 Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliVTrack& part)
124 float myphi = part.Phi();
125 myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
126 return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
127 && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
128 && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
131 Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
134 cluster.GetPosition(pos);
136 float myphi = cp.Phi();
137 myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
138 return TMath::Abs(cp.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
139 && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
140 && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
145 AliAnalysisEtSelectorEmcal::CalcTrackClusterDistance(const Float_t clsPos[3],Int_t *trkMatchId) const
146 { // calculate distance between cluster and closest track
148 Double_t trkPos[3] = {0,0,0};
150 Int_t bestTrkMatchId = -1;
151 Double_t distance = 9999; // init to a big number
154 Double_t distX = 0, distY = 0, distZ = 0;
156 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
157 AliESDtrack *track = static_cast<AliESDtrack*>( fEvent->GetTrack(iTrack) );
159 AliError(Form("ERROR: Could not get track %d", iTrack));
163 // check for approx. eta and phi range before we propagate..
166 AliEMCALTrack emctrack(*track);
167 if (!emctrack.PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
170 emctrack.GetXYZ(trkPos);
172 distX = clsPos[0]-trkPos[0];
173 distY = clsPos[1]-trkPos[1];
174 distZ = clsPos[2]-trkPos[2];
175 dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);
177 if (dist < distance) {
179 bestTrkMatchId = iTrack;
183 // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
184 *trkMatchId = bestTrkMatchId;
188 UInt_t AliAnalysisEtSelectorEmcal::GetLabel(const AliESDCaloCluster *cluster, AliStack *stack){
189 UInt_t correctLabel = TMath::Abs(cluster->GetLabel());
190 correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
192 //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
193 // if(!stack){return 0;}
194 // if( correctLabel > 0 && stack->IsSecondaryFromMaterial(correctLabel)){//if this is flagged as a secondary then we look to see where it really came from
195 // TParticle *hitParticle = stack->Particle(correctLabel);
197 // Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >420);
198 // if(partVtxSecondary){//at this point we have something which converted near the detector. Let's find the mother particle
199 // UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
201 // TParticle *mother = stack->Particle(mothIdx);
203 // if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
204 // if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
205 // cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
206 // cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
207 // //so now we know that the particle originated near the cover and within the acceptance of the detector
217 return correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?