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 "AliMCEventHandler.h"
  52 //#include "AliStack.h"
  53 #include "AliVParticle.h"
  54
  55 ////////////////////////////////////////////////
  56 //---------------------------------------------
  57 // Class for transverse regions analysis
  58 //---------------------------------------------
  59 ////////////////////////////////////////////////
  60
  61
  62 using namespace std;
  63
  64 ClassImp(AliAnalyseLeadingTrackUE)
  65
  66 //-------------------------------------------------------------------
  67 AliAnalyseLeadingTrackUE::AliAnalyseLeadingTrackUE() :
  68   TObject(),
  69   fDebug(0),
  70   fFilterBit(16),
  71   fOnlyHadrons(kFALSE),
  72   fTrackEtaCut(0.8)
  73 {
  74   // constructor
  75 }
  76
  77
  78 //-------------------------------------------------------------------
  79 AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
  80 {
  81   // assignment operator
  82   return *this;
  83 }
  84
  85
  86 //-------------------------------------------------------------------
  87 AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
  88 {
  89
  90   //clear memory
  91
  92 }
  93
  94
  95 //____________________________________________________________________
  96 Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track, Int_t filterbit)
  97 {
  98   // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
  99   // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
 100
 101   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
 102   switch (filterbit){
 103    /*case 16:
 104    // tighter cuts on primary particles for high pT tracks
 105    // needed as input for jetfinder
 106    esdTrackCuts->SetMinNClustersTPC(50);
 107    esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 108    esdTrackCuts->SetRequireTPCRefit(kTRUE);
 109    esdTrackCuts->SetMaxDCAToVertexXY(2.4);
 110    esdTrackCuts->SetMaxDCAToVertexZ(3.2);
 111    esdTrackCuts->SetDCAToVertex2D(kTRUE);
 112    esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 113    esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 114    esdTrackCuts->SetRequireITSRefit(kTRUE); // additional cut
 115    break;
 116    */
 117
 118    case 16:
 119    esdTrackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);
 120    break;
 121
 122    default:
 123    if (fDebug > 1) AliFatal("Set of cuts not defined");
 124    break;
 125    }
 126
 127   // select track according to set of cuts
 128   if (! esdTrackCuts->IsSelected(track) )return kFALSE;
 129
 130   return kTRUE;
 131 }
 132
 133
 134
 135
 136
 137 //____________________________________________________________________
 138 TObjArray*  AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
 139 {
 140
 141   // Returns an array of charged particles (or jets) ordered according to their pT.
 142
 143   Int_t nTracks = NParticles(obj);
 144
 145
 146   if( !nTracks ) return 0;
 147
 148   // Define array of AliVParticle objects
 149   TObjArray* tracks = new TObjArray(nTracks);
 150
 151   // Loop over tracks or jets
 152   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 153         AliVParticle* part = ParticleWithCuts( obj, ipart );
 154         if (!part) continue;
 155         // Accept leading-tracks in a limited pseudo-rapidity range
 156         if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
 157         tracks->AddLast( part );
 158         }
 159   // Order tracks by pT
 160   QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
 161
 162   nTracks = tracks->GetEntriesFast();
 163   if( !nTracks ) return 0;
 164
 165   return tracks;
 166   }
 167
 168
 169 //-------------------------------------------------------------------
 170 TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
 171 {
 172
 173  // Returns two lists of particles, one for MIN and one for MAX region
 174   Double_t sumpT1 = 0.;
 175   Double_t sumpT2 = 0.;
 176
 177   Int_t particles1 = transv1->GetEntries();
 178   Int_t particles2 = transv2->GetEntries();
 179
 180   // Loop on transverse region 1
 181   for (Int_t i=0; i<particles1; i++){
 182         AliVParticle *part = (AliVParticle*)transv1->At(i);
 183         sumpT1 +=  part->Pt();
 184         }
 185
 186   // Loop on transverse region 2
 187   for (Int_t i=0; i<particles2; i++){
 188         AliVParticle *part = (AliVParticle*)transv2->At(i);
 189         sumpT2 +=  part->Pt();
 190         }
 191
 192   TObjArray *regionParticles = new TObjArray;
 193   if ( sumpT2 >= sumpT1 ){
 194         regionParticles->AddLast(transv1); // MIN
 195         regionParticles->AddLast(transv2); // MAX
 196   }else {
 197         regionParticles->AddLast(transv2); // MIN
 198         regionParticles->AddLast(transv1); // MAX
 199         }
 200
 201  return regionParticles;
 202 }
 203
 204 //-------------------------------------------------------------------
 205 Int_t  AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
 206 {
 207
 208   //Returns the number of particles in AliAODMCParticle array  or AliAODTracks or AliESDTracks
 209
 210   Int_t nTracks;
 211
 212   if (obj->InheritsFrom("TClonesArray")){ // MC particles
 213         TClonesArray *arrayMC = dynamic_cast<TClonesArray*>(obj);
 214         nTracks = arrayMC->GetEntriesFast();
 215   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 216         TObjArray *array = dynamic_cast<TObjArray*>(obj);
 217         nTracks = array->GetEntriesFast();
 218   }else if (obj->InheritsFrom("AliAODEvent")){  // RECO AOD tracks
 219         AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(obj);
 220         nTracks = aodEvent->GetNTracks();
 221   }else if (obj->InheritsFrom("AliESDEvent")){  // RECO ESD tracks
 222         AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*>(obj);
 223         nTracks = esdEvent->GetNumberOfTracks();
 224   }else {
 225         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 226         return 0;
 227         }
 228
 229   return nTracks;
 230 }
 231
 232
 233 //-------------------------------------------------------------------
 234 AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries)
 235 {
 236   // Returns track or MC particle at position "ipart" if passes selection criteria
 237   AliVParticle *part=0;
 238
 239   if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
 240         TClonesArray *arrayMC = dynamic_cast<TClonesArray*>(obj);
 241         part = (AliVParticle*)arrayMC->At( ipart );
 242         if (!part)return 0;
 243         // eventually only primaries
 244         if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
 245         // eventually only hadrons
 246         if (fOnlyHadrons){
 247                 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
 248                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 249                                   TMath::Abs(pdgCode)==2212 || // Proton
 250                                   TMath::Abs(pdgCode)==321;    // Kaon
 251                 if (!isHadron) return 0;
 252                 }
 253
 254   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 255         TObjArray *array = dynamic_cast<TObjArray*>(obj);
 256         part = (AliVParticle*)array->At( ipart );
 257         if (!part)return 0;
 258   }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
 259         AliMCEvent* mcEvent =  dynamic_cast<AliMCEvent*>(obj);
 260         part = mcEvent->GetTrack( ipart );
 261         if (!part) return 0;
 262         // eventually only primaries
 263         if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
 264         // eventually only hadrons
 265         //TO-DO
 266         /*if (fOnlyHadrons){
 267                 Int_t pdgCode =  part->GetPdgCode();
 268                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 269                                   TMath::Abs(pdgCode)==2212 || // Proton
 270                                   TMath::Abs(pdgCode)==321;    // Kaon
 271                 if (!isHadron) return 0;
 272                 }
 273         */
 274
 275   }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
 276         AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(obj);
 277         part = aodEvent->GetTrack(ipart);
 278         // track selection cuts
 279         if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
 280         //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0;
 281         // only primary candidates
 282         //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
 283         // eventually only hadrons
 284         if (fOnlyHadrons){
 285                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 286                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 287                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 288                 if (!isHadron) return 0;
 289                 }
 290
 291   }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
 292         AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*>(obj);
 293         part = esdEvent->GetTrack(ipart);
 294         if (!part)return 0;
 295         // track selection cuts
 296         if (!( ApplyCuts(part, fFilterBit)) )return 0;
 297
 298         // only primary candidates (does not exist for ESD tracks??????)
 299         //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
 300
 301         // eventually only hadrons
 302         //TO-DO
 303         /*if (fOnlyHadrons){
 304                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 305                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 306                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 307                 if (!isHadron) return 0;
 308                 }
 309         */
 310   }else {
 311         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 312         return 0;
 313         }
 314
 315   // only charged
 316   if (!part->Charge())return 0;
 317
 318   return part;
 319 }
 320
 321
 322 //-------------------------------------------------------------------
 323 void  AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
 324 {
 325   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
 326
 327   static TObject *tmp;
 328   static int i;           // "static" to save stack space
 329   int j;
 330
 331   while (last - first > 1) {
 332     i = first;
 333     j = last;
 334     for (;;) {
 335       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
 336         ;
 337       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
 338         ;
 339       if (i >= j)
 340         break;
 341
 342       tmp  = a[i];
 343       a[i] = a[j];
 344       a[j] = tmp;
 345     }
 346     if (j == first) {
 347       ++first;
 348       continue;
 349     }
 350     tmp = a[first];
 351     a[first] = a[j];
 352     a[j] = tmp;
 353     if (j - first < last - (j + 1)) {
 354       QSortTracks(a, first, j);
 355       first = j + 1;   // QSortTracks(j + 1, last);
 356     } else {
 357       QSortTracks(a, j + 1, last);
 358       last = j;        // QSortTracks(first, j);
 359     }
 360   }
 361 }
 362
 363 //____________________________________________________________________
 364 TObjArray*  AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
 365 {
 366
 367   // Assign particles to towards, away or transverse regions.
 368   // Returns a lists of particles for each region.
 369
 370   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 371   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 372
 373   // Define output lists of particles
 374   TList *toward = new TList();
 375   TList *away = new TList();
 376   // Two transverse regions, for the moment those are not yet MIN and MAX!!!
 377   // MIN and MAX can be sorted in GetMinMaxRegion function
 378   TList *transverse1 = new TList();
 379   TList *transverse2 = new TList();
 380
 381   TObjArray *regionParticles = new TObjArray;
 382   regionParticles->AddLast(toward);
 383   regionParticles->AddLast(away);
 384   regionParticles->AddLast(transverse1);
 385   regionParticles->AddLast(transverse2);
 386
 387   if (!leading)
 388     return regionParticles;
 389
 390   // Switch to vector for leading particle
 391   TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
 392
 393   Int_t nTracks = NParticles(obj);
 394   if( !nTracks ) return 0;
 395   // Loop over tracks
 396   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 397         AliVParticle* part = ParticleWithCuts(obj, ipart);
 398         if (!part)continue;
 399         //Switch to vectors for particles
 400         TVector3 partVect(part->Px(), part->Py(), part->Pz());
 401
 402         Int_t region = 0;
 403         if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
 404         // transverse regions
 405         if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
 406         if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1;   //right
 407
 408         if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2;    //forward
 409         if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2;  //backward
 410
 411         // skip leading particle
 412         if (leading == part)
 413           continue;
 414
 415         if (!region)continue;
 416         if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
 417                 Int_t label = ((AliAODTrack*)part)->GetLabel();
 418                 // re-define part as the matched MC particle
 419                 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
 420                 if (!part)continue;
 421                 // skip leading particle
 422                 if (leading == part)
 423                   continue;
 424                 }
 425         if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
 426                 Int_t label = ((AliESDtrack*)part)->GetLabel();
 427                 // look for the matched MC particle (but do not re-define part)
 428                 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
 429                 }
 430
 431         if ( region == 1 ) transverse1->Add(part);
 432         if ( region == -1 ) transverse2->Add(part);
 433         if ( region == 2 ) toward->Add(part);
 434         if ( region == -2 ) away->Add(part);
 435
 436         }//end loop on tracks
 437
 438   return regionParticles;
 439
 440 }
 441
 442
 443 //____________________________________________________________________
 444 Bool_t  AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
 445 {
 446
 447   //Use AliPhysicsSelection to select good events
 448   if( !obj->InheritsFrom("AliAODInputHandler") ) { // not needed for AOD input
 449         if (!(((AliInputEventHandler*)obj)->IsEventSelected()&AliVEvent::kMB))return kFALSE;
 450         }
 451
 452   return kTRUE;
 453
 454 }
 455
 456 //____________________________________________________________________
 457 Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
 458 {
 459
 460   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
 461
 462   if (obj->InheritsFrom("AliAODEvent")){
 463         Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
 464         if( nVertex > 0 ) {
 465                 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
 466                 Int_t nTracksPrim = vertex->GetNContributors();
 467                 Double_t zVertex = vertex->GetZ();
 468                 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 469                 // Reject TPC only vertex
 470                 TString name(vertex->GetName());
 471                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 472
 473                 // Select a quality vertex by number of tracks?
 474                 if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 475                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 476                         return kFALSE;
 477                         }
 478                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 479                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 480                 //  return kFALSE;
 481                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 482                 } else {
 483                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 484                         return kFALSE;
 485                         }
 486         }
 487
 488   if (obj->InheritsFrom("AliMCEvent"))
 489   {
 490     if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
 491     {
 492       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 493       return kFALSE;
 494     }
 495   }
 496
 497   // ESD case for DCA studies
 498   if (obj->InheritsFrom("AliESDEvent")){
 499        AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
 500        if ( vertex){
 501                Int_t nTracksPrim = vertex->GetNContributors();
 502                Double_t zVertex = vertex->GetZ();
 503                if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 504                // Reject SPD or TPC only vertex
 505                TString name(vertex->GetName());
 506                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 507
 508                // Select a quality vertex by number of tracks?
 509                if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 510                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 511                        return kFALSE;
 512                        }
 513                // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 514                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 515                 //  return kFALSE;
 516                if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 517                } else {
 518                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 519                        return kFALSE;
 520                        }
 521        }
 522
 523   return kTRUE;
 524 }