PWGCF/Correlations/Base/AliAnalyseLeadingTrackUE.cxx

   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  **************************************************************************/
  15
  16 #include <TList.h>
  17 #include <TMath.h>
  18 #include <TObjArray.h>
  19 #include <TObject.h>
  20 #include <TVector3.h>
  21
  22 #include "AliAnalyseLeadingTrackUE.h"
  23
  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"
  34 #include "TFormula.h"
  35
  36 #include "AliAnalysisManager.h"
  37 #include "AliMCEventHandler.h"
  38 #include "AliStack.h"
  39 #include "AliPIDResponse.h"
  40 #include "AliHelperPID.h"
  41
  42
  43 ////////////////////////////////////////////////
  44 //---------------------------------------------
  45 // Class for transverse regions analysis
  46 //---------------------------------------------
  47 ////////////////////////////////////////////////
  48
  49
  50 using namespace std;
  51
  52 ClassImp(AliAnalyseLeadingTrackUE)
  53
  54 //-------------------------------------------------------------------
  55 AliAnalyseLeadingTrackUE::AliAnalyseLeadingTrackUE() :
  56   TObject(),
  57   fDebug(0),
  58   fFilterBit(16),
  59   fTrackStatus(0),
  60   fOnlyHadrons(kFALSE),
  61   fCheckMotherPDG(kTRUE),
  62   fTrackEtaCut(0.8),
  63   fTrackEtaCutMin(-1.),
  64   fTrackPhiCutEvPlMin(0.),
  65   fTrackPhiCutEvPlMax(0.),
  66   fTrackPtMin(0),
  67   fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
  68   fDCAXYCut(0),
  69   fSharedClusterCut(-1),
  70   fCrossedRowsCut(-1),
  71   fFoundFractionCut(-1),
  72   fEsdTrackCuts(0x0),
  73   fEsdTrackCutsExtra1(0x0),
  74   fEsdTrackCutsExtra2(0x0),
  75   fHelperPID(0x0),
  76   fEventCounter(0)
  77 {
  78   // constructor
  79 }
  80
  81
  82 //-------------------------------------------------------------------
  83 AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
  84 {
  85   // assignment operator
  86   return *this;
  87 }
  88
  89
  90 //-------------------------------------------------------------------
  91 AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
  92 {
  93
  94   //clear memory
  95
  96 }
  97
  98
  99 //____________________________________________________________________
 100 Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
 101 {
 102
 103   // select track according to set of cuts
 104   if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
 105   if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->IsSelected(track)) return kFALSE;
 106
 107   return kTRUE;
 108 }
 109
 110
 111 //____________________________________________________________________
 112 void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
 113
 114   // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
 115   // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
 116
 117   if (filterbit == -1)
 118     filterbit = fFilterBit;
 119
 120   if (filterbit == 128)
 121   {
 122     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 123     fEsdTrackCuts->SetMinNClustersTPC(70);
 124   }
 125   else if (filterbit == 256)
 126   {
 127     // syst study
 128     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 129     fEsdTrackCuts->SetMinNClustersTPC(80);
 130     fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
 131     fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
 132     fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
 133   }
 134   else if (filterbit == 512)
 135   {
 136     // syst study
 137     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 138     fEsdTrackCuts->SetMinNClustersTPC(60);
 139     fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
 140     fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
 141     fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
 142   }
 143   else if (filterbit == 1024)
 144   {
 145     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 146     fEsdTrackCuts->SetMinNClustersTPC(-1);
 147     fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
 148     fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 149   }
 150   else if (filterbit == 2048) // mimic hybrid tracks
 151   {
 152     // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
 153     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 154     fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
 155     fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
 156     fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
 157     fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 158     fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 159     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 160     fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
 161
 162     // Add SPD requirement
 163     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 164     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 165     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
 166
 167     fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
 168     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
 169     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
 170   }
 171   else if (filterbit == 4096) // require TOF matching
 172   {
 173     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
 174     // FIXME: TOF REQUIREMENTS ARE IN GetParticleSpecies FOR THE MOMENT
 175   }
 176   else
 177   {
 178     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 179     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 180
 181     // Add SPD requirement
 182     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 183     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 184
 185     // Add SDD requirement
 186     fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
 187     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
 188   }
 189 }
 190
 191 //____________________________________________________________________
 192 TObjArray*  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 || TMath::Abs(part->Eta()) < fTrackEtaCutMin ) 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
 223 void 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     mcObj->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 //     Printf("%d %d %d", i, part->GetLabel(), mother->GetLabel());
 292     if (mother->GetLabel() >= maxLabel)
 293     {
 294 //       Printf("Removing %d with label %d", i, part->GetLabel()); ((AliMCParticle*)part)->Particle()->Print(); ((AliMCParticle*)mother)->Particle()->Print();
 295       TObject* object = tracks->RemoveAt(i);
 296       if (tracks->IsOwner())
 297         delete object;
 298     }
 299   }
 300
 301   tracks->Compress();
 302
 303   AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
 304 }
 305
 306 //-------------------------------------------------------------------
 307 void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
 308 {
 309   // remove particles from weak decays
 310   // <tracks> can be the following cases:
 311   // a. tracks: in this case the label is taken and then case b.
 312   // b. particles: it is checked if IsSecondaryFromWeakDecay is true
 313   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
 314
 315   TClonesArray* arrayMC = 0;
 316   AliMCEvent* mcEvent = 0;
 317   if (mcObj->InheritsFrom("AliMCEvent"))
 318     mcEvent = static_cast<AliMCEvent*>(mcObj);
 319   else if (mcObj->InheritsFrom("TClonesArray"))
 320     arrayMC = static_cast<TClonesArray*>(mcObj);
 321   else
 322   {
 323     mcObj->Dump();
 324     AliFatal("Invalid object passed");
 325   }
 326
 327   Int_t before = tracks->GetEntriesFast();
 328
 329   for (Int_t i=0; i<before; ++i)
 330   {
 331     AliVParticle* part = (AliVParticle*) tracks->At(i);
 332
 333     if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
 334       part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
 335
 336     if (part->InheritsFrom("AliAODMCParticle"))
 337     {
 338       if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
 339         continue;
 340     }
 341     else if (part->InheritsFrom("AliMCParticle") && mcEvent)
 342     {
 343       if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
 344         continue;
 345     }
 346     else
 347     {
 348       part->Dump();
 349       AliFatal("Unknown particle");
 350     }
 351
 352 //     Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
 353     TObject* object = tracks->RemoveAt(i);
 354     if (tracks->IsOwner())
 355       delete object;
 356   }
 357
 358   tracks->Compress();
 359
 360   if (before > tracks->GetEntriesFast())
 361     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
 362 }
 363
 364 //-------------------------------------------------------------------
 365 TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts, Bool_t speciesOnTracks, Double_t evtPlane)
 366 {
 367   // Returns an array of particles that pass the cuts, if arrayMC is given each reconstructed particle is replaced by its corresponding MC particles, depending on the parameter onlyprimaries only for primaries
 368   // particleSpecies: -1 all particles are returned
 369   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 370   // speciesOnTracks if kFALSE, particleSpecies is only applied on the matched MC particle (not on the track itself)
 371   // Passing down the Double_t* evtPlane (range [-pi/2,pi/2]) will apply a phi cut with respect to the eventplane between fTrackPhiCutEvPlMin and fTrackPhiCutEvPlMax. For values outside [-pi/2,pi/2], this will be ignored.
 372
 373   Int_t nTracks = NParticles(obj);
 374   TObjArray* tracks = new TObjArray;
 375
 376   // for TPC only tracks
 377   Bool_t hasOwnership = kFALSE;
 378   if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
 379     hasOwnership = kTRUE;
 380
 381   if (!arrayMC)
 382     tracks->SetOwner(hasOwnership);
 383
 384   // Loop over tracks or jets
 385   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 386     AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, (speciesOnTracks) ? particleSpecies : -1);
 387     if (!part) continue;
 388
 389     if (TMath::Abs(evtPlane)<=TMath::Pi()/2) { //evtPlane range: (-pi/2,pi/2)
 390       Double_t phiPart = part->Phi(); //range: [0,2*pi)
 391       if(phiPart>TMath::Pi()) phiPart-=2*TMath::Pi();
 392
 393       Double_t dPhi = 0; //range: [0,pi/2], i.e. the difference over the shortest angle.
 394       Double_t diff = TMath::Abs(phiPart-evtPlane);
 395       if(diff<=TMath::Pi()/2) dPhi = diff;
 396       else if(diff<=TMath::Pi()) dPhi = TMath::Pi()-diff;
 397       else dPhi = diff-TMath::Pi();
 398
 399       if(dPhi<fTrackPhiCutEvPlMin || dPhi>fTrackPhiCutEvPlMax) {
 400         if (hasOwnership)
 401           delete part;
 402         continue;
 403       }
 404     }
 405
 406     if (useEtaPtCuts)
 407       if (TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin || part->Pt() < fTrackPtMin)
 408       {
 409         if (hasOwnership)
 410           delete part;
 411         continue;
 412       }
 413
 414 //     Printf("%p %p %d Accepted %d %f %f %f", obj, arrayMC, particleSpecies, ipart, part->Eta(), part->Phi(), part->Pt());
 415
 416     if (arrayMC) {
 417       Int_t label = part->GetLabel();
 418       if (hasOwnership)
 419         delete part;
 420       // re-define part as the matched MC particle
 421       part = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
 422       if (!part)continue;
 423     }
 424
 425     tracks->Add(part);
 426   }
 427
 428   return tracks;
 429 }
 430
 431 //-------------------------------------------------------------------
 432 TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
 433 {
 434   // particleSpecies: -1 all particles are returned
 435   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 436
 437   Int_t nTracks = NParticles(obj);
 438   TObjArray* tracksReconstructed = new TObjArray;
 439   TObjArray* tracksOriginal = new TObjArray;
 440   TObjArray* tracksFake = new TObjArray;
 441
 442   // for TPC only tracks
 443   Bool_t hasOwnership = kFALSE;
 444   if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
 445     hasOwnership = kTRUE;
 446
 447   tracksReconstructed->SetOwner(hasOwnership);
 448   tracksFake->SetOwner(hasOwnership);
 449
 450   // Loop over tracks or jets
 451   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 452     AliVParticle* partReconstructed = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
 453     if (!partReconstructed) continue;
 454
 455     if (useEtaPtCuts)
 456       if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || TMath::Abs(partReconstructed->Eta()) < fTrackEtaCutMin || partReconstructed->Pt() < fTrackPtMin)
 457       {
 458         if (hasOwnership)
 459           delete partReconstructed;
 460         continue;
 461       }
 462
 463     Int_t label = partReconstructed->GetLabel();
 464     if (label == 0)
 465     {
 466       /*
 467       Printf(">>> TPC only track:");
 468       partReconstructed->Print();
 469       partReconstructed->Dump();
 470       Printf(">>> Global track:");
 471       ((AliESDEvent*) obj)->GetTrack(ipart)->Dump();
 472       Printf("Fake (TPC only): eta = %f, phi = %f, pT = %f, ncl = %d, dedx = %f", partReconstructed->Eta(), partReconstructed->Phi(), partReconstructed->Pt(), ((AliESDtrack*) partReconstructed)->GetTPCclusters(0), ((AliESDtrack*) partReconstructed)->GetTPCsignal());
 473       Printf("Fake (global  ): eta = %f, phi = %f, pT = %f, ncl = %d, dedx = %f", ((AliESDEvent*) obj)->GetTrack(ipart)->Eta(), ((AliESDEvent*) obj)->GetTrack(ipart)->Phi(), ((AliESDEvent*) obj)->GetTrack(ipart)->Pt(), ((AliESDEvent*) obj)->GetTrack(ipart)->GetTPCclusters(0), ((AliESDEvent*) obj)->GetTrack(ipart)->GetTPCsignal());
 474       */
 475       tracksFake->AddLast(partReconstructed);
 476       continue;
 477     }
 478
 479     AliVParticle* partOriginal = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
 480     if (!partOriginal)
 481     {
 482       if (hasOwnership)
 483         delete partReconstructed;
 484       continue;
 485     }
 486
 487     tracksReconstructed->AddLast(partReconstructed);
 488     tracksOriginal->AddLast(partOriginal);
 489   }
 490   TObjArray* pairs = new TObjArray;
 491   pairs->SetOwner(kTRUE);
 492   pairs->Add(tracksReconstructed);
 493   pairs->Add(tracksOriginal);
 494   pairs->Add(tracksFake);
 495   return pairs;
 496 }
 497
 498 //-------------------------------------------------------------------
 499 TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
 500 {
 501
 502  // Returns two lists of particles, one for MIN and one for MAX region
 503   Double_t sumpT1 = 0.;
 504   Double_t sumpT2 = 0.;
 505
 506   Int_t particles1 = transv1->GetEntries();
 507   Int_t particles2 = transv2->GetEntries();
 508
 509   // Loop on transverse region 1
 510   for (Int_t i=0; i<particles1; i++){
 511         AliVParticle *part = (AliVParticle*)transv1->At(i);
 512         sumpT1 +=  part->Pt();
 513         }
 514
 515   // Loop on transverse region 2
 516   for (Int_t i=0; i<particles2; i++){
 517         AliVParticle *part = (AliVParticle*)transv2->At(i);
 518         sumpT2 +=  part->Pt();
 519         }
 520
 521   TObjArray *regionParticles = new TObjArray;
 522   if ( sumpT2 >= sumpT1 ){
 523         regionParticles->AddLast(transv1); // MIN
 524         regionParticles->AddLast(transv2); // MAX
 525   }else {
 526         regionParticles->AddLast(transv2); // MIN
 527         regionParticles->AddLast(transv1); // MAX
 528         }
 529
 530  return regionParticles;
 531 }
 532
 533 //-------------------------------------------------------------------
 534 Int_t  AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
 535 {
 536
 537   //Returns the number of particles in AliAODMCParticle array  or AliAODTracks or AliESDTracks
 538
 539   Int_t nTracks;
 540
 541   if (obj->InheritsFrom("TClonesArray")){ // MC particles
 542         TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
 543         nTracks = arrayMC->GetEntriesFast();
 544   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 545         TObjArray *array = static_cast<TObjArray*>(obj);
 546         nTracks = array->GetEntriesFast();
 547   }else if (obj->InheritsFrom("AliAODEvent")){  // RECO AOD tracks
 548         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 549         nTracks = aodEvent->GetNumberOfTracks();
 550   }else if (obj->InheritsFrom("AliESDEvent")){  // RECO ESD tracks
 551         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 552         nTracks = esdEvent->GetNumberOfTracks();
 553   }else if (obj->InheritsFrom("AliMCEvent")){  // RECO ESD tracks
 554         AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
 555         nTracks = mcEvent->GetNumberOfTracks();
 556   }else {
 557         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 558         return 0;
 559         }
 560
 561   return nTracks;
 562 }
 563
 564 //-------------------------------------------------------------------
 565 AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
 566 {
 567   // Returns track or MC particle at position "ipart" if passes selection criteria
 568   // particleSpecies: -1 all particles are returned
 569   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 570   AliVParticle *part=0;
 571
 572   if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
 573         TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
 574         part = (AliVParticle*)arrayMC->At( ipart );
 575         if (!part)return 0;
 576         // eventually only primaries
 577         if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
 578         // eventually only hadrons
 579         if (fOnlyHadrons){
 580                 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
 581                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 582                                   TMath::Abs(pdgCode)==2212 || // Proton
 583                                   TMath::Abs(pdgCode)==321;    // Kaon
 584                 if (!isHadron) return 0;
 585                 }
 586         if (particleSpecies != -1) {
 587                 // find the primary mother
 588                 AliVParticle* mother = part;
 589                 if(fCheckMotherPDG) {
 590                   while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 591                   {
 592                   if (((AliAODMCParticle*)mother)->GetMother() < 0)
 593                   {
 594                     mother = 0;
 595                     break;
 596                   }
 597
 598                   mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 599                   if (!mother)
 600                     break;
 601                  }
 602                 }
 603                 if (mother)
 604                 {
 605                   Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
 606                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 607                           return 0;
 608                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 609                           return 0;
 610                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 611                           return 0;
 612                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 613                           return 0;
 614                 }
 615                 else
 616                 {
 617                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 618                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 619                   part->Print();
 620
 621                   /*
 622                   // this code prints the details of the mother that is missing in the AOD
 623                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 624
 625                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 626
 627                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 628                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 629                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 630                   */
 631
 632                   if (particleSpecies != 3)
 633                     return 0;
 634                 }
 635         }
 636
 637   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 638         TObjArray *array = static_cast<TObjArray*>(obj);
 639         part = (AliVParticle*)array->At( ipart );
 640         if (!part)return 0;
 641   }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
 642         AliMCEvent* mcEvent =  static_cast<AliMCEvent*>(obj);
 643         part = mcEvent->GetTrack( ipart );
 644         if (!part) return 0;
 645         // eventually only primaries
 646         if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
 647         // eventually only hadrons
 648         //TO-DO
 649         /*if (fOnlyHadrons){
 650                 Int_t pdgCode =  part->GetPdgCode();
 651                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 652                                   TMath::Abs(pdgCode)==2212 || // Proton
 653                                   TMath::Abs(pdgCode)==321;    // Kaon
 654                 if (!isHadron) return 0;
 655                 }
 656         */
 657        if (particleSpecies != -1) {
 658                 // find the primary mother
 659                 AliMCParticle* mother = (AliMCParticle*) part;
 660 //              Printf("");
 661                 if(fCheckMotherPDG) {
 662                   while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
 663                   {
 664 //                Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
 665                   if (mother->GetMother() < 0)
 666                   {
 667                     mother = 0;
 668                     break;
 669                   }
 670
 671                   mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
 672                   if (!mother)
 673                     break;
 674                   }
 675                 }
 676                 if (mother)
 677                 {
 678                   Int_t pdgCode = mother->PdgCode();
 679                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 680                           return 0;
 681                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 682                           return 0;
 683                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 684                           return 0;
 685                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 686                           return 0;
 687                 }
 688                 else
 689                 {
 690                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 691                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 692                   //part->Dump();
 693                   //part->Print();
 694
 695                   /*
 696                   // this code prints the details of the mother that is missing in the AOD
 697                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 698
 699                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 700
 701                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 702                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 703                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 704                   */
 705
 706                   if (particleSpecies != 3)
 707                     return 0;
 708                 }
 709         }
 710   }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
 711         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 712         part = aodEvent->GetTrack(ipart);
 713
 714         // track selection cuts
 715         if (fFilterBit != 0 && !(((AliAODTrack*)part)->TestFilterBit(fFilterBit))) return 0;
 716         if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
 717
 718         // DCA XY
 719         if (fDCAXYCut)
 720         {
 721           const AliVVertex* vertex = aodEvent->GetPrimaryVertex();
 722           if (!vertex)
 723             return 0;
 724
 725           Double_t pos[2];
 726           Double_t covar[3];
 727           AliAODTrack* clone = (AliAODTrack*) part->Clone();
 728           Bool_t success = clone->PropagateToDCA(vertex, ((AliVAODHeader*)aodEvent->GetHeader())->GetMagneticField(), 3, pos, covar);
 729           delete clone;
 730           if (!success)
 731             return 0;
 732
 733 //        Printf("%f", ((AliAODTrack*)part)->DCA());
 734 //        Printf("%f", pos[0]);
 735           if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(part->Pt()))
 736             return 0;
 737         }
 738
 739         if (fSharedClusterCut >= 0)
 740         {
 741           Double_t frac = Double_t(((AliAODTrack*)part)->GetTPCnclsS()) / Double_t(((AliAODTrack*)part)->GetTPCncls());
 742           if (frac > fSharedClusterCut)
 743             return 0;
 744         }
 745
 746         if (fCrossedRowsCut >= 0)
 747         {
 748           if (((AliAODTrack*) part)->GetTPCNCrossedRows() < fCrossedRowsCut)
 749             return 0;
 750         }
 751
 752         if (fFoundFractionCut >= 0)
 753         {
 754           UInt_t findableClusters = ((AliAODTrack*) part)->GetTPCNclsF();
 755           if (findableClusters == 0)
 756             return 0;
 757           if (((Double_t) ((AliAODTrack*) part)->GetTPCNCrossedRows() / findableClusters) < fFoundFractionCut)
 758             return 0;
 759         }
 760
 761         // eventually only hadrons
 762         if (fOnlyHadrons){
 763                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 764                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 765                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 766                 if (!isHadron) return 0;
 767                 }
 768
 769         if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0;
 770
 771   }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
 772         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 773         part = esdEvent->GetTrack(ipart);
 774         if (!part)return 0;
 775
 776         // track selection cuts
 777         if (!( ApplyCuts(part)) )
 778           return 0;
 779
 780         if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
 781
 782         if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
 783         {
 784           // create TPC only tracks constrained to the SPD vertex
 785
 786           const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
 787
 788           AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
 789           if(!track) return 0;
 790
 791 //        Printf(">%f %f %f", track->Eta(), track->Phi(), track->Pt());
 792
 793           if(track->Pt()>0.){
 794             // only constrain tracks above threshold
 795             AliExternalTrackParam exParam;
 796             // take the B-feild from the ESD, no 3D fieldMap available at this point
 797             Bool_t relate = kFALSE;
 798             relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
 799             if(!relate)
 800             {
 801 //                 Printf("relating failed");
 802               delete track;
 803               return 0;
 804             }
 805             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 806           }
 807
 808 //        Printf(">%f %f %f\n", track->Eta(), track->Phi(), track->Pt());
 809
 810           part = track;
 811         }
 812         else if (fFilterBit == 2048)
 813         {
 814           // hybrid tracks
 815
 816           // clone
 817           AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
 818 //        Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
 819
 820           if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
 821           {
 822 //          Float_t ptBefore = esdTrack->Pt();
 823             // set constrained pT as default pT
 824             if (!esdTrack->GetConstrainedParam())
 825               return 0;
 826             esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
 827 //          Printf("%f %f", ptBefore, esdTrack->Pt());
 828           }
 829           part = esdTrack;
 830         }
 831
 832         // eventually only hadrons
 833         //TO-DO
 834         /*if (fOnlyHadrons){
 835                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 836                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 837                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 838                 if (!isHadron) return 0;
 839                 }
 840         */
 841         if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0; // If it is -1 you take all the particles
 842
 843   }else {
 844         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 845         return 0;
 846         }
 847
 848   // only charged
 849   if (!part->Charge())return 0;
 850
 851   part->SetUniqueID(fEventCounter * 100000 + ipart);
 852   return part;
 853 }
 854
 855
 856 //-------------------------------------------------------------------
 857 void  AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
 858 {
 859   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
 860
 861   static TObject *tmp;
 862   static int i;           // "static" to save stack space
 863   int j;
 864
 865   while (last - first > 1) {
 866     i = first;
 867     j = last;
 868     for (;;) {
 869       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
 870         ;
 871       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
 872         ;
 873       if (i >= j)
 874         break;
 875
 876       tmp  = a[i];
 877       a[i] = a[j];
 878       a[j] = tmp;
 879     }
 880     if (j == first) {
 881       ++first;
 882       continue;
 883     }
 884     tmp = a[first];
 885     a[first] = a[j];
 886     a[j] = tmp;
 887     if (j - first < last - (j + 1)) {
 888       QSortTracks(a, first, j);
 889       first = j + 1;   // QSortTracks(j + 1, last);
 890     } else {
 891       QSortTracks(a, j + 1, last);
 892       last = j;        // QSortTracks(first, j);
 893     }
 894   }
 895 }
 896
 897 //____________________________________________________________________
 898 TObjArray*  AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
 899 {
 900
 901   // Assign particles to towards, away or transverse regions.
 902   // Returns a lists of particles for each region.
 903
 904   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 905   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 906
 907   // Define output lists of particles
 908   TList *toward = new TList();
 909   TList *away = new TList();
 910   // Two transverse regions, for the moment those are not yet MIN and MAX!!!
 911   // MIN and MAX can be sorted in GetMinMaxRegion function
 912   TList *transverse1 = new TList();
 913   TList *transverse2 = new TList();
 914
 915   TObjArray *regionParticles = new TObjArray;
 916   regionParticles->SetOwner(kTRUE);
 917
 918   regionParticles->AddLast(toward);
 919   regionParticles->AddLast(away);
 920   regionParticles->AddLast(transverse1);
 921   regionParticles->AddLast(transverse2);
 922
 923   if (!leading)
 924     return regionParticles;
 925
 926   // Switch to vector for leading particle
 927   TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
 928
 929   Int_t nTracks = NParticles(obj);
 930   if( !nTracks ) return 0;
 931   // Loop over tracks
 932   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 933         AliVParticle* part = ParticleWithCuts(obj, ipart);
 934         if (!part)continue;
 935         //Switch to vectors for particles
 936         TVector3 partVect(part->Px(), part->Py(), part->Pz());
 937
 938         Int_t region = 0;
 939         if( TMath::Abs(partVect.Eta()) > fTrackEtaCut || TMath::Abs(partVect.Eta()) < fTrackEtaCutMin) continue;
 940         // transverse regions
 941         if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
 942         if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1;   //right
 943
 944         if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2;    //forward
 945         if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2;  //backward
 946
 947         // skip leading particle
 948         if (leading == part)
 949           continue;
 950
 951         if (!region)continue;
 952         if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
 953                 Int_t label = ((AliAODTrack*)part)->GetLabel();
 954                 // re-define part as the matched MC particle
 955                 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
 956                 if (!part)continue;
 957                 // skip leading particle
 958                 if (leading == part)
 959                   continue;
 960                 }
 961         if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
 962                 Int_t label = ((AliESDtrack*)part)->GetLabel();
 963                 // look for the matched MC particle (but do not re-define part)
 964                 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
 965                 }
 966
 967         if ( region == 1 ) transverse1->Add(part);
 968         if ( region == -1 ) transverse2->Add(part);
 969         if ( region == 2 ) toward->Add(part);
 970         if ( region == -2 ) away->Add(part);
 971
 972         }//end loop on tracks
 973
 974   return regionParticles;
 975
 976 }
 977
 978
 979 //____________________________________________________________________
 980 Bool_t  AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
 981 {
 982   if (!obj) // MC
 983     return kFALSE;
 984
 985   // Use AliPhysicsSelection to select good events, works for ESD and AOD
 986   if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
 987     return kFALSE;
 988
 989   return kTRUE;
 990 }
 991
 992 //____________________________________________________________________
 993 Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
 994 {
 995
 996   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
 997
 998   if (obj->InheritsFrom("AliAODEvent")){
 999         Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
1000         if( nVertex > 0 ) {
1001                 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
1002                 Int_t nTracksPrim = vertex->GetNContributors();
1003                 Double_t zVertex = vertex->GetZ();
1004                 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1005                 // Reject TPC only vertex
1006                 TString name(vertex->GetName());
1007                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
1008
1009                 // Select a quality vertex by number of tracks?
1010                 if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
1011                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1012                         return kFALSE;
1013                         }
1014                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1015                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1016                 //  return kFALSE;
1017                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1018                 } else {
1019                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1020                         return kFALSE;
1021                         }
1022         }
1023
1024   if (obj->InheritsFrom("AliMCEvent"))
1025   {
1026     if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) >= zed)
1027     {
1028       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1029       return kFALSE;
1030     }
1031   }
1032
1033   if (obj->InheritsFrom("AliAODMCHeader"))
1034   {
1035     if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= zed)
1036     {
1037       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1038       return kFALSE;
1039     }
1040   }
1041
1042   // ESD case for DCA studies
1043   if (obj->InheritsFrom("AliESDEvent")){
1044        AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
1045        if ( vertex){
1046                Int_t nTracksPrim = vertex->GetNContributors();
1047                Double_t zVertex = vertex->GetZ();
1048                if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1049                // Reject SPD or TPC only vertex
1050                TString name(vertex->GetName());
1051                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
1052
1053                // Select a quality vertex by number of tracks?
1054                if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
1055                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1056                        return kFALSE;
1057                        }
1058                // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1059                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1060                 //  return kFALSE;
1061                if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1062                } else {
1063                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1064                        return kFALSE;
1065                        }
1066        }
1067
1068   return kTRUE;
1069 }
1070
1071 //____________________________________________________________________
1072
1073 Bool_t AliAnalyseLeadingTrackUE::CheckTrack(AliVParticle * part)
1074 {
1075   // check if the track status flags are set
1076
1077   UInt_t status=((AliVTrack*)part)->GetStatus();
1078   if ((status & fTrackStatus) == fTrackStatus)
1079     return kTRUE;
1080   return kFALSE;
1081 }