2 #include "AliESDEvent.h"
3 #include "AliAODEvent.h"
4 #include "AliVVertex.h"
6 #include "AliAODVertex.h"
9 #include <TMatrixDSym.h>
12 #include "AliAnalysisUtils.h"
14 ClassImp(AliAnalysisUtils)
16 //______________________________________________________________________
17 AliAnalysisUtils::AliAnalysisUtils():TObject(),
21 fCutOnZVertexSPD(kTRUE),
22 fUseMVPlpSelection(kFALSE),
23 fUseOutOfBunchPileUp(kFALSE),
27 fCheckPlpFromDifferentBCMV(kFALSE),
30 fnSigmaPlpZdistSPD(3.),
31 fnSigmaPlpDiamXYSPD(2.),
32 fnSigmaPlpDiamZSPD(5.),
33 fUseSPDCutInMultBins(kFALSE)
38 //______________________________________________________________________
39 Bool_t AliAnalysisUtils::IsVertexSelected2013pA(AliVEvent *event)
41 Bool_t accept = kFALSE;
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);
54 AliFatal("Event is neither of AOD nor ESD type");
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){
67 TString vtxTtl = trkVtx->GetTitle();
68 if (!vtxTtl.Contains("VertexerTracks")) return accept;
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;
76 TString vtxTyp = spdVtx->GetTitle();
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;
83 if (TMath::Abs(zvtx) > fMaxVtxZ) return accept;
88 //______________________________________________________________________
89 Bool_t AliAnalysisUtils::IsFirstEventInChunk(AliVEvent *event)
92 Bool_t accept = kFALSE;
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);
105 AliFatal("Event is neither of AOD nor ESD type");
111 AliAODHeader *aodheader = 0x0;
112 aodheader = aod->GetHeader();
114 AliFatal("AOD header not there ?!");
117 if(aodheader->GetEventNumberESDFile()==0) accept = kTRUE;
119 if(esd->GetEventNumberInFile()==0) accept = kTRUE;
125 //______________________________________________________________________
126 Bool_t AliAnalysisUtils::IsPileUpEvent(AliVEvent *event)
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;
138 //______________________________________________________________________
139 Bool_t AliAnalysisUtils::IsPileUpMV(AliVEvent *event)
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);
146 AliFatal("Event is neither of AOD nor ESD type");
150 const AliVVertex* vtPrm = 0;
151 const AliVVertex* vtPlp = 0;
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
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
164 Int_t bcPrim = vtPrm->GetBC();
166 for (Int_t ipl=0;ipl<nPlp;ipl++) {
167 vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
169 if (vtPlp->GetNContributors() < fMinPlpContribMV) continue;
170 if (vtPlp->GetChi2perNDF() > fMaxPlpChi2MV) continue;
171 if(fCheckPlpFromDifferentBCMV)
173 Int_t bcPlp = vtPlp->GetBC();
174 if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
177 Double_t wDst = GetWDist(vtPrm,vtPlp);
178 if (wDst<fMinWDistMV) continue;
180 return kTRUE; // pile-up: well separated vertices
187 //______________________________________________________________________
188 Bool_t AliAnalysisUtils::IsPileUpSPD(AliVEvent *event)
190 // check for SPD pile-up
191 const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
192 const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
195 AliFatal("Event is neither of AOD nor ESD type");
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);
203 //______________________________________________________________________
204 Bool_t AliAnalysisUtils::IsOutOfBunchPileUp(AliVEvent *event)
206 // check for SPD pile-up
207 const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
208 const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
211 AliFatal("Event is neither of AOD nor ESD type");
214 Int_t bc2 = (aod)?aod->GetHeader()->GetIRInt2ClosestInteractionMap():esd->GetHeader()->GetIRInt2ClosestInteractionMap();
218 Int_t bc1 = (aod)?aod->GetHeader()->GetIRInt1ClosestInteractionMap():esd->GetHeader()->GetIRInt1ClosestInteractionMap();
225 //______________________________________________________________________
226 Double_t AliAnalysisUtils::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
228 // calculate sqrt of weighted distance to other vertex
230 printf("One of vertices is not valid\n");
233 static TMatrixDSym vVb(3);
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.;
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;