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