1 /*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: A.Abrahantes, E.Lopez, S.Vallero *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 #include <TObjArray.h>
22 #include "AliAnalyseLeadingTrackUE.h"
24 #include "AliAODEvent.h"
25 #include "AliAODMCParticle.h"
26 #include "AliAODTrack.h"
27 #include "AliESDEvent.h"
28 #include "AliESDtrack.h"
29 #include "AliESDtrackCuts.h"
30 #include "AliInputEventHandler.h"
31 #include "AliMCEvent.h"
32 #include "AliVParticle.h"
33 #include "AliAODMCHeader.h"
36 #include "AliAnalysisManager.h"
37 #include "AliMCEventHandler.h"
39 #include "AliPIDResponse.h"
40 #include "AliHelperPID.h"
43 ////////////////////////////////////////////////
44 //---------------------------------------------
45 // Class for transverse regions analysis
46 //---------------------------------------------
47 ////////////////////////////////////////////////
52 ClassImp(AliAnalyseLeadingTrackUE)
54 //-------------------------------------------------------------------
55 AliAnalyseLeadingTrackUE::AliAnalyseLeadingTrackUE() :
61 fCheckMotherPDG(kTRUE),
64 fTrackPhiCutEvPlMin(0.),
65 fTrackPhiCutEvPlMax(0.),
67 fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
69 fSharedClusterCut(-1),
71 fFoundFractionCut(-1),
73 fEsdTrackCutsExtra1(0x0),
74 fEsdTrackCutsExtra2(0x0),
82 //-------------------------------------------------------------------
83 AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
85 // assignment operator
90 //-------------------------------------------------------------------
91 AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
99 //____________________________________________________________________
100 Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
103 // select track according to set of cuts
104 if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
105 if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->IsSelected(track)) return kFALSE;
111 //____________________________________________________________________
112 void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
114 // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
115 // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
118 filterbit = fFilterBit;
120 if (filterbit == 128)
122 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
123 fEsdTrackCuts->SetMinNClustersTPC(70);
125 else if (filterbit == 256)
128 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
129 fEsdTrackCuts->SetMinNClustersTPC(80);
130 fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
131 fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
132 fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
134 else if (filterbit == 512)
137 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
138 fEsdTrackCuts->SetMinNClustersTPC(60);
139 fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
140 fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
141 fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
143 else if (filterbit == 1024)
145 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
146 fEsdTrackCuts->SetMinNClustersTPC(-1);
147 fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
148 fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
150 else if (filterbit == 2048) // mimic hybrid tracks
152 // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
153 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
154 fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
155 fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
156 fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
157 fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
158 fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
159 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
160 fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
162 // Add SPD requirement
163 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
164 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
165 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
167 fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
168 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
169 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
171 else if (filterbit == 4096) // require TOF matching
173 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
174 // FIXME: TOF REQUIREMENTS ARE IN GetParticleSpecies FOR THE MOMENT
178 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
179 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
181 // Add SPD requirement
182 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
183 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
185 // Add SDD requirement
186 fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
187 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
191 //____________________________________________________________________
192 TObjArray* AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
195 // Returns an array of charged particles (or jets) ordered according to their pT.
197 Int_t nTracks = NParticles(obj);
200 if( !nTracks ) return 0;
202 // Define array of AliVParticle objects
203 TObjArray* tracks = new TObjArray(nTracks);
205 // Loop over tracks or jets
206 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
207 AliVParticle* part = ParticleWithCuts( obj, ipart );
209 // Accept leading-tracks in a limited pseudo-rapidity range
210 if( TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin ) continue;
211 tracks->AddLast( part );
213 // Order tracks by pT
214 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
216 nTracks = tracks->GetEntriesFast();
217 if( !nTracks ) return 0;
223 void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
225 // remove injected signals (primaries above <maxLabel>)
226 // <tracks> can be the following cases:
227 // a. tracks: in this case the label is taken and then case b.
228 // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
229 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
231 TClonesArray* arrayMC = 0;
232 AliMCEvent* mcEvent = 0;
233 if (mcObj->InheritsFrom("AliMCEvent"))
234 mcEvent = static_cast<AliMCEvent*>(mcObj);
235 else if (mcObj->InheritsFrom("TClonesArray"))
236 arrayMC = static_cast<TClonesArray*>(mcObj);
240 AliFatal("Invalid object passed");
243 Int_t before = tracks->GetEntriesFast();
245 for (Int_t i=0; i<before; ++i)
247 AliVParticle* part = (AliVParticle*) tracks->At(i);
249 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
250 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
252 AliVParticle* mother = part;
255 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
257 if (((AliMCParticle*)mother)->GetMother() < 0)
263 mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
270 // find the primary mother
271 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
273 if (((AliAODMCParticle*)mother)->GetMother() < 0)
279 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
287 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
291 // Printf("%d %d %d", i, part->GetLabel(), mother->GetLabel());
292 if (mother->GetLabel() >= maxLabel)
294 // Printf("Removing %d with label %d", i, part->GetLabel()); ((AliMCParticle*)part)->Particle()->Print(); ((AliMCParticle*)mother)->Particle()->Print();
295 TObject* object = tracks->RemoveAt(i);
296 if (tracks->IsOwner())
303 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
306 //-------------------------------------------------------------------
307 void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
309 // remove particles from weak decays
310 // <tracks> can be the following cases:
311 // a. tracks: in this case the label is taken and then case b.
312 // b. particles: it is checked if IsSecondaryFromWeakDecay is true
313 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
315 TClonesArray* arrayMC = 0;
316 AliMCEvent* mcEvent = 0;
317 if (mcObj->InheritsFrom("AliMCEvent"))
318 mcEvent = static_cast<AliMCEvent*>(mcObj);
319 else if (mcObj->InheritsFrom("TClonesArray"))
320 arrayMC = static_cast<TClonesArray*>(mcObj);
324 AliFatal("Invalid object passed");
327 Int_t before = tracks->GetEntriesFast();
329 for (Int_t i=0; i<before; ++i)
331 AliVParticle* part = (AliVParticle*) tracks->At(i);
333 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
334 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
336 if (part->InheritsFrom("AliAODMCParticle"))
338 if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
341 else if (part->InheritsFrom("AliMCParticle") && mcEvent)
343 if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
349 AliFatal("Unknown particle");
352 // Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
353 TObject* object = tracks->RemoveAt(i);
354 if (tracks->IsOwner())
360 if (before > tracks->GetEntriesFast())
361 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
364 //-------------------------------------------------------------------
365 TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts, Bool_t speciesOnTracks, Double_t evtPlane)
367 // Returns an array of particles that pass the cuts, if arrayMC is given each reconstructed particle is replaced by its corresponding MC particles, depending on the parameter onlyprimaries only for primaries
368 // particleSpecies: -1 all particles are returned
369 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
370 // speciesOnTracks if kFALSE, particleSpecies is only applied on the matched MC particle (not on the track itself)
371 // Passing down the Double_t* evtPlane (range [-pi/2,pi/2]) will apply a phi cut with respect to the eventplane between fTrackPhiCutEvPlMin and fTrackPhiCutEvPlMax. For values outside [-pi/2,pi/2], this will be ignored.
373 Int_t nTracks = NParticles(obj);
374 TObjArray* tracks = new TObjArray;
376 // for TPC only tracks
377 Bool_t hasOwnership = kFALSE;
378 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
379 hasOwnership = kTRUE;
382 tracks->SetOwner(hasOwnership);
384 // Loop over tracks or jets
385 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
386 AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, (speciesOnTracks) ? particleSpecies : -1);
389 if (TMath::Abs(evtPlane)<=TMath::Pi()/2) { //evtPlane range: (-pi/2,pi/2)
390 Double_t phiPart = part->Phi(); //range: [0,2*pi)
391 if(phiPart>TMath::Pi()) phiPart-=2*TMath::Pi();
393 Double_t dPhi = 0; //range: [0,pi/2], i.e. the difference over the shortest angle.
394 Double_t diff = TMath::Abs(phiPart-evtPlane);
395 if(diff<=TMath::Pi()/2) dPhi = diff;
396 else if(diff<=TMath::Pi()) dPhi = TMath::Pi()-diff;
397 else dPhi = diff-TMath::Pi();
399 if(dPhi<fTrackPhiCutEvPlMin || dPhi>fTrackPhiCutEvPlMax) {
407 if (TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin || part->Pt() < fTrackPtMin)
414 // Printf("%p %p %d Accepted %d %f %f %f", obj, arrayMC, particleSpecies, ipart, part->Eta(), part->Phi(), part->Pt());
417 Int_t label = part->GetLabel();
420 // re-define part as the matched MC particle
421 part = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
431 //-------------------------------------------------------------------
432 TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
434 // particleSpecies: -1 all particles are returned
435 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
437 Int_t nTracks = NParticles(obj);
438 TObjArray* tracksReconstructed = new TObjArray;
439 TObjArray* tracksOriginal = new TObjArray;
440 TObjArray* tracksFake = new TObjArray;
442 // for TPC only tracks
443 Bool_t hasOwnership = kFALSE;
444 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
445 hasOwnership = kTRUE;
447 tracksReconstructed->SetOwner(hasOwnership);
448 tracksFake->SetOwner(hasOwnership);
450 // Loop over tracks or jets
451 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
452 AliVParticle* partReconstructed = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
453 if (!partReconstructed) continue;
456 if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || TMath::Abs(partReconstructed->Eta()) < fTrackEtaCutMin || partReconstructed->Pt() < fTrackPtMin)
459 delete partReconstructed;
463 Int_t label = partReconstructed->GetLabel();
467 Printf(">>> TPC only track:");
468 partReconstructed->Print();
469 partReconstructed->Dump();
470 Printf(">>> Global track:");
471 ((AliESDEvent*) obj)->GetTrack(ipart)->Dump();
472 Printf("Fake (TPC only): eta = %f, phi = %f, pT = %f, ncl = %d, dedx = %f", partReconstructed->Eta(), partReconstructed->Phi(), partReconstructed->Pt(), ((AliESDtrack*) partReconstructed)->GetTPCclusters(0), ((AliESDtrack*) partReconstructed)->GetTPCsignal());
473 Printf("Fake (global ): eta = %f, phi = %f, pT = %f, ncl = %d, dedx = %f", ((AliESDEvent*) obj)->GetTrack(ipart)->Eta(), ((AliESDEvent*) obj)->GetTrack(ipart)->Phi(), ((AliESDEvent*) obj)->GetTrack(ipart)->Pt(), ((AliESDEvent*) obj)->GetTrack(ipart)->GetTPCclusters(0), ((AliESDEvent*) obj)->GetTrack(ipart)->GetTPCsignal());
475 tracksFake->AddLast(partReconstructed);
479 AliVParticle* partOriginal = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
483 delete partReconstructed;
487 tracksReconstructed->AddLast(partReconstructed);
488 tracksOriginal->AddLast(partOriginal);
490 TObjArray* pairs = new TObjArray;
491 pairs->SetOwner(kTRUE);
492 pairs->Add(tracksReconstructed);
493 pairs->Add(tracksOriginal);
494 pairs->Add(tracksFake);
498 //-------------------------------------------------------------------
499 TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
502 // Returns two lists of particles, one for MIN and one for MAX region
503 Double_t sumpT1 = 0.;
504 Double_t sumpT2 = 0.;
506 Int_t particles1 = transv1->GetEntries();
507 Int_t particles2 = transv2->GetEntries();
509 // Loop on transverse region 1
510 for (Int_t i=0; i<particles1; i++){
511 AliVParticle *part = (AliVParticle*)transv1->At(i);
512 sumpT1 += part->Pt();
515 // Loop on transverse region 2
516 for (Int_t i=0; i<particles2; i++){
517 AliVParticle *part = (AliVParticle*)transv2->At(i);
518 sumpT2 += part->Pt();
521 TObjArray *regionParticles = new TObjArray;
522 if ( sumpT2 >= sumpT1 ){
523 regionParticles->AddLast(transv1); // MIN
524 regionParticles->AddLast(transv2); // MAX
526 regionParticles->AddLast(transv2); // MIN
527 regionParticles->AddLast(transv1); // MAX
530 return regionParticles;
533 //-------------------------------------------------------------------
534 Int_t AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
537 //Returns the number of particles in AliAODMCParticle array or AliAODTracks or AliESDTracks
541 if (obj->InheritsFrom("TClonesArray")){ // MC particles
542 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
543 nTracks = arrayMC->GetEntriesFast();
544 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
545 TObjArray *array = static_cast<TObjArray*>(obj);
546 nTracks = array->GetEntriesFast();
547 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD tracks
548 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
549 nTracks = aodEvent->GetNTracks();
550 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD tracks
551 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
552 nTracks = esdEvent->GetNumberOfTracks();
553 }else if (obj->InheritsFrom("AliMCEvent")){ // RECO ESD tracks
554 AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
555 nTracks = mcEvent->GetNumberOfTracks();
557 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
564 //-------------------------------------------------------------------
565 AliVParticle* AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
567 // Returns track or MC particle at position "ipart" if passes selection criteria
568 // particleSpecies: -1 all particles are returned
569 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
570 AliVParticle *part=0;
572 if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
573 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
574 part = (AliVParticle*)arrayMC->At( ipart );
576 // eventually only primaries
577 if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
578 // eventually only hadrons
580 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
581 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
582 TMath::Abs(pdgCode)==2212 || // Proton
583 TMath::Abs(pdgCode)==321; // Kaon
584 if (!isHadron) return 0;
586 if (particleSpecies != -1) {
587 // find the primary mother
588 AliVParticle* mother = part;
589 if(fCheckMotherPDG) {
590 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
592 if (((AliAODMCParticle*)mother)->GetMother() < 0)
598 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
605 Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
606 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
608 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
610 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
612 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
617 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
618 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
622 // this code prints the details of the mother that is missing in the AOD
623 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
625 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
627 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
628 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
629 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
632 if (particleSpecies != 3)
637 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
638 TObjArray *array = static_cast<TObjArray*>(obj);
639 part = (AliVParticle*)array->At( ipart );
641 }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
642 AliMCEvent* mcEvent = static_cast<AliMCEvent*>(obj);
643 part = mcEvent->GetTrack( ipart );
645 // eventually only primaries
646 if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
647 // eventually only hadrons
650 Int_t pdgCode = part->GetPdgCode();
651 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
652 TMath::Abs(pdgCode)==2212 || // Proton
653 TMath::Abs(pdgCode)==321; // Kaon
654 if (!isHadron) return 0;
657 if (particleSpecies != -1) {
658 // find the primary mother
659 AliMCParticle* mother = (AliMCParticle*) part;
661 if(fCheckMotherPDG) {
662 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
664 // Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
665 if (mother->GetMother() < 0)
671 mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
678 Int_t pdgCode = mother->PdgCode();
679 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
681 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
683 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
685 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
690 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
691 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
696 // this code prints the details of the mother that is missing in the AOD
697 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
699 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
701 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
702 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
703 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
706 if (particleSpecies != 3)
710 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
711 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
712 part = aodEvent->GetTrack(ipart);
714 // track selection cuts
715 if (fFilterBit != 0 && !(((AliAODTrack*)part)->TestFilterBit(fFilterBit))) return 0;
716 if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
721 const AliVVertex* vertex = aodEvent->GetPrimaryVertex();
727 AliAODTrack* clone = (AliAODTrack*) part->Clone();
728 Bool_t success = clone->PropagateToDCA(vertex, aodEvent->GetHeader()->GetMagneticField(), 3, pos, covar);
733 // Printf("%f", ((AliAODTrack*)part)->DCA());
734 // Printf("%f", pos[0]);
735 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(part->Pt()))
739 if (fSharedClusterCut >= 0)
741 Double_t frac = Double_t(((AliAODTrack*)part)->GetTPCnclsS()) / Double_t(((AliAODTrack*)part)->GetTPCncls());
742 if (frac > fSharedClusterCut)
746 if (fCrossedRowsCut >= 0)
748 if (((AliAODTrack*) part)->GetTPCNCrossedRows() < fCrossedRowsCut)
752 if (fFoundFractionCut >= 0)
754 UInt_t findableClusters = ((AliAODTrack*) part)->GetTPCNclsF();
755 if (findableClusters == 0)
757 if (((Double_t) ((AliAODTrack*) part)->GetTPCNCrossedRows() / findableClusters) < fFoundFractionCut)
761 // eventually only hadrons
763 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
764 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
765 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
766 if (!isHadron) return 0;
769 if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0;
771 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
772 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
773 part = esdEvent->GetTrack(ipart);
776 // track selection cuts
777 if (!( ApplyCuts(part)) )
780 if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
782 if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
784 // create TPC only tracks constrained to the SPD vertex
786 const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
788 AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
791 // Printf(">%f %f %f", track->Eta(), track->Phi(), track->Pt());
794 // only constrain tracks above threshold
795 AliExternalTrackParam exParam;
796 // take the B-feild from the ESD, no 3D fieldMap available at this point
797 Bool_t relate = kFALSE;
798 relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
801 // Printf("relating failed");
805 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
808 // Printf(">%f %f %f\n", track->Eta(), track->Phi(), track->Pt());
812 else if (fFilterBit == 2048)
817 AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
818 // Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
820 if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
822 // Float_t ptBefore = esdTrack->Pt();
823 // set constrained pT as default pT
824 if (!esdTrack->GetConstrainedParam())
826 esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
827 // Printf("%f %f", ptBefore, esdTrack->Pt());
832 // eventually only hadrons
835 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
836 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
837 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
838 if (!isHadron) return 0;
841 if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0; // If it is -1 you take all the particles
844 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
849 if (!part->Charge())return 0;
851 part->SetUniqueID(fEventCounter * 100000 + ipart);
856 //-------------------------------------------------------------------
857 void AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
859 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
862 static int i; // "static" to save stack space
865 while (last - first > 1) {
869 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
871 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
887 if (j - first < last - (j + 1)) {
888 QSortTracks(a, first, j);
889 first = j + 1; // QSortTracks(j + 1, last);
891 QSortTracks(a, j + 1, last);
892 last = j; // QSortTracks(first, j);
897 //____________________________________________________________________
898 TObjArray* AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
901 // Assign particles to towards, away or transverse regions.
902 // Returns a lists of particles for each region.
904 static const Double_t k60rad = 60.*TMath::Pi()/180.;
905 static const Double_t k120rad = 120.*TMath::Pi()/180.;
907 // Define output lists of particles
908 TList *toward = new TList();
909 TList *away = new TList();
910 // Two transverse regions, for the moment those are not yet MIN and MAX!!!
911 // MIN and MAX can be sorted in GetMinMaxRegion function
912 TList *transverse1 = new TList();
913 TList *transverse2 = new TList();
915 TObjArray *regionParticles = new TObjArray;
916 regionParticles->SetOwner(kTRUE);
918 regionParticles->AddLast(toward);
919 regionParticles->AddLast(away);
920 regionParticles->AddLast(transverse1);
921 regionParticles->AddLast(transverse2);
924 return regionParticles;
926 // Switch to vector for leading particle
927 TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
929 Int_t nTracks = NParticles(obj);
930 if( !nTracks ) return 0;
932 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
933 AliVParticle* part = ParticleWithCuts(obj, ipart);
935 //Switch to vectors for particles
936 TVector3 partVect(part->Px(), part->Py(), part->Pz());
939 if( TMath::Abs(partVect.Eta()) > fTrackEtaCut || TMath::Abs(partVect.Eta()) < fTrackEtaCutMin) continue;
940 // transverse regions
941 if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
942 if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1; //right
944 if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2; //forward
945 if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2; //backward
947 // skip leading particle
951 if (!region)continue;
952 if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
953 Int_t label = ((AliAODTrack*)part)->GetLabel();
954 // re-define part as the matched MC particle
955 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
957 // skip leading particle
961 if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
962 Int_t label = ((AliESDtrack*)part)->GetLabel();
963 // look for the matched MC particle (but do not re-define part)
964 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
967 if ( region == 1 ) transverse1->Add(part);
968 if ( region == -1 ) transverse2->Add(part);
969 if ( region == 2 ) toward->Add(part);
970 if ( region == -2 ) away->Add(part);
972 }//end loop on tracks
974 return regionParticles;
979 //____________________________________________________________________
980 Bool_t AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
985 // Use AliPhysicsSelection to select good events, works for ESD and AOD
986 if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
992 //____________________________________________________________________
993 Bool_t AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
996 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
998 if (obj->InheritsFrom("AliAODEvent")){
999 Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
1001 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
1002 Int_t nTracksPrim = vertex->GetNContributors();
1003 Double_t zVertex = vertex->GetZ();
1004 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1005 // Reject TPC only vertex
1006 TString name(vertex->GetName());
1007 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
1009 // Select a quality vertex by number of tracks?
1010 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
1011 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1014 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1015 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1017 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1019 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1024 if (obj->InheritsFrom("AliMCEvent"))
1026 if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) >= zed)
1028 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1033 if (obj->InheritsFrom("AliAODMCHeader"))
1035 if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= zed)
1037 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1042 // ESD case for DCA studies
1043 if (obj->InheritsFrom("AliESDEvent")){
1044 AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
1046 Int_t nTracksPrim = vertex->GetNContributors();
1047 Double_t zVertex = vertex->GetZ();
1048 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1049 // Reject SPD or TPC only vertex
1050 TString name(vertex->GetName());
1051 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
1053 // Select a quality vertex by number of tracks?
1054 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
1055 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1058 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1059 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1061 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1063 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1071 //____________________________________________________________________
1073 Bool_t AliAnalyseLeadingTrackUE::CheckTrack(AliVParticle * part)
1075 // check if the track status flags are set
1077 UInt_t status=((AliVTrack*)part)->GetStatus();
1078 if ((status & fTrackStatus) == fTrackStatus)