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::GetMinMaxRegion(TList *transv1, TList *transv2)
 332 {
 333
 334  // Returns two lists of particles, one for MIN and one for MAX region
 335   Double_t sumpT1 = 0.;
 336   Double_t sumpT2 = 0.;
 337
 338   Int_t particles1 = transv1->GetEntries();
 339   Int_t particles2 = transv2->GetEntries();
 340
 341   // Loop on transverse region 1
 342   for (Int_t i=0; i<particles1; i++){
 343         AliVParticle *part = (AliVParticle*)transv1->At(i);
 344         sumpT1 +=  part->Pt();
 345         }
 346
 347   // Loop on transverse region 2
 348   for (Int_t i=0; i<particles2; i++){
 349         AliVParticle *part = (AliVParticle*)transv2->At(i);
 350         sumpT2 +=  part->Pt();
 351         }
 352
 353   TObjArray *regionParticles = new TObjArray;
 354   if ( sumpT2 >= sumpT1 ){
 355         regionParticles->AddLast(transv1); // MIN
 356         regionParticles->AddLast(transv2); // MAX
 357   }else {
 358         regionParticles->AddLast(transv2); // MIN
 359         regionParticles->AddLast(transv1); // MAX
 360         }
 361
 362  return regionParticles;
 363 }
 364
 365 //-------------------------------------------------------------------
 366 Int_t  AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
 367 {
 368
 369   //Returns the number of particles in AliAODMCParticle array  or AliAODTracks or AliESDTracks
 370
 371   Int_t nTracks;
 372
 373   if (obj->InheritsFrom("TClonesArray")){ // MC particles
 374         TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
 375         nTracks = arrayMC->GetEntriesFast();
 376   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 377         TObjArray *array = static_cast<TObjArray*>(obj);
 378         nTracks = array->GetEntriesFast();
 379   }else if (obj->InheritsFrom("AliAODEvent")){  // RECO AOD tracks
 380         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 381         nTracks = aodEvent->GetNTracks();
 382   }else if (obj->InheritsFrom("AliESDEvent")){  // RECO ESD tracks
 383         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 384         nTracks = esdEvent->GetNumberOfTracks();
 385   }else if (obj->InheritsFrom("AliMCEvent")){  // RECO ESD tracks
 386         AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
 387         nTracks = mcEvent->GetNumberOfTracks();
 388   }else {
 389         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 390         return 0;
 391         }
 392
 393   return nTracks;
 394 }
 395
 396
 397 //-------------------------------------------------------------------
 398 AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
 399 {
 400   // Returns track or MC particle at position "ipart" if passes selection criteria
 401   // particleSpecies: -1 all particles are returned
 402   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 403   AliVParticle *part=0;
 404
 405   if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
 406         TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
 407         part = (AliVParticle*)arrayMC->At( ipart );
 408         if (!part)return 0;
 409         // eventually only primaries
 410         if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
 411         // eventually only hadrons
 412         if (fOnlyHadrons){
 413                 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
 414                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 415                                   TMath::Abs(pdgCode)==2212 || // Proton
 416                                   TMath::Abs(pdgCode)==321;    // Kaon
 417                 if (!isHadron) return 0;
 418                 }
 419         if (particleSpecies != -1) {
 420                 // find the primary mother
 421                 AliVParticle* mother = part;
 422                 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 423                 {
 424                   if (((AliAODMCParticle*)mother)->GetMother() < 0)
 425                   {
 426                     mother = 0;
 427                     break;
 428                   }
 429
 430                   mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 431                   if (!mother)
 432                     break;
 433                 }
 434
 435                 if (mother)
 436                 {
 437                   Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
 438                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 439                           return 0;
 440                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 441                           return 0;
 442                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 443                           return 0;
 444                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 445                           return 0;
 446                 }
 447                 else
 448                 {
 449                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 450                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 451                   part->Print();
 452
 453                   /*
 454                   // this code prints the details of the mother that is missing in the AOD
 455                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 456
 457                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 458
 459                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 460                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 461                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 462                   */
 463
 464                   if (particleSpecies != 3)
 465                     return 0;
 466                 }
 467         }
 468
 469   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 470         TObjArray *array = static_cast<TObjArray*>(obj);
 471         part = (AliVParticle*)array->At( ipart );
 472         if (!part)return 0;
 473   }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
 474         AliMCEvent* mcEvent =  static_cast<AliMCEvent*>(obj);
 475         part = mcEvent->GetTrack( ipart );
 476         if (!part) return 0;
 477         // eventually only primaries
 478         if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
 479         // eventually only hadrons
 480         //TO-DO
 481         /*if (fOnlyHadrons){
 482                 Int_t pdgCode =  part->GetPdgCode();
 483                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 484                                   TMath::Abs(pdgCode)==2212 || // Proton
 485                                   TMath::Abs(pdgCode)==321;    // Kaon
 486                 if (!isHadron) return 0;
 487                 }
 488         */
 489        if (particleSpecies != -1) {
 490                 // find the primary mother
 491                 AliMCParticle* mother = (AliMCParticle*) part;
 492 //              Printf("");
 493                 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
 494                 {
 495 //                Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
 496                   if (mother->GetMother() < 0)
 497                   {
 498                     mother = 0;
 499                     break;
 500                   }
 501
 502                   mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
 503                   if (!mother)
 504                     break;
 505                 }
 506
 507                 if (mother)
 508                 {
 509                   Int_t pdgCode = mother->PdgCode();
 510                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 511                           return 0;
 512                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 513                           return 0;
 514                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 515                           return 0;
 516                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 517                           return 0;
 518                 }
 519                 else
 520                 {
 521                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 522                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 523                   //part->Dump();
 524                   //part->Print();
 525
 526                   /*
 527                   // this code prints the details of the mother that is missing in the AOD
 528                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 529
 530                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 531
 532                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 533                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 534                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 535                   */
 536
 537                   if (particleSpecies != 3)
 538                     return 0;
 539                 }
 540         }
 541   }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
 542         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 543         part = aodEvent->GetTrack(ipart);
 544         // track selection cuts
 545         if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
 546         //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0;
 547         // only primary candidates
 548         //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
 549         // eventually only hadrons
 550         if (fOnlyHadrons){
 551                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 552                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 553                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 554                 if (!isHadron) return 0;
 555                 }
 556
 557   }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
 558         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 559         part = esdEvent->GetTrack(ipart);
 560         if (!part)return 0;
 561         // track selection cuts
 562
 563         if (!( ApplyCuts(part)) )
 564          return 0;
 565
 566         if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
 567         {
 568           // create TPC only tracks constrained to the SPD vertex
 569
 570           const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
 571
 572           AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
 573           if(!track) return 0;
 574
 575           if(track->Pt()>0.){
 576             // only constrain tracks above threshold
 577             AliExternalTrackParam exParam;
 578             // take the B-feild from the ESD, no 3D fieldMap available at this point
 579             Bool_t relate = kFALSE;
 580             relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
 581             if(!relate)
 582             {
 583 //                 Printf("relating failed");
 584               delete track;
 585               return 0;
 586             }
 587             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 588           }
 589
 590           part = track;
 591         }
 592
 593         // eventually only hadrons
 594         //TO-DO
 595         /*if (fOnlyHadrons){
 596                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 597                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 598                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 599                 if (!isHadron) return 0;
 600                 }
 601         */
 602   }else {
 603         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 604         return 0;
 605         }
 606
 607   // only charged
 608   if (!part->Charge())return 0;
 609
 610   return part;
 611 }
 612
 613
 614 //-------------------------------------------------------------------
 615 void  AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
 616 {
 617   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
 618
 619   static TObject *tmp;
 620   static int i;           // "static" to save stack space
 621   int j;
 622
 623   while (last - first > 1) {
 624     i = first;
 625     j = last;
 626     for (;;) {
 627       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
 628         ;
 629       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
 630         ;
 631       if (i >= j)
 632         break;
 633
 634       tmp  = a[i];
 635       a[i] = a[j];
 636       a[j] = tmp;
 637     }
 638     if (j == first) {
 639       ++first;
 640       continue;
 641     }
 642     tmp = a[first];
 643     a[first] = a[j];
 644     a[j] = tmp;
 645     if (j - first < last - (j + 1)) {
 646       QSortTracks(a, first, j);
 647       first = j + 1;   // QSortTracks(j + 1, last);
 648     } else {
 649       QSortTracks(a, j + 1, last);
 650       last = j;        // QSortTracks(first, j);
 651     }
 652   }
 653 }
 654
 655 //____________________________________________________________________
 656 TObjArray*  AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
 657 {
 658
 659   // Assign particles to towards, away or transverse regions.
 660   // Returns a lists of particles for each region.
 661
 662   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 663   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 664
 665   // Define output lists of particles
 666   TList *toward = new TList();
 667   TList *away = new TList();
 668   // Two transverse regions, for the moment those are not yet MIN and MAX!!!
 669   // MIN and MAX can be sorted in GetMinMaxRegion function
 670   TList *transverse1 = new TList();
 671   TList *transverse2 = new TList();
 672
 673   TObjArray *regionParticles = new TObjArray;
 674   regionParticles->SetOwner(kTRUE);
 675
 676   regionParticles->AddLast(toward);
 677   regionParticles->AddLast(away);
 678   regionParticles->AddLast(transverse1);
 679   regionParticles->AddLast(transverse2);
 680
 681   if (!leading)
 682     return regionParticles;
 683
 684   // Switch to vector for leading particle
 685   TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
 686
 687   Int_t nTracks = NParticles(obj);
 688   if( !nTracks ) return 0;
 689   // Loop over tracks
 690   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 691         AliVParticle* part = ParticleWithCuts(obj, ipart);
 692         if (!part)continue;
 693         //Switch to vectors for particles
 694         TVector3 partVect(part->Px(), part->Py(), part->Pz());
 695
 696         Int_t region = 0;
 697         if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
 698         // transverse regions
 699         if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
 700         if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1;   //right
 701
 702         if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2;    //forward
 703         if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2;  //backward
 704
 705         // skip leading particle
 706         if (leading == part)
 707           continue;
 708
 709         if (!region)continue;
 710         if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
 711                 Int_t label = ((AliAODTrack*)part)->GetLabel();
 712                 // re-define part as the matched MC particle
 713                 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
 714                 if (!part)continue;
 715                 // skip leading particle
 716                 if (leading == part)
 717                   continue;
 718                 }
 719         if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
 720                 Int_t label = ((AliESDtrack*)part)->GetLabel();
 721                 // look for the matched MC particle (but do not re-define part)
 722                 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
 723                 }
 724
 725         if ( region == 1 ) transverse1->Add(part);
 726         if ( region == -1 ) transverse2->Add(part);
 727         if ( region == 2 ) toward->Add(part);
 728         if ( region == -2 ) away->Add(part);
 729
 730         }//end loop on tracks
 731
 732   return regionParticles;
 733
 734 }
 735
 736
 737 //____________________________________________________________________
 738 Bool_t  AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
 739 {
 740   if (!obj) // MC
 741     return kFALSE;
 742
 743   // Use AliPhysicsSelection to select good events, works for ESD and AOD
 744   if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
 745     return kFALSE;
 746
 747   return kTRUE;
 748 }
 749
 750 //____________________________________________________________________
 751 Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
 752 {
 753
 754   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
 755
 756   if (obj->InheritsFrom("AliAODEvent")){
 757         Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
 758         if( nVertex > 0 ) {
 759                 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
 760                 Int_t nTracksPrim = vertex->GetNContributors();
 761                 Double_t zVertex = vertex->GetZ();
 762                 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 763                 // Reject TPC only vertex
 764                 TString name(vertex->GetName());
 765                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 766
 767                 // Select a quality vertex by number of tracks?
 768                 if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 769                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 770                         return kFALSE;
 771                         }
 772                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 773                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 774                 //  return kFALSE;
 775                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 776                 } else {
 777                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 778                         return kFALSE;
 779                         }
 780         }
 781
 782   if (obj->InheritsFrom("AliMCEvent"))
 783   {
 784     if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
 785     {
 786       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 787       return kFALSE;
 788     }
 789   }
 790
 791   if (obj->InheritsFrom("AliAODMCHeader"))
 792   {
 793     if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) > zed)
 794     {
 795       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 796       return kFALSE;
 797     }
 798   }
 799
 800   // ESD case for DCA studies
 801   if (obj->InheritsFrom("AliESDEvent")){
 802        AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
 803        if ( vertex){
 804                Int_t nTracksPrim = vertex->GetNContributors();
 805                Double_t zVertex = vertex->GetZ();
 806                if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 807                // Reject SPD or TPC only vertex
 808                TString name(vertex->GetName());
 809                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 810
 811                // Select a quality vertex by number of tracks?
 812                if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 813                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 814                        return kFALSE;
 815                        }
 816                // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 817                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 818                 //  return kFALSE;
 819                if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 820                } else {
 821                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 822                        return kFALSE;
 823                        }
 824        }
 825
 826   return kTRUE;
 827 }