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