Changing to centrality dependent corrections

[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisUtils.cxx

diff --git a/ANALYSIS/AliAnalysisUtils.cxx b/ANALYSIS/AliAnalysisUtils.cxx
index 988c9607cc0df1c6d5b9a0286f3910e6b6489c8e..225c2b0504a7dfd776a9cc1d03ed9c347d5f4738 100644 (file)
--- a/ANALYSIS/AliAnalysisUtils.cxx
+++ b/ANALYSIS/AliAnalysisUtils.cxx
@@ -4,6 +4,10 @@
 #include "AliVVertex.h"
 #include "AliLog.h"
 #include "AliAODVertex.h"
+#include "AliVTrack.h"
+#include "AliVEvent.h"
+#include <TMatrixDSym.h>
+#include <TMath.h>
 
 #include "AliAnalysisUtils.h"
 
@@ -14,7 +18,19 @@ AliAnalysisUtils::AliAnalysisUtils():TObject(),
   fisAOD(kTRUE),
   fMinVtxContr(0),
   fMaxVtxZ(10.),
-  fCutOnZVertexSPD(kTRUE)
+  fCutOnZVertexSPD(kTRUE),
+  fUseMVPlpSelection(kFALSE),
+  fUseOutOfBunchPileUp(kFALSE),
+  fMinPlpContribMV(5),
+  fMaxPlpChi2MV(5.),
+  fMinWDistMV(15.),
+  fCheckPlpFromDifferentBCMV(kFALSE),
+  fMinPlpContribSPD(5),
+  fMinPlpZdistSPD(0.8),
+  fnSigmaPlpZdistSPD(3.),
+  fnSigmaPlpDiamXYSPD(2.),
+  fnSigmaPlpDiamZSPD(5.),
+  fUseSPDCutInMultBins(kFALSE)
 {
   // Default contructor
 }
@@ -23,7 +39,7 @@ AliAnalysisUtils::AliAnalysisUtils():TObject(),
 Bool_t AliAnalysisUtils::IsVertexSelected2013pA(AliVEvent *event)
 {
   Bool_t accept = kFALSE;
-
+  
   // Check whether the event is of AOD or ESD type
   AliAODEvent *aod = 0x0;
   AliESDEvent *esd = 0x0;
@@ -102,6 +118,135 @@ Bool_t AliAnalysisUtils::IsFirstEventInChunk(AliVEvent *event)
   } else {
     if(esd->GetEventNumberInFile()==0) accept = kTRUE;
   }
-
+  
   return accept;
 }
+
+//______________________________________________________________________
+Bool_t AliAnalysisUtils::IsPileUpEvent(AliVEvent *event)
+{
+  Bool_t isPileUp=kFALSE;
+  //check for multiple vertices
+  if(fUseMVPlpSelection)isPileUp=IsPileUpMV(event);
+  else isPileUp=IsPileUpSPD(event);
+  //check for different BC 
+  if(fUseOutOfBunchPileUp && IsOutOfBunchPileUp(event))isPileUp=kTRUE;
+  
+  return isPileUp;
+}
+
+//______________________________________________________________________
+Bool_t AliAnalysisUtils::IsPileUpMV(AliVEvent *event)
+{
+  // check for multi-vertexer pile-up
+  const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
+  const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
+  //
+  if (!aod && !esd) {
+    AliFatal("Event is neither of AOD nor ESD type");
+    return kFALSE;
+  }
+  //
+  const AliVVertex* vtPrm = 0;
+  const AliVVertex* vtPlp = 0;
+  Int_t nPlp = 0;
+  //
+  if (aod) {
+    if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
+    vtPrm = aod->GetPrimaryVertex();
+    if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
+  }
+  else {
+    if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
+    vtPrm = esd->GetPrimaryVertexTracks();
+    if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
+  }
+  Int_t bcPrim = vtPrm->GetBC();
+  //
+  for (Int_t ipl=0;ipl<nPlp;ipl++) {
+    vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
+    //
+    if (vtPlp->GetNContributors() < fMinPlpContribMV) continue;
+    if (vtPlp->GetChi2perNDF() > fMaxPlpChi2MV) continue;
+    if(fCheckPlpFromDifferentBCMV)
+      {
+       Int_t bcPlp = vtPlp->GetBC();
+       if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
+      }
+    //
+    Double_t wDst = GetWDist(vtPrm,vtPlp);
+    if (wDst<fMinWDistMV) continue;
+    //
+    return kTRUE; // pile-up: well separated vertices
+  }
+  //
+  return kFALSE;
+  //
+}
+
+//______________________________________________________________________
+Bool_t AliAnalysisUtils::IsPileUpSPD(AliVEvent *event)
+{
+  // check for SPD pile-up
+  const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
+  const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
+  //
+  if (!aod && !esd) {
+    AliFatal("Event is neither of AOD nor ESD type");
+    return kFALSE;
+  }
+  //
+  if (aod) return (fUseSPDCutInMultBins)?aod->IsPileupFromSPDInMultBins():aod->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
+  else return (fUseSPDCutInMultBins)?esd->IsPileupFromSPDInMultBins():esd->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
+}
+
+//______________________________________________________________________
+Bool_t AliAnalysisUtils::IsOutOfBunchPileUp(AliVEvent *event)
+{
+  // check for SPD pile-up
+  const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
+  const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
+  //
+  if (!aod && !esd) {
+    AliFatal("Event is neither of AOD nor ESD type");
+    return kFALSE;
+  }
+  Int_t bc2 = (aod)?aod->GetHeader()->GetIRInt2ClosestInteractionMap():esd->GetHeader()->GetIRInt2ClosestInteractionMap();
+  if (bc2 != 0)
+    return kTRUE;
+  
+  Int_t bc1 = (aod)?aod->GetHeader()->GetIRInt1ClosestInteractionMap():esd->GetHeader()->GetIRInt1ClosestInteractionMap();
+  if (bc1 != 0)
+    return kTRUE;
+  
+  return kFALSE;
+}
+
+//______________________________________________________________________
+Double_t AliAnalysisUtils::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
+{
+  // calculate sqrt of weighted distance to other vertex
+  if (!v0 || !v1) {
+    printf("One of vertices is not valid\n");
+    return 0;
+  }
+  static TMatrixDSym vVb(3);
+  Double_t dist = -1;
+  Double_t dx = v0->GetX()-v1->GetX();
+  Double_t dy = v0->GetY()-v1->GetY();
+  Double_t dz = v0->GetZ()-v1->GetZ();
+  Double_t cov0[6],cov1[6];
+  v0->GetCovarianceMatrix(cov0);
+  v1->GetCovarianceMatrix(cov1);
+  vVb(0,0) = cov0[0]+cov1[0];
+  vVb(1,1) = cov0[2]+cov1[2];
+  vVb(2,2) = cov0[5]+cov1[5];
+  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
+  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
+  vVb.InvertFast();
+  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
+  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
+    +    2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
+  return dist>0 ? TMath::Sqrt(dist) : -1; 
+
+}