]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/Base/AliAnalyseLeadingTrackUE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
including histo to check tracks inside the cone
[u/mrichter/AliRoot.git] / PWGCF / Correlations / Base / AliAnalyseLeadingTrackUE.cxx
CommitLineData
a75aacd6 1/*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: A.Abrahantes, E.Lopez, S.Vallero *
5 * Version 1.0 *
6 * *
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 **************************************************************************/
a75aacd6 15//#include <TBranch.h>
16//#include <TCanvas.h>
17//#include <TChain.h>
18//#include <TFile.h>
19//#include <TH1F.h>
20//#include <TH1I.h>
21//#include <TH2F.h>
22#include <TList.h>
23//#include <TLorentzVector.h>
24#include <TMath.h>
25#include <TObjArray.h>
26#include <TObject.h>
27//#include <TProfile.h>
28//#include <TRandom.h>
29//#include <TSystem.h>
30//#include <TTree.h>
31#include <TVector3.h>
32
33#include "AliAnalyseLeadingTrackUE.h"
34//#include "AliAnalysisTask.h"
35
36//#include "AliAnalysisHelperJetTasks.h"
37//#include "AliAnalysisManager.h"
38#include "AliAODEvent.h"
39//#include "AliAODHandler.h"
a75aacd6 40//#include "AliAODJet.h"
41#include "AliAODMCParticle.h"
42#include "AliAODTrack.h"
43#include "AliESDEvent.h"
44#include "AliESDtrack.h"
45#include "AliESDtrackCuts.h"
46//#include "AliGenPythiaEventHeader.h"
47#include "AliInputEventHandler.h"
48//#include "AliKFVertex.h"
144bd037 49//#include "AliLog.h"
a75aacd6 50#include "AliMCEvent.h"
a75aacd6 51#include "AliVParticle.h"
7dd4dec4 52#include "AliAODMCHeader.h"
a75aacd6 53
ea1919ac 54#include "AliAnalysisManager.h"
55#include "AliMCEventHandler.h"
56#include "AliStack.h"
57
58
a75aacd6 59////////////////////////////////////////////////
60//---------------------------------------------
61// Class for transverse regions analysis
62//---------------------------------------------
63////////////////////////////////////////////////
64
65
66using namespace std;
67
68ClassImp(AliAnalyseLeadingTrackUE)
69
70//-------------------------------------------------------------------
71AliAnalyseLeadingTrackUE::AliAnalyseLeadingTrackUE() :
72 TObject(),
73 fDebug(0),
74 fFilterBit(16),
75 fOnlyHadrons(kFALSE),
3712ba26 76 fTrackEtaCut(0.8),
e0331fd9 77 fTrackPtMin(0),
5cabd2cd 78 fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
3712ba26 79 fEsdTrackCuts(0x0),
8b43ff78 80 fEsdTrackCutsExtra1(0x0),
81 fEsdTrackCutsExtra2(0x0)
a75aacd6 82{
83 // constructor
84}
85
86
87//-------------------------------------------------------------------
88AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
89{
90 // assignment operator
91 return *this;
92}
93
94
95//-------------------------------------------------------------------
96AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
97{
98
99 //clear memory
100
101}
102
103
104//____________________________________________________________________
3712ba26 105Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
a75aacd6 106{
a75aacd6 107
108 // select track according to set of cuts
85bfac17 109 if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
8b43ff78 110 if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->IsSelected(track)) return kFALSE;
3712ba26 111
a75aacd6 112 return kTRUE;
113}
114
115
3712ba26 116//____________________________________________________________________
1be4eb33 117void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
3712ba26 118
85bfac17 119 // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
120 // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
1be4eb33 121
122 if (filterbit == -1)
123 filterbit = fFilterBit;
a75aacd6 124
1be4eb33 125 if (filterbit == 128)
85bfac17 126 {
127 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
128 fEsdTrackCuts->SetMinNClustersTPC(70);
129 }
1be4eb33 130 else if (filterbit == 256)
131 {
132 // syst study
133 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
134 fEsdTrackCuts->SetMinNClustersTPC(80);
135 fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
136 fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
137 fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
138 }
139 else if (filterbit == 512)
140 {
141 // syst study
142 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
143 fEsdTrackCuts->SetMinNClustersTPC(60);
144 fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
145 fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
146 fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
147 }
4bbcfff9 148 else if (filterbit == 1024)
149 {
150 fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
151 fEsdTrackCuts->SetMinNClustersTPC(-1);
152 fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
153 fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
154 }
8b43ff78 155 else if (filterbit == 2048) // mimic hybrid tracks
156 {
157 // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
158 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
159 fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
160 fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
161 fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
162 fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
163 fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
164 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
165
166 // Add SPD requirement
167 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
168 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
169 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
170
171 // empty object
172 fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
173 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
174 // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
175 }
85bfac17 176 else
177 {
178 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
8b43ff78 179 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
85bfac17 180
181 // Add SPD requirement
8b43ff78 182 fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
183 fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
85bfac17 184
185 // Add SDD requirement
8b43ff78 186 fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
187 fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
85bfac17 188 }
3712ba26 189}
a75aacd6 190
a75aacd6 191//____________________________________________________________________
192TObjArray* AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
193{
194
195 // Returns an array of charged particles (or jets) ordered according to their pT.
196
197 Int_t nTracks = NParticles(obj);
198
199
200 if( !nTracks ) return 0;
201
202 // Define array of AliVParticle objects
203 TObjArray* tracks = new TObjArray(nTracks);
204
205 // Loop over tracks or jets
206 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
207 AliVParticle* part = ParticleWithCuts( obj, ipart );
208 if (!part) continue;
209 // Accept leading-tracks in a limited pseudo-rapidity range
210 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
211 tracks->AddLast( part );
212 }
213 // Order tracks by pT
214 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
215
216 nTracks = tracks->GetEntriesFast();
217 if( !nTracks ) return 0;
218
219 return tracks;
220 }
221
222
1ccd8a0a 223void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
224{
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)
230
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);
237 else
238 {
239 arrayMC->Dump();
240 AliFatal("Invalid object passed");
241 }
242
243 Int_t before = tracks->GetEntriesFast();
244
245 for (Int_t i=0; i<before; ++i)
246 {
247 AliVParticle* part = (AliVParticle*) tracks->At(i);
248
249 if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
250 part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
251
252 AliVParticle* mother = part;
253 if (mcEvent)
254 {
255 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
256 {
257 if (((AliMCParticle*)mother)->GetMother() < 0)
258 {
259 mother = 0;
260 break;
261 }
262
263 mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
264 if (!mother)
265 break;
266 }
267 }
268 else
269 {
270 // find the primary mother
271 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
272 {
273 if (((AliAODMCParticle*)mother)->GetMother() < 0)
274 {
275 mother = 0;
276 break;
277 }
278
279 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
280 if (!mother)
281 break;
282 }
283 }
284
285 if (!mother)
286 {
287 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
288 continue;
289 }
290
291 if (mother->GetLabel() > maxLabel)
292 {
293// Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
294 TObject* object = tracks->RemoveAt(i);
295 if (tracks->IsOwner())
296 delete object;
297 }
298 }
299
300 tracks->Compress();
301
302 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
303}
304
b1831bcb 305//-------------------------------------------------------------------
e0331fd9 306TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
b1831bcb 307{
308 // 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
309 // particleSpecies: -1 all particles are returned
310 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
311
312 Int_t nTracks = NParticles(obj);
313 TObjArray* tracks = new TObjArray;
85bfac17 314
315 // for TPC only tracks
eed401dc 316 Bool_t hasOwnership = kFALSE;
8b43ff78 317 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
eed401dc 318 hasOwnership = kTRUE;
319
320 if (!arrayMC)
321 tracks->SetOwner(hasOwnership);
b1831bcb 322
323 // Loop over tracks or jets
324 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
325 AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
326 if (!part) continue;
327
e0331fd9 328 if (useEtaPtCuts)
329 if (TMath::Abs(part->Eta()) > fTrackEtaCut || part->Pt() < fTrackPtMin)
92834c86 330 {
331 if (hasOwnership)
332 delete part;
e0331fd9 333 continue;
92834c86 334 }
e0331fd9 335
eed401dc 336 if (arrayMC) {
337 Int_t label = part->GetLabel();
338 if (hasOwnership)
339 delete part;
b1831bcb 340 // re-define part as the matched MC particle
eed401dc 341 part = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
b1831bcb 342 if (!part)continue;
343 }
344
345 tracks->Add(part);
346 }
347
348 return tracks;
349}
350
b591fb9c 351//-------------------------------------------------------------------
352TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
353{
354 // particleSpecies: -1 all particles are returned
355 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
356
357 Int_t nTracks = NParticles(obj);
358 TObjArray* tracksReconstructed = new TObjArray;
359 TObjArray* tracksOriginal = new TObjArray;
360 TObjArray* tracksFake = new TObjArray;
361
362 // for TPC only tracks
363 Bool_t hasOwnership = kFALSE;
8b43ff78 364 if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
b591fb9c 365 hasOwnership = kTRUE;
366
367 tracksReconstructed->SetOwner(hasOwnership);
368 tracksFake->SetOwner(hasOwnership);
369
370 // Loop over tracks or jets
371 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
372 AliVParticle* partReconstructed = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
373 if (!partReconstructed) continue;
374
375 if (useEtaPtCuts)
376 if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || partReconstructed->Pt() < fTrackPtMin)
377 {
378 if (hasOwnership)
379 delete partReconstructed;
380 continue;
381 }
382
383 Int_t label = partReconstructed->GetLabel();
384 if (label == 0)
385 {
d728f490 386 /*
387 Printf(">>> TPC only track:");
388 partReconstructed->Print();
389 partReconstructed->Dump();
390 Printf(">>> Global track:");
391 ((AliESDEvent*) obj)->GetTrack(ipart)->Dump();
392 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());
393 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());
394 */
b591fb9c 395 tracksFake->AddLast(partReconstructed);
396 continue;
397 }
d728f490 398
b591fb9c 399 AliVParticle* partOriginal = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
d728f490 400 if (!partOriginal)
401 {
402 if (hasOwnership)
403 delete partReconstructed;
404 continue;
405 }
b591fb9c 406
407 tracksReconstructed->AddLast(partReconstructed);
408 tracksOriginal->AddLast(partOriginal);
409 }
410 TObjArray* pairs = new TObjArray;
411 pairs->SetOwner(kTRUE);
412 pairs->Add(tracksReconstructed);
413 pairs->Add(tracksOriginal);
414 pairs->Add(tracksFake);
415 return pairs;
416}
417
a75aacd6 418//-------------------------------------------------------------------
419TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
420{
421
422 // Returns two lists of particles, one for MIN and one for MAX region
423 Double_t sumpT1 = 0.;
424 Double_t sumpT2 = 0.;
425
426 Int_t particles1 = transv1->GetEntries();
427 Int_t particles2 = transv2->GetEntries();
428
429 // Loop on transverse region 1
430 for (Int_t i=0; i<particles1; i++){
431 AliVParticle *part = (AliVParticle*)transv1->At(i);
432 sumpT1 += part->Pt();
433 }
434
435 // Loop on transverse region 2
436 for (Int_t i=0; i<particles2; i++){
437 AliVParticle *part = (AliVParticle*)transv2->At(i);
438 sumpT2 += part->Pt();
439 }
440
441 TObjArray *regionParticles = new TObjArray;
442 if ( sumpT2 >= sumpT1 ){
443 regionParticles->AddLast(transv1); // MIN
444 regionParticles->AddLast(transv2); // MAX
445 }else {
446 regionParticles->AddLast(transv2); // MIN
447 regionParticles->AddLast(transv1); // MAX
448 }
449
450 return regionParticles;
451}
452
453//-------------------------------------------------------------------
454Int_t AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
455{
456
457 //Returns the number of particles in AliAODMCParticle array or AliAODTracks or AliESDTracks
458
459 Int_t nTracks;
460
461 if (obj->InheritsFrom("TClonesArray")){ // MC particles
ea1919ac 462 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
a75aacd6 463 nTracks = arrayMC->GetEntriesFast();
464 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
ea1919ac 465 TObjArray *array = static_cast<TObjArray*>(obj);
a75aacd6 466 nTracks = array->GetEntriesFast();
467 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD tracks
ea1919ac 468 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
a75aacd6 469 nTracks = aodEvent->GetNTracks();
470 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD tracks
ea1919ac 471 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
a75aacd6 472 nTracks = esdEvent->GetNumberOfTracks();
3712ba26 473 }else if (obj->InheritsFrom("AliMCEvent")){ // RECO ESD tracks
ea1919ac 474 AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
3712ba26 475 nTracks = mcEvent->GetNumberOfTracks();
a75aacd6 476 }else {
477 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
478 return 0;
479 }
480
481 return nTracks;
482}
483
484
485//-------------------------------------------------------------------
b1831bcb 486AliVParticle* AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
a75aacd6 487{
488 // Returns track or MC particle at position "ipart" if passes selection criteria
b1831bcb 489 // particleSpecies: -1 all particles are returned
490 // 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
a75aacd6 491 AliVParticle *part=0;
492
493 if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
ea1919ac 494 TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
a75aacd6 495 part = (AliVParticle*)arrayMC->At( ipart );
496 if (!part)return 0;
497 // eventually only primaries
498 if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
499 // eventually only hadrons
500 if (fOnlyHadrons){
501 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
502 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
503 TMath::Abs(pdgCode)==2212 || // Proton
504 TMath::Abs(pdgCode)==321; // Kaon
505 if (!isHadron) return 0;
506 }
b1831bcb 507 if (particleSpecies != -1) {
508 // find the primary mother
509 AliVParticle* mother = part;
510 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
511 {
512 if (((AliAODMCParticle*)mother)->GetMother() < 0)
513 {
514 mother = 0;
515 break;
516 }
517
518 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
519 if (!mother)
520 break;
521 }
522
523 if (mother)
524 {
525 Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
526 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
527 return 0;
528 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
529 return 0;
530 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
531 return 0;
532 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
533 return 0;
534 }
1c6b6c97 535 else
536 {
537 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
ea1919ac 538 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
539 part->Print();
540
541 /*
542 // this code prints the details of the mother that is missing in the AOD
543 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
544
545 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
546
547 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
548 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
549 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
550 */
551
1c6b6c97 552 if (particleSpecies != 3)
553 return 0;
554 }
b1831bcb 555 }
a75aacd6 556
557 }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
ea1919ac 558 TObjArray *array = static_cast<TObjArray*>(obj);
a75aacd6 559 part = (AliVParticle*)array->At( ipart );
560 if (!part)return 0;
561 }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
ea1919ac 562 AliMCEvent* mcEvent = static_cast<AliMCEvent*>(obj);
a75aacd6 563 part = mcEvent->GetTrack( ipart );
564 if (!part) return 0;
565 // eventually only primaries
566 if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
567 // eventually only hadrons
568 //TO-DO
569 /*if (fOnlyHadrons){
570 Int_t pdgCode = part->GetPdgCode();
571 Bool_t isHadron = TMath::Abs(pdgCode)==211 || // Pion
572 TMath::Abs(pdgCode)==2212 || // Proton
573 TMath::Abs(pdgCode)==321; // Kaon
574 if (!isHadron) return 0;
575 }
576 */
eed401dc 577 if (particleSpecies != -1) {
578 // find the primary mother
579 AliMCParticle* mother = (AliMCParticle*) part;
580// Printf("");
581 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
582 {
583// Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
584 if (mother->GetMother() < 0)
585 {
586 mother = 0;
587 break;
588 }
589
590 mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
591 if (!mother)
592 break;
593 }
594
595 if (mother)
596 {
597 Int_t pdgCode = mother->PdgCode();
598 if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
599 return 0;
600 if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
601 return 0;
602 if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
603 return 0;
604 if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
605 return 0;
606 }
607 else
608 {
609 // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
610 Printf("WARNING: No mother found for particle %d:", part->GetLabel());
611 //part->Dump();
612 //part->Print();
613
614 /*
615 // this code prints the details of the mother that is missing in the AOD
616 AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
617
618 AliMCEvent* fMcEvent = fMcHandler->MCEvent();
619
620 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
621 fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
622 Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
623 */
624
625 if (particleSpecies != 3)
626 return 0;
627 }
628 }
a75aacd6 629 }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
ea1919ac 630 AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
a75aacd6 631 part = aodEvent->GetTrack(ipart);
632 // track selection cuts
144bd037 633 if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
634 //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0;
a75aacd6 635 // only primary candidates
636 //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
637 // eventually only hadrons
638 if (fOnlyHadrons){
639 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
640 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
641 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
642 if (!isHadron) return 0;
643 }
644
645 }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
ea1919ac 646 AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
a75aacd6 647 part = esdEvent->GetTrack(ipart);
648 if (!part)return 0;
649 // track selection cuts
e0331fd9 650
651 if (!( ApplyCuts(part)) )
652 return 0;
a75aacd6 653
1bba939a 654 if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
85bfac17 655 {
656 // create TPC only tracks constrained to the SPD vertex
657
658 const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
659
660 AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
661 if(!track) return 0;
662
85bfac17 663 if(track->Pt()>0.){
664 // only constrain tracks above threshold
665 AliExternalTrackParam exParam;
666 // take the B-feild from the ESD, no 3D fieldMap available at this point
667 Bool_t relate = kFALSE;
caa8c7b2 668 relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
85bfac17 669 if(!relate)
670 {
671// Printf("relating failed");
672 delete track;
673 return 0;
674 }
675 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
676 }
677
678 part = track;
679 }
8b43ff78 680 else if (fFilterBit == 2048)
681 {
682 // hybrid tracks
683
684 // clone
685 AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
686// Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
687
688 if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
689 {
690// Float_t ptBefore = esdTrack->Pt();
691 // set constrained pT as default pT
692 if (!esdTrack->GetConstrainedParam())
693 return 0;
694 esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
695// Printf("%f %f", ptBefore, esdTrack->Pt());
696 }
697 part = esdTrack;
698 }
a75aacd6 699
700 // eventually only hadrons
701 //TO-DO
702 /*if (fOnlyHadrons){
703 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
704 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
705 ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
706 if (!isHadron) return 0;
707 }
708 */
709 }else {
710 if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
711 return 0;
712 }
713
714 // only charged
715 if (!part->Charge())return 0;
716
717 return part;
718}
719
720
721//-------------------------------------------------------------------
722void AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
723{
724 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
725
726 static TObject *tmp;
727 static int i; // "static" to save stack space
728 int j;
729
730 while (last - first > 1) {
731 i = first;
732 j = last;
733 for (;;) {
734 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
735 ;
736 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
737 ;
738 if (i >= j)
739 break;
740
741 tmp = a[i];
742 a[i] = a[j];
743 a[j] = tmp;
744 }
745 if (j == first) {
746 ++first;
747 continue;
748 }
749 tmp = a[first];
750 a[first] = a[j];
751 a[j] = tmp;
752 if (j - first < last - (j + 1)) {
753 QSortTracks(a, first, j);
754 first = j + 1; // QSortTracks(j + 1, last);
755 } else {
756 QSortTracks(a, j + 1, last);
757 last = j; // QSortTracks(first, j);
758 }
759 }
760}
761
762//____________________________________________________________________
763TObjArray* AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
764{
765
766 // Assign particles to towards, away or transverse regions.
767 // Returns a lists of particles for each region.
768
769 static const Double_t k60rad = 60.*TMath::Pi()/180.;
770 static const Double_t k120rad = 120.*TMath::Pi()/180.;
771
772 // Define output lists of particles
773 TList *toward = new TList();
774 TList *away = new TList();
775 // Two transverse regions, for the moment those are not yet MIN and MAX!!!
776 // MIN and MAX can be sorted in GetMinMaxRegion function
777 TList *transverse1 = new TList();
778 TList *transverse2 = new TList();
779
780 TObjArray *regionParticles = new TObjArray;
05c47aff 781 regionParticles->SetOwner(kTRUE);
782
a75aacd6 783 regionParticles->AddLast(toward);
784 regionParticles->AddLast(away);
785 regionParticles->AddLast(transverse1);
786 regionParticles->AddLast(transverse2);
787
788 if (!leading)
789 return regionParticles;
790
791 // Switch to vector for leading particle
792 TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
793
794 Int_t nTracks = NParticles(obj);
795 if( !nTracks ) return 0;
796 // Loop over tracks
797 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
798 AliVParticle* part = ParticleWithCuts(obj, ipart);
799 if (!part)continue;
800 //Switch to vectors for particles
801 TVector3 partVect(part->Px(), part->Py(), part->Pz());
802
803 Int_t region = 0;
804 if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
805 // transverse regions
806 if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
807 if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1; //right
808
809 if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2; //forward
810 if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2; //backward
811
812 // skip leading particle
b1831bcb 813 if (leading == part)
a75aacd6 814 continue;
815
816 if (!region)continue;
817 if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
818 Int_t label = ((AliAODTrack*)part)->GetLabel();
819 // re-define part as the matched MC particle
820 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
821 if (!part)continue;
822 // skip leading particle
823 if (leading == part)
824 continue;
825 }
826 if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
827 Int_t label = ((AliESDtrack*)part)->GetLabel();
828 // look for the matched MC particle (but do not re-define part)
829 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
830 }
831
832 if ( region == 1 ) transverse1->Add(part);
833 if ( region == -1 ) transverse2->Add(part);
834 if ( region == 2 ) toward->Add(part);
835 if ( region == -2 ) away->Add(part);
836
837 }//end loop on tracks
838
839 return regionParticles;
840
841}
842
843
844//____________________________________________________________________
845Bool_t AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
846{
1c6b6c97 847 if (!obj) // MC
848 return kFALSE;
a75aacd6 849
85bfac17 850 // Use AliPhysicsSelection to select good events, works for ESD and AOD
5cabd2cd 851 if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
85bfac17 852 return kFALSE;
05c47aff 853
a75aacd6 854 return kTRUE;
a75aacd6 855}
856
857//____________________________________________________________________
858Bool_t AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
859{
860
861 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
862
863 if (obj->InheritsFrom("AliAODEvent")){
864 Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
865 if( nVertex > 0 ) {
866 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
144bd037 867 Int_t nTracksPrim = vertex->GetNContributors();
a75aacd6 868 Double_t zVertex = vertex->GetZ();
869 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
144bd037 870 // Reject TPC only vertex
871 TString name(vertex->GetName());
872 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
873
874 // Select a quality vertex by number of tracks?
a75aacd6 875 if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
876 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
877 return kFALSE;
878 }
879 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
880 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
881 // return kFALSE;
882 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
883 } else {
884 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
885 return kFALSE;
886 }
887 }
888
889 if (obj->InheritsFrom("AliMCEvent"))
890 {
891 if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
892 {
893 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
894 return kFALSE;
895 }
896 }
144bd037 897
7dd4dec4 898 if (obj->InheritsFrom("AliAODMCHeader"))
899 {
900 if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) > zed)
901 {
902 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
903 return kFALSE;
904 }
905 }
906
144bd037 907 // ESD case for DCA studies
908 if (obj->InheritsFrom("AliESDEvent")){
909 AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
910 if ( vertex){
911 Int_t nTracksPrim = vertex->GetNContributors();
912 Double_t zVertex = vertex->GetZ();
913 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
914 // Reject SPD or TPC only vertex
915 TString name(vertex->GetName());
916 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
917
918 // Select a quality vertex by number of tracks?
919 if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
920 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
921 return kFALSE;
922 }
923 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
924 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
925 // return kFALSE;
926 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
927 } else {
928 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
929 return kFALSE;
930 }
931 }
a75aacd6 932
933 return kTRUE;
934}