Disable caching for async prefetching
[u/mrichter/AliRoot.git] / 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 }