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