]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - CORRFW/AliCFTrackIsPrimaryCuts.cxx
Updates of trending macro for QA outputs from AliAnalysisTaskTOFqa:
[u/mrichter/AliRoot.git] / CORRFW / AliCFTrackIsPrimaryCuts.cxx
index 3d4534e66f780bb15e884c467d7ccf17934e67ec..a9675355a75e8b5e16cc19bf44269f8db7418444 100644 (file)
 // - accept or not accept daughter tracks of kink decays
 //
 // By default, the distance to 'vertex calculated from tracks' is used.
-// Optionally the TPC (TPC only tracks based) vertex can be used.
+// Optionally the SPD (tracklet based) or TPC (TPC only tracks based) vertex
+// can be used.
+// Note: the distance to the TPC-vertex is already stored in the ESD,
+// the distance to the SPD-vertex has to be re-calculated by propagating each
+// track while executing this cut.
 //
 // The cut values for these cuts are set with the corresponding set functions.
 // All cut classes provided by the correction framework are supposed to be
 #include <TBits.h>
 
 #include <AliESDtrack.h>
+#include <AliAODTrack.h>
+#include <AliESDEvent.h>
+#include <AliAODEvent.h>
 #include <AliLog.h>
-#include <AliExternalTrackParam.h>
 #include "AliCFTrackIsPrimaryCuts.h"
 
 ClassImp(AliCFTrackIsPrimaryCuts)
@@ -59,6 +65,8 @@ ClassImp(AliCFTrackIsPrimaryCuts)
 //__________________________________________________________________________________
 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts() :
   AliCFCutBase(),
+  fEvt(0x0),
+  fUseSPDvertex(0),
   fUseTPCvertex(0),
   fMinDCAToVertexXY(0),
   fMinDCAToVertexZ(0),
@@ -103,6 +111,8 @@ AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts() :
 //__________________________________________________________________________________
 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(Char_t* name, Char_t* title) :
   AliCFCutBase(name,title),
+  fEvt(0x0),
+  fUseSPDvertex(0),
   fUseTPCvertex(0),
   fMinDCAToVertexXY(0),
   fMinDCAToVertexZ(0),
@@ -147,6 +157,8 @@ AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(Char_t* name, Char_t* title) :
 //__________________________________________________________________________________
 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(const AliCFTrackIsPrimaryCuts& c) :
   AliCFCutBase(c),
+  fEvt(c.fEvt),
+  fUseSPDvertex(c.fUseSPDvertex),
   fUseTPCvertex(c.fUseTPCvertex),
   fMinDCAToVertexXY(c.fMinDCAToVertexXY),
   fMinDCAToVertexZ(c.fMinDCAToVertexZ),
@@ -196,6 +208,8 @@ AliCFTrackIsPrimaryCuts& AliCFTrackIsPrimaryCuts::operator=(const AliCFTrackIsPr
   //
   if (this != &c) {
     AliCFCutBase::operator=(c) ;
+    fEvt = c.fEvt;
+    fUseSPDvertex = c.fUseSPDvertex;
     fUseTPCvertex = c.fUseTPCvertex;
     fMinDCAToVertexXY = c.fMinDCAToVertexXY;
     fMinDCAToVertexZ = c.fMinDCAToVertexZ;
@@ -257,6 +271,7 @@ AliCFTrackIsPrimaryCuts::~AliCFTrackIsPrimaryCuts()
     for (Int_t i=0; i<kNHist; i++)
       if(fhQA[i][j])           delete fhQA[i][j];
   }
+  if(fEvt)                     delete fEvt;
   if(fBitmap)                  delete fBitmap;
   if(fhBinLimNSigma)           delete fhBinLimNSigma;
   if(fhBinLimRequireSigma)     delete fhBinLimRequireSigma;
@@ -274,6 +289,7 @@ void AliCFTrackIsPrimaryCuts::Initialise()
   //
   // sets everything to zero
   //
+  fUseSPDvertex = 0;
   fUseTPCvertex = 0;
   fMinDCAToVertexXY = 0;
   fMinDCAToVertexZ = 0;
@@ -346,14 +362,58 @@ void AliCFTrackIsPrimaryCuts::Copy(TObject &c) const
   TNamed::Copy(c);
 }
 //__________________________________________________________________________________
+void AliCFTrackIsPrimaryCuts::SetRecEventInfo(const TObject* evt) {
+  //
+  // Sets pointer to event information (AliESDEvent or AliAODEvent)
+  //
+  if (!evt) {
+    AliError("Pointer to AliVEvent !");
+    return;
+  }
+  TString className(evt->ClassName());
+  if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
+    AliError("argument must point to an AliESDEvent or AliAODEvent !");
+    return ;
+  }
+  fEvt = (AliVEvent*) evt;
+}
+//__________________________________________________________________________________
+void AliCFTrackIsPrimaryCuts::UseSPDvertex(Bool_t b) {
+  fUseSPDvertex = b;
+  if(fUseTPCvertex && fUseSPDvertex) {
+       fUseSPDvertex = kFALSE;
+       AliError("SPD and TPC vertex chosen. TPC vertex is preferred.");
+  }
+}
+//__________________________________________________________________________________
+void AliCFTrackIsPrimaryCuts::UseTPCvertex(Bool_t b) {
+  fUseTPCvertex = b;
+  if(fUseTPCvertex) fUseSPDvertex = kFALSE;
+}
+//__________________________________________________________________________________
 void AliCFTrackIsPrimaryCuts::GetDCA(AliESDtrack* esdTrack)
 {
   if (!esdTrack) return;
 
   Float_t b[2] = {0.,0.};
   Float_t bCov[3] = {0.,0.,0.};
-  if(!fUseTPCvertex) esdTrack->GetImpactParameters(b,bCov);
-  if( fUseTPCvertex) esdTrack->GetImpactParametersTPC(b,bCov);
+  if(!fUseSPDvertex && !fUseTPCvertex) esdTrack->GetImpactParameters(b,bCov);
+  if( fUseTPCvertex)                   esdTrack->GetImpactParametersTPC(b,bCov);
+
+  if( fUseSPDvertex) {
+       if (!fEvt) return;
+       AliESDEvent * evt = 0x0 ; 
+       evt = dynamic_cast<AliESDEvent*>(fEvt);
+       if (!evt) {
+         AliError("event not found");
+         return;
+       }
+       const AliESDVertex *vtx = evt->GetVertex();
+       const Double_t Bz = evt->GetMagneticField();
+       AliExternalTrackParam *cParam = 0;
+       Bool_t success = esdTrack->RelateToVertex(vtx, Bz, kVeryBig, cParam);
+       if (success) esdTrack->GetImpactParameters(b,bCov);
+  }
 
   if (bCov[0]<=0 || bCov[2]<=0) {
       bCov[0]=0; bCov[2]=0;
@@ -369,7 +429,7 @@ void AliCFTrackIsPrimaryCuts::GetDCA(AliESDtrack* esdTrack)
   }
 
   // get n_sigma
-  if(!fUseTPCvertex)
+  if(!fUseSPDvertex && !fUseTPCvertex)
        fDCA[5] = AliESDtrackCuts::GetSigmaToVertex(esdTrack);
 
   if(fUseTPCvertex) {
@@ -385,6 +445,71 @@ void AliCFTrackIsPrimaryCuts::GetDCA(AliESDtrack* esdTrack)
   return;
 }
 //__________________________________________________________________________________
+void AliCFTrackIsPrimaryCuts::GetDCA(AliAODTrack* aodTrack)
+{
+  if (!aodTrack) return;
+
+  Double_t p[3] = {0.,0.,0.};
+  Double_t cov[21] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
+
+  aodTrack->XYZAtDCA(p); // position at DCA
+  aodTrack->GetCovarianceXYZPxPyPz(cov);
+
+
+  fDCA[5] = -1; // n_sigma
+  if (p[0]==-999. || p[1]==-999. ||  p[2]==-999.) {
+    AliError("dca info not available !");
+    fDCA[0] = -999.; // impact parameter xy
+    fDCA[1] = -999.; // impact parameter z
+    return;
+  }
+
+  AliAODEvent * evt = 0x0;
+  evt = dynamic_cast<AliAODEvent*>(fEvt);
+  if (!evt) return;
+
+  // primary vertex is the "best": tracks, SPD or TPC vertex
+  AliAODVertex * primaryVertex = evt->GetVertex(0);
+  // dca = track postion - primary vertex position
+  p[0] = p[0] - primaryVertex->GetX();
+  p[1] = p[1] - primaryVertex->GetY();
+  p[2] = p[2] - primaryVertex->GetZ();
+  // calculate dca in transverse plane
+  Float_t b[2] = {0.,0.};
+  b[0] = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
+  b[1] = p[2];
+  // resolution
+  Double_t bCov[3] = {0.,0.,0.};
+  // how to calculate the resoultion in the transverse plane ?
+  bCov[0] = 0.; // to do: calculate cov in transverse plane
+  bCov[2] = 0.; // from cov in x and y, need to know correlation
+
+  if (bCov[0]<=0 || bCov[2]<=0) {
+      bCov[0]=0; bCov[2]=0;
+  }
+  fDCA[0] = b[0]; // impact parameter xy
+  fDCA[1] = b[1]; // impact parameter z
+  fDCA[2] = TMath::Sqrt(bCov[0]); // resolution xy
+  fDCA[3] = TMath::Sqrt(bCov[2]); // resolution z
+
+  if (!fAbsDCAToVertex) {
+       AliError("resolution of impact parameter not available, use absolute dca cut instead !");
+       if (fDCA[2] > 0) fDCA[0] = fDCA[0]/fDCA[2]; // normalised impact parameter xy
+       if (fDCA[3] > 0) fDCA[1] = fDCA[1]/fDCA[3]; // normalised impact parameter z
+  }
+
+  // get n_sigma
+  if (fDCA[2]==0 || fDCA[3]==0)
+       return;
+  fDCA[5] = 1000.;
+  Float_t d = TMath::Sqrt(TMath::Power(b[0]/fDCA[2],2) + TMath::Power(b[1]/fDCA[3],2));
+  if (TMath::Exp(-d * d / 2) < 1e-15)
+       return;
+  fDCA[5] = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
+
+  return;
+}
+//__________________________________________________________________________________
 void AliCFTrackIsPrimaryCuts::SelectionBitMap(TObject* obj)
 {
   //
@@ -401,113 +526,82 @@ void AliCFTrackIsPrimaryCuts::SelectionBitMap(TObject* obj)
     AliError("object must derived from AliVParticle !");
     return;
   }
+  
+  AliESDtrack * esdTrack = dynamic_cast<AliESDtrack*>(obj);
+  AliAODTrack * aodTrack = dynamic_cast<AliAODTrack*>(obj);
 
-  Bool_t isESDTrack = strcmp(obj->ClassName(),"AliESDtrack") == 0 ? kTRUE : kFALSE ;
-  Bool_t isAODTrack = strcmp(obj->ClassName(),"AliAODTrack") == 0 ? kTRUE : kFALSE ;
+  if (!(esdTrack || aodTrack)) {
+    AliError("object must be an ESDtrack or an AODtrack !");
+    return;
+  }
 
-  AliESDtrack * esdTrack = 0x0 ;
-  AliAODTrack * aodTrack = 0x0 ; 
-  if (isESDTrack) esdTrack = dynamic_cast<AliESDtrack*>(obj);
-  if (isAODTrack) aodTrack = dynamic_cast<AliAODTrack*>(obj);
+  Bool_t isESDTrack = kFALSE;
+  Bool_t isAODTrack = kFALSE;
+
+  if (esdTrack) isESDTrack = strcmp(obj->ClassName(),"AliESDtrack") == 0 ? kTRUE : kFALSE ;
+  if (aodTrack) isAODTrack = strcmp(obj->ClassName(),"AliAODTrack") == 0 ? kTRUE : kFALSE ;
 
   // get the track to vertex parameter for ESD track
   if (isESDTrack) GetDCA(esdTrack);
+  if (isAODTrack) GetDCA(aodTrack);
+
+  // check whether dca info is filled
+  Bool_t dcaInfo = 0;
+  if (fDCA[0]>-990. && fDCA[1]>-990.) dcaInfo = 1;
 
   Float_t bxy = 0, bz = 0;
-  if (isESDTrack) {
-       bxy = TMath::Abs(fDCA[0]);
-       bz  = TMath::Abs(fDCA[1]);
- }
+  bxy = TMath::Abs(fDCA[0]);
+  bz  = TMath::Abs(fDCA[1]);
+
   Float_t b2Dmin = 0, b2Dmax = 0;
-  if (isESDTrack) {
-    if (fMinDCAToVertexXY>0 && fMinDCAToVertexZ>0)
-       b2Dmin = fDCA[0]*fDCA[0]/fMinDCAToVertexXY/fMinDCAToVertexXY + fDCA[1]*fDCA[1]/fMinDCAToVertexZ/fMinDCAToVertexZ;
-    if (fMaxDCAToVertexXY>0 && fMaxDCAToVertexZ>0) 
-       b2Dmax = fDCA[0]*fDCA[0]/fMaxDCAToVertexXY/fMaxDCAToVertexXY + fDCA[1]*fDCA[1]/fMaxDCAToVertexZ/fMaxDCAToVertexZ;
-  }
+  if (fMinDCAToVertexXY>0 && fMinDCAToVertexZ>0)
+    b2Dmin = fDCA[0]*fDCA[0]/fMinDCAToVertexXY/fMinDCAToVertexXY + fDCA[1]*fDCA[1]/fMinDCAToVertexZ/fMinDCAToVertexZ;
+  if (fMaxDCAToVertexXY>0 && fMaxDCAToVertexZ>0) 
+    b2Dmax = fDCA[0]*fDCA[0]/fMaxDCAToVertexXY/fMaxDCAToVertexXY + fDCA[1]*fDCA[1]/fMaxDCAToVertexZ/fMaxDCAToVertexZ;
+
 
   // fill the bitmap
   Int_t iCutBit = 0;
 
-  if (isESDTrack) {
-    if (fDCAToVertex2D || (!fDCAToVertex2D && bxy >= fMinDCAToVertexXY && bxy <= fMaxDCAToVertexXY)) {
-      fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
+  if (!dcaInfo || fDCAToVertex2D || (!fDCAToVertex2D && bxy >= fMinDCAToVertexXY && bxy <= fMaxDCAToVertexXY))
+       fBitmap->SetBitNumber(iCutBit,kTRUE);
   iCutBit++;
 
-  if (isESDTrack) {
-    if (fDCAToVertex2D || (!fDCAToVertex2D && bz  >= fMinDCAToVertexZ && bz  <= fMaxDCAToVertexZ)) {
-      fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
+  if (!dcaInfo || fDCAToVertex2D || (!fDCAToVertex2D && bz  >= fMinDCAToVertexZ && bz  <= fMaxDCAToVertexZ))
+       fBitmap->SetBitNumber(iCutBit,kTRUE);
   iCutBit++;
 
-  if (isESDTrack) {
-    if (!fDCAToVertex2D || (fDCAToVertex2D && b2Dmin > 1  && b2Dmax < 1)) {
+  if (!dcaInfo || !fDCAToVertex2D || (fDCAToVertex2D && TMath::Sqrt(b2Dmin) > 1  && TMath::Sqrt(b2Dmax) < 1))
       fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
   iCutBit++;
 
-  if (isESDTrack) {
-    if (fDCA[5] >= fNSigmaToVertexMin && fDCA[5] <= fNSigmaToVertexMax) {
+  if (!dcaInfo || (fDCA[5] >= fNSigmaToVertexMin && fDCA[5] <= fNSigmaToVertexMax))
       fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
   iCutBit++;
 
-  if (isESDTrack) {
-    if (fDCA[2] < fSigmaDCAxy) {
+  if (!dcaInfo || fDCA[2] < fSigmaDCAxy)
       fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
   iCutBit++;
 
-  if (isESDTrack) {
-    if (fDCA[3] < fSigmaDCAz) {
+  if (!dcaInfo || fDCA[3] < fSigmaDCAz)
       fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
   iCutBit++;
 
-  if (isESDTrack) {
-    if (!fRequireSigmaToVertex || (fDCA[5]>=0 && fRequireSigmaToVertex)) {
+  if (!dcaInfo || !fRequireSigmaToVertex || (fDCA[5]>=0 && fRequireSigmaToVertex))
       fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-
   iCutBit++;
 
-  if (esdTrack) {
-    if (fAcceptKinkDaughters || (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)<=0)) {
+  if (!dcaInfo || fAcceptKinkDaughters || (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)<=0))
       fBitmap->SetBitNumber(iCutBit,kTRUE);
-    }
-  }
-  else fBitmap->SetBitNumber(iCutBit,kTRUE);
-  
   iCutBit++;
-  
+
   if (isAODTrack) {
     if (fAODType==AliAODTrack::kUndef || fAODType == aodTrack->GetType()) {
       fBitmap->SetBitNumber(iCutBit,kTRUE);
     }
   }
   else fBitmap->SetBitNumber(iCutBit,kTRUE);
-  
+
   return;
 }
 //__________________________________________________________________________________
@@ -526,6 +620,7 @@ Bool_t AliCFTrackIsPrimaryCuts::IsSelected(TObject* obj) {
        break;
     }
   }
+
   if (!isSelected) return kFALSE ;
   if (fIsQAOn) FillHistograms(obj,1);
   return kTRUE;
@@ -699,10 +794,10 @@ void AliCFTrackIsPrimaryCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_
   fhCutCorrelation->GetYaxis()->SetBinLabel(9,"AOD type");
 
   // book QA histograms
-  Char_t str[256];
+  Char_t str[5];
   for (Int_t i=0; i<kNStepQA; i++) {
-    if (i==0) sprintf(str," ");
-    else sprintf(str,"_cut");
+    if (i==0) snprintf(str,5," ");
+    else snprintf(str,5,"_cut");
 
     fhDcaXYvsDcaZ[i]            = new  TH2F(Form("%s_dcaXYvsDcaZ%s",GetName(),str),"",200,-10,10,200,-10,10);
     fhQA[kCutNSigmaToVertex][i]        = new TH1F(Form("%s_nSigmaToVertex%s",GetName(),str),"",fhNBinsNSigma-1,fhBinLimNSigma);
@@ -738,7 +833,6 @@ void AliCFTrackIsPrimaryCuts::FillHistograms(TObject* obj, Bool_t f)
   //
 
   if (!obj) return;
-
   Bool_t isESDTrack = strcmp(obj->ClassName(),"AliESDtrack") == 0 ? kTRUE : kFALSE ;
   Bool_t isAODTrack = strcmp(obj->ClassName(),"AliAODTrack") == 0 ? kTRUE : kFALSE ;
 
@@ -763,18 +857,21 @@ void AliCFTrackIsPrimaryCuts::FillHistograms(TObject* obj, Bool_t f)
       fhQA[kDcaXYnorm][f]->Fill(fDCA[0]);
 
     fhQA[kCutNSigmaToVertex][f]->Fill(fDCA[5]);
-    if (fDCA[5]<0 && fRequireSigmaToVertex) fhQA[kCutRequireSigmaToVertex][f]->Fill(0.);
-    if (!(fDCA[5]<0 && fRequireSigmaToVertex)) fhQA[kCutRequireSigmaToVertex][f]->Fill(1.);
-
-    if (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) fhQA[kCutAcceptKinkDaughters][f]->Fill(0.);
-    if (!(!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)) fhQA[kCutAcceptKinkDaughters][f]->Fill(0.);
+    if (fDCA[5]<0 && fRequireSigmaToVertex) 
+      fhQA[kCutRequireSigmaToVertex][f]->Fill(0.);
+    else 
+      fhQA[kCutRequireSigmaToVertex][f]->Fill(1.);
+
+    if (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
+      fhQA[kCutAcceptKinkDaughters][f]->Fill(0.);
+    else 
+      fhQA[kCutAcceptKinkDaughters][f]->Fill(1.);
   }
 
   // fill cut statistics and cut correlation histograms with information from the bitmap
   if (f) return;
 
   // Get the bitmap of the single cuts
-  if ( !obj ) return;
   SelectionBitMap(obj);
 
   // Number of single cuts in this class