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