ANALYSIS/AliAnalysisUtils.cxx

   1 #include "AliVEvent.h"
   2 #include "AliESDEvent.h"
   3 #include "AliAODEvent.h"
   4 #include "AliVVertex.h"
   5 #include "AliLog.h"
   6 #include "AliAODVertex.h"
   7 #include "AliVTrack.h"
   8 #include "AliVEvent.h"
   9 #include <TMatrixDSym.h>
  10 #include <TMath.h>
  11 #include "AliVMultiplicity.h"
  12
  13 #include "AliAnalysisUtils.h"
  14
  15 ClassImp(AliAnalysisUtils)
  16
  17 //______________________________________________________________________
  18 AliAnalysisUtils::AliAnalysisUtils():TObject(),
  19   fisAOD(kTRUE),
  20   fMinVtxContr(0),
  21   fMaxVtxZ(10.),
  22   fCutOnZVertexSPD(kTRUE),
  23   fUseMVPlpSelection(kFALSE),
  24   fUseOutOfBunchPileUp(kFALSE),
  25   fMinPlpContribMV(5),
  26   fMaxPlpChi2MV(5.),
  27   fMinWDistMV(15.),
  28   fCheckPlpFromDifferentBCMV(kFALSE),
  29   fMinPlpContribSPD(5),
  30   fMinPlpZdistSPD(0.8),
  31   fnSigmaPlpZdistSPD(3.),
  32   fnSigmaPlpDiamXYSPD(2.),
  33   fnSigmaPlpDiamZSPD(5.),
  34   fUseSPDCutInMultBins(kFALSE),
  35   fASPDCvsTCut(65.),
  36   fBSPDCvsTCut(4.)
  37 {
  38   // Default contructor
  39 }
  40
  41 //______________________________________________________________________
  42 Bool_t AliAnalysisUtils::IsVertexSelected2013pA(AliVEvent *event)
  43 {
  44   Bool_t accept = kFALSE;
  45
  46   // Check whether the event is of AOD or ESD type
  47   AliAODEvent *aod = 0x0;
  48   AliESDEvent *esd = 0x0;
  49   aod = dynamic_cast<AliAODEvent*>(event);
  50   esd = dynamic_cast<AliESDEvent*>(event);
  51
  52   if(aod) {
  53     fisAOD = kTRUE;
  54   } else {
  55     fisAOD = kFALSE;
  56     if(!esd) {
  57       AliFatal("Event is neither of AOD nor ESD type");
  58       return accept;
  59     }
  60   }
  61
  62   const AliVVertex *trkVtx = fisAOD ?
  63     dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
  64     dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex()) ;
  65   if(!trkVtx || trkVtx->GetNContributors()<=fMinVtxContr){
  66     accept = kFALSE;
  67     return accept;
  68   }
  69
  70   TString vtxTtl = trkVtx->GetTitle();
  71   if (!vtxTtl.Contains("VertexerTracks")) return accept;
  72
  73   Float_t zvtx = trkVtx->GetZ();
  74   const AliVVertex* spdVtx = fisAOD ?
  75     dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) :
  76     dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD()) ;
  77   if (spdVtx->GetNContributors()<=fMinVtxContr) return accept;
  78
  79   TString vtxTyp = spdVtx->GetTitle();
  80   Double_t cov[6]={0};
  81   spdVtx->GetCovarianceMatrix(cov);
  82   Double_t zRes = TMath::Sqrt(cov[5]);
  83   if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return accept;
  84   if (fCutOnZVertexSPD && TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return accept;
  85
  86   if (TMath::Abs(zvtx) > fMaxVtxZ) return accept;
  87
  88   return kTRUE;
  89 }
  90
  91 //______________________________________________________________________
  92 Bool_t AliAnalysisUtils::IsFirstEventInChunk(AliVEvent *event)
  93 {
  94
  95   Bool_t accept = kFALSE;
  96
  97   // Check whether the event is of AOD or ESD type
  98   AliAODEvent *aod = 0x0;
  99   AliESDEvent *esd = 0x0;
 100   aod = dynamic_cast<AliAODEvent*>(event);
 101   esd = dynamic_cast<AliESDEvent*>(event);
 102
 103   if(aod) {
 104     fisAOD = kTRUE;
 105   } else {
 106     fisAOD = kFALSE;
 107     if(!esd) {
 108       AliFatal("Event is neither of AOD nor ESD type");
 109       return accept;
 110     }
 111   }
 112
 113   if(fisAOD){
 114     AliAODHeader *aodheader = 0x0;
 115     aodheader = aod->GetHeader();
 116     if(!aodheader){
 117       AliFatal("AOD header not there ?!");
 118       return kFALSE;
 119     }
 120     if(aodheader->GetEventNumberESDFile()==0) accept = kTRUE;
 121   } else {
 122     if(esd->GetEventNumberInFile()==0) accept = kTRUE;
 123   }
 124
 125   return accept;
 126 }
 127
 128 //______________________________________________________________________
 129 Bool_t AliAnalysisUtils::IsPileUpEvent(AliVEvent *event)
 130 {
 131   Bool_t isPileUp=kFALSE;
 132   //check for multiple vertices
 133   if(fUseMVPlpSelection)isPileUp=IsPileUpMV(event);
 134   else isPileUp=IsPileUpSPD(event);
 135   //check for different BC
 136   if(fUseOutOfBunchPileUp && IsOutOfBunchPileUp(event))isPileUp=kTRUE;
 137
 138   return isPileUp;
 139 }
 140
 141 //______________________________________________________________________
 142 Bool_t AliAnalysisUtils::IsPileUpMV(AliVEvent *event)
 143 {
 144   // check for multi-vertexer pile-up
 145   const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
 146   const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
 147   //
 148   if (!aod && !esd) {
 149     AliFatal("Event is neither of AOD nor ESD type");
 150     return kFALSE;
 151   }
 152   //
 153   const AliVVertex* vtPrm = 0;
 154   const AliVVertex* vtPlp = 0;
 155   Int_t nPlp = 0;
 156   //
 157   if (aod) {
 158     if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
 159     vtPrm = aod->GetPrimaryVertex();
 160     if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
 161   }
 162   else {
 163     if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
 164     vtPrm = esd->GetPrimaryVertexTracks();
 165     if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
 166   }
 167   Int_t bcPrim = vtPrm->GetBC();
 168   //
 169   for (Int_t ipl=0;ipl<nPlp;ipl++) {
 170     vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
 171     //
 172     if (vtPlp->GetNContributors() < fMinPlpContribMV) continue;
 173     if (vtPlp->GetChi2perNDF() > fMaxPlpChi2MV) continue;
 174     if(fCheckPlpFromDifferentBCMV)
 175       {
 176         Int_t bcPlp = vtPlp->GetBC();
 177         if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
 178       }
 179     //
 180     Double_t wDst = GetWDist(vtPrm,vtPlp);
 181     if (wDst<fMinWDistMV) continue;
 182     //
 183     return kTRUE; // pile-up: well separated vertices
 184   }
 185   //
 186   return kFALSE;
 187   //
 188 }
 189
 190 //______________________________________________________________________
 191 Bool_t AliAnalysisUtils::IsPileUpSPD(AliVEvent *event)
 192 {
 193   // check for SPD pile-up
 194   const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
 195   const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
 196   //
 197   if (!aod && !esd) {
 198     AliFatal("Event is neither of AOD nor ESD type");
 199     return kFALSE;
 200   }
 201   //
 202   if (aod) return (fUseSPDCutInMultBins)?aod->IsPileupFromSPDInMultBins():aod->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
 203   else return (fUseSPDCutInMultBins)?esd->IsPileupFromSPDInMultBins():esd->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
 204 }
 205
 206 //______________________________________________________________________
 207 Bool_t AliAnalysisUtils::IsOutOfBunchPileUp(AliVEvent *event)
 208 {
 209   // check for SPD pile-up
 210   const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
 211   const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
 212   //
 213   if (!aod && !esd) {
 214     AliFatal("Event is neither of AOD nor ESD type");
 215     return kFALSE;
 216   }
 217   Int_t bc2 = (aod)?aod->GetHeader()->GetIRInt2ClosestInteractionMap():esd->GetHeader()->GetIRInt2ClosestInteractionMap();
 218   if (bc2 != 0)
 219     return kTRUE;
 220
 221   Int_t bc1 = (aod)?aod->GetHeader()->GetIRInt1ClosestInteractionMap():esd->GetHeader()->GetIRInt1ClosestInteractionMap();
 222   if (bc1 != 0)
 223     return kTRUE;
 224
 225   return kFALSE;
 226 }
 227
 228
 229 //______________________________________________________________________
 230 Bool_t AliAnalysisUtils::IsSPDClusterVsTrackletBG(AliVEvent *event){
 231   Int_t nClustersLayer0 = event->GetNumberOfITSClusters(0);
 232   Int_t nClustersLayer1 = event->GetNumberOfITSClusters(1);
 233   Int_t nTracklets      = event->GetMultiplicity()->GetNumberOfTracklets();
 234   if (nClustersLayer0 + nClustersLayer1 > fASPDCvsTCut + nTracklets*fBSPDCvsTCut) return kTRUE;
 235   return kFALSE;
 236 }
 237
 238
 239 //______________________________________________________________________
 240 Double_t AliAnalysisUtils::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
 241 {
 242   // calculate sqrt of weighted distance to other vertex
 243   if (!v0 || !v1) {
 244     printf("One of vertices is not valid\n");
 245     return 0;
 246   }
 247   static TMatrixDSym vVb(3);
 248   Double_t dist = -1;
 249   Double_t dx = v0->GetX()-v1->GetX();
 250   Double_t dy = v0->GetY()-v1->GetY();
 251   Double_t dz = v0->GetZ()-v1->GetZ();
 252   Double_t cov0[6],cov1[6];
 253   v0->GetCovarianceMatrix(cov0);
 254   v1->GetCovarianceMatrix(cov1);
 255   vVb(0,0) = cov0[0]+cov1[0];
 256   vVb(1,1) = cov0[2]+cov1[2];
 257   vVb(2,2) = cov0[5]+cov1[5];
 258   vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
 259   vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
 260   vVb.InvertFast();
 261   if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
 262   dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
 263     +    2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
 264   return dist>0 ? TMath::Sqrt(dist) : -1;
 265
 266 }
 267