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