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 fFoundFractionCut(-1),
71 fEsdTrackCutsExtra1(0x0),
72 fEsdTrackCutsExtra2(0x0),
80 //-------------------------------------------------------------------
81 AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
83 // assignment operator
88 //-------------------------------------------------------------------
89 AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
97 //____________________________________________________________________
98 Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
101 // select track according to set of cuts
102 if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
103 if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->IsSelected(track)) return kFALSE;
109 //____________________________________________________________________
110 void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
112 // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
113 // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
116 filterbit = fFilterBit;
118 if (filterbit == 128)
120 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
121 fEsdTrackCuts->SetMinNClustersTPC(70);
123 else if (filterbit == 256)
126 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
127 fEsdTrackCuts->SetMinNClustersTPC(80);
128 fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
129 fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
130 fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
132 else if (filterbit == 512)
135 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
136 fEsdTrackCuts->SetMinNClustersTPC(60);
137 fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
138 fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
139 fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
141 else if (filterbit == 1024)
143 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
144 fEsdTrackCuts->SetMinNClustersTPC(-1);
145 fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
146 fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
148 else if (filterbit == 2048) // mimic hybrid tracks
150 // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
151 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
152 fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
153 fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
154 fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
155 fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
156 fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
157 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
158 fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
160 // Add SPD requirement
161 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
162 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
163 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
165 fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
166 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
167 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
169 else if (filterbit == 4096) // require TOF matching
171 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
172 // FIXME: TOF REQUIREMENTS ARE IN GetParticleSpecies FOR THE MOMENT
176 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
177 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
179 // Add SPD requirement
180 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
181 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
183 // Add SDD requirement
184 fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
185 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
189 //____________________________________________________________________
190 TObjArray* AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
193 // Returns an array of charged particles (or jets) ordered according to their pT.
195 Int_t nTracks = NParticles(obj);
198 if( !nTracks ) return 0;
200 // Define array of AliVParticle objects
201 TObjArray* tracks = new TObjArray(nTracks);
203 // Loop over tracks or jets
204 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
205 AliVParticle* part = ParticleWithCuts( obj, ipart );
207 // Accept leading-tracks in a limited pseudo-rapidity range
208 if( TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin ) continue;
209 tracks->AddLast( part );
211 // Order tracks by pT
212 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
214 nTracks = tracks->GetEntriesFast();
215 if( !nTracks ) return 0;
221 void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
223 // remove injected signals (primaries above <maxLabel>)
224 // <tracks> can be the following cases:
225 // a. tracks: in this case the label is taken and then case b.
226 // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
227 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
229 TClonesArray* arrayMC = 0;
230 AliMCEvent* mcEvent = 0;
231 if (mcObj->InheritsFrom("AliMCEvent"))
232 mcEvent = static_cast<AliMCEvent*>(mcObj);
233 else if (mcObj->InheritsFrom("TClonesArray"))
234 arrayMC = static_cast<TClonesArray*>(mcObj);
238 AliFatal("Invalid object passed");
241 Int_t before = tracks->GetEntriesFast();
243 for (Int_t i=0; i<before; ++i)
245 AliVParticle* part = (AliVParticle*) tracks->At(i);
247 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
248 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
250 AliVParticle* mother = part;
253 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
255 if (((AliMCParticle*)mother)->GetMother() < 0)
261 mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
268 // find the primary mother
269 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
271 if (((AliAODMCParticle*)mother)->GetMother() < 0)
277 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
285 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
289 // Printf("%d %d %d", i, part->GetLabel(), mother->GetLabel());
290 if (mother->GetLabel() >= maxLabel)
292 // Printf("Removing %d with label %d", i, part->GetLabel()); ((AliMCParticle*)part)->Particle()->Print(); ((AliMCParticle*)mother)->Particle()->Print();
293 TObject* object = tracks->RemoveAt(i);
294 if (tracks->IsOwner())
301 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
304 //-------------------------------------------------------------------
305 void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
307 // remove particles from weak decays
308 // <tracks> can be the following cases:
309 // a. tracks: in this case the label is taken and then case b.
310 // b. particles: it is checked if IsSecondaryFromWeakDecay is true
311 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
313 TClonesArray* arrayMC = 0;
314 AliMCEvent* mcEvent = 0;
315 if (mcObj->InheritsFrom("AliMCEvent"))
316 mcEvent = static_cast<AliMCEvent*>(mcObj);
317 else if (mcObj->InheritsFrom("TClonesArray"))
318 arrayMC = static_cast<TClonesArray*>(mcObj);
322 AliFatal("Invalid object passed");
325 Int_t before = tracks->GetEntriesFast();
327 for (Int_t i=0; i<before; ++i)
329 AliVParticle* part = (AliVParticle*) tracks->At(i);
331 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
332 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
334 if (part->InheritsFrom("AliAODMCParticle"))
336 if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
339 else if (part->InheritsFrom("AliMCParticle") && mcEvent)
341 if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
347 AliFatal("Unknown particle");
350 // Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
351 TObject* object = tracks->RemoveAt(i);
352 if (tracks->IsOwner())
358 if (before > tracks->GetEntriesFast())
359 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
362 //-------------------------------------------------------------------
363 TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts, Bool_t speciesOnTracks)
365 // 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
366 // particleSpecies: -1 all particles are returned
367 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
368 // speciesOnTracks if kFALSE, particleSpecies is only applied on the matched MC particle (not on the track itself)
370 Int_t nTracks = NParticles(obj);
371 TObjArray* tracks = new TObjArray;
373 // for TPC only tracks
374 Bool_t hasOwnership = kFALSE;
375 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
376 hasOwnership = kTRUE;
379 tracks->SetOwner(hasOwnership);
381 // Loop over tracks or jets
382 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
383 AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, (speciesOnTracks) ? particleSpecies : -1);
387 if (TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin || part->Pt() < fTrackPtMin)
394 // Printf("%p %p %d Accepted %d %f %f %f", obj, arrayMC, particleSpecies, ipart, part->Eta(), part->Phi(), part->Pt());
397 Int_t label = part->GetLabel();
400 // re-define part as the matched MC particle
401 part = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
411 //-------------------------------------------------------------------
412 TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
414 // particleSpecies: -1 all particles are returned
415 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
417 Int_t nTracks = NParticles(obj);
418 TObjArray* tracksReconstructed = new TObjArray;
419 TObjArray* tracksOriginal = new TObjArray;
420 TObjArray* tracksFake = new TObjArray;
422 // for TPC only tracks
423 Bool_t hasOwnership = kFALSE;
424 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
425 hasOwnership = kTRUE;
427 tracksReconstructed->SetOwner(hasOwnership);
428 tracksFake->SetOwner(hasOwnership);
430 // Loop over tracks or jets
431 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
432 AliVParticle* partReconstructed = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
433 if (!partReconstructed) continue;
436 if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || TMath::Abs(partReconstructed->Eta()) < fTrackEtaCutMin || partReconstructed->Pt() < fTrackPtMin)
439 delete partReconstructed;
443 Int_t label = partReconstructed->GetLabel();
447 Printf(">>> TPC only track:");
448 partReconstructed->Print();
449 partReconstructed->Dump();
450 Printf(">>> Global track:");
451 ((AliESDEvent*) obj)->GetTrack(ipart)->Dump();
452 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());
453 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());
455 tracksFake->AddLast(partReconstructed);
459 AliVParticle* partOriginal = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
463 delete partReconstructed;
467 tracksReconstructed->AddLast(partReconstructed);
468 tracksOriginal->AddLast(partOriginal);
470 TObjArray* pairs = new TObjArray;
471 pairs->SetOwner(kTRUE);
472 pairs->Add(tracksReconstructed);
473 pairs->Add(tracksOriginal);
474 pairs->Add(tracksFake);
478 //-------------------------------------------------------------------
479 TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
482 // Returns two lists of particles, one for MIN and one for MAX region
483 Double_t sumpT1 = 0.;
484 Double_t sumpT2 = 0.;
486 Int_t particles1 = transv1->GetEntries();
487 Int_t particles2 = transv2->GetEntries();
489 // Loop on transverse region 1
490 for (Int_t i=0; i<particles1; i++){
491 AliVParticle *part = (AliVParticle*)transv1->At(i);
492 sumpT1 += part->Pt();
495 // Loop on transverse region 2
496 for (Int_t i=0; i<particles2; i++){
497 AliVParticle *part = (AliVParticle*)transv2->At(i);
498 sumpT2 += part->Pt();
501 TObjArray *regionParticles = new TObjArray;
502 if ( sumpT2 >= sumpT1 ){
503 regionParticles->AddLast(transv1); // MIN
504 regionParticles->AddLast(transv2); // MAX
506 regionParticles->AddLast(transv2); // MIN
507 regionParticles->AddLast(transv1); // MAX
510 return regionParticles;
513 //-------------------------------------------------------------------
514 Int_t AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
517 //Returns the number of particles in AliAODMCParticle array or AliAODTracks or AliESDTracks
521 if (obj->InheritsFrom("TClonesArray")){ // MC particles
522 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
523 nTracks = arrayMC->GetEntriesFast();
524 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
525 TObjArray *array = static_cast<TObjArray*>(obj);
526 nTracks = array->GetEntriesFast();
527 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD tracks
528 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
529 nTracks = aodEvent->GetNTracks();
530 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD tracks
531 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
532 nTracks = esdEvent->GetNumberOfTracks();
533 }else if (obj->InheritsFrom("AliMCEvent")){ // RECO ESD tracks
534 AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
535 nTracks = mcEvent->GetNumberOfTracks();
537 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
544 //-------------------------------------------------------------------
545 AliVParticle* AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
547 // Returns track or MC particle at position "ipart" if passes selection criteria
548 // particleSpecies: -1 all particles are returned
549 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
550 AliVParticle *part=0;
552 if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
553 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
554 part = (AliVParticle*)arrayMC->At( ipart );
556 // eventually only primaries
557 if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
558 // eventually only hadrons
560 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
561 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
562 TMath::Abs(pdgCode)==2212 || // Proton
563 TMath::Abs(pdgCode)==321; // Kaon
564 if (!isHadron) return 0;
566 if (particleSpecies != -1) {
567 // find the primary mother
568 AliVParticle* mother = part;
569 if(fCheckMotherPDG) {
570 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
572 if (((AliAODMCParticle*)mother)->GetMother() < 0)
578 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
585 Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
586 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
588 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
590 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
592 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
597 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
598 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
602 // this code prints the details of the mother that is missing in the AOD
603 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
605 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
607 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
608 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
609 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
612 if (particleSpecies != 3)
617 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
618 TObjArray *array = static_cast<TObjArray*>(obj);
619 part = (AliVParticle*)array->At( ipart );
621 }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
622 AliMCEvent* mcEvent = static_cast<AliMCEvent*>(obj);
623 part = mcEvent->GetTrack( ipart );
625 // eventually only primaries
626 if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
627 // eventually only hadrons
630 Int_t pdgCode = part->GetPdgCode();
631 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
632 TMath::Abs(pdgCode)==2212 || // Proton
633 TMath::Abs(pdgCode)==321; // Kaon
634 if (!isHadron) return 0;
637 if (particleSpecies != -1) {
638 // find the primary mother
639 AliMCParticle* mother = (AliMCParticle*) part;
641 if(fCheckMotherPDG) {
642 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
644 // Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
645 if (mother->GetMother() < 0)
651 mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
658 Int_t pdgCode = mother->PdgCode();
659 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
661 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
663 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
665 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
670 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
671 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
676 // this code prints the details of the mother that is missing in the AOD
677 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
679 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
681 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
682 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
683 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
686 if (particleSpecies != 3)
690 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
691 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
692 part = aodEvent->GetTrack(ipart);
694 // track selection cuts
695 if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
696 if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
701 const AliVVertex* vertex = aodEvent->GetPrimaryVertex();
707 AliAODTrack* clone = (AliAODTrack*) part->Clone();
708 Bool_t success = clone->PropagateToDCA(vertex, aodEvent->GetHeader()->GetMagneticField(), 3, pos, covar);
713 // Printf("%f", ((AliAODTrack*)part)->DCA());
714 // Printf("%f", pos[0]);
715 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(part->Pt()))
719 if (fSharedClusterCut >= 0)
721 Double_t frac = Double_t(((AliAODTrack*)part)->GetTPCnclsS()) / Double_t(((AliAODTrack*)part)->GetTPCncls());
722 if (frac > fSharedClusterCut)
726 if (fCrossedRowsCut >= 0)
728 if (((AliAODTrack*) part)->GetTPCNCrossedRows() < fCrossedRowsCut)
732 if (fFoundFractionCut >= 0)
734 UInt_t findableClusters = ((AliAODTrack*) part)->GetTPCNclsF();
735 if (findableClusters == 0)
737 if (((Double_t) ((AliAODTrack*) part)->GetTPCNCrossedRows() / findableClusters) < fFoundFractionCut)
741 // eventually only hadrons
743 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
744 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
745 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
746 if (!isHadron) return 0;
749 if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0;
751 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
752 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
753 part = esdEvent->GetTrack(ipart);
756 // track selection cuts
757 if (!( ApplyCuts(part)) )
760 if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
762 if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
764 // create TPC only tracks constrained to the SPD vertex
766 const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
768 AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
771 // Printf(">%f %f %f", track->Eta(), track->Phi(), track->Pt());
774 // only constrain tracks above threshold
775 AliExternalTrackParam exParam;
776 // take the B-feild from the ESD, no 3D fieldMap available at this point
777 Bool_t relate = kFALSE;
778 relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
781 // Printf("relating failed");
785 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
788 // Printf(">%f %f %f\n", track->Eta(), track->Phi(), track->Pt());
792 else if (fFilterBit == 2048)
797 AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
798 // Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
800 if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
802 // Float_t ptBefore = esdTrack->Pt();
803 // set constrained pT as default pT
804 if (!esdTrack->GetConstrainedParam())
806 esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
807 // Printf("%f %f", ptBefore, esdTrack->Pt());
812 // eventually only hadrons
815 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
816 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
817 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
818 if (!isHadron) return 0;
821 if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0; // If it is -1 you take all the particles
824 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
829 if (!part->Charge())return 0;
831 part->SetUniqueID(fEventCounter * 100000 + ipart);
836 //-------------------------------------------------------------------
837 void AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
839 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
842 static int i; // "static" to save stack space
845 while (last - first > 1) {
849 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
851 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
867 if (j - first < last - (j + 1)) {
868 QSortTracks(a, first, j);
869 first = j + 1; // QSortTracks(j + 1, last);
871 QSortTracks(a, j + 1, last);
872 last = j; // QSortTracks(first, j);
877 //____________________________________________________________________
878 TObjArray* AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
881 // Assign particles to towards, away or transverse regions.
882 // Returns a lists of particles for each region.
884 static const Double_t k60rad = 60.*TMath::Pi()/180.;
885 static const Double_t k120rad = 120.*TMath::Pi()/180.;
887 // Define output lists of particles
888 TList *toward = new TList();
889 TList *away = new TList();
890 // Two transverse regions, for the moment those are not yet MIN and MAX!!!
891 // MIN and MAX can be sorted in GetMinMaxRegion function
892 TList *transverse1 = new TList();
893 TList *transverse2 = new TList();
895 TObjArray *regionParticles = new TObjArray;
896 regionParticles->SetOwner(kTRUE);
898 regionParticles->AddLast(toward);
899 regionParticles->AddLast(away);
900 regionParticles->AddLast(transverse1);
901 regionParticles->AddLast(transverse2);
904 return regionParticles;
906 // Switch to vector for leading particle
907 TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
909 Int_t nTracks = NParticles(obj);
910 if( !nTracks ) return 0;
912 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
913 AliVParticle* part = ParticleWithCuts(obj, ipart);
915 //Switch to vectors for particles
916 TVector3 partVect(part->Px(), part->Py(), part->Pz());
919 if( TMath::Abs(partVect.Eta()) > fTrackEtaCut || TMath::Abs(partVect.Eta()) < fTrackEtaCutMin) continue;
920 // transverse regions
921 if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
922 if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1; //right
924 if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2; //forward
925 if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2; //backward
927 // skip leading particle
931 if (!region)continue;
932 if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
933 Int_t label = ((AliAODTrack*)part)->GetLabel();
934 // re-define part as the matched MC particle
935 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
937 // skip leading particle
941 if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
942 Int_t label = ((AliESDtrack*)part)->GetLabel();
943 // look for the matched MC particle (but do not re-define part)
944 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
947 if ( region == 1 ) transverse1->Add(part);
948 if ( region == -1 ) transverse2->Add(part);
949 if ( region == 2 ) toward->Add(part);
950 if ( region == -2 ) away->Add(part);
952 }//end loop on tracks
954 return regionParticles;
959 //____________________________________________________________________
960 Bool_t AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
965 // Use AliPhysicsSelection to select good events, works for ESD and AOD
966 if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
972 //____________________________________________________________________
973 Bool_t AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
976 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
978 if (obj->InheritsFrom("AliAODEvent")){
979 Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
981 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
982 Int_t nTracksPrim = vertex->GetNContributors();
983 Double_t zVertex = vertex->GetZ();
984 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
985 // Reject TPC only vertex
986 TString name(vertex->GetName());
987 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
989 // Select a quality vertex by number of tracks?
990 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
991 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
994 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
995 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
997 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
999 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1004 if (obj->InheritsFrom("AliMCEvent"))
1006 if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) >= zed)
1008 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1013 if (obj->InheritsFrom("AliAODMCHeader"))
1015 if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= zed)
1017 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1022 // ESD case for DCA studies
1023 if (obj->InheritsFrom("AliESDEvent")){
1024 AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
1026 Int_t nTracksPrim = vertex->GetNContributors();
1027 Double_t zVertex = vertex->GetZ();
1028 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1029 // Reject SPD or TPC only vertex
1030 TString name(vertex->GetName());
1031 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
1033 // Select a quality vertex by number of tracks?
1034 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
1035 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1038 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1039 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1041 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1043 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1051 //____________________________________________________________________
1053 Bool_t AliAnalyseLeadingTrackUE::CheckTrack(AliVParticle * part)
1055 // check if the track status flags are set
1057 UInt_t status=((AliVTrack*)part)->GetStatus();
1058 if ((status & fTrackStatus) == fTrackStatus)