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