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