]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/Correlations/Base/AliAnalyseLeadingTrackUE.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / Correlations / Base / AliAnalyseLeadingTrackUE.cxx
index 5eee46b3f08116a4f06af381c081561d461cfa13..d156f4e14aa633f8df2a58346bf9000a7bae15cc 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-//#include <TBranch.h>
-//#include <TCanvas.h>
-//#include <TChain.h>
-//#include <TFile.h>
-//#include <TH1F.h>
-//#include <TH1I.h>
-//#include <TH2F.h>
+
 #include <TList.h>
-//#include <TLorentzVector.h>
 #include <TMath.h>
 #include <TObjArray.h>
 #include <TObject.h>
-//#include <TProfile.h>
-//#include <TRandom.h>
-//#include <TSystem.h>
-//#include <TTree.h>
 #include <TVector3.h>
 
 #include "AliAnalyseLeadingTrackUE.h"
-//#include "AliAnalysisTask.h"
 
-//#include "AliAnalysisHelperJetTasks.h"
-//#include "AliAnalysisManager.h"
 #include "AliAODEvent.h"
-//#include "AliAODHandler.h"
-//#include "AliAODJet.h"
 #include "AliAODMCParticle.h"
 #include "AliAODTrack.h"
 #include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliESDtrackCuts.h"
-//#include "AliGenPythiaEventHeader.h"
 #include "AliInputEventHandler.h"
-//#include "AliKFVertex.h"
-//#include "AliLog.h"
 #include "AliMCEvent.h"
 #include "AliVParticle.h"
 #include "AliAODMCHeader.h"
+#include "TFormula.h"
 
 #include "AliAnalysisManager.h"
 #include "AliMCEventHandler.h"
 #include "AliStack.h"
+#include "AliPIDResponse.h"
+#include "AliHelperPID.h"
 
 
 ////////////////////////////////////////////////
@@ -72,13 +56,24 @@ AliAnalyseLeadingTrackUE::AliAnalyseLeadingTrackUE() :
   TObject(),
   fDebug(0),
   fFilterBit(16),
+  fTrackStatus(0),
   fOnlyHadrons(kFALSE),
+  fCheckMotherPDG(kTRUE),
   fTrackEtaCut(0.8),
+  fTrackEtaCutMin(-1.),
+  fTrackPhiCutEvPlMin(0.),
+  fTrackPhiCutEvPlMax(0.),
   fTrackPtMin(0),
   fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
+  fDCAXYCut(0),
+  fSharedClusterCut(-1),
+  fCrossedRowsCut(-1),
+  fFoundFractionCut(-1),
   fEsdTrackCuts(0x0), 
   fEsdTrackCutsExtra1(0x0), 
-  fEsdTrackCutsExtra2(0x0) 
+  fEsdTrackCutsExtra2(0x0), 
+  fHelperPID(0x0),
+  fEventCounter(0)
 {
   // constructor
 }
@@ -162,6 +157,7 @@ void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
     fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
     fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
+    fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
 
     // Add SPD requirement 
     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
@@ -172,9 +168,14 @@ void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
   }
+  else if (filterbit == 4096) // require TOF matching
+  {
+    fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
+    // FIXME: TOF REQUIREMENTS ARE IN GetParticleSpecies FOR THE MOMENT
+  }
   else
   {
-    fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
+    fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 
     // Add SPD requirement 
@@ -206,7 +207,7 @@ TObjArray*  AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
        AliVParticle* part = ParticleWithCuts( obj, ipart );
         if (!part) continue;
        // Accept leading-tracks in a limited pseudo-rapidity range     
-       if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
+       if( TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin ) continue;
        tracks->AddLast( part );
        }
   // Order tracks by pT        
@@ -235,7 +236,7 @@ void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject*
     arrayMC = static_cast<TClonesArray*>(mcObj);
   else
   {
-    arrayMC->Dump();
+    mcObj->Dump();
     AliFatal("Invalid object passed");
   }
   
@@ -286,10 +287,11 @@ void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject*
       Printf("WARNING: No mother found for particle %d:", part->GetLabel());
       continue;
     }
-   
-    if (mother->GetLabel() > maxLabel)
+
+//     Printf("%d %d %d", i, part->GetLabel(), mother->GetLabel());
+    if (mother->GetLabel() >= maxLabel)
     {
-//       Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
+//       Printf("Removing %d with label %d", i, part->GetLabel()); ((AliMCParticle*)part)->Particle()->Print(); ((AliMCParticle*)mother)->Particle()->Print();
       TObject* object = tracks->RemoveAt(i);
       if (tracks->IsOwner())
        delete object;
@@ -318,7 +320,7 @@ void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcOb
     arrayMC = static_cast<TClonesArray*>(mcObj);
   else
   {
-    arrayMC->Dump();
+    mcObj->Dump();
     AliFatal("Invalid object passed");
   }
   
@@ -338,7 +340,7 @@ void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcOb
     }
     else if (part->InheritsFrom("AliMCParticle") && mcEvent)
     {
-      if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(((AliMCParticle*) part)->Label())))
+      if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
        continue;
     }
     else
@@ -360,12 +362,14 @@ void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcOb
 }
 
 //-------------------------------------------------------------------
-TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
+TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts, Bool_t speciesOnTracks, Double_t evtPlane)
 {
   // 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 
   // particleSpecies: -1 all particles are returned
   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
-
+  // speciesOnTracks if kFALSE, particleSpecies is only applied on the matched MC particle (not on the track itself)
+  // Passing down the Double_t* evtPlane (range [-pi/2,pi/2]) will apply a phi cut with respect to the eventplane between fTrackPhiCutEvPlMin and fTrackPhiCutEvPlMax. For values outside [-pi/2,pi/2], this will be ignored.
+  
   Int_t nTracks = NParticles(obj);
   TObjArray* tracks = new TObjArray;
   
@@ -379,16 +383,35 @@ TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject*
  
   // Loop over tracks or jets
   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
-    AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
+    AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, (speciesOnTracks) ? particleSpecies : -1);
     if (!part) continue;
     
+    if (TMath::Abs(evtPlane)<=TMath::Pi()/2) { //evtPlane range: (-pi/2,pi/2)
+      Double_t phiPart = part->Phi(); //range: [0,2*pi)
+      if(phiPart>TMath::Pi()) phiPart-=2*TMath::Pi();
+
+      Double_t dPhi = 0; //range: [0,pi/2], i.e. the difference over the shortest angle.
+      Double_t diff = TMath::Abs(phiPart-evtPlane);
+      if(diff<=TMath::Pi()/2) dPhi = diff;
+      else if(diff<=TMath::Pi()) dPhi = TMath::Pi()-diff;
+      else dPhi = diff-TMath::Pi();
+      
+      if(dPhi<fTrackPhiCutEvPlMin || dPhi>fTrackPhiCutEvPlMax) {
+        if (hasOwnership)
+          delete part;
+        continue;
+      }
+    }
+    
     if (useEtaPtCuts)
-      if (TMath::Abs(part->Eta()) > fTrackEtaCut || part->Pt() < fTrackPtMin)
+      if (TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin || part->Pt() < fTrackPtMin)
       {
        if (hasOwnership)
          delete part;
         continue;
       }
+      
+//     Printf("%p %p %d Accepted %d %f %f %f", obj, arrayMC, particleSpecies, ipart, part->Eta(), part->Phi(), part->Pt());
     
     if (arrayMC) {
       Int_t label = part->GetLabel();
@@ -430,7 +453,7 @@ TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arr
     if (!partReconstructed) continue;
 
     if (useEtaPtCuts)
-      if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || partReconstructed->Pt() < fTrackPtMin)
+      if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || TMath::Abs(partReconstructed->Eta()) < fTrackEtaCutMin || partReconstructed->Pt() < fTrackPtMin)
       {
         if (hasOwnership)
           delete partReconstructed;
@@ -523,7 +546,7 @@ Int_t  AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
         nTracks = array->GetEntriesFast();
   }else if (obj->InheritsFrom("AliAODEvent")){  // RECO AOD tracks
        AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
-        nTracks = aodEvent->GetNTracks();
+        nTracks = aodEvent->GetNumberOfTracks();
   }else if (obj->InheritsFrom("AliESDEvent")){  // RECO ESD tracks
        AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
         nTracks = esdEvent->GetNumberOfTracks();
@@ -563,8 +586,9 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
         if (particleSpecies != -1) {
                 // find the primary mother
                 AliVParticle* mother = part;
-                while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
-                {
+                if(fCheckMotherPDG) {
+                 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
+                 {
                   if (((AliAODMCParticle*)mother)->GetMother() < 0)
                   {
                     mother = 0;
@@ -574,8 +598,8 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
                   mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
                   if (!mother)
                     break;
+                 }
                 }
-                
                 if (mother)
                 {
                   Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
@@ -634,8 +658,9 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
                 // find the primary mother
                 AliMCParticle* mother = (AliMCParticle*) part;
 //             Printf("");
-                while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
-                {
+                if(fCheckMotherPDG) {
+                 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
+                  {
 //               Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
                   if (mother->GetMother() < 0)
                   {
@@ -646,8 +671,8 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
                   mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
                   if (!mother)
                     break;
+                  }
                 }
-                
                 if (mother)
                 {
                   Int_t pdgCode = mother->PdgCode();
@@ -685,11 +710,54 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
   }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
        AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
         part = aodEvent->GetTrack(ipart);
+       
        // track selection cuts
-       if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0; 
-       //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0; 
-       // only primary candidates
-       //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
+       if (fFilterBit != 0 && !(((AliAODTrack*)part)->TestFilterBit(fFilterBit))) return 0; 
+       if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
+       
+       // DCA XY
+       if (fDCAXYCut)
+       {
+         const AliVVertex* vertex = aodEvent->GetPrimaryVertex();
+         if (!vertex)
+           return 0;
+         
+         Double_t pos[2];
+         Double_t covar[3];
+         AliAODTrack* clone = (AliAODTrack*) part->Clone();
+         Bool_t success = clone->PropagateToDCA(vertex, ((AliVAODHeader*)aodEvent->GetHeader())->GetMagneticField(), 3, pos, covar);
+         delete clone;
+         if (!success)
+           return 0;
+
+//       Printf("%f", ((AliAODTrack*)part)->DCA());
+//       Printf("%f", pos[0]);
+         if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(part->Pt()))
+           return 0;
+       }
+       
+       if (fSharedClusterCut >= 0)
+       {
+         Double_t frac = Double_t(((AliAODTrack*)part)->GetTPCnclsS()) / Double_t(((AliAODTrack*)part)->GetTPCncls());
+         if (frac > fSharedClusterCut)
+           return 0;
+       }
+       
+       if (fCrossedRowsCut >= 0)
+       {
+         if (((AliAODTrack*) part)->GetTPCNCrossedRows() < fCrossedRowsCut)
+           return 0;
+       }
+       
+       if (fFoundFractionCut >= 0)
+       {
+         UInt_t findableClusters = ((AliAODTrack*) part)->GetTPCNclsF();
+         if (findableClusters == 0)
+           return 0;
+         if (((Double_t) ((AliAODTrack*) part)->GetTPCNCrossedRows() / findableClusters) < fFoundFractionCut)
+           return 0;
+       }
+
        // eventually only hadrons
        if (fOnlyHadrons){
                Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
@@ -697,16 +765,20 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
                                  ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
                if (!isHadron) return 0;                                  
                }
+               
+       if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0;
   
   }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
        AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
         part = esdEvent->GetTrack(ipart);
        if (!part)return 0;
+
        // track selection cuts
-       
        if (!( ApplyCuts(part)) )
-        return 0; 
+         return 0; 
        
+       if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
+
        if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
        {
          // create TPC only tracks constrained to the SPD vertex
@@ -766,6 +838,8 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
                if (!isHadron) return 0;                                  
                }
        */      
+       if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0; // If it is -1 you take all the particles
+
   }else {
        if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
        return 0;
@@ -774,6 +848,7 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
   // only charged
   if (!part->Charge())return 0;
   
+  part->SetUniqueID(fEventCounter * 100000 + ipart);
   return part;
 }
 
@@ -861,7 +936,7 @@ TObjArray*  AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, T
        TVector3 partVect(part->Px(), part->Py(), part->Pz());
  
        Int_t region = 0;
-       if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
+       if( TMath::Abs(partVect.Eta()) > fTrackEtaCut || TMath::Abs(partVect.Eta()) < fTrackEtaCutMin) continue;
        // transverse regions
        if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
        if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1;   //right
@@ -932,7 +1007,7 @@ Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntra
                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 
                // Select a quality vertex by number of tracks?
-               if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
+               if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
                        return kFALSE;
                        }
@@ -948,7 +1023,7 @@ Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntra
 
   if (obj->InheritsFrom("AliMCEvent"))
   { 
-    if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
+    if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) >= zed)
     {
       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
       return kFALSE;
@@ -957,7 +1032,7 @@ Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntra
 
   if (obj->InheritsFrom("AliAODMCHeader"))
   { 
-    if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) > zed)
+    if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= zed)
     {
       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
       return kFALSE;
@@ -976,7 +1051,7 @@ Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntra
                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 
                // Select a quality vertex by number of tracks?
-               if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
+               if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
                        return kFALSE;
                        }
@@ -992,3 +1067,15 @@ Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntra
        
   return kTRUE;
 }
+
+//____________________________________________________________________
+
+Bool_t AliAnalyseLeadingTrackUE::CheckTrack(AliVParticle * part)
+{
+  // check if the track status flags are set
+  
+  UInt_t status=((AliVTrack*)part)->GetStatus();
+  if ((status & fTrackStatus) == fTrackStatus)
+    return kTRUE;
+  return kFALSE;
+}