]>
Commit | Line | Data |
---|---|---|
6a213b59 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //---------------------------------------------------------------------------- | |
17 | // Implementation of the heavy-flavour vertexing analysis class | |
18 | // Candidates are store as objects deriving from AliAODRecoDecay. | |
19 | // An example of usage can be found in the macro AliVertexingHFTest.C. | |
20 | // Can be used as a task of AliAnalysisManager by means of the interface | |
b056c5e3 | 21 | // classes AliAnalysisTaskSEVertexingHF or AliAnalysisTaskVertexingHF. |
6a213b59 | 22 | // |
23 | // Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita | |
b056c5e3 | 24 | // Contact: andrea.dainese@lnl.infn.it |
6a213b59 | 25 | //---------------------------------------------------------------------------- |
26 | #include <TTree.h> | |
b82f6d67 | 27 | #include <TFile.h> |
6a213b59 | 28 | #include <TDatabasePDG.h> |
29 | #include <TString.h> | |
5e6f195c | 30 | #include "AliESDEvent.h" |
6a213b59 | 31 | #include "AliVertexerTracks.h" |
b82f6d67 | 32 | #include "AliKFParticle.h" |
6a213b59 | 33 | #include "AliESDVertex.h" |
b056c5e3 | 34 | #include "AliAODEvent.h" |
6a213b59 | 35 | #include "AliAODRecoDecay.h" |
36 | #include "AliAODRecoDecayHF.h" | |
37 | #include "AliAODRecoDecayHF2Prong.h" | |
38 | #include "AliAODRecoDecayHF3Prong.h" | |
39 | #include "AliAODRecoDecayHF4Prong.h" | |
40 | #include "AliAnalysisVertexingHF.h" | |
41 | ||
42 | ClassImp(AliAnalysisVertexingHF) | |
43 | ||
44 | //---------------------------------------------------------------------------- | |
45 | AliAnalysisVertexingHF::AliAnalysisVertexingHF(): | |
b82f6d67 | 46 | fBzkG(0.), |
47 | fSecVtxWithKF(kFALSE), | |
b056c5e3 | 48 | fUseTRef(kFALSE), |
6a213b59 | 49 | fRecoPrimVtxSkippingTrks(kFALSE), |
50 | fRmTrksFromPrimVtx(kFALSE), | |
51 | fV1(0x0), | |
52 | fDebug(0), | |
b82f6d67 | 53 | fD0toKpi(kTRUE), |
54 | fJPSItoEle(kTRUE), | |
55 | f3Prong(kTRUE), | |
56 | f4Prong(kTRUE), | |
57 | fITSrefit(kFALSE), | |
58 | fBothSPD(kTRUE), | |
59 | fMinITSCls(5), | |
60 | fMinPtCut(0.), | |
61 | fMind0rphiCut(0.) | |
6a213b59 | 62 | { |
63 | // Default constructor | |
64 | ||
65 | SetD0toKpiCuts(); | |
66 | SetBtoJPSICuts(); | |
67 | SetDplusCuts(); | |
b82f6d67 | 68 | SetDsCuts(); |
69 | SetLcCuts(); | |
6a213b59 | 70 | } |
13977a79 | 71 | //-------------------------------------------------------------------------- |
72 | AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) : | |
73 | TNamed(source), | |
b82f6d67 | 74 | fBzkG(source.fBzkG), |
75 | fSecVtxWithKF(source.fSecVtxWithKF), | |
b056c5e3 | 76 | fUseTRef(source.fUseTRef), |
13977a79 | 77 | fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks), |
78 | fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx), | |
79 | fV1(source.fV1), | |
80 | fDebug(source.fDebug), | |
81 | fD0toKpi(source.fD0toKpi), | |
82 | fJPSItoEle(source.fJPSItoEle), | |
83 | f3Prong(source.f3Prong), | |
84 | f4Prong(source.f4Prong), | |
85 | fITSrefit(source.fITSrefit), | |
86 | fBothSPD(source.fBothSPD), | |
87 | fMinITSCls(source.fMinITSCls), | |
88 | fMinPtCut(source.fMinPtCut), | |
89 | fMind0rphiCut(source.fMind0rphiCut) | |
90 | { | |
91 | // | |
92 | // Copy constructor | |
93 | // | |
94 | for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i]; | |
95 | for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i]; | |
96 | for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i]; | |
b82f6d67 | 97 | for(Int_t i=0; i<1; i++) fDsCuts[i]=source.fDsCuts[i]; |
98 | for(Int_t i=0; i<1; i++) fLcCuts[i]=source.fLcCuts[i]; | |
13977a79 | 99 | } |
100 | //-------------------------------------------------------------------------- | |
101 | AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source) | |
102 | { | |
103 | // | |
104 | // assignment operator | |
105 | // | |
106 | if(&source == this) return *this; | |
b82f6d67 | 107 | fBzkG = source.fBzkG; |
108 | fSecVtxWithKF = source.fSecVtxWithKF; | |
b056c5e3 | 109 | fUseTRef = source.fUseTRef; |
13977a79 | 110 | fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks; |
111 | fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx; | |
112 | fV1 = source.fV1; | |
113 | fDebug = source.fDebug; | |
114 | fD0toKpi = source.fD0toKpi; | |
115 | fJPSItoEle = source.fJPSItoEle; | |
116 | f3Prong = source.f3Prong; | |
117 | f4Prong = source.f4Prong; | |
118 | fITSrefit = source.fITSrefit; | |
119 | fBothSPD = source.fBothSPD; | |
120 | fMinITSCls = source.fMinITSCls; | |
121 | fMinPtCut = source.fMinPtCut; | |
122 | fMind0rphiCut = source.fMind0rphiCut; | |
123 | ||
124 | for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i]; | |
125 | for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i]; | |
126 | for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i]; | |
b82f6d67 | 127 | for(Int_t i=0; i<1; i++) fDsCuts[i]=source.fDsCuts[i]; |
128 | for(Int_t i=0; i<1; i++) fLcCuts[i]=source.fLcCuts[i]; | |
13977a79 | 129 | |
130 | return *this; | |
131 | } | |
6a213b59 | 132 | //---------------------------------------------------------------------------- |
133 | AliAnalysisVertexingHF::~AliAnalysisVertexingHF() { | |
134 | // Destructor | |
13977a79 | 135 | if(fV1) { delete fV1; fV1=0; } |
6a213b59 | 136 | } |
137 | //---------------------------------------------------------------------------- | |
b056c5e3 | 138 | void AliAnalysisVertexingHF::FindCandidatesESDtoAOD(AliESDEvent *esd, |
139 | TClonesArray *aodVerticesHFTClArr, | |
140 | TClonesArray *aodD0toKpiTClArr, | |
141 | TClonesArray *aodJPSItoEleTClArr, | |
142 | TClonesArray *aodCharm3ProngTClArr, | |
143 | TClonesArray *aodCharm4ProngTClArr) | |
144 | { | |
145 | // Find heavy-flavour vertex candidates | |
146 | // Input: ESD | |
147 | // Output: AOD (additional branches added) | |
148 | ||
149 | if(!aodVerticesHFTClArr) { | |
150 | printf("ERROR: no aodVerticesHFTClArr"); | |
151 | return; | |
152 | } | |
153 | if(fD0toKpi && !aodD0toKpiTClArr) { | |
154 | printf("ERROR: no aodD0toKpiTClArr"); | |
155 | return; | |
156 | } | |
157 | if(fJPSItoEle && !aodJPSItoEleTClArr) { | |
158 | printf("ERROR: no aodJPSItoEleTClArr"); | |
159 | return; | |
160 | } | |
161 | if(f3Prong && !aodCharm3ProngTClArr) { | |
162 | printf("ERROR: no aodCharm3ProngTClArr"); | |
163 | return; | |
164 | } | |
165 | if(f4Prong && !aodCharm4ProngTClArr) { | |
166 | printf("ERROR: no aodCharm4ProngTClArr"); | |
167 | return; | |
168 | } | |
169 | ||
170 | // delete candidates from previous event and create references | |
171 | Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0; | |
172 | aodVerticesHFTClArr->Delete(); | |
173 | iVerticesHF = aodVerticesHFTClArr->GetEntriesFast(); | |
174 | TClonesArray &verticesHFRef = *aodVerticesHFTClArr; | |
175 | if(fD0toKpi) { | |
176 | aodD0toKpiTClArr->Delete(); | |
177 | iD0toKpi = aodD0toKpiTClArr->GetEntriesFast(); | |
178 | } | |
179 | if(fJPSItoEle) { | |
180 | aodJPSItoEleTClArr->Delete(); | |
181 | iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast(); | |
182 | } | |
183 | if(f3Prong) { | |
184 | aodCharm3ProngTClArr->Delete(); | |
185 | i3Prong = aodCharm3ProngTClArr->GetEntriesFast(); | |
186 | } | |
187 | if(f4Prong) { | |
188 | aodCharm4ProngTClArr->Delete(); | |
189 | i4Prong = aodCharm4ProngTClArr->GetEntriesFast(); | |
190 | } | |
191 | TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr; | |
192 | TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr; | |
193 | TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr; | |
194 | TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr; | |
195 | ||
196 | ||
197 | AliAODRecoDecayHF2Prong *io2Prong = 0; | |
198 | AliAODRecoDecayHF3Prong *io3Prong = 0; | |
199 | AliAODRecoDecayHF4Prong *io4Prong = 0; | |
200 | ||
201 | Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries; | |
202 | Int_t nTrksP=0,nTrksN=0; | |
203 | Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2; | |
204 | Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE; | |
205 | AliESDtrack *postrack1 = 0; | |
206 | AliESDtrack *postrack2 = 0; | |
207 | AliESDtrack *negtrack1 = 0; | |
208 | AliESDtrack *negtrack2 = 0; | |
209 | Double_t dcaMax = fD0toKpiCuts[1]; | |
210 | if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1]; | |
211 | if(dcaMax < fDplusCuts[11]) dcaMax=fDplusCuts[11]; | |
212 | if(fDebug) printf(" dca cut set to %f cm\n",dcaMax); | |
213 | ||
214 | //------- SINGLE EVENT ANALYSIS --------------------------------- | |
215 | // | |
216 | Int_t ev = (Int_t)esd->GetEventNumberInFile(); | |
217 | printf("--- Finding candidates in event %d\n",ev); | |
218 | ||
219 | ||
220 | // get Bz from ESD | |
221 | fBzkG = (Double_t)esd->GetMagneticField(); | |
222 | ||
223 | trkEntries = (Int_t)esd->GetNumberOfTracks(); | |
224 | printf(" Number of tracks: %d\n",trkEntries); | |
225 | ||
226 | if(trkEntries<2 || !esd->GetPrimaryVertex()) { | |
227 | return; | |
228 | } | |
229 | ||
230 | // retrieve primary vertex from the AliESDEvent | |
231 | AliESDVertex copy(*(esd->GetPrimaryVertex())); | |
232 | SetPrimaryVertex(©); | |
233 | ||
234 | // call function which applies sigle-track selection and | |
235 | // separates positives and negatives | |
236 | TObjArray trksP(trkEntries/2); | |
237 | TObjArray trksN(trkEntries/2); | |
238 | SelectTracks(esd,trksP,nTrksP,trksN,nTrksN); | |
239 | ||
240 | printf(" Pos. tracks: %d Neg. tracks: %d\n",nTrksP,nTrksN); | |
241 | ||
242 | TObjArray *twoTrackArray1 = new TObjArray(2); | |
243 | TObjArray *twoTrackArray2 = new TObjArray(2); | |
244 | TObjArray *threeTrackArray = new TObjArray(3); | |
245 | TObjArray *fourTrackArray = new TObjArray(4); | |
246 | ||
247 | // LOOP ON POSITIVE TRACKS | |
248 | for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) { | |
249 | if(fDebug && iTrkP1%1==0) printf(" Processing positive track number %d of %d\n",iTrkP1,nTrksP); | |
250 | // get track from tracks array | |
251 | postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1); | |
252 | ||
253 | // LOOP ON NEGATIVE TRACKS | |
254 | for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) { | |
255 | if(fDebug && iTrkN1%1==0) printf(" Processing negative track number %d of %d\n",iTrkN1,nTrksN); | |
256 | // get track from tracks array | |
257 | negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1); | |
258 | // DCA between the two tracks | |
259 | dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy); | |
260 | if(dcap1n1>dcaMax) { negtrack1=0; continue; } | |
261 | // Vertexing | |
262 | twoTrackArray1->AddAt(postrack1,0); | |
263 | twoTrackArray1->AddAt(negtrack1,1); | |
264 | AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1); | |
265 | if(!vertexp1n1) { | |
266 | twoTrackArray1->Clear(); | |
267 | negtrack1=0; | |
268 | continue; | |
269 | } | |
270 | if(fD0toKpi || fJPSItoEle) { | |
271 | io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI); | |
272 | if(okD0 || okJPSI) { | |
273 | if(fUseTRef) { | |
274 | AliAODVertex *v = new(verticesHFRef[iVerticesHF++]) | |
275 | AliAODVertex(*(io2Prong->GetOwnSecondaryVtx())); | |
276 | v->SetType(AliAODVertex::kUndef); // to be changed | |
277 | Double_t px[2]={io2Prong->PxProng(0),io2Prong->PxProng(1)}; | |
278 | Double_t py[2]={io2Prong->PyProng(0),io2Prong->PyProng(1)}; | |
279 | Double_t pz[2]={io2Prong->PzProng(0),io2Prong->PzProng(1)}; | |
280 | Double_t d0[2]={io2Prong->Getd0Prong(0),io2Prong->Getd0Prong(1)}; | |
281 | Double_t d0err[2]={io2Prong->Getd0errProng(0),io2Prong->Getd0errProng(1)}; | |
282 | UShort_t id[2]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID()}; | |
283 | if(okD0) { | |
284 | AliAODRecoDecayHF2Prong *rd=new(aodD0toKpiRef[iD0toKpi++]) | |
285 | AliAODRecoDecayHF2Prong(v,px,py,pz,d0,d0err,dcap1n1); | |
286 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io2Prong->GetOwnPrimaryVtx()); | |
287 | rd->SetProngIDs(2,id); | |
b0ebe389 | 288 | v->SetParent(rd); |
b056c5e3 | 289 | } |
290 | if(okJPSI) { | |
291 | AliAODRecoDecayHF2Prong *rd=new(aodJPSItoEleRef[iJPSItoEle++]) | |
292 | AliAODRecoDecayHF2Prong(v,px,py,pz,d0,d0err,dcap1n1); | |
293 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io2Prong->GetOwnPrimaryVtx()); | |
294 | rd->SetProngIDs(2,id); | |
b0ebe389 | 295 | if(!okD0) v->SetParent(rd); // do something better here... |
b056c5e3 | 296 | } |
297 | //printf("DCA: %f\n",rd->GetDCA()); | |
298 | } else { | |
299 | if(okD0) new(aodD0toKpiRef[iD0toKpi++]) AliAODRecoDecayHF2Prong(*io2Prong); | |
300 | if(okJPSI) new(aodJPSItoEleRef[iJPSItoEle++]) AliAODRecoDecayHF2Prong(*io2Prong); | |
301 | } | |
302 | } | |
303 | //delete io2Prong; | |
304 | io2Prong=NULL; | |
305 | } | |
306 | ||
307 | twoTrackArray1->Clear(); | |
308 | if(!f3Prong && !f4Prong) { | |
309 | negtrack1=0; | |
310 | delete vertexp1n1; | |
311 | continue; | |
312 | } | |
313 | ||
314 | ||
315 | // 2nd LOOP ON POSITIVE TRACKS | |
316 | for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) { | |
317 | // get track from tracks array | |
318 | postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2); | |
319 | dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy); | |
320 | if(dcap2n1>dcaMax) { postrack2=0; continue; } | |
321 | dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy); | |
322 | if(dcap1p2>dcaMax) { postrack2=0; continue; } | |
323 | ||
324 | // Vertexing | |
325 | twoTrackArray2->AddAt(postrack2,0); | |
326 | twoTrackArray2->AddAt(negtrack1,1); | |
327 | AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2); | |
328 | if(!vertexp2n1) { | |
329 | twoTrackArray2->Clear(); | |
330 | postrack2=0; | |
331 | continue; | |
332 | } | |
333 | if(f3Prong) { | |
334 | threeTrackArray->AddAt(postrack1,0); | |
335 | threeTrackArray->AddAt(negtrack1,1); | |
336 | threeTrackArray->AddAt(postrack2,2); | |
337 | io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong); | |
338 | if(ok3Prong) { | |
339 | if(fUseTRef) { | |
340 | AliAODVertex *v = new(verticesHFRef[iVerticesHF++]) | |
341 | AliAODVertex(*(io3Prong->GetOwnSecondaryVtx())); | |
342 | v->SetType(AliAODVertex::kUndef); | |
343 | Double_t px[3]={io3Prong->PxProng(0),io3Prong->PxProng(1),io3Prong->PxProng(2)}; | |
344 | Double_t py[3]={io3Prong->PyProng(0),io3Prong->PyProng(1),io3Prong->PyProng(2)}; | |
345 | Double_t pz[3]={io3Prong->PzProng(0),io3Prong->PzProng(1),io3Prong->PzProng(2)}; | |
346 | Double_t d0[3]={io3Prong->Getd0Prong(0),io3Prong->Getd0Prong(1),io3Prong->Getd0Prong(2)}; | |
347 | Double_t d0err[3]={io3Prong->Getd0errProng(0),io3Prong->Getd0errProng(1),io3Prong->Getd0errProng(2)}; | |
348 | Double_t dcas[3]={io3Prong->GetDCA(0),io3Prong->GetDCA(1),io3Prong->GetDCA(2)}; | |
349 | UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID()}; | |
350 | AliAODRecoDecayHF3Prong *rd=new(aodCharm3ProngRef[i3Prong++]) | |
351 | AliAODRecoDecayHF3Prong(v,px,py,pz,d0,d0err,dcas,io3Prong->GetSigmaVert(),io3Prong->GetDist12toPrim(),io3Prong->GetDist23toPrim(),io3Prong->GetCharge()); | |
352 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io3Prong->GetOwnPrimaryVtx()); | |
353 | rd->SetProngIDs(3,id); | |
b0ebe389 | 354 | v->SetParent(rd); |
b056c5e3 | 355 | } else { |
356 | new(aodD0toKpiRef[i3Prong++]) AliAODRecoDecayHF3Prong(*io3Prong); | |
357 | } | |
358 | } | |
359 | if(io3Prong) { /*delete io3Prong;*/ io3Prong=NULL; } | |
360 | } | |
361 | if(f4Prong) { | |
362 | // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong) | |
363 | for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) { | |
364 | // get track from tracks array | |
365 | negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2); | |
366 | dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy); | |
367 | if(dcap1n2>dcaMax) { negtrack2=0; continue; } | |
368 | // Vertexing | |
369 | fourTrackArray->AddAt(postrack1,0); | |
370 | fourTrackArray->AddAt(negtrack1,1); | |
371 | fourTrackArray->AddAt(postrack2,2); | |
372 | fourTrackArray->AddAt(negtrack2,3); | |
373 | io4Prong = Make4Prong(fourTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong); | |
374 | if(ok4Prong) { | |
375 | if(fUseTRef) { | |
376 | AliAODVertex *v = new(verticesHFRef[iVerticesHF++]) | |
377 | AliAODVertex(*(io4Prong->GetOwnSecondaryVtx())); | |
378 | v->SetType(AliAODVertex::kUndef); | |
379 | Double_t px[4]={io4Prong->PxProng(0),io4Prong->PxProng(1), | |
380 | io4Prong->PxProng(2),io4Prong->PxProng(3)}; | |
381 | Double_t py[4]={io4Prong->PyProng(0),io4Prong->PyProng(1), | |
382 | io4Prong->PyProng(2),io4Prong->PyProng(3)}; | |
383 | Double_t pz[4]={io4Prong->PzProng(0),io4Prong->PzProng(1), | |
384 | io4Prong->PzProng(2),io4Prong->PzProng(3)}; | |
385 | Double_t d0[4]={io4Prong->Getd0Prong(0),io4Prong->Getd0Prong(1), | |
386 | io4Prong->Getd0Prong(2),io4Prong->Getd0Prong(3)}; | |
387 | Double_t d0err[4]={io4Prong->Getd0errProng(0),io4Prong->Getd0errProng(1), | |
388 | io4Prong->Getd0errProng(2),io4Prong->Getd0errProng(3)}; | |
389 | Double_t dcas[6]; io4Prong->GetDCAs(dcas); | |
390 | UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()}; | |
391 | AliAODRecoDecayHF4Prong *rd=new(aodCharm4ProngRef[i4Prong++]) | |
392 | AliAODRecoDecayHF4Prong(v,px,py,pz,d0,d0err,dcas,io4Prong->GetDist12toPrim(),io4Prong->GetDist23toPrim(),io4Prong->GetDist14toPrim(),io4Prong->GetDist34toPrim(),io4Prong->GetCharge()); | |
393 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io4Prong->GetOwnPrimaryVtx()); | |
394 | rd->SetProngIDs(4,id); | |
b0ebe389 | 395 | v->SetParent(rd); |
b056c5e3 | 396 | } else { |
397 | new(aodD0toKpiRef[i4Prong++]) AliAODRecoDecayHF4Prong(*io4Prong); | |
398 | } | |
399 | } | |
400 | if(io4Prong) { /*delete io4Prong;*/ io4Prong=NULL; } | |
401 | fourTrackArray->Clear(); | |
402 | negtrack2 = 0; | |
403 | } // end loop on negative tracks | |
404 | } | |
405 | postrack2 = 0; | |
406 | delete vertexp2n1; | |
407 | } // end 2nd loop on positive tracks | |
408 | twoTrackArray2->Clear(); | |
409 | ||
410 | // 2nd LOOP ON NEGATIVE TRACKS | |
411 | for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) { | |
412 | // get track from tracks array | |
413 | negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2); | |
414 | dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy); | |
415 | if(dcap1n2>dcaMax) { negtrack2=0; continue; } | |
416 | dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy); | |
417 | if(dcan1n2>dcaMax) { negtrack2=0; continue; } | |
418 | ||
419 | // Vertexing | |
420 | twoTrackArray2->AddAt(postrack1,0); | |
421 | twoTrackArray2->AddAt(negtrack2,1); | |
422 | AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2); | |
423 | if(!vertexp1n2) { | |
424 | twoTrackArray2->Clear(); | |
425 | negtrack2=0; | |
426 | continue; | |
427 | } | |
428 | if(f3Prong) { | |
429 | threeTrackArray->AddAt(negtrack1,0); | |
430 | threeTrackArray->AddAt(postrack1,1); | |
431 | threeTrackArray->AddAt(negtrack2,2); | |
432 | io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong); | |
433 | if(ok3Prong) { | |
434 | if(fUseTRef) { | |
435 | AliAODVertex *v = new(verticesHFRef[iVerticesHF++]) | |
436 | AliAODVertex(*(io3Prong->GetOwnSecondaryVtx())); | |
437 | v->SetType(AliAODVertex::kUndef); | |
438 | Double_t px[3]={io3Prong->PxProng(0),io3Prong->PxProng(1),io3Prong->PxProng(2)}; | |
439 | Double_t py[3]={io3Prong->PyProng(0),io3Prong->PyProng(1),io3Prong->PyProng(2)}; | |
440 | Double_t pz[3]={io3Prong->PzProng(0),io3Prong->PzProng(1),io3Prong->PzProng(2)}; | |
441 | Double_t d0[3]={io3Prong->Getd0Prong(0),io3Prong->Getd0Prong(1),io3Prong->Getd0Prong(2)}; | |
442 | Double_t d0err[3]={io3Prong->Getd0errProng(0),io3Prong->Getd0errProng(1),io3Prong->Getd0errProng(2)}; | |
443 | Double_t dcas[3]={io3Prong->GetDCA(0),io3Prong->GetDCA(1),io3Prong->GetDCA(2)}; | |
444 | UShort_t id[3]={(UShort_t)negtrack1->GetID(),(UShort_t)postrack1->GetID(),(UShort_t)negtrack2->GetID()}; | |
445 | AliAODRecoDecayHF3Prong *rd=new(aodCharm3ProngRef[i3Prong++]) | |
446 | AliAODRecoDecayHF3Prong(v,px,py,pz,d0,d0err,dcas,io3Prong->GetSigmaVert(),io3Prong->GetDist12toPrim(),io3Prong->GetDist23toPrim(),io3Prong->GetCharge()); | |
447 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io3Prong->GetOwnPrimaryVtx()); | |
448 | rd->SetProngIDs(3,id); | |
b0ebe389 | 449 | v->SetParent(rd); |
b056c5e3 | 450 | } else { |
451 | new(aodD0toKpiRef[i3Prong++]) AliAODRecoDecayHF3Prong(*io3Prong); | |
452 | } | |
453 | } | |
454 | if(io3Prong) { /*delete io3Prong;*/ io3Prong=NULL; } | |
455 | } | |
456 | negtrack2 = 0; | |
457 | delete vertexp1n2; | |
458 | } // end 2nd loop on negative tracks | |
459 | twoTrackArray2->Clear(); | |
460 | ||
461 | negtrack1 = 0; | |
462 | delete vertexp1n1; | |
463 | } // end 1st loop on negative tracks | |
464 | ||
465 | postrack1 = 0; | |
466 | } // end 1st loop on positive tracks | |
467 | ||
468 | ||
469 | if(fD0toKpi) { | |
470 | printf(" D0->Kpi in event %d = %d;\n", | |
471 | (Int_t)esd->GetEventNumberInFile(), | |
472 | (Int_t)aodD0toKpiTClArr->GetEntriesFast()); | |
473 | } | |
474 | if(fJPSItoEle) { | |
475 | printf(" JPSI->ee in event %d = %d;\n", | |
476 | (Int_t)esd->GetEventNumberInFile(), | |
477 | (Int_t)aodJPSItoEleTClArr->GetEntriesFast()); | |
478 | } | |
479 | if(f3Prong) { | |
480 | printf(" Charm->3Prong in event %d = %d;\n", | |
481 | (Int_t)esd->GetEventNumberInFile(), | |
482 | (Int_t)aodCharm3ProngTClArr->GetEntriesFast()); | |
483 | } | |
484 | if(f4Prong) { | |
485 | printf(" Charm->4Prong in event %d = %d;\n", | |
486 | (Int_t)esd->GetEventNumberInFile(), | |
487 | (Int_t)aodCharm4ProngTClArr->GetEntriesFast()); | |
488 | } | |
489 | ||
490 | ||
491 | //printf("delete twoTr 1\n"); | |
492 | twoTrackArray1->Delete(); delete twoTrackArray1; | |
493 | //printf("delete twoTr 2\n"); | |
494 | twoTrackArray2->Delete(); delete twoTrackArray2; | |
495 | //printf("delete threeTr 1\n"); | |
496 | threeTrackArray->Clear(); | |
497 | threeTrackArray->Delete(); delete threeTrackArray; | |
498 | //printf("delete fourTr 1\n"); | |
499 | fourTrackArray->Delete(); delete fourTrackArray; | |
500 | ||
501 | //------- END SINGLE EVENT ANALYSIS -------------------------------- | |
502 | ||
503 | return; | |
504 | } | |
505 | //---------------------------------------------------------------------------- | |
260836a9 | 506 | void AliAnalysisVertexingHF::FindCandidates(AliESDEvent *esd,TTree *treeout[]) |
6a213b59 | 507 | { |
508 | // Find heavy-flavour vertex candidates | |
b056c5e3 | 509 | // |
510 | // DEPRECATED: use FindCandidatesESDtoAOD! | |
6a213b59 | 511 | |
b056c5e3 | 512 | fUseTRef=kFALSE; // cannot use TRefs outside AOD |
513 | ||
b82f6d67 | 514 | AliAODRecoDecayHF2Prong *io2Prong = new AliAODRecoDecayHF2Prong(); |
515 | AliAODRecoDecayHF3Prong *io3Prong = new AliAODRecoDecayHF3Prong(); | |
516 | AliAODRecoDecayHF4Prong *io4Prong = new AliAODRecoDecayHF4Prong(); | |
6a213b59 | 517 | Int_t itree=0; |
518 | Int_t itreeD0toKpi=-1,itreeJPSItoEle=-1,itree3Prong=-1,itree4Prong=-1; | |
519 | Int_t initEntriesD0toKpi=0,initEntriesJPSItoEle=0,initEntries3Prong=0,initEntries4Prong=0; | |
520 | if(fD0toKpi) { | |
521 | itreeD0toKpi=itree; | |
260836a9 | 522 | treeout[itree]->SetBranchAddress("D0toKpi",&io2Prong); |
6a213b59 | 523 | itree++; |
260836a9 | 524 | initEntriesD0toKpi = treeout[itreeD0toKpi]->GetEntries(); |
6a213b59 | 525 | } |
526 | if(fJPSItoEle) { | |
527 | itreeJPSItoEle=itree; | |
260836a9 | 528 | treeout[itree]->SetBranchAddress("JPSItoEle",&io2Prong); |
6a213b59 | 529 | itree++; |
260836a9 | 530 | initEntriesJPSItoEle = treeout[itreeJPSItoEle]->GetEntries(); |
6a213b59 | 531 | } |
532 | if(f3Prong) { | |
533 | itree3Prong=itree; | |
260836a9 | 534 | treeout[itree]->SetBranchAddress("Charmto3Prong",&io3Prong); |
6a213b59 | 535 | itree++; |
260836a9 | 536 | initEntries3Prong = treeout[itree3Prong]->GetEntries(); |
6a213b59 | 537 | } |
538 | if(f4Prong) { | |
539 | itree4Prong=itree; | |
260836a9 | 540 | treeout[itree]->SetBranchAddress("D0to4Prong",&io4Prong); |
6a213b59 | 541 | itree++; |
260836a9 | 542 | initEntries4Prong = treeout[itree4Prong]->GetEntries(); |
6a213b59 | 543 | } |
544 | delete io2Prong; io2Prong = NULL; | |
545 | delete io3Prong; io3Prong = NULL; | |
546 | delete io4Prong; io4Prong = NULL; | |
547 | ||
548 | Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries; | |
549 | Int_t nTrksP=0,nTrksN=0; | |
550 | Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2; | |
551 | Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE; | |
552 | AliESDtrack *postrack1 = 0; | |
553 | AliESDtrack *postrack2 = 0; | |
554 | AliESDtrack *negtrack1 = 0; | |
555 | AliESDtrack *negtrack2 = 0; | |
556 | Double_t dcaMax = fD0toKpiCuts[1]; | |
557 | if(dcaMax<fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1]; | |
558 | if(dcaMax<fDplusCuts[11]) dcaMax=fDplusCuts[11]; | |
559 | if(fDebug) printf(" dca cut set to %f cm\n",dcaMax); | |
560 | ||
6a213b59 | 561 | Int_t ev = (Int_t)esd->GetEventNumberInFile(); |
562 | printf("--- Finding candidates in event %d\n",ev); | |
563 | ||
b82f6d67 | 564 | fBzkG = (Double_t)esd->GetMagneticField(); |
6a213b59 | 565 | |
566 | trkEntries = (Int_t)esd->GetNumberOfTracks(); | |
567 | printf(" Number of tracks: %d\n",trkEntries); | |
568 | if(trkEntries<2) return; | |
569 | ||
b82f6d67 | 570 | // retrieve primary vertex from the AliESDEvent |
6a213b59 | 571 | if(!esd->GetPrimaryVertex()) { |
572 | printf(" No vertex in AliESD\n"); | |
573 | return; | |
574 | } | |
575 | AliESDVertex copy(*(esd->GetPrimaryVertex())); | |
576 | SetPrimaryVertex(©); | |
6a213b59 | 577 | |
578 | // call function which applies sigle-track selection and | |
579 | // separetes positives and negatives | |
b056c5e3 | 580 | TObjArray trksP(trkEntries/2); |
581 | TObjArray trksN(trkEntries/2); | |
6a213b59 | 582 | SelectTracks(esd,trksP,nTrksP, |
583 | trksN,nTrksN); | |
584 | ||
585 | printf(" Pos. tracks: %d Neg. tracks: %d\n",nTrksP,nTrksN); | |
586 | ||
587 | TObjArray *twoTrackArray1 = new TObjArray(2); | |
588 | TObjArray *twoTrackArray2 = new TObjArray(2); | |
589 | TObjArray *threeTrackArray = new TObjArray(3); | |
590 | TObjArray *fourTrackArray = new TObjArray(4); | |
591 | ||
592 | // LOOP ON POSITIVE TRACKS | |
593 | for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) { | |
594 | if(fDebug) if(iTrkP1%1==0) printf(" Processing positive track number %d of %d\n",iTrkP1,nTrksP); | |
595 | // get track from track array | |
596 | postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1); | |
597 | ||
598 | // LOOP ON NEGATIVE TRACKS | |
599 | for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) { | |
600 | if(fDebug) if(iTrkN1%1==0) printf(" Processing negative track number %d of %d\n",iTrkN1,nTrksN); | |
601 | // get track from tracks array | |
602 | negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1); | |
603 | // DCA between the two tracks | |
b82f6d67 | 604 | dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy); |
6a213b59 | 605 | if(dcap1n1>dcaMax) { negtrack1=0; continue; } |
b82f6d67 | 606 | // Vertexing |
6a213b59 | 607 | twoTrackArray1->AddAt(postrack1,0); |
608 | twoTrackArray1->AddAt(negtrack1,1); | |
b82f6d67 | 609 | AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1); |
b056c5e3 | 610 | if(!vertexp1n1) { |
806873d4 | 611 | twoTrackArray1->Clear(); |
806873d4 | 612 | negtrack1=0; |
613 | continue; | |
6a213b59 | 614 | } |
615 | if(fD0toKpi || fJPSItoEle) { | |
616 | io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI); | |
260836a9 | 617 | if(okD0) treeout[itreeD0toKpi]->Fill(); |
618 | if(okJPSI) treeout[itreeJPSItoEle]->Fill(); | |
6a213b59 | 619 | delete io2Prong; io2Prong=NULL; |
620 | } | |
621 | ||
622 | twoTrackArray1->Clear(); | |
623 | if(!f3Prong && !f4Prong) { | |
624 | negtrack1=0; | |
625 | delete vertexp1n1; | |
626 | continue; | |
627 | } | |
628 | ||
629 | // 2nd LOOP ON POSITIVE TRACKS | |
630 | for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) { | |
631 | // get track from tracks array | |
632 | postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2); | |
b82f6d67 | 633 | dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy); |
6a213b59 | 634 | if(dcap2n1>dcaMax) { postrack2=0; continue; } |
b82f6d67 | 635 | dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy); |
6a213b59 | 636 | if(dcap1p2>dcaMax) { postrack2=0; continue; } |
637 | ||
b82f6d67 | 638 | // Vertexing |
6a213b59 | 639 | twoTrackArray2->AddAt(postrack2,0); |
640 | twoTrackArray2->AddAt(negtrack1,1); | |
b82f6d67 | 641 | AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2); |
b056c5e3 | 642 | if(!vertexp2n1) { |
643 | twoTrackArray2->Clear(); | |
644 | postrack2=0; | |
645 | continue; | |
646 | } | |
6a213b59 | 647 | if(f3Prong) { |
648 | threeTrackArray->AddAt(postrack1,0); | |
649 | threeTrackArray->AddAt(negtrack1,1); | |
650 | threeTrackArray->AddAt(postrack2,2); | |
651 | io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong); | |
260836a9 | 652 | if(ok3Prong) treeout[itree3Prong]->Fill(); |
6a213b59 | 653 | if(io3Prong) delete io3Prong; io3Prong=NULL; |
654 | } | |
655 | if(f4Prong) { | |
656 | // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong) | |
657 | for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) { | |
658 | // get track from tracks array | |
659 | negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2); | |
b82f6d67 | 660 | dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy); |
6a213b59 | 661 | if(dcap1n2>dcaMax) { negtrack2=0; continue; } |
b82f6d67 | 662 | // Vertexing |
6a213b59 | 663 | fourTrackArray->AddAt(postrack1,0); |
664 | fourTrackArray->AddAt(negtrack1,1); | |
665 | fourTrackArray->AddAt(postrack2,2); | |
666 | fourTrackArray->AddAt(negtrack2,3); | |
667 | io4Prong = Make4Prong(fourTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong); | |
260836a9 | 668 | if(ok4Prong) treeout[itree4Prong]->Fill(); |
6a213b59 | 669 | delete io4Prong; io4Prong=NULL; |
670 | fourTrackArray->Clear(); | |
671 | negtrack2 = 0; | |
672 | } // end loop on negative tracks | |
673 | } | |
674 | postrack2 = 0; | |
675 | delete vertexp2n1; | |
676 | } // end 2nd loop on positive tracks | |
677 | twoTrackArray2->Clear(); | |
678 | ||
679 | // 2nd LOOP ON NEGATIVE TRACKS | |
680 | for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) { | |
681 | // get track from tracks array | |
682 | negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2); | |
b82f6d67 | 683 | dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy); |
6a213b59 | 684 | if(dcap1n2>dcaMax) { negtrack2=0; continue; } |
b82f6d67 | 685 | dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy); |
6a213b59 | 686 | if(dcan1n2>dcaMax) { negtrack2=0; continue; } |
687 | ||
b82f6d67 | 688 | // Vertexing |
6a213b59 | 689 | twoTrackArray2->AddAt(postrack1,0); |
690 | twoTrackArray2->AddAt(negtrack2,1); | |
b82f6d67 | 691 | AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2); |
b056c5e3 | 692 | if(!vertexp1n2) { |
693 | twoTrackArray2->Clear(); | |
694 | negtrack2=0; | |
695 | continue; | |
696 | } | |
6a213b59 | 697 | if(f3Prong) { |
698 | threeTrackArray->AddAt(negtrack1,0); | |
699 | threeTrackArray->AddAt(postrack1,1); | |
700 | threeTrackArray->AddAt(negtrack2,2); | |
701 | io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong); | |
260836a9 | 702 | if(ok3Prong) treeout[itree3Prong]->Fill(); |
6a213b59 | 703 | if(io3Prong) delete io3Prong; io3Prong=NULL; |
704 | } | |
705 | negtrack2 = 0; | |
706 | delete vertexp1n2; | |
707 | } // end 2nd loop on negative tracks | |
708 | twoTrackArray2->Clear(); | |
709 | ||
710 | negtrack1 = 0; | |
711 | delete vertexp1n1; | |
712 | } // end 1st loop on negative tracks | |
713 | ||
714 | postrack1 = 0; | |
715 | } // end 1st loop on positive tracks | |
716 | ||
717 | ||
718 | ||
6a213b59 | 719 | //printf("delete twoTr 1\n"); |
720 | twoTrackArray1->Delete(); delete twoTrackArray1; | |
721 | //printf("delete twoTr 2\n"); | |
722 | twoTrackArray2->Delete(); delete twoTrackArray2; | |
723 | //printf("delete threeTr 1\n"); | |
724 | threeTrackArray->Clear(); | |
725 | threeTrackArray->Delete(); delete threeTrackArray; | |
726 | //printf("delete fourTr 1\n"); | |
727 | fourTrackArray->Delete(); delete fourTrackArray; | |
728 | ||
729 | ||
730 | // create a copy of this class to be written to output file | |
731 | //AliAnalysisVertexingHF *copy = (AliAnalysisVertexingHF*)this->Clone("AnalysisVertexingHF"); | |
732 | ||
733 | // print statistics | |
734 | if(fD0toKpi) { | |
735 | printf(" D0->Kpi: event %d = %d; total = %d;\n", | |
736 | (Int_t)esd->GetEventNumberInFile(), | |
260836a9 | 737 | (Int_t)treeout[itreeD0toKpi]->GetEntries()-initEntriesD0toKpi, |
738 | (Int_t)treeout[itreeD0toKpi]->GetEntries()); | |
6a213b59 | 739 | } |
740 | if(fJPSItoEle) { | |
741 | printf(" JPSI->ee: event %d = %d; total = %d;\n", | |
742 | (Int_t)esd->GetEventNumberInFile(), | |
260836a9 | 743 | (Int_t)treeout[itreeJPSItoEle]->GetEntries()-initEntriesJPSItoEle, |
744 | (Int_t)treeout[itreeJPSItoEle]->GetEntries()); | |
6a213b59 | 745 | } |
746 | if(f3Prong) { | |
747 | printf(" Charm->3Prong: event %d = %d; total = %d;\n", | |
748 | (Int_t)esd->GetEventNumberInFile(), | |
260836a9 | 749 | (Int_t)treeout[itree3Prong]->GetEntries()-initEntries3Prong, |
750 | (Int_t)treeout[itree3Prong]->GetEntries()); | |
6a213b59 | 751 | } |
752 | if(f4Prong) { | |
753 | printf(" Charm->4Prong: event %d = %d; total = %d;\n", | |
754 | (Int_t)esd->GetEventNumberInFile(), | |
260836a9 | 755 | (Int_t)treeout[itree4Prong]->GetEntries()-initEntries4Prong, |
756 | (Int_t)treeout[itree4Prong]->GetEntries()); | |
6a213b59 | 757 | } |
758 | ||
759 | ||
760 | return; | |
761 | } | |
762 | //---------------------------------------------------------------------------- | |
763 | AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong( | |
5e6f195c | 764 | TObjArray *twoTrackArray1,AliESDEvent *esd, |
6a213b59 | 765 | AliESDVertex *secVertexESD,Double_t dca, |
766 | Bool_t &okD0,Bool_t &okJPSI) const | |
767 | { | |
768 | // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI | |
769 | // reconstruction cuts | |
770 | // G.E.Bruno (J/psi), A.Dainese (D0->Kpi) | |
771 | ||
772 | okD0=kFALSE; okJPSI=kFALSE; | |
773 | ||
774 | Double_t px[2],py[2],pz[2],d0[2],d0err[2]; | |
6a213b59 | 775 | |
776 | AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0); | |
777 | AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1); | |
778 | ||
779 | // propagate tracks to secondary vertex, to compute inv. mass | |
b82f6d67 | 780 | postrack->RelateToVertex(secVertexESD,fBzkG,10.); |
781 | negtrack->RelateToVertex(secVertexESD,fBzkG,10.); | |
6a213b59 | 782 | |
783 | Double_t momentum[3]; | |
784 | postrack->GetPxPyPz(momentum); | |
785 | px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; | |
786 | negtrack->GetPxPyPz(momentum); | |
787 | px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; | |
788 | ||
789 | ||
790 | // invariant mass cut (try to improve coding here..) | |
791 | Bool_t okMassCut=kFALSE; | |
792 | if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE; | |
793 | if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE; | |
794 | if(!okMassCut) { | |
795 | if(fDebug) printf(" candidate didn't pass mass cut\n"); | |
796 | return 0x0; | |
797 | } | |
798 | ||
799 | ||
800 | AliESDVertex *primVertex = fV1; | |
801 | AliESDVertex *ownPrimVertex=0; | |
802 | ||
803 | // primary vertex from *other* tracks in the event | |
804 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) { | |
805 | ownPrimVertex = OwnPrimaryVertex(2,twoTrackArray1,esd); | |
806 | if(!ownPrimVertex) { | |
807 | return 0x0; | |
808 | } else { | |
809 | if(ownPrimVertex->GetNContributors()<2) { | |
810 | delete ownPrimVertex; | |
811 | return 0x0; | |
812 | } else { | |
813 | primVertex = ownPrimVertex; | |
814 | } | |
815 | } | |
816 | } | |
817 | ||
818 | Float_t d0z0[2],covd0z0[3]; | |
b82f6d67 | 819 | postrack->RelateToVertex(primVertex,fBzkG,10.); |
6a213b59 | 820 | postrack->GetImpactParameters(d0z0,covd0z0); |
821 | d0[0] = d0z0[0]; | |
822 | d0err[0] = TMath::Sqrt(covd0z0[0]); | |
b82f6d67 | 823 | negtrack->RelateToVertex(primVertex,fBzkG,10.); |
6a213b59 | 824 | negtrack->GetImpactParameters(d0z0,covd0z0); |
825 | d0[1] = d0z0[0]; | |
826 | d0err[1] = TMath::Sqrt(covd0z0[0]); | |
827 | ||
828 | // create the object AliAODRecoDecayHF2Prong | |
829 | Double_t pos[3],cov[6]; | |
830 | secVertexESD->GetXYZ(pos); // position | |
831 | secVertexESD->GetCovMatrix(cov); //covariance matrix | |
832 | AliAODVertex *secVertexAOD = new AliAODVertex(pos,cov,secVertexESD->GetChi2toNDF()); | |
b056c5e3 | 833 | AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVertexAOD,px,py,pz,d0,d0err,dca); |
834 | the2Prong->SetOwnSecondaryVtx(secVertexAOD); | |
6a213b59 | 835 | primVertex->GetXYZ(pos); // position |
836 | primVertex->GetCovMatrix(cov); //covariance matrix | |
837 | AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF()); | |
6a213b59 | 838 | the2Prong->SetOwnPrimaryVtx(primVertexAOD); |
839 | ||
b056c5e3 | 840 | |
6a213b59 | 841 | // select D0->Kpi |
842 | Int_t checkD0,checkD0bar; | |
843 | if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar); | |
844 | //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0); | |
845 | // select J/psi from B | |
846 | Int_t checkJPSI; | |
847 | if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI); | |
848 | //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI); | |
849 | ||
850 | ||
851 | if(okD0 || okJPSI) { | |
852 | // get PID info from ESD | |
853 | Double_t esdpid0[5]; | |
854 | postrack->GetESDpid(esdpid0); | |
855 | Double_t esdpid1[5]; | |
856 | negtrack->GetESDpid(esdpid1); | |
857 | Double_t esdpid[10]; | |
858 | for(Int_t i=0;i<5;i++) { | |
859 | esdpid[i] = esdpid0[i]; | |
860 | esdpid[5+i] = esdpid1[i]; | |
861 | } | |
862 | the2Prong->SetPID(2,esdpid); | |
863 | } | |
864 | ||
865 | if(ownPrimVertex) delete ownPrimVertex; | |
b82f6d67 | 866 | |
6a213b59 | 867 | return the2Prong; |
868 | } | |
869 | //---------------------------------------------------------------------------- | |
870 | AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong( | |
5e6f195c | 871 | TObjArray *threeTrackArray,AliESDEvent *esd, |
6a213b59 | 872 | AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1, |
873 | Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2, | |
874 | Bool_t &ok3Prong) const | |
875 | { | |
876 | // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac | |
877 | // reconstruction cuts | |
878 | // E.Bruna, F.Prino | |
879 | ||
880 | ok3Prong=kFALSE; | |
881 | Double_t px[3],py[3],pz[3],d0[3],d0err[3];//d0z[3]; | |
882 | Float_t d0z0[2],covd0z0[3]; | |
883 | ||
6a213b59 | 884 | |
885 | AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0); | |
886 | AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1); | |
887 | AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2); | |
888 | ||
889 | AliESDVertex *primVertex = fV1; | |
890 | ||
b82f6d67 | 891 | postrack1->RelateToVertex(primVertex,fBzkG,10.); |
892 | negtrack->RelateToVertex(primVertex,fBzkG,10.); | |
893 | postrack2->RelateToVertex(primVertex,fBzkG,10.); | |
6a213b59 | 894 | |
895 | Double_t momentum[3]; | |
896 | postrack1->GetPxPyPz(momentum); | |
897 | px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; | |
898 | negtrack->GetPxPyPz(momentum); | |
899 | px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; | |
900 | postrack2->GetPxPyPz(momentum); | |
901 | px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2]; | |
902 | ||
903 | postrack1->GetImpactParameters(d0z0,covd0z0); | |
904 | d0[0]=d0z0[0]; | |
905 | d0err[0] = TMath::Sqrt(covd0z0[0]); | |
906 | negtrack->GetImpactParameters(d0z0,covd0z0); | |
907 | d0[1]=d0z0[0]; | |
908 | d0err[1] = TMath::Sqrt(covd0z0[0]); | |
909 | postrack2->GetImpactParameters(d0z0,covd0z0); | |
910 | d0[2]=d0z0[0]; | |
911 | d0err[2] = TMath::Sqrt(covd0z0[0]); | |
912 | ||
913 | ||
b82f6d67 | 914 | // invariant mass cut for D+, Ds, Lc |
6a213b59 | 915 | Bool_t okMassCut=kFALSE; |
916 | if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE; | |
917 | if(!okMassCut) { | |
918 | if(fDebug) printf(" candidate didn't pass mass cut\n"); | |
919 | return 0x0; | |
920 | } | |
921 | ||
922 | //charge | |
923 | Short_t charge=(Short_t)(postrack1->GetSign()*postrack2->GetSign()*negtrack->GetSign()); | |
924 | ||
925 | AliESDVertex *ownPrimVertex = 0; | |
926 | // primary vertex from *other* tracks in the event | |
927 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) { | |
928 | ownPrimVertex = OwnPrimaryVertex(3,threeTrackArray,esd); | |
929 | if(!ownPrimVertex) { | |
930 | return 0x0; | |
931 | } else { | |
932 | if(ownPrimVertex->GetNContributors()<2) { | |
933 | delete ownPrimVertex; | |
934 | return 0x0; | |
935 | } else { | |
936 | primVertex = ownPrimVertex; | |
937 | } | |
938 | } | |
939 | } | |
940 | ||
941 | // create the object AliAODRecoDecayHF3Prong | |
b82f6d67 | 942 | AliESDVertex* secVert3Prong = ReconstructSecondaryVertex(threeTrackArray); |
b056c5e3 | 943 | if(!secVert3Prong) { |
944 | if(ownPrimVertex) delete ownPrimVertex; | |
945 | return 0x0; | |
946 | } | |
6a213b59 | 947 | Double_t pos[3],cov[6],sigmavert; |
948 | secVert3Prong->GetXYZ(pos); // position | |
949 | secVert3Prong->GetCovMatrix(cov); //covariance matrix | |
950 | sigmavert=secVert3Prong->GetDispersion(); | |
951 | ||
952 | AliAODVertex *secVert3PrAOD = new AliAODVertex(pos,cov,secVert3Prong->GetChi2toNDF()); | |
953 | primVertex->GetXYZ(pos); // position | |
954 | primVertex->GetCovMatrix(cov); //covariance matrix | |
955 | AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF()); | |
956 | Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2}; | |
957 | ||
958 | Double_t dist12=TMath::Sqrt((vertexp1n1->GetXv()-pos[0])*(vertexp1n1->GetXv()-pos[0])+(vertexp1n1->GetYv()-pos[1])*(vertexp1n1->GetYv()-pos[1])+(vertexp1n1->GetZv()-pos[2])*(vertexp1n1->GetZv()-pos[2])); | |
959 | Double_t dist23=TMath::Sqrt((vertexp2n1->GetXv()-pos[0])*(vertexp2n1->GetXv()-pos[0])+(vertexp2n1->GetYv()-pos[1])*(vertexp2n1->GetYv()-pos[1])+(vertexp2n1->GetZv()-pos[2])*(vertexp2n1->GetZv()-pos[2])); | |
960 | ||
961 | AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert3PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge); | |
b056c5e3 | 962 | the3Prong->SetOwnSecondaryVtx(secVert3PrAOD); |
6a213b59 | 963 | the3Prong->SetOwnPrimaryVtx(primVertexAOD); |
964 | ||
965 | ||
b82f6d67 | 966 | // select D+->Kpipi, Ds->KKpi, Lc->pKpi |
967 | if(f3Prong) { | |
968 | ok3Prong = kFALSE; | |
969 | Int_t ok1,ok2; | |
970 | if(the3Prong->SelectDplus(fDplusCuts)) ok3Prong = kTRUE; | |
971 | if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE; | |
972 | if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE; | |
973 | } | |
6a213b59 | 974 | //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong); |
975 | if(ok3Prong) { | |
976 | // get PID info from ESD | |
977 | Double_t esdpid0[5]; | |
978 | postrack1->GetESDpid(esdpid0); | |
979 | Double_t esdpid1[5]; | |
980 | negtrack->GetESDpid(esdpid1); | |
981 | Double_t esdpid2[5]; | |
982 | postrack2->GetESDpid(esdpid2); | |
983 | ||
984 | ||
985 | Double_t esdpid[15]; | |
986 | for(Int_t i=0;i<5;i++) { | |
987 | esdpid[i] = esdpid0[i]; | |
988 | esdpid[5+i] = esdpid1[i]; | |
989 | esdpid[10+i] = esdpid2[i]; | |
990 | } | |
991 | the3Prong->SetPID(3,esdpid); | |
992 | } | |
993 | ||
994 | if(ownPrimVertex) delete ownPrimVertex; | |
995 | ||
996 | return the3Prong; | |
997 | } | |
998 | //---------------------------------------------------------------------------- | |
999 | AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong( | |
5e6f195c | 1000 | TObjArray *fourTrackArray,AliESDEvent *esd, |
6a213b59 | 1001 | AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1, |
1002 | Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1, | |
1003 | Bool_t &ok4Prong) const | |
1004 | { | |
1005 | // Make 4Prong candidates and check if they pass D0toKpipipi | |
1006 | // reconstruction cuts | |
1007 | // G.E.Bruno, R.Romita | |
1008 | ||
1009 | ok4Prong=kFALSE; | |
1010 | ||
1011 | Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3]; | |
1012 | //Float_t d0z0[2],covd0z0[3]; | |
1013 | ||
b82f6d67 | 1014 | px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used) |
6a213b59 | 1015 | |
1016 | //charge | |
1017 | Short_t charge=0; | |
1018 | ||
1019 | AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0); | |
1020 | AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1); | |
1021 | AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2); | |
1022 | AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3); | |
1023 | ||
1024 | AliESDVertex *primVertex = fV1; | |
1025 | ||
b82f6d67 | 1026 | postrack1->RelateToVertex(primVertex,fBzkG,10.); |
1027 | negtrack1->RelateToVertex(primVertex,fBzkG,10.); | |
1028 | postrack2->RelateToVertex(primVertex,fBzkG,10.); | |
1029 | negtrack2->RelateToVertex(primVertex,fBzkG,10.); | |
6a213b59 | 1030 | |
1031 | Double_t momentum[3]; | |
1032 | postrack1->GetPxPyPz(momentum); | |
1033 | px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; | |
1034 | negtrack1->GetPxPyPz(momentum); | |
1035 | px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; | |
1036 | postrack2->GetPxPyPz(momentum); | |
1037 | px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2]; | |
1038 | negtrack2->GetPxPyPz(momentum); | |
1039 | px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2]; | |
1040 | ||
1041 | // invariant mass cut for rho or D0 (try to improve coding here..) | |
1042 | //Bool_t okMassCut=kFALSE; | |
1043 | //if(!okMassCut) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE; | |
1044 | //if(!okMassCut) { | |
1045 | // if(fDebug) printf(" candidate didn't pass mass cut\n"); | |
1046 | // return 0x0; | |
1047 | //} | |
1048 | ||
1049 | AliESDVertex *ownPrimVertex = 0; | |
1050 | // primary vertex from *other* tracks in the event | |
1051 | if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) { | |
1052 | ownPrimVertex = OwnPrimaryVertex(4,fourTrackArray,esd); | |
1053 | if(!ownPrimVertex) { | |
1054 | return 0x0; | |
1055 | } else { | |
1056 | if(ownPrimVertex->GetNContributors()<2) { | |
1057 | delete ownPrimVertex; | |
1058 | return 0x0; | |
1059 | } else { | |
1060 | primVertex = ownPrimVertex; | |
1061 | } | |
1062 | } | |
1063 | } | |
1064 | ||
1065 | // create the object AliAODRecoDecayHF4Prong | |
b82f6d67 | 1066 | AliESDVertex* secVert4Prong = ReconstructSecondaryVertex(fourTrackArray); |
b056c5e3 | 1067 | if(!secVert4Prong) { |
1068 | if(ownPrimVertex) delete ownPrimVertex; | |
1069 | return 0x0; | |
1070 | } | |
6a213b59 | 1071 | Double_t pos[3],cov[6],sigmavert; |
1072 | secVert4Prong->GetXYZ(pos); // position | |
1073 | secVert4Prong->GetCovMatrix(cov); //covariance matrix | |
1074 | sigmavert=secVert4Prong->GetDispersion(); | |
1075 | ||
1076 | AliAODVertex *secVert4PrAOD = new AliAODVertex(pos,cov,secVert4Prong->GetChi2toNDF()); | |
1077 | primVertex->GetXYZ(pos); // position | |
1078 | primVertex->GetCovMatrix(cov); //covariance matrix | |
1079 | AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF()); | |
1080 | //Double_t dca[6]={dcap1n1,dcap2n1,dcap1p2,0.,0.,0.}; // | |
1081 | Double_t dca[6]={0.,0.,0.,0.,0.,0.}; // modify it | |
1082 | ||
1083 | Double_t dist12=TMath::Sqrt((vertexp1n1->GetXv()-pos[0])*(vertexp1n1->GetXv()-pos[0])+(vertexp1n1->GetYv()-pos[1])*(vertexp1n1->GetYv()-pos[1])+(vertexp1n1->GetZv()-pos[2])*(vertexp1n1->GetZv()-pos[2])); | |
1084 | Double_t dist23=TMath::Sqrt((vertexp2n1->GetXv()-pos[0])*(vertexp2n1->GetXv()-pos[0])+(vertexp2n1->GetYv()-pos[1])*(vertexp2n1->GetYv()-pos[1])+(vertexp2n1->GetZv()-pos[2])*(vertexp2n1->GetZv()-pos[2])); | |
1085 | Double_t dist14=0.; // to be implemented | |
1086 | Double_t dist34=0.; // to be implemented | |
1087 | ||
1088 | //AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge); | |
1089 | AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge); | |
1090 | the4Prong->SetOwnPrimaryVtx(primVertexAOD); | |
b056c5e3 | 1091 | the4Prong->SetOwnSecondaryVtx(secVert4PrAOD); |
6a213b59 | 1092 | |
1093 | ||
1094 | // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available | |
1095 | // select D0->Kpipipi | |
1096 | //Int_t checkD0,checkD0bar; | |
1097 | // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar); | |
1098 | ok4Prong=kFALSE; //for the time being ... | |
1099 | ||
1100 | ||
1101 | // get PID info from ESD | |
1102 | Double_t esdpid0[5]; | |
1103 | postrack1->GetESDpid(esdpid0); | |
1104 | Double_t esdpid1[5]; | |
1105 | negtrack1->GetESDpid(esdpid1); | |
1106 | Double_t esdpid2[5]; | |
1107 | postrack2->GetESDpid(esdpid2); | |
1108 | Double_t esdpid3[5]; | |
1109 | negtrack2->GetESDpid(esdpid3); | |
1110 | ||
1111 | Double_t esdpid[20]; | |
1112 | for(Int_t i=0;i<5;i++) { | |
1113 | esdpid[i] = esdpid0[i]; | |
1114 | esdpid[5+i] = esdpid1[i]; | |
1115 | esdpid[10+i] = esdpid2[i]; | |
1116 | esdpid[15+i] = esdpid3[i]; | |
1117 | } | |
1118 | the4Prong->SetPID(4,esdpid); | |
1119 | ||
1120 | if(ownPrimVertex) delete ownPrimVertex; | |
1121 | ||
1122 | return the4Prong; | |
1123 | } | |
1124 | //----------------------------------------------------------------------------- | |
1125 | AliESDVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(Int_t ntrks, | |
1126 | TObjArray *trkArray, | |
5e6f195c | 1127 | AliESDEvent *esd) const |
6a213b59 | 1128 | { |
1129 | // Returns primary vertex specific to this candidate | |
1130 | ||
1131 | AliVertexerTracks *vertexer1 = new AliVertexerTracks(esd->GetMagneticField()); | |
1132 | AliESDVertex *ownPrimVertex = 0; | |
1133 | ||
1134 | // recalculating the vertex | |
1135 | if(fRecoPrimVtxSkippingTrks) { | |
1136 | if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) { | |
1137 | Float_t diamondcovxy[3]; | |
1138 | esd->GetDiamondCovXY(diamondcovxy); | |
1139 | Double_t pos[3]={esd->GetDiamondX(),esd->GetDiamondY(),0.}; | |
1140 | Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.}; | |
1141 | AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1); | |
1142 | vertexer1->SetVtxStart(diamond); | |
1143 | delete diamond; diamond=NULL; | |
1144 | if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) | |
1145 | vertexer1->SetOnlyFitter(); | |
1146 | } | |
d70e06c7 | 1147 | Int_t skipped[10]; |
6a213b59 | 1148 | AliESDtrack *t = 0; |
1149 | for(Int_t i=0; i<ntrks; i++) { | |
1150 | t = (AliESDtrack*)trkArray->UncheckedAt(i); | |
d70e06c7 | 1151 | skipped[i] = (Int_t)t->GetID(); |
6a213b59 | 1152 | } |
1153 | vertexer1->SetSkipTracks(ntrks,skipped); | |
1154 | ownPrimVertex = (AliESDVertex*)vertexer1->FindPrimaryVertex(esd); | |
1155 | } | |
1156 | ||
1157 | // removing the prongs tracks | |
1158 | if(fRmTrksFromPrimVtx) { | |
f65f93f3 | 1159 | TObjArray rmArray(ntrks); |
1160 | UShort_t *rmId = new UShort_t[ntrks]; | |
6a213b59 | 1161 | AliESDtrack *esdTrack = 0; |
6a213b59 | 1162 | AliESDtrack *t = 0; |
1163 | for(Int_t i=0; i<ntrks; i++) { | |
1164 | t = (AliESDtrack*)trkArray->UncheckedAt(i); | |
1165 | esdTrack = new AliESDtrack(*t); | |
f65f93f3 | 1166 | rmArray.AddLast(esdTrack); |
1167 | rmId[i]=(UShort_t)esdTrack->GetID(); | |
6a213b59 | 1168 | delete esdTrack; |
1169 | } | |
1170 | Float_t diamondxy[2]={esd->GetDiamondX(),esd->GetDiamondY()}; | |
f65f93f3 | 1171 | ownPrimVertex = vertexer1->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy); |
1172 | delete [] rmId; rmId=NULL; | |
1173 | rmArray.Delete(); | |
6a213b59 | 1174 | } |
1175 | ||
1176 | delete vertexer1; vertexer1=NULL; | |
1177 | ||
1178 | return ownPrimVertex; | |
1179 | } | |
1180 | //----------------------------------------------------------------------------- | |
1181 | void AliAnalysisVertexingHF::PrintStatus() const { | |
1182 | // Print parameters being used | |
1183 | ||
1184 | printf("Preselections:\n"); | |
1185 | printf(" fITSrefit = %d\n",(Int_t)fITSrefit); | |
1186 | printf(" fBothSPD = %d\n",(Int_t)fBothSPD); | |
1187 | printf(" fMinITSCls = %d\n",fMinITSCls); | |
1188 | printf(" fMinPtCut = %f GeV/c\n",fMinPtCut); | |
1189 | printf(" fMind0rphiCut = %f cm\n",fMind0rphiCut); | |
4d5c9633 | 1190 | if(fSecVtxWithKF) { |
1191 | printf("Secondary vertex with Kalman filter package (AliKFParticle)\n"); | |
1192 | } else { | |
1193 | printf("Secondary vertex with AliVertexerTracks\n"); | |
1194 | } | |
6a213b59 | 1195 | if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n"); |
1196 | if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n"); | |
1197 | if(fD0toKpi) { | |
1198 | printf("Reconstruct D0->Kpi candidates with cuts:\n"); | |
1199 | printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]); | |
1200 | printf(" dca [cm] < %f\n",fD0toKpiCuts[1]); | |
1201 | printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]); | |
1202 | printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]); | |
1203 | printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]); | |
1204 | printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]); | |
1205 | printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]); | |
1206 | printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]); | |
1207 | printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]); | |
1208 | } | |
1209 | if(fJPSItoEle) { | |
1210 | printf("Reconstruct J/psi from B candidates with cuts:\n"); | |
1211 | printf(" |M-MJPSI| [GeV] < %f\n",fBtoJPSICuts[0]); | |
1212 | printf(" dca [cm] < %f\n",fBtoJPSICuts[1]); | |
1213 | printf(" cosThetaStar < %f\n",fBtoJPSICuts[2]); | |
1214 | printf(" pTP [GeV/c] > %f\n",fBtoJPSICuts[3]); | |
1215 | printf(" pTN [GeV/c] > %f\n",fBtoJPSICuts[4]); | |
1216 | printf(" |d0P| [cm] < %f\n",fBtoJPSICuts[5]); | |
1217 | printf(" |d0N| [cm] < %f\n",fBtoJPSICuts[6]); | |
1218 | printf(" d0d0 [cm^2] < %f\n",fBtoJPSICuts[7]); | |
1219 | printf(" cosThetaPoint > %f\n",fBtoJPSICuts[8]); | |
1220 | } | |
1221 | if(f3Prong) { | |
1222 | printf("Reconstruct 3 prong candidates.\n"); | |
4d5c9633 | 1223 | printf(" D+->Kpipi cuts:\n"); |
6a213b59 | 1224 | printf(" |M-MD+| [GeV] < %f\n",fDplusCuts[0]); |
1225 | printf(" pTK [GeV/c] > %f\n",fDplusCuts[1]); | |
1226 | printf(" pTPi [GeV/c] > %f\n",fDplusCuts[2]); | |
1227 | printf(" |d0K| [cm] > %f\n",fDplusCuts[3]); | |
1228 | printf(" |d0Pi| [cm] > %f\n",fDplusCuts[4]); | |
1229 | printf(" dist12 [cm] < %f\n",fDplusCuts[5]); | |
1230 | printf(" sigmavert [cm] < %f\n",fDplusCuts[6]); | |
1231 | printf(" dist prim-sec [cm] > %f\n",fDplusCuts[7]); | |
1232 | printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]); | |
1233 | printf(" cosThetaPoint > %f\n",fDplusCuts[9]); | |
1234 | printf(" Sum d0^2 [cm^2] > %f\n",fDplusCuts[10]); | |
1235 | printf(" dca cut [cm] < %f\n",fDplusCuts[11]); | |
4d5c9633 | 1236 | printf(" Ds->KKpi cuts:\n"); |
1237 | printf(" |M-MDs| [GeV] < %f\n",fDsCuts[0]); | |
1238 | printf(" Lc->pKpi cuts:\n"); | |
1239 | printf(" |M-MD+| [GeV] < %f\n",fLcCuts[0]); | |
6a213b59 | 1240 | } |
1241 | ||
1242 | return; | |
1243 | } | |
1244 | //----------------------------------------------------------------------------- | |
1245 | Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay, | |
1246 | Int_t nprongs, | |
1247 | Double_t *px, | |
1248 | Double_t *py, | |
1249 | Double_t *pz) const { | |
1250 | // Check invariant mass cut | |
1251 | ||
1252 | Short_t dummycharge=0; | |
1253 | Double_t *dummyd0 = new Double_t[nprongs]; | |
b056c5e3 | 1254 | for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.; |
6a213b59 | 1255 | AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0); |
1256 | delete [] dummyd0; | |
1257 | ||
1258 | UInt_t pdg2[2],pdg3[3]; | |
1259 | Double_t mPDG,minv; | |
1260 | ||
1261 | Bool_t retval=kFALSE; | |
1262 | switch (decay) | |
1263 | { | |
1264 | case 0: // D0->Kpi | |
1265 | pdg2[0]=211; pdg2[1]=321; | |
1266 | mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
1267 | minv = rd->InvMass(nprongs,pdg2); | |
1268 | if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE; | |
1269 | pdg2[0]=321; pdg2[1]=211; | |
1270 | minv = rd->InvMass(nprongs,pdg2); | |
1271 | if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE; | |
1272 | break; | |
1273 | case 1: // JPSI->ee | |
1274 | pdg2[0]=11; pdg2[1]=11; | |
1275 | mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass(); | |
1276 | minv = rd->InvMass(nprongs,pdg2); | |
1277 | if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE; | |
1278 | break; | |
1279 | case 2: // D+->Kpipi | |
1280 | pdg3[0]=211; pdg3[1]=321; pdg3[2]=211; | |
1281 | mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass(); | |
1282 | minv = rd->InvMass(nprongs,pdg3); | |
1283 | if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE; | |
b82f6d67 | 1284 | // Ds+->KKpi |
1285 | pdg3[0]=321; pdg3[1]=321; pdg3[2]=211; | |
1286 | mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass(); | |
1287 | minv = rd->InvMass(nprongs,pdg3); | |
1288 | if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE; | |
1289 | pdg3[0]=211; pdg3[1]=321; pdg3[2]=321; | |
1290 | minv = rd->InvMass(nprongs,pdg3); | |
1291 | if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE; | |
1292 | // Lc->pKpi | |
1293 | pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211; | |
1294 | mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass(); | |
1295 | minv = rd->InvMass(nprongs,pdg3); | |
1296 | if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE; | |
1297 | pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212; | |
1298 | minv = rd->InvMass(nprongs,pdg3); | |
1299 | if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE; | |
6a213b59 | 1300 | break; |
1301 | default: | |
1302 | break; | |
1303 | } | |
1304 | ||
1305 | delete rd; | |
1306 | ||
1307 | return retval; | |
1308 | } | |
1309 | //----------------------------------------------------------------------------- | |
5e6f195c | 1310 | void AliAnalysisVertexingHF::SelectTracks(AliESDEvent *esd, |
6a213b59 | 1311 | TObjArray &trksP,Int_t &nTrksP, |
1312 | TObjArray &trksN,Int_t &nTrksN) const | |
1313 | { | |
1314 | // Fill two TObjArrays with positive and negative tracks and | |
1315 | // apply single-track preselection | |
1316 | ||
1317 | nTrksP=0,nTrksN=0; | |
1318 | ||
1319 | Int_t entries = (Int_t)esd->GetNumberOfTracks(); | |
1320 | ||
1321 | // transfer ITS tracks from ESD to arrays and to a tree | |
1322 | for(Int_t i=0; i<entries; i++) { | |
1323 | ||
1324 | AliESDtrack *esdtrack = esd->GetTrack(i); | |
1325 | UInt_t status = esdtrack->GetStatus(); | |
1326 | ||
1327 | // require refit in ITS | |
1328 | if(fITSrefit && !(status&AliESDtrack::kITSrefit)) { | |
1329 | if(fDebug) printf("track %d is not kITSrefit\n",i); | |
1330 | continue; | |
1331 | } | |
1332 | ||
1333 | // require minimum # of ITS points | |
1334 | if(esdtrack->GetNcls(0)<fMinITSCls) { | |
1335 | if(fDebug) printf("track %d has %d ITS cls\n",i,esdtrack->GetNcls(0)); | |
1336 | continue; | |
1337 | } | |
1338 | // require points on the 2 pixel layers | |
1339 | if(fBothSPD) | |
1340 | if(!TESTBIT(esdtrack->GetITSClusterMap(),0) || | |
1341 | !TESTBIT(esdtrack->GetITSClusterMap(),1)) continue; | |
1342 | ||
1343 | // single track selection | |
b82f6d67 | 1344 | if(!SingleTrkCuts(*esdtrack)) continue; |
6a213b59 | 1345 | |
1346 | if(esdtrack->GetSign()<0) { // negative track | |
1347 | trksN.AddLast(esdtrack); | |
1348 | nTrksN++; | |
1349 | } else { // positive track | |
1350 | trksP.AddLast(esdtrack); | |
1351 | nTrksP++; | |
1352 | } | |
1353 | ||
1354 | } // loop on ESD tracks | |
1355 | ||
1356 | return; | |
1357 | } | |
1358 | //----------------------------------------------------------------------------- | |
1359 | void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1, | |
1360 | Double_t cut2,Double_t cut3,Double_t cut4, | |
1361 | Double_t cut5,Double_t cut6, | |
1362 | Double_t cut7,Double_t cut8) | |
1363 | { | |
1364 | // Set the cuts for D0 selection | |
1365 | fD0toKpiCuts[0] = cut0; | |
1366 | fD0toKpiCuts[1] = cut1; | |
1367 | fD0toKpiCuts[2] = cut2; | |
1368 | fD0toKpiCuts[3] = cut3; | |
1369 | fD0toKpiCuts[4] = cut4; | |
1370 | fD0toKpiCuts[5] = cut5; | |
1371 | fD0toKpiCuts[6] = cut6; | |
1372 | fD0toKpiCuts[7] = cut7; | |
1373 | fD0toKpiCuts[8] = cut8; | |
1374 | ||
1375 | return; | |
1376 | } | |
1377 | //----------------------------------------------------------------------------- | |
1378 | void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9]) | |
1379 | { | |
1380 | // Set the cuts for D0 selection | |
1381 | ||
1382 | for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i]; | |
1383 | ||
1384 | return; | |
1385 | } | |
1386 | //----------------------------------------------------------------------------- | |
1387 | void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1, | |
1388 | Double_t cut2,Double_t cut3,Double_t cut4, | |
1389 | Double_t cut5,Double_t cut6, | |
1390 | Double_t cut7,Double_t cut8) | |
1391 | { | |
1392 | // Set the cuts for J/psi from B selection | |
1393 | fBtoJPSICuts[0] = cut0; | |
1394 | fBtoJPSICuts[1] = cut1; | |
1395 | fBtoJPSICuts[2] = cut2; | |
1396 | fBtoJPSICuts[3] = cut3; | |
1397 | fBtoJPSICuts[4] = cut4; | |
1398 | fBtoJPSICuts[5] = cut5; | |
1399 | fBtoJPSICuts[6] = cut6; | |
1400 | fBtoJPSICuts[7] = cut7; | |
1401 | fBtoJPSICuts[8] = cut8; | |
1402 | ||
1403 | return; | |
1404 | } | |
1405 | //----------------------------------------------------------------------------- | |
1406 | void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9]) | |
1407 | { | |
1408 | // Set the cuts for J/psi from B selection | |
1409 | ||
1410 | for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i]; | |
1411 | ||
1412 | return; | |
1413 | } | |
1414 | //----------------------------------------------------------------------------- | |
1415 | void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1, | |
1416 | Double_t cut2,Double_t cut3,Double_t cut4, | |
1417 | Double_t cut5,Double_t cut6, | |
1418 | Double_t cut7,Double_t cut8, | |
1419 | Double_t cut9,Double_t cut10,Double_t cut11) | |
1420 | { | |
1421 | // Set the cuts for Dplus->Kpipi selection | |
1422 | fDplusCuts[0] = cut0; | |
1423 | fDplusCuts[1] = cut1; | |
1424 | fDplusCuts[2] = cut2; | |
1425 | fDplusCuts[3] = cut3; | |
1426 | fDplusCuts[4] = cut4; | |
1427 | fDplusCuts[5] = cut5; | |
1428 | fDplusCuts[6] = cut6; | |
1429 | fDplusCuts[7] = cut7; | |
1430 | fDplusCuts[8] = cut8; | |
1431 | fDplusCuts[9] = cut9; | |
1432 | fDplusCuts[10] = cut10; | |
1433 | fDplusCuts[11] = cut11; | |
1434 | ||
1435 | return; | |
1436 | } | |
1437 | //----------------------------------------------------------------------------- | |
1438 | void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12]) | |
1439 | { | |
b82f6d67 | 1440 | // Set the cuts for Dplus->Kpipi selection |
6a213b59 | 1441 | |
1442 | for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i]; | |
1443 | ||
1444 | return; | |
1445 | } | |
1446 | //----------------------------------------------------------------------------- | |
b82f6d67 | 1447 | void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0) |
1448 | { | |
1449 | // Set the cuts for Ds->KKpi selection | |
1450 | fDsCuts[0] = cut0; | |
1451 | ||
1452 | return; | |
1453 | } | |
1454 | //----------------------------------------------------------------------------- | |
1455 | void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[1]) | |
1456 | { | |
1457 | // Set the cuts for Ds->KKpi selection | |
1458 | ||
1459 | for(Int_t i=0; i<1; i++) fDsCuts[i] = cuts[i]; | |
1460 | ||
1461 | return; | |
1462 | } | |
1463 | //----------------------------------------------------------------------------- | |
1464 | void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0) | |
1465 | { | |
1466 | // Set the cuts for Lc->pKpi selection | |
1467 | fLcCuts[0] = cut0; | |
1468 | ||
1469 | return; | |
1470 | } | |
1471 | //----------------------------------------------------------------------------- | |
1472 | void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[1]) | |
1473 | { | |
1474 | // Set the cuts for Lc->pKpi selection | |
1475 | ||
1476 | for(Int_t i=0; i<1; i++) fLcCuts[i] = cuts[i]; | |
1477 | ||
1478 | return; | |
1479 | } | |
1480 | //----------------------------------------------------------------------------- | |
1481 | Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk) const | |
6a213b59 | 1482 | { |
1483 | // Check if track passes some kinematical cuts | |
6a213b59 | 1484 | |
b056c5e3 | 1485 | if(TMath::Abs(trk.Pt()) < fMinPtCut) { |
1486 | //printf("pt %f\n",1./trk.GetParameter()[4]); | |
6a213b59 | 1487 | return kFALSE; |
1488 | } | |
b82f6d67 | 1489 | trk.RelateToVertex(fV1,fBzkG,10.); |
6a213b59 | 1490 | Float_t d0z0[2],covd0z0[3]; |
1491 | trk.GetImpactParameters(d0z0,covd0z0); | |
1492 | if(TMath::Abs(d0z0[0]) < fMind0rphiCut) { | |
1493 | printf("d0rphi %f\n",TMath::Abs(d0z0[0])); | |
1494 | return kFALSE; | |
1495 | } | |
1496 | ||
1497 | return kTRUE; | |
1498 | } | |
1499 | //----------------------------------------------------------------------------- | |
b82f6d67 | 1500 | AliESDVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray) const |
1501 | { | |
1502 | // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle | |
1503 | ||
1504 | AliESDVertex *vertex = 0; | |
1505 | ||
1506 | if(!fSecVtxWithKF) { // AliVertexerTracks | |
1507 | ||
1508 | AliVertexerTracks *vertexer2 = new AliVertexerTracks(fBzkG); | |
1509 | vertexer2->SetDebug(0); | |
1510 | vertexer2->SetVtxStart(fV1); | |
f65f93f3 | 1511 | vertex = (AliESDVertex*)vertexer2->VertexForSelectedESDTracks(trkArray); |
b82f6d67 | 1512 | delete vertexer2; |
1513 | ||
b056c5e3 | 1514 | if(vertex->GetNContributors()!=trkArray->GetEntriesFast()) { |
1515 | if(fDebug) printf("vertexing failed\n"); | |
1516 | delete vertex; vertex=0; | |
1517 | } | |
1518 | ||
b82f6d67 | 1519 | } else { // Kalman Filter vertexer (AliKFParticle) |
1520 | ||
1521 | AliKFParticle::SetField(fBzkG); | |
1522 | ||
1523 | AliKFParticle vertexKF; | |
1524 | ||
1525 | Int_t nTrks = trkArray->GetEntriesFast(); | |
1526 | for(Int_t i=0; i<nTrks; i++) { | |
1527 | AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i); | |
1528 | AliKFParticle daughterKF(*esdTrack,211); | |
1529 | vertexKF.AddDaughter(daughterKF); | |
1530 | } | |
1531 | vertex = new AliESDVertex(); | |
1532 | vertexKF.CopyToESDVertex(*vertex); | |
1533 | ||
1534 | } | |
1535 | ||
1536 | return vertex; | |
1537 | } | |
1538 | //----------------------------------------------------------------------------- | |
6a213b59 | 1539 | |
1540 | ||
1541 |