]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ANALYSIS/AliAnalysisUtils.cxx
prettier logs
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisUtils.cxx
index 988c9607cc0df1c6d5b9a0286f3910e6b6489c8e..15fef1d74d790c569be83fc2d1e99c70e2a99469 100644 (file)
@@ -4,6 +4,12 @@
 #include "AliVVertex.h"
 #include "AliLog.h"
 #include "AliAODVertex.h"
+#include "AliVTrack.h"
+#include "AliVEvent.h"
+#include <TMatrixDSym.h>
+#include <TMath.h>
+#include "AliVMultiplicity.h"
+#include "AliPPVsMultUtils.h"
 
 #include "AliAnalysisUtils.h"
 
@@ -14,7 +20,22 @@ 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),
+  fASPDCvsTCut(65.),
+  fBSPDCvsTCut(4.),
+  fPPVsMultUtils(0x0)
 {
   // Default contructor
 }
@@ -23,7 +44,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;
@@ -93,7 +114,8 @@ Bool_t AliAnalysisUtils::IsFirstEventInChunk(AliVEvent *event)
 
   if(fisAOD){
     AliAODHeader *aodheader = 0x0;
-    aodheader = aod->GetHeader();
+    aodheader = dynamic_cast<AliAODHeader*>(aod->GetHeader());
+    if(!aodheader) AliFatal("Not a standard AOD");
     if(!aodheader){
       AliFatal("AOD header not there ?!");
       return kFALSE;
@@ -102,6 +124,157 @@ 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)?((AliVAODHeader*)aod->GetHeader())->GetIRInt2ClosestInteractionMap():esd->GetHeader()->GetIRInt2ClosestInteractionMap();
+  if (bc2 != 0)
+    return kTRUE;
+  
+  Int_t bc1 = (aod)?((AliVAODHeader*)aod->GetHeader())->GetIRInt1ClosestInteractionMap():esd->GetHeader()->GetIRInt1ClosestInteractionMap();
+  if (bc1 != 0)
+    return kTRUE;
+  
+  return kFALSE;
+}
+
+
+//______________________________________________________________________
+Bool_t AliAnalysisUtils::IsSPDClusterVsTrackletBG(AliVEvent *event){
+  Int_t nClustersLayer0 = event->GetNumberOfITSClusters(0);
+  Int_t nClustersLayer1 = event->GetNumberOfITSClusters(1);
+  Int_t nTracklets      = event->GetMultiplicity()->GetNumberOfTracklets();
+  if (nClustersLayer0 + nClustersLayer1 > fASPDCvsTCut + nTracklets*fBSPDCvsTCut) 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; 
+
+}
+
+//______________________________________________________________________
+Float_t AliAnalysisUtils::GetMultiplicityPercentile(AliVEvent *event, TString lMethod ){
+  if(!fPPVsMultUtils)
+    fPPVsMultUtils=new AliPPVsMultUtils();
+  if( (event->InheritsFrom("AliAODEvent")) || (event->InheritsFrom("AliESDEvent")) ) return fPPVsMultUtils->GetMultiplicityPercentile(event,lMethod);
+  else {
+    AliFatal("Event is neither of AOD nor ESD type"); 
+    return -999.;
+  }
+}