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),
65 fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
67 fSharedClusterCut(-1),
69 fEsdTrackCutsExtra1(0x0),
70 fEsdTrackCutsExtra2(0x0),
78 //-------------------------------------------------------------------
79 AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
81 // assignment operator
86 //-------------------------------------------------------------------
87 AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
95 //____________________________________________________________________
96 Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
99 // select track according to set of cuts
100 if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
101 if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->IsSelected(track)) return kFALSE;
107 //____________________________________________________________________
108 void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
110 // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
111 // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
114 filterbit = fFilterBit;
116 if (filterbit == 128)
118 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
119 fEsdTrackCuts->SetMinNClustersTPC(70);
121 else if (filterbit == 256)
124 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
125 fEsdTrackCuts->SetMinNClustersTPC(80);
126 fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
127 fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
128 fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
130 else if (filterbit == 512)
133 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
134 fEsdTrackCuts->SetMinNClustersTPC(60);
135 fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
136 fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
137 fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
139 else if (filterbit == 1024)
141 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
142 fEsdTrackCuts->SetMinNClustersTPC(-1);
143 fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
144 fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
146 else if (filterbit == 2048) // mimic hybrid tracks
148 // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
149 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
150 fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
151 fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
152 fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
153 fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
154 fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
155 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
156 fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
158 // Add SPD requirement
159 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
160 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
161 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
163 fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
164 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
165 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
167 else if (filterbit == 4096) // require TOF matching
169 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
170 // FIXME: TOF REQUIREMENTS ARE IN GetParticleSpecies FOR THE MOMENT
174 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
175 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
177 // Add SPD requirement
178 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
179 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
181 // Add SDD requirement
182 fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
183 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
187 //____________________________________________________________________
188 TObjArray* AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
191 // Returns an array of charged particles (or jets) ordered according to their pT.
193 Int_t nTracks = NParticles(obj);
196 if( !nTracks ) return 0;
198 // Define array of AliVParticle objects
199 TObjArray* tracks = new TObjArray(nTracks);
201 // Loop over tracks or jets
202 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
203 AliVParticle* part = ParticleWithCuts( obj, ipart );
205 // Accept leading-tracks in a limited pseudo-rapidity range
206 if( TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin ) continue;
207 tracks->AddLast( part );
209 // Order tracks by pT
210 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
212 nTracks = tracks->GetEntriesFast();
213 if( !nTracks ) return 0;
219 void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
221 // remove injected signals (primaries above <maxLabel>)
222 // <tracks> can be the following cases:
223 // a. tracks: in this case the label is taken and then case b.
224 // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
225 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
227 TClonesArray* arrayMC = 0;
228 AliMCEvent* mcEvent = 0;
229 if (mcObj->InheritsFrom("AliMCEvent"))
230 mcEvent = static_cast<AliMCEvent*>(mcObj);
231 else if (mcObj->InheritsFrom("TClonesArray"))
232 arrayMC = static_cast<TClonesArray*>(mcObj);
236 AliFatal("Invalid object passed");
239 Int_t before = tracks->GetEntriesFast();
241 for (Int_t i=0; i<before; ++i)
243 AliVParticle* part = (AliVParticle*) tracks->At(i);
245 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
246 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
248 AliVParticle* mother = part;
251 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
253 if (((AliMCParticle*)mother)->GetMother() < 0)
259 mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
266 // find the primary mother
267 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
269 if (((AliAODMCParticle*)mother)->GetMother() < 0)
275 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
283 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
287 // Printf("%d %d %d", i, part->GetLabel(), mother->GetLabel());
288 if (mother->GetLabel() >= maxLabel)
290 // Printf("Removing %d with label %d", i, part->GetLabel()); ((AliMCParticle*)part)->Particle()->Print(); ((AliMCParticle*)mother)->Particle()->Print();
291 TObject* object = tracks->RemoveAt(i);
292 if (tracks->IsOwner())
299 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
302 //-------------------------------------------------------------------
303 void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
305 // remove particles from weak decays
306 // <tracks> can be the following cases:
307 // a. tracks: in this case the label is taken and then case b.
308 // b. particles: it is checked if IsSecondaryFromWeakDecay is true
309 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
311 TClonesArray* arrayMC = 0;
312 AliMCEvent* mcEvent = 0;
313 if (mcObj->InheritsFrom("AliMCEvent"))
314 mcEvent = static_cast<AliMCEvent*>(mcObj);
315 else if (mcObj->InheritsFrom("TClonesArray"))
316 arrayMC = static_cast<TClonesArray*>(mcObj);
320 AliFatal("Invalid object passed");
323 Int_t before = tracks->GetEntriesFast();
325 for (Int_t i=0; i<before; ++i)
327 AliVParticle* part = (AliVParticle*) tracks->At(i);
329 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
330 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
332 if (part->InheritsFrom("AliAODMCParticle"))
334 if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
337 else if (part->InheritsFrom("AliMCParticle") && mcEvent)
339 if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
345 AliFatal("Unknown particle");
348 // Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
349 TObject* object = tracks->RemoveAt(i);
350 if (tracks->IsOwner())
356 if (before > tracks->GetEntriesFast())
357 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
360 //-------------------------------------------------------------------
361 TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts, Bool_t speciesOnTracks)
363 // 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
364 // particleSpecies: -1 all particles are returned
365 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
366 // speciesOnTracks if kFALSE, particleSpecies is only applied on the matched MC particle (not on the track itself)
368 Int_t nTracks = NParticles(obj);
369 TObjArray* tracks = new TObjArray;
371 // for TPC only tracks
372 Bool_t hasOwnership = kFALSE;
373 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
374 hasOwnership = kTRUE;
377 tracks->SetOwner(hasOwnership);
379 // Loop over tracks or jets
380 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
381 AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, (speciesOnTracks) ? particleSpecies : -1);
385 if (TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin || part->Pt() < fTrackPtMin)
392 // Printf("%p %p %d Accepted %d %f %f %f", obj, arrayMC, particleSpecies, ipart, part->Eta(), part->Phi(), part->Pt());
395 Int_t label = part->GetLabel();
398 // re-define part as the matched MC particle
399 part = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
409 //-------------------------------------------------------------------
410 TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
412 // particleSpecies: -1 all particles are returned
413 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
415 Int_t nTracks = NParticles(obj);
416 TObjArray* tracksReconstructed = new TObjArray;
417 TObjArray* tracksOriginal = new TObjArray;
418 TObjArray* tracksFake = new TObjArray;
420 // for TPC only tracks
421 Bool_t hasOwnership = kFALSE;
422 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
423 hasOwnership = kTRUE;
425 tracksReconstructed->SetOwner(hasOwnership);
426 tracksFake->SetOwner(hasOwnership);
428 // Loop over tracks or jets
429 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
430 AliVParticle* partReconstructed = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
431 if (!partReconstructed) continue;
434 if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || TMath::Abs(partReconstructed->Eta()) < fTrackEtaCutMin || partReconstructed->Pt() < fTrackPtMin)
437 delete partReconstructed;
441 Int_t label = partReconstructed->GetLabel();
445 Printf(">>> TPC only track:");
446 partReconstructed->Print();
447 partReconstructed->Dump();
448 Printf(">>> Global track:");
449 ((AliESDEvent*) obj)->GetTrack(ipart)->Dump();
450 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());
451 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());
453 tracksFake->AddLast(partReconstructed);
457 AliVParticle* partOriginal = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
461 delete partReconstructed;
465 tracksReconstructed->AddLast(partReconstructed);
466 tracksOriginal->AddLast(partOriginal);
468 TObjArray* pairs = new TObjArray;
469 pairs->SetOwner(kTRUE);
470 pairs->Add(tracksReconstructed);
471 pairs->Add(tracksOriginal);
472 pairs->Add(tracksFake);
476 //-------------------------------------------------------------------
477 TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
480 // Returns two lists of particles, one for MIN and one for MAX region
481 Double_t sumpT1 = 0.;
482 Double_t sumpT2 = 0.;
484 Int_t particles1 = transv1->GetEntries();
485 Int_t particles2 = transv2->GetEntries();
487 // Loop on transverse region 1
488 for (Int_t i=0; i<particles1; i++){
489 AliVParticle *part = (AliVParticle*)transv1->At(i);
490 sumpT1 += part->Pt();
493 // Loop on transverse region 2
494 for (Int_t i=0; i<particles2; i++){
495 AliVParticle *part = (AliVParticle*)transv2->At(i);
496 sumpT2 += part->Pt();
499 TObjArray *regionParticles = new TObjArray;
500 if ( sumpT2 >= sumpT1 ){
501 regionParticles->AddLast(transv1); // MIN
502 regionParticles->AddLast(transv2); // MAX
504 regionParticles->AddLast(transv2); // MIN
505 regionParticles->AddLast(transv1); // MAX
508 return regionParticles;
511 //-------------------------------------------------------------------
512 Int_t AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
515 //Returns the number of particles in AliAODMCParticle array or AliAODTracks or AliESDTracks
519 if (obj->InheritsFrom("TClonesArray")){ // MC particles
520 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
521 nTracks = arrayMC->GetEntriesFast();
522 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
523 TObjArray *array = static_cast<TObjArray*>(obj);
524 nTracks = array->GetEntriesFast();
525 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD tracks
526 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
527 nTracks = aodEvent->GetNTracks();
528 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD tracks
529 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
530 nTracks = esdEvent->GetNumberOfTracks();
531 }else if (obj->InheritsFrom("AliMCEvent")){ // RECO ESD tracks
532 AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
533 nTracks = mcEvent->GetNumberOfTracks();
535 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
542 //-------------------------------------------------------------------
543 AliVParticle* AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
545 // Returns track or MC particle at position "ipart" if passes selection criteria
546 // particleSpecies: -1 all particles are returned
547 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
548 AliVParticle *part=0;
550 if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
551 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
552 part = (AliVParticle*)arrayMC->At( ipart );
554 // eventually only primaries
555 if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
556 // eventually only hadrons
558 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
559 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
560 TMath::Abs(pdgCode)==2212 || // Proton
561 TMath::Abs(pdgCode)==321; // Kaon
562 if (!isHadron) return 0;
564 if (particleSpecies != -1) {
565 // find the primary mother
566 AliVParticle* mother = part;
567 if(fCheckMotherPDG) {
568 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
570 if (((AliAODMCParticle*)mother)->GetMother() < 0)
576 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
583 Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
584 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
586 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
588 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
590 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
595 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
596 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
600 // this code prints the details of the mother that is missing in the AOD
601 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
603 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
605 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
606 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
607 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
610 if (particleSpecies != 3)
615 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
616 TObjArray *array = static_cast<TObjArray*>(obj);
617 part = (AliVParticle*)array->At( ipart );
619 }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
620 AliMCEvent* mcEvent = static_cast<AliMCEvent*>(obj);
621 part = mcEvent->GetTrack( ipart );
623 // eventually only primaries
624 if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
625 // eventually only hadrons
628 Int_t pdgCode = part->GetPdgCode();
629 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
630 TMath::Abs(pdgCode)==2212 || // Proton
631 TMath::Abs(pdgCode)==321; // Kaon
632 if (!isHadron) return 0;
635 if (particleSpecies != -1) {
636 // find the primary mother
637 AliMCParticle* mother = (AliMCParticle*) part;
639 if(fCheckMotherPDG) {
640 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
642 // Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
643 if (mother->GetMother() < 0)
649 mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
656 Int_t pdgCode = mother->PdgCode();
657 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
659 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
661 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
663 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
668 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
669 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
674 // this code prints the details of the mother that is missing in the AOD
675 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
677 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
679 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
680 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
681 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
684 if (particleSpecies != 3)
688 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
689 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
690 part = aodEvent->GetTrack(ipart);
692 // track selection cuts
693 if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
694 if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
699 const AliVVertex* vertex = aodEvent->GetPrimaryVertex();
705 AliAODTrack* clone = (AliAODTrack*) part->Clone();
706 Bool_t success = clone->PropagateToDCA(vertex, aodEvent->GetHeader()->GetMagneticField(), 3, pos, covar);
711 // Printf("%f", ((AliAODTrack*)part)->DCA());
712 // Printf("%f", pos[0]);
713 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(part->Pt()))
717 if (fSharedClusterCut >= 0)
719 Double_t frac = Double_t(((AliAODTrack*)part)->GetTPCnclsS()) / Double_t(((AliAODTrack*)part)->GetTPCncls());
720 if (frac > fSharedClusterCut)
724 // eventually only hadrons
726 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
727 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
728 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
729 if (!isHadron) return 0;
732 if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0;
734 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
735 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
736 part = esdEvent->GetTrack(ipart);
739 // track selection cuts
740 if (!( ApplyCuts(part)) )
743 if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
745 if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
747 // create TPC only tracks constrained to the SPD vertex
749 const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
751 AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
754 // Printf(">%f %f %f", track->Eta(), track->Phi(), track->Pt());
757 // only constrain tracks above threshold
758 AliExternalTrackParam exParam;
759 // take the B-feild from the ESD, no 3D fieldMap available at this point
760 Bool_t relate = kFALSE;
761 relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
764 // Printf("relating failed");
768 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
771 // Printf(">%f %f %f\n", track->Eta(), track->Phi(), track->Pt());
775 else if (fFilterBit == 2048)
780 AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
781 // Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
783 if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
785 // Float_t ptBefore = esdTrack->Pt();
786 // set constrained pT as default pT
787 if (!esdTrack->GetConstrainedParam())
789 esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
790 // Printf("%f %f", ptBefore, esdTrack->Pt());
795 // eventually only hadrons
798 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
799 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
800 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
801 if (!isHadron) return 0;
804 if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0; // If it is -1 you take all the particles
807 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
812 if (!part->Charge())return 0;
814 part->SetUniqueID(fEventCounter * 100000 + ipart);
819 //-------------------------------------------------------------------
820 void AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
822 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
825 static int i; // "static" to save stack space
828 while (last - first > 1) {
832 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
834 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
850 if (j - first < last - (j + 1)) {
851 QSortTracks(a, first, j);
852 first = j + 1; // QSortTracks(j + 1, last);
854 QSortTracks(a, j + 1, last);
855 last = j; // QSortTracks(first, j);
860 //____________________________________________________________________
861 TObjArray* AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
864 // Assign particles to towards, away or transverse regions.
865 // Returns a lists of particles for each region.
867 static const Double_t k60rad = 60.*TMath::Pi()/180.;
868 static const Double_t k120rad = 120.*TMath::Pi()/180.;
870 // Define output lists of particles
871 TList *toward = new TList();
872 TList *away = new TList();
873 // Two transverse regions, for the moment those are not yet MIN and MAX!!!
874 // MIN and MAX can be sorted in GetMinMaxRegion function
875 TList *transverse1 = new TList();
876 TList *transverse2 = new TList();
878 TObjArray *regionParticles = new TObjArray;
879 regionParticles->SetOwner(kTRUE);
881 regionParticles->AddLast(toward);
882 regionParticles->AddLast(away);
883 regionParticles->AddLast(transverse1);
884 regionParticles->AddLast(transverse2);
887 return regionParticles;
889 // Switch to vector for leading particle
890 TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
892 Int_t nTracks = NParticles(obj);
893 if( !nTracks ) return 0;
895 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
896 AliVParticle* part = ParticleWithCuts(obj, ipart);
898 //Switch to vectors for particles
899 TVector3 partVect(part->Px(), part->Py(), part->Pz());
902 if( TMath::Abs(partVect.Eta()) > fTrackEtaCut || TMath::Abs(partVect.Eta()) < fTrackEtaCutMin) continue;
903 // transverse regions
904 if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
905 if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1; //right
907 if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2; //forward
908 if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2; //backward
910 // skip leading particle
914 if (!region)continue;
915 if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
916 Int_t label = ((AliAODTrack*)part)->GetLabel();
917 // re-define part as the matched MC particle
918 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
920 // skip leading particle
924 if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
925 Int_t label = ((AliESDtrack*)part)->GetLabel();
926 // look for the matched MC particle (but do not re-define part)
927 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
930 if ( region == 1 ) transverse1->Add(part);
931 if ( region == -1 ) transverse2->Add(part);
932 if ( region == 2 ) toward->Add(part);
933 if ( region == -2 ) away->Add(part);
935 }//end loop on tracks
937 return regionParticles;
942 //____________________________________________________________________
943 Bool_t AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
948 // Use AliPhysicsSelection to select good events, works for ESD and AOD
949 if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
955 //____________________________________________________________________
956 Bool_t AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
959 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
961 if (obj->InheritsFrom("AliAODEvent")){
962 Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
964 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
965 Int_t nTracksPrim = vertex->GetNContributors();
966 Double_t zVertex = vertex->GetZ();
967 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
968 // Reject TPC only vertex
969 TString name(vertex->GetName());
970 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
972 // Select a quality vertex by number of tracks?
973 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
974 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
977 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
978 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
980 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
982 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
987 if (obj->InheritsFrom("AliMCEvent"))
989 if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) >= zed)
991 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
996 if (obj->InheritsFrom("AliAODMCHeader"))
998 if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= zed)
1000 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1005 // ESD case for DCA studies
1006 if (obj->InheritsFrom("AliESDEvent")){
1007 AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
1009 Int_t nTracksPrim = vertex->GetNContributors();
1010 Double_t zVertex = vertex->GetZ();
1011 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1012 // Reject SPD or TPC only vertex
1013 TString name(vertex->GetName());
1014 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
1016 // Select a quality vertex by number of tracks?
1017 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
1018 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1021 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1022 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1024 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1026 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1034 //____________________________________________________________________
1036 Bool_t AliAnalyseLeadingTrackUE::CheckTrack(AliVParticle * part)
1038 // check if the track status flags are set
1040 UInt_t status=((AliVTrack*)part)->GetStatus();
1041 if ((status & fTrackStatus) == fTrackStatus)