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