Adding pile-up removal cuts for p-Pb data (Leonardo Milano).
authorakalweit <akalweit@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 30 Apr 2013 16:03:29 +0000 (16:03 +0000)
committerakalweit <akalweit@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 30 Apr 2013 16:03:29 +0000 (16:03 +0000)
ANALYSIS/AliAnalysisUtils.cxx
ANALYSIS/AliAnalysisUtils.h

index 988c960..225c2b0 100644 (file)
@@ -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; 
+
+}
index 6f3480d..8218dc7 100644 (file)
@@ -1,18 +1,20 @@
 #ifndef ALIANALYSISUTILS_H
 #define ALIANALYSISUTILS_H
 
-//////////////////////////////////////////////////////////////////////////////
-//                                                                          //
-// Class with functions useful for different analyses                       //
+////////////////////////////////////////////////////////////////////
+//                                                                           //
+// Class with functions useful for different analyses                      //
 // - vertex selection                                                       //
-//    * 2013 pA default cuts                                                //
-// - identification of the fist event of the chunk                          //
+//    * 2013 pA default cuts                                             //
+// - identification of the fist event of the chunk                         //
+// - identification pileup events                                           //
 //                                                                          //
-//////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////
 
 #include "TObject.h"
 
 class AliVEvent;
+class AliVVertex;
 
 class AliAnalysisUtils : public TObject {
 
@@ -20,24 +22,64 @@ class AliAnalysisUtils : public TObject {
 
   AliAnalysisUtils();
   virtual ~AliAnalysisUtils(){};
-
+  
   Bool_t IsVertexSelected2013pA(AliVEvent *event);
   Bool_t IsFirstEventInChunk(AliVEvent *event);
-
+  
+  Bool_t IsPileUpEvent(AliVEvent *event); //to be used in the analysis
+  Bool_t IsPileUpMV(AliVEvent *event); //MV pileup selection implemented here
+  Bool_t IsPileUpSPD(AliVEvent *event); //this calls IsPileUpFromSPD
+  Bool_t IsOutOfBunchPileUp(AliVEvent *event); //out-of-bunch pileup rejection using trigger information
+  
+  Double_t GetWDist(const AliVVertex* v0, const AliVVertex* v1);
+  
   void SetMinVtxContr(Int_t contr=1) {fMinVtxContr=contr;}
   void SetMaxVtxZ(Float_t z=1e6) {fMaxVtxZ=z;}
   void SetCutOnZVertexSPD(Bool_t iscut=true) { fCutOnZVertexSPD = iscut; }
-
+  
+  //general pileup selection settings
+  void SetUseMVPlpSelection(Bool_t useMVPlpSelection) { fUseMVPlpSelection = useMVPlpSelection;}
+  void SetUseOutOfBunchPileUp(Bool_t useOutOfBunchPileUp) { fUseOutOfBunchPileUp = useOutOfBunchPileUp;}
+  
+  //Multi Vertex pileup selection
+  void SetMinPlpContribMV(Int_t minPlpContribMV) { fMinPlpContribMV = minPlpContribMV;}
+  void SetMaxPlpChi2MV(Float_t maxPlpChi2MV) { fMaxPlpChi2MV = maxPlpChi2MV;}
+  void SetMinWDistMV(Float_t minWDistMV) { fMinWDistMV = minWDistMV;}
+  void SetCheckPlpFromDifferentBCMV(Bool_t checkPlpFromDifferentBCMV) { fCheckPlpFromDifferentBCMV = checkPlpFromDifferentBCMV;}
+  //SPD Pileup slection
+  void SetMinPlpContribSPD(Int_t minPlpContribSPD) { fMinPlpContribSPD = minPlpContribSPD;}
+  void SetMinPlpZdistSPD(Float_t minPlpZdistSPD) { fMinPlpZdistSPD = minPlpZdistSPD;}
+  void SetnSigmaPlpZdistSPD(Float_t nSigmaPlpZdistSPD) { fnSigmaPlpZdistSPD = nSigmaPlpZdistSPD;}
+  void SetnSigmaPlpDiamXYSPD(Float_t nSigmaPlpDiamXYSPD) { fnSigmaPlpDiamXYSPD = nSigmaPlpDiamXYSPD;}
+  void SetnSigmaPlpDiamZSPD(Float_t nSigmaPlpDiamZSPD) { fnSigmaPlpDiamZSPD = nSigmaPlpDiamZSPD;}
+  void SetUseSPDCutInMultBins(Bool_t useSPDCutInMultBins) { fUseSPDCutInMultBins = useSPDCutInMultBins;}
+  
+  
  private:
-
+  
   Bool_t fisAOD; // flag for AOD:1 or ESD:0
-
+  
   Int_t fMinVtxContr; // minimum vertex contributors
   Float_t fMaxVtxZ;   // maximum |z| of primary vertex
-
+  
   Bool_t fCutOnZVertexSPD; // 0: no cut, 1: |zvtx-SPD - zvtx-TPC|<0.5cm
-
-  ClassDef(AliAnalysisUtils,0) // base helper class
+  
+  Bool_t  fUseMVPlpSelection; //switch between SPD and MV pileup
+  Bool_t  fUseOutOfBunchPileUp; //BC information from the trigger is used to tag out-of-bunch pileup
+  
+  Int_t    fMinPlpContribMV; //minimum contributors to the pilup vertices, multi-vertex
+  Float_t  fMaxPlpChi2MV; //minimum value of Chi2perNDF of the pileup vertex, multi-vertex
+  Float_t  fMinWDistMV; //minimum of the sqrt of weighted distance between the primary and the pilup vertex, multi-vertex
+  Bool_t  fCheckPlpFromDifferentBCMV; //pileup from different BC (Bunch Crossings)
+  
+  Int_t    fMinPlpContribSPD; //minimum contributors to the pilup vertices, SPD
+  Float_t  fMinPlpZdistSPD;  //minimum distance for the SPD pileup vertex
+  Float_t  fnSigmaPlpZdistSPD;  //cut on Z for SPD pileup
+  Float_t  fnSigmaPlpDiamXYSPD;  //cut on nsigma diamond XY for SPD pileup
+  Float_t  fnSigmaPlpDiamZSPD;  //cut on nsigma diamond Z for SPD pileup
+  Bool_t  fUseSPDCutInMultBins;  //use IsPileupFromSPDInMultBins instead of IsPileupFromSPD
+  
+  ClassDef(AliAnalysisUtils,1) // base helper class
 };
 #endif