ANALYSIS and ANALYSIScalib
[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 #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 }