]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FLOW/AliFlowTasks/AliFlowTrackCuts.cxx
new parameter mixing option: kTrackWithPtFromFirstMother
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliFlowTrackCuts.cxx
index 4e948d7e26ed73a2165b1cba303ee89b35db7a8d..3e92f87a6f4b15ea1aac8f690a9ba5e01b83e11f 100644 (file)
 
 #include <limits.h>
 #include <float.h>
+#include "TParticle.h"
+#include "AliStack.h"
 #include "AliMCEvent.h"
+#include "AliESDEvent.h"
 #include "AliVParticle.h"
 #include "AliMCParticle.h"
 #include "AliESDtrack.h"
@@ -53,81 +56,107 @@ AliFlowTrackCuts::AliFlowTrackCuts():
   AliFlowTrackSimpleCuts(),
   fAliESDtrackCuts(new AliESDtrackCuts()),
   fQA(kFALSE),
+  fCutMC(kFALSE),
   fCutMCprocessType(kFALSE),
   fMCprocessType(kPNoProcess),
   fCutMCPID(kFALSE),
   fMCPID(0),
+  fIgnoreSignInPID(kFALSE),
   fCutMCisPrimary(kFALSE),
   fMCisPrimary(kFALSE),
   fRequireCharge(kFALSE),
   fFakesAreOK(kTRUE),
+  fCutSPDtrackletDeltaPhi(kFALSE),
+  fSPDtrackletDeltaPhiMax(FLT_MAX),
+  fSPDtrackletDeltaPhiMin(-FLT_MAX),
+  fIgnoreTPCzRange(kFALSE),
+  fIgnoreTPCzRangeMax(FLT_MAX),
+  fIgnoreTPCzRangeMin(-FLT_MAX),
   fParamType(kGlobal),
   fParamMix(kPure),
-  fCleanupTrack(kFALSE),
   fTrack(NULL),
   fTrackPhi(0.),
   fTrackEta(0.),
   fTrackWeight(0.),
   fTrackLabel(INT_MIN),
   fMCevent(NULL),
-  fMCparticle(NULL)
+  fMCparticle(NULL),
+  fEvent(NULL),
+  fTPCtrack()
 {
   //constructor 
 }
 
 //-----------------------------------------------------------------------
-AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& someCuts):
-  AliFlowTrackSimpleCuts(someCuts),
-  fAliESDtrackCuts(new AliESDtrackCuts(*(someCuts.fAliESDtrackCuts))),
-  fQA(someCuts.fQA),
-  fCutMCprocessType(someCuts.fCutMCprocessType),
-  fMCprocessType(someCuts.fMCprocessType),
-  fCutMCPID(someCuts.fCutMCPID),
-  fMCPID(someCuts.fMCPID),
-  fCutMCisPrimary(someCuts.fCutMCisPrimary),
-  fMCisPrimary(someCuts.fMCisPrimary),
-  fRequireCharge(someCuts.fRequireCharge),
-  fFakesAreOK(someCuts.fFakesAreOK),
-  fParamType(someCuts.fParamType),
-  fParamMix(someCuts.fParamMix),
-  fCleanupTrack(kFALSE),
+AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
+  AliFlowTrackSimpleCuts(that),
+  fAliESDtrackCuts(new AliESDtrackCuts(*(that.fAliESDtrackCuts))),
+  fQA(that.fQA),
+  fCutMC(that.fCutMC),
+  fCutMCprocessType(that.fCutMCprocessType),
+  fMCprocessType(that.fMCprocessType),
+  fCutMCPID(that.fCutMCPID),
+  fMCPID(that.fMCPID),
+  fIgnoreSignInPID(that.fIgnoreSignInPID),
+  fCutMCisPrimary(that.fCutMCisPrimary),
+  fMCisPrimary(that.fMCisPrimary),
+  fRequireCharge(that.fRequireCharge),
+  fFakesAreOK(that.fFakesAreOK),
+  fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
+  fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
+  fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
+  fIgnoreTPCzRange(that.fIgnoreTPCzRange),
+  fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
+  fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
+  fParamType(that.fParamType),
+  fParamMix(that.fParamMix),
   fTrack(NULL),
-  fTrackPhi(someCuts.fTrackPhi),
-  fTrackEta(someCuts.fTrackEta),
-  fTrackWeight(someCuts.fTrackWeight),
+  fTrackPhi(0.),
+  fTrackEta(0.),
+  fTrackWeight(0.),
   fTrackLabel(INT_MIN),
   fMCevent(NULL),
-  fMCparticle(NULL)
+  fMCparticle(NULL),
+  fEvent(NULL),
+  fTPCtrack(that.fTPCtrack)
 {
   //copy constructor
 }
 
 //-----------------------------------------------------------------------
-AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& someCuts)
+AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
 {
   //assignment
-  AliFlowTrackSimpleCuts::operator=(someCuts);
-  *fAliESDtrackCuts=*(someCuts.fAliESDtrackCuts);
-  fQA=someCuts.fQA;
-  fCutMCprocessType=someCuts.fCutMCprocessType;
-  fMCprocessType=someCuts.fMCprocessType;
-  fCutMCPID=someCuts.fCutMCPID;
-  fMCPID=someCuts.fMCPID;
-  fCutMCisPrimary=someCuts.fCutMCisPrimary;
-  fMCisPrimary=someCuts.fMCisPrimary;
-  fRequireCharge=someCuts.fRequireCharge;
-  fFakesAreOK=someCuts.fFakesAreOK;
-  fParamType=someCuts.fParamType;
-  fParamMix=someCuts.fParamMix;
-
-  fCleanupTrack=kFALSE;
+  AliFlowTrackSimpleCuts::operator=(that);
+  *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
+  fQA=that.fQA;
+  fCutMC=that.fCutMC;
+  fCutMCprocessType=that.fCutMCprocessType;
+  fMCprocessType=that.fMCprocessType;
+  fCutMCPID=that.fCutMCPID;
+  fMCPID=that.fMCPID;
+  fIgnoreSignInPID=that.fIgnoreSignInPID,
+  fCutMCisPrimary=that.fCutMCisPrimary;
+  fMCisPrimary=that.fMCisPrimary;
+  fRequireCharge=that.fRequireCharge;
+  fFakesAreOK=that.fFakesAreOK;
+  fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
+  fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
+  fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
+  fIgnoreTPCzRange=that.fIgnoreTPCzRange;
+  fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
+  fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
+  fParamType=that.fParamType;
+  fParamMix=that.fParamMix;
+
   fTrack=NULL;
-  fTrackPhi=someCuts.fTrackPhi;
-  fTrackPhi=someCuts.fTrackPhi;
-  fTrackWeight=someCuts.fTrackWeight;
+  fTrackPhi=0.;
+  fTrackPhi=0.;
+  fTrackWeight=0.;
   fTrackLabel=INT_MIN;
   fMCevent=NULL;
   fMCparticle=NULL;
+  fEvent=NULL;
 
   return *this;
 }
@@ -136,7 +165,6 @@ AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& someCuts)
 AliFlowTrackCuts::~AliFlowTrackCuts()
 {
   //dtor
-  if (fCleanupTrack) delete fTrack;
   delete fAliESDtrackCuts;
 }
 
@@ -153,13 +181,33 @@ Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
   return kFALSE;  //default when passed wrong type of object
 }
 
+//-----------------------------------------------------------------------
+Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
+{
+  //check cuts
+  AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
+  if (vparticle) 
+  {
+    return PassesMCcuts(fMCevent,vparticle->GetLabel());
+  }
+  AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
+  if (tracklets)
+  {
+    Int_t label0 = tracklets->GetLabel(id,0);
+    Int_t label1 = tracklets->GetLabel(id,1);
+    Int_t label = (label0==label1)?tracklets->GetLabel(id,1):-666;
+    return PassesMCcuts(fMCevent,label);
+  }
+  return kFALSE;  //default when passed wrong type of object
+}
+
 //-----------------------------------------------------------------------
 Bool_t AliFlowTrackCuts::PassesCuts(AliFlowTrackSimple* track)
 {
   //check cuts on a flowtracksimple
 
   //clean up from last iteration
-  if (fCleanupTrack) delete fTrack; fTrack = NULL;
+  fTrack = NULL;
   return AliFlowTrackSimpleCuts::PassesCuts(track);
 }
 
@@ -168,9 +216,10 @@ Bool_t AliFlowTrackCuts::PassesCuts(AliMultiplicity* tracklet, Int_t id)
 {
   //check cuts on a tracklets
 
-  //clean up from last iteration
-  if (fCleanupTrack) delete fTrack; fTrack = NULL;
+  //clean up from last iteration, and init label
+  fTrack = NULL;
   fMCparticle=NULL;
+  fTrackLabel=-1;
 
   fTrackPhi = tracklet->GetPhi(id);
   fTrackEta = tracklet->GetEta(id);
@@ -179,36 +228,59 @@ Bool_t AliFlowTrackCuts::PassesCuts(AliMultiplicity* tracklet, Int_t id)
   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
 
   //check MC info if available
-  fTrackLabel = tracklet->GetLabel(id,1); //TODO: this can be improved
-  if (!PassesMCcuts()) return kFALSE;
+  //if the 2 clusters have different label track cannot be good
+  //and should therefore not pass the mc cuts
+  Int_t label0 = tracklet->GetLabel(id,0);
+  Int_t label1 = tracklet->GetLabel(id,1);
+  //if possible get label and mcparticle
+  fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
+  if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
+  //check MC cuts
+  if (fCutMC && !PassesMCcuts()) return kFALSE;
   return kTRUE;
 }
 
 //-----------------------------------------------------------------------
-Bool_t AliFlowTrackCuts::PassesMCcuts()
+Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
 {
   //check the MC info
-  if (!fMCevent) {AliError("no MC info"); return kFALSE;}
-  fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
-  if (!fMCparticle) {AliError("no MC info"); return kFALSE;}
+  if (!mcEvent) return kFALSE;
+  if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
+  AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
+  if (!mcparticle) {AliError("no MC track"); return kFALSE;}
 
   if (fCutMCisPrimary)
   {
-    if (IsPhysicalPrimary() != fMCisPrimary) return kFALSE;
+    if (IsPhysicalPrimary(mcEvent,label) != fMCisPrimary) return kFALSE;
   }
   if (fCutMCPID)
   {
-    Int_t pdgCode = fMCparticle->PdgCode();
-    if (fMCPID != pdgCode) return kFALSE;
+    Int_t pdgCode = mcparticle->PdgCode();
+    if (fIgnoreSignInPID) 
+    {
+      if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
+    }
+    else 
+    {
+      if (fMCPID != pdgCode) return kFALSE;
+    }
   }
   if ( fCutMCprocessType )
   {
-    TParticle* particle = fMCparticle->Particle();
+    TParticle* particle = mcparticle->Particle();
     Int_t processID = particle->GetUniqueID();
     if (processID != fMCprocessType ) return kFALSE;
   }
   return kTRUE;
 }
+//-----------------------------------------------------------------------
+Bool_t AliFlowTrackCuts::PassesMCcuts()
+{
+  if (!fMCevent) return kFALSE;
+  if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
+  fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
+  return PassesMCcuts(fMCevent,fTrackLabel);
+}
 
 //-----------------------------------------------------------------------
 Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
@@ -219,7 +291,7 @@ Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
   //  start by preparing the track parameters to cut on //////////
   ////////////////////////////////////////////////////////////////
   //clean up from last iteration
-  if (fCleanupTrack) delete fTrack; fTrack=NULL; 
+  fTrack=NULL; 
 
   //get the label and the mc particle
   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
@@ -239,8 +311,9 @@ Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
   ////////////////////////////////////////////////////////////////
   ////////////////////////////////////////////////////////////////
 
+  if (!fTrack) return kFALSE;
   Bool_t pass=kTRUE;
-  //check the common cuts for the current particle (MC,AOD,ESD)
+  //check the common cuts for the current particle fTrack (MC,AOD,ESD)
   if (fCutPt) {if (fTrack->Pt() < fPtMin || fTrack->Pt() >= fPtMax ) pass=kFALSE;}
   if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
   if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
@@ -255,11 +328,26 @@ Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
   //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
 
   //when additionally MC info is required
-  if (!PassesMCcuts()) pass=kFALSE;
+  if (fCutMC && !PassesMCcuts()) pass=kFALSE;
 
   //check all else for ESDs using aliesdtrackcuts
   if (esdTrack && (fParamType!=kMC) ) 
+  {
+    if (fIgnoreTPCzRange)
+    {
+      const AliExternalTrackParam* pin = esdTrack->GetOuterParam();
+      const AliExternalTrackParam* pout = esdTrack->GetInnerParam();
+      if (pin&&pout)
+      {
+        Double_t zin = pin->GetZ();
+        Double_t zout = pout->GetZ();
+        if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
+        if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
+        if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
+      }
+    }
     if (!fAliESDtrackCuts->IsSelected(static_cast<AliESDtrack*>(fTrack))) pass=kFALSE;
+  }
 
   return pass; //true by default, if we didn't set any cuts
 }
@@ -271,7 +359,6 @@ void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
   switch (fParamType)
   {
     default:
-      fCleanupTrack = kFALSE;
       fTrack = track;
   }
 }
@@ -284,12 +371,16 @@ void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
   {
     case kGlobal:
       fTrack = track;
-      fCleanupTrack = kFALSE;
       break;
     case kESD_TPConly:
-      fTrack = new AliESDtrack();
-      track->FillTPCOnlyTrack(*(static_cast<AliESDtrack*>(fTrack)));
-      fCleanupTrack = kTRUE;
+      if (!track->FillTPCOnlyTrack(fTPCtrack))
+      {
+        fTrack=NULL;
+        fMCparticle=NULL;
+        fTrackLabel=-1;
+        return;
+      }
+      fTrack = &fTPCtrack;
       //recalculate the label and mc particle, they may differ as TPClabel != global label
       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
@@ -297,7 +388,6 @@ void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
       break;
     default:
       fTrack = track;
-      fCleanupTrack = kFALSE;
   }
 }
 
@@ -330,16 +420,51 @@ AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
 {
   //get a flow track constructed from whatever we applied cuts on
   //caller is resposible for deletion
+  //if construction fails return NULL
   AliFlowTrack* flowtrack=NULL;
+  TParticle *tmpTParticle=NULL;
+  AliMCParticle* tmpAliMCParticle=NULL;
   if (fParamType==kESD_SPDtracklet)
   {
-    flowtrack = new AliFlowTrack();
-    flowtrack->SetPhi(fTrackPhi);
-    flowtrack->SetEta(fTrackEta);
+    switch (fParamMix)
+    {
+      case kPure:
+        flowtrack = new AliFlowTrack();
+        flowtrack->SetPhi(fTrackPhi);
+        flowtrack->SetEta(fTrackEta);
+        break;
+      case kTrackWithMCkine:
+        if (!fMCparticle) return NULL;
+        flowtrack = new AliFlowTrack();
+        flowtrack->SetPhi( fMCparticle->Phi() );
+        flowtrack->SetEta( fMCparticle->Eta() );
+        flowtrack->SetPt( fMCparticle->Pt() );
+        break;
+      case kTrackWithMCpt:
+        if (!fMCparticle) return NULL;
+        flowtrack = new AliFlowTrack();
+        flowtrack->SetPhi(fTrackPhi);
+        flowtrack->SetEta(fTrackEta);
+        flowtrack->SetPt(fMCparticle->Pt());
+        break;
+      case kTrackWithPtFromFirstMother:
+        if (!fMCparticle) return NULL;
+        flowtrack = new AliFlowTrack();
+        flowtrack->SetPhi(fTrackPhi);
+        flowtrack->SetEta(fTrackEta);
+        tmpTParticle = fMCparticle->Particle();
+        tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
+        flowtrack->SetPt(tmpAliMCParticle->Pt());
+      default:
+        flowtrack = new AliFlowTrack();
+        flowtrack->SetPhi(fTrackPhi);
+        flowtrack->SetEta(fTrackEta);
+    }
     flowtrack->SetSource(AliFlowTrack::kFromTracklet);
   }
   else
   {
+    if (!fTrack) return NULL;
     switch(fParamMix)
     {
       case kPure:
@@ -350,7 +475,18 @@ AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
         break;
       case kTrackWithMCPID:
         flowtrack = new AliFlowTrack(fTrack);
+        //flowtrack->setPID(...) from mc, when implemented
         break;
+      case kTrackWithMCpt:
+        if (!fMCparticle) return NULL;
+        flowtrack = new AliFlowTrack(fTrack);
+        flowtrack->SetPt(fMCparticle->Pt());
+      case kTrackWithPtFromFirstMother:
+        if (!fMCparticle) return NULL;
+        flowtrack = new AliFlowTrack(fTrack);
+        tmpTParticle = fMCparticle->Particle();
+        tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
+        flowtrack->SetPt(tmpAliMCParticle->Pt());
       default:
         flowtrack = new AliFlowTrack(fTrack);
     }
@@ -366,7 +502,23 @@ AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
 {
   //check if current particle is a physical primary
-  return fMCevent->IsPhysicalPrimary(fTrackLabel);
+  if (!fMCevent) return kFALSE;
+  if (fTrackLabel<0) return kFALSE;
+  return IsPhysicalPrimary(fMCevent, fTrackLabel);
+}
+
+//-----------------------------------------------------------------------
+Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label)
+{
+  //check if current particle is a physical primary
+  Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
+  if (!physprim) return kFALSE;
+  AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
+  if (!track) return kFALSE;
+  TParticle* particle = track->Particle();
+  Bool_t transported = particle->TestBit(kTransportBit);
+  //printf("prim: %s, transp: %s\n",(physprim)?"YES":"NO ",(transported)?"YES":"NO ");
+  return (physprim && transported);
 }
 
 //-----------------------------------------------------------------------
@@ -393,3 +545,46 @@ const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
 void AliFlowTrackCuts::DefineHistograms()
 {
 }
+
+//-----------------------------------------------------------------------
+Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
+{
+  //get the number of tracks in the input event according source
+  //selection (ESD tracks, tracklets, MC particles etc.)
+  AliESDEvent* esd=NULL;
+  switch (fParamType)
+  {
+    case kESD_SPDtracklet:
+      esd = dynamic_cast<AliESDEvent*>(fEvent);
+      if (!esd) return 0;
+      return esd->GetMultiplicity()->GetNumberOfTracklets();
+    case kMC:
+      if (!fMCevent) return 0;
+      return fMCevent->GetNumberOfTracks();
+    default:
+      if (!fEvent) return 0;
+      return fEvent->GetNumberOfTracks();
+  }
+  return 0;
+}
+
+//-----------------------------------------------------------------------
+TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
+{
+  //get the input object according the data source selection:
+  //(esd tracks, traclets, mc particles,etc...)
+  AliESDEvent* esd=NULL;
+  switch (fParamType)
+  {
+    case kESD_SPDtracklet:
+      esd = dynamic_cast<AliESDEvent*>(fEvent);
+      if (!esd) return NULL;
+      return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
+    case kMC:
+      if (!fMCevent) return NULL;
+      return fMCevent->GetTrack(i);
+    default:
+      if (!fEvent) return NULL;
+      return fEvent->GetTrack(i);
+  }
+}