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