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   fEsdTrackCutsExtra1(0x0),
  81   fEsdTrackCutsExtra2(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 (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->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 if (filterbit == 2048) // mimic hybrid tracks
 156   {
 157     // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
 158     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 159     fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
 160     fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
 161     fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
 162     fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 163     fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 164     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 165
 166     // Add SPD requirement
 167     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 168     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 169     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
 170
 171     fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
 172     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
 173     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
 174   }
 175   else
 176   {
 177     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 178     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 179
 180     // Add SPD requirement
 181     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 182     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 183
 184     // Add SDD requirement
 185     fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
 186     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
 187   }
 188 }
 189
 190 //____________________________________________________________________
 191 TObjArray*  AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
 192 {
 193
 194   // Returns an array of charged particles (or jets) ordered according to their pT.
 195
 196   Int_t nTracks = NParticles(obj);
 197
 198
 199   if( !nTracks ) return 0;
 200
 201   // Define array of AliVParticle objects
 202   TObjArray* tracks = new TObjArray(nTracks);
 203
 204   // Loop over tracks or jets
 205   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 206         AliVParticle* part = ParticleWithCuts( obj, ipart );
 207         if (!part) continue;
 208         // Accept leading-tracks in a limited pseudo-rapidity range
 209         if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
 210         tracks->AddLast( part );
 211         }
 212   // Order tracks by pT
 213   QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
 214
 215   nTracks = tracks->GetEntriesFast();
 216   if( !nTracks ) return 0;
 217
 218   return tracks;
 219   }
 220
 221
 222 void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
 223 {
 224   // remove injected signals (primaries above <maxLabel>)
 225   // <tracks> can be the following cases:
 226   // a. tracks: in this case the label is taken and then case b.
 227   // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
 228   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
 229
 230   TClonesArray* arrayMC = 0;
 231   AliMCEvent* mcEvent = 0;
 232   if (mcObj->InheritsFrom("AliMCEvent"))
 233     mcEvent = static_cast<AliMCEvent*>(mcObj);
 234   else if (mcObj->InheritsFrom("TClonesArray"))
 235     arrayMC = static_cast<TClonesArray*>(mcObj);
 236   else
 237   {
 238     arrayMC->Dump();
 239     AliFatal("Invalid object passed");
 240   }
 241
 242   Int_t before = tracks->GetEntriesFast();
 243
 244   for (Int_t i=0; i<before; ++i)
 245   {
 246     AliVParticle* part = (AliVParticle*) tracks->At(i);
 247
 248     if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
 249       part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
 250
 251     AliVParticle* mother = part;
 252     if (mcEvent)
 253     {
 254       while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
 255       {
 256         if (((AliMCParticle*)mother)->GetMother() < 0)
 257         {
 258           mother = 0;
 259           break;
 260         }
 261
 262         mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
 263         if (!mother)
 264           break;
 265       }
 266     }
 267     else
 268     {
 269       // find the primary mother
 270       while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 271       {
 272         if (((AliAODMCParticle*)mother)->GetMother() < 0)
 273         {
 274           mother = 0;
 275           break;
 276         }
 277
 278         mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 279         if (!mother)
 280           break;
 281       }
 282     }
 283
 284     if (!mother)
 285     {
 286       Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 287       continue;
 288     }
 289
 290     if (mother->GetLabel() > maxLabel)
 291     {
 292 //       Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
 293       TObject* object = tracks->RemoveAt(i);
 294       if (tracks->IsOwner())
 295         delete object;
 296     }
 297   }
 298
 299   tracks->Compress();
 300
 301   AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
 302 }
 303
 304 //-------------------------------------------------------------------
 305 void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
 306 {
 307   // remove particles from weak decays
 308   // <tracks> can be the following cases:
 309   // a. tracks: in this case the label is taken and then case b.
 310   // b. particles: it is checked if IsSecondaryFromWeakDecay is true
 311   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
 312
 313   TClonesArray* arrayMC = 0;
 314   AliMCEvent* mcEvent = 0;
 315   if (mcObj->InheritsFrom("AliMCEvent"))
 316     mcEvent = static_cast<AliMCEvent*>(mcObj);
 317   else if (mcObj->InheritsFrom("TClonesArray"))
 318     arrayMC = static_cast<TClonesArray*>(mcObj);
 319   else
 320   {
 321     arrayMC->Dump();
 322     AliFatal("Invalid object passed");
 323   }
 324
 325   Int_t before = tracks->GetEntriesFast();
 326
 327   for (Int_t i=0; i<before; ++i)
 328   {
 329     AliVParticle* part = (AliVParticle*) tracks->At(i);
 330
 331     if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
 332       part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
 333
 334     if (part->InheritsFrom("AliAODMCParticle"))
 335     {
 336       if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
 337         continue;
 338     }
 339     else if (part->InheritsFrom("AliMCParticle") && mcEvent)
 340     {
 341       if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(((AliMCParticle*) part)->Label())))
 342         continue;
 343     }
 344     else
 345     {
 346       part->Dump();
 347       AliFatal("Unknown particle");
 348     }
 349
 350 //     Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
 351     TObject* object = tracks->RemoveAt(i);
 352     if (tracks->IsOwner())
 353       delete object;
 354   }
 355
 356   tracks->Compress();
 357
 358   if (before > tracks->GetEntriesFast())
 359     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
 360 }
 361
 362 //-------------------------------------------------------------------
 363 TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
 364 {
 365   // 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
 366   // particleSpecies: -1 all particles are returned
 367   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 368
 369   Int_t nTracks = NParticles(obj);
 370   TObjArray* tracks = new TObjArray;
 371
 372   // for TPC only tracks
 373   Bool_t hasOwnership = kFALSE;
 374   if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
 375     hasOwnership = kTRUE;
 376
 377   if (!arrayMC)
 378     tracks->SetOwner(hasOwnership);
 379
 380   // Loop over tracks or jets
 381   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 382     AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
 383     if (!part) continue;
 384
 385     if (useEtaPtCuts)
 386       if (TMath::Abs(part->Eta()) > fTrackEtaCut || part->Pt() < fTrackPtMin)
 387       {
 388         if (hasOwnership)
 389           delete part;
 390         continue;
 391       }
 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 || 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                 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 567                 {
 568                   if (((AliAODMCParticle*)mother)->GetMother() < 0)
 569                   {
 570                     mother = 0;
 571                     break;
 572                   }
 573
 574                   mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 575                   if (!mother)
 576                     break;
 577                 }
 578
 579                 if (mother)
 580                 {
 581                   Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
 582                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 583                           return 0;
 584                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 585                           return 0;
 586                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 587                           return 0;
 588                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 589                           return 0;
 590                 }
 591                 else
 592                 {
 593                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 594                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 595                   part->Print();
 596
 597                   /*
 598                   // this code prints the details of the mother that is missing in the AOD
 599                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 600
 601                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 602
 603                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 604                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 605                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 606                   */
 607
 608                   if (particleSpecies != 3)
 609                     return 0;
 610                 }
 611         }
 612
 613   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 614         TObjArray *array = static_cast<TObjArray*>(obj);
 615         part = (AliVParticle*)array->At( ipart );
 616         if (!part)return 0;
 617   }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
 618         AliMCEvent* mcEvent =  static_cast<AliMCEvent*>(obj);
 619         part = mcEvent->GetTrack( ipart );
 620         if (!part) return 0;
 621         // eventually only primaries
 622         if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
 623         // eventually only hadrons
 624         //TO-DO
 625         /*if (fOnlyHadrons){
 626                 Int_t pdgCode =  part->GetPdgCode();
 627                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 628                                   TMath::Abs(pdgCode)==2212 || // Proton
 629                                   TMath::Abs(pdgCode)==321;    // Kaon
 630                 if (!isHadron) return 0;
 631                 }
 632         */
 633        if (particleSpecies != -1) {
 634                 // find the primary mother
 635                 AliMCParticle* mother = (AliMCParticle*) part;
 636 //              Printf("");
 637                 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
 638                 {
 639 //                Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
 640                   if (mother->GetMother() < 0)
 641                   {
 642                     mother = 0;
 643                     break;
 644                   }
 645
 646                   mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
 647                   if (!mother)
 648                     break;
 649                 }
 650
 651                 if (mother)
 652                 {
 653                   Int_t pdgCode = mother->PdgCode();
 654                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 655                           return 0;
 656                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 657                           return 0;
 658                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 659                           return 0;
 660                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 661                           return 0;
 662                 }
 663                 else
 664                 {
 665                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 666                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 667                   //part->Dump();
 668                   //part->Print();
 669
 670                   /*
 671                   // this code prints the details of the mother that is missing in the AOD
 672                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 673
 674                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 675
 676                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 677                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 678                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 679                   */
 680
 681                   if (particleSpecies != 3)
 682                     return 0;
 683                 }
 684         }
 685   }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
 686         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 687         part = aodEvent->GetTrack(ipart);
 688         // track selection cuts
 689         if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
 690         //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0;
 691         // only primary candidates
 692         //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
 693         // eventually only hadrons
 694         if (fOnlyHadrons){
 695                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 696                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 697                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 698                 if (!isHadron) return 0;
 699                 }
 700
 701   }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
 702         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 703         part = esdEvent->GetTrack(ipart);
 704         if (!part)return 0;
 705         // track selection cuts
 706
 707         if (!( ApplyCuts(part)) )
 708          return 0;
 709
 710         if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
 711         {
 712           // create TPC only tracks constrained to the SPD vertex
 713
 714           const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
 715
 716           AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
 717           if(!track) return 0;
 718
 719 //        Printf(">%f %f %f", track->Eta(), track->Phi(), track->Pt());
 720
 721           if(track->Pt()>0.){
 722             // only constrain tracks above threshold
 723             AliExternalTrackParam exParam;
 724             // take the B-feild from the ESD, no 3D fieldMap available at this point
 725             Bool_t relate = kFALSE;
 726             relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
 727             if(!relate)
 728             {
 729 //                 Printf("relating failed");
 730               delete track;
 731               return 0;
 732             }
 733             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 734           }
 735
 736 //        Printf(">%f %f %f\n", track->Eta(), track->Phi(), track->Pt());
 737
 738           part = track;
 739         }
 740         else if (fFilterBit == 2048)
 741         {
 742           // hybrid tracks
 743
 744           // clone
 745           AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
 746 //        Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
 747
 748           if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
 749           {
 750 //          Float_t ptBefore = esdTrack->Pt();
 751             // set constrained pT as default pT
 752             if (!esdTrack->GetConstrainedParam())
 753               return 0;
 754             esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
 755 //          Printf("%f %f", ptBefore, esdTrack->Pt());
 756           }
 757           part = esdTrack;
 758         }
 759
 760         // eventually only hadrons
 761         //TO-DO
 762         /*if (fOnlyHadrons){
 763                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 764                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 765                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 766                 if (!isHadron) return 0;
 767                 }
 768         */
 769   }else {
 770         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 771         return 0;
 772         }
 773
 774   // only charged
 775   if (!part->Charge())return 0;
 776
 777   return part;
 778 }
 779
 780
 781 //-------------------------------------------------------------------
 782 void  AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
 783 {
 784   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
 785
 786   static TObject *tmp;
 787   static int i;           // "static" to save stack space
 788   int j;
 789
 790   while (last - first > 1) {
 791     i = first;
 792     j = last;
 793     for (;;) {
 794       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
 795         ;
 796       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
 797         ;
 798       if (i >= j)
 799         break;
 800
 801       tmp  = a[i];
 802       a[i] = a[j];
 803       a[j] = tmp;
 804     }
 805     if (j == first) {
 806       ++first;
 807       continue;
 808     }
 809     tmp = a[first];
 810     a[first] = a[j];
 811     a[j] = tmp;
 812     if (j - first < last - (j + 1)) {
 813       QSortTracks(a, first, j);
 814       first = j + 1;   // QSortTracks(j + 1, last);
 815     } else {
 816       QSortTracks(a, j + 1, last);
 817       last = j;        // QSortTracks(first, j);
 818     }
 819   }
 820 }
 821
 822 //____________________________________________________________________
 823 TObjArray*  AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
 824 {
 825
 826   // Assign particles to towards, away or transverse regions.
 827   // Returns a lists of particles for each region.
 828
 829   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 830   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 831
 832   // Define output lists of particles
 833   TList *toward = new TList();
 834   TList *away = new TList();
 835   // Two transverse regions, for the moment those are not yet MIN and MAX!!!
 836   // MIN and MAX can be sorted in GetMinMaxRegion function
 837   TList *transverse1 = new TList();
 838   TList *transverse2 = new TList();
 839
 840   TObjArray *regionParticles = new TObjArray;
 841   regionParticles->SetOwner(kTRUE);
 842
 843   regionParticles->AddLast(toward);
 844   regionParticles->AddLast(away);
 845   regionParticles->AddLast(transverse1);
 846   regionParticles->AddLast(transverse2);
 847
 848   if (!leading)
 849     return regionParticles;
 850
 851   // Switch to vector for leading particle
 852   TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
 853
 854   Int_t nTracks = NParticles(obj);
 855   if( !nTracks ) return 0;
 856   // Loop over tracks
 857   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 858         AliVParticle* part = ParticleWithCuts(obj, ipart);
 859         if (!part)continue;
 860         //Switch to vectors for particles
 861         TVector3 partVect(part->Px(), part->Py(), part->Pz());
 862
 863         Int_t region = 0;
 864         if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
 865         // transverse regions
 866         if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
 867         if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1;   //right
 868
 869         if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2;    //forward
 870         if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2;  //backward
 871
 872         // skip leading particle
 873         if (leading == part)
 874           continue;
 875
 876         if (!region)continue;
 877         if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
 878                 Int_t label = ((AliAODTrack*)part)->GetLabel();
 879                 // re-define part as the matched MC particle
 880                 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
 881                 if (!part)continue;
 882                 // skip leading particle
 883                 if (leading == part)
 884                   continue;
 885                 }
 886         if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
 887                 Int_t label = ((AliESDtrack*)part)->GetLabel();
 888                 // look for the matched MC particle (but do not re-define part)
 889                 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
 890                 }
 891
 892         if ( region == 1 ) transverse1->Add(part);
 893         if ( region == -1 ) transverse2->Add(part);
 894         if ( region == 2 ) toward->Add(part);
 895         if ( region == -2 ) away->Add(part);
 896
 897         }//end loop on tracks
 898
 899   return regionParticles;
 900
 901 }
 902
 903
 904 //____________________________________________________________________
 905 Bool_t  AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
 906 {
 907   if (!obj) // MC
 908     return kFALSE;
 909
 910   // Use AliPhysicsSelection to select good events, works for ESD and AOD
 911   if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
 912     return kFALSE;
 913
 914   return kTRUE;
 915 }
 916
 917 //____________________________________________________________________
 918 Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
 919 {
 920
 921   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
 922
 923   if (obj->InheritsFrom("AliAODEvent")){
 924         Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
 925         if( nVertex > 0 ) {
 926                 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
 927                 Int_t nTracksPrim = vertex->GetNContributors();
 928                 Double_t zVertex = vertex->GetZ();
 929                 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 930                 // Reject TPC only vertex
 931                 TString name(vertex->GetName());
 932                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 933
 934                 // Select a quality vertex by number of tracks?
 935                 if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 936                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 937                         return kFALSE;
 938                         }
 939                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 940                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 941                 //  return kFALSE;
 942                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 943                 } else {
 944                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 945                         return kFALSE;
 946                         }
 947         }
 948
 949   if (obj->InheritsFrom("AliMCEvent"))
 950   {
 951     if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
 952     {
 953       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 954       return kFALSE;
 955     }
 956   }
 957
 958   if (obj->InheritsFrom("AliAODMCHeader"))
 959   {
 960     if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) > zed)
 961     {
 962       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 963       return kFALSE;
 964     }
 965   }
 966
 967   // ESD case for DCA studies
 968   if (obj->InheritsFrom("AliESDEvent")){
 969        AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
 970        if ( vertex){
 971                Int_t nTracksPrim = vertex->GetNContributors();
 972                Double_t zVertex = vertex->GetZ();
 973                if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 974                // Reject SPD or TPC only vertex
 975                TString name(vertex->GetName());
 976                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 977
 978                // Select a quality vertex by number of tracks?
 979                if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 980                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 981                        return kFALSE;
 982                        }
 983                // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 984                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 985                 //  return kFALSE;
 986                if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 987                } else {
 988                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 989                        return kFALSE;
 990                        }
 991        }
 992
 993   return kTRUE;
 994 }