1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
21 // classes AliAnalysisTaskSEVertexingHF or AliAnalysisTaskVertexingHF.
23 // Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita
24 // Contact: andrea.dainese@lnl.infn.it
25 //----------------------------------------------------------------------------
28 #include <TDatabasePDG.h>
30 #include "AliESDEvent.h"
31 #include "AliVertexerTracks.h"
32 #include "AliKFParticle.h"
33 #include "AliESDVertex.h"
34 #include "AliAODEvent.h"
35 #include "AliAODRecoDecay.h"
36 #include "AliAODRecoDecayHF.h"
37 #include "AliAODRecoDecayHF2Prong.h"
38 #include "AliAODRecoDecayHF3Prong.h"
39 #include "AliAODRecoDecayHF4Prong.h"
40 #include "AliAnalysisVertexingHF.h"
42 ClassImp(AliAnalysisVertexingHF)
44 //----------------------------------------------------------------------------
45 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
47 fSecVtxWithKF(kFALSE),
49 fRecoPrimVtxSkippingTrks(kFALSE),
50 fRmTrksFromPrimVtx(kFALSE),
63 // Default constructor
71 //--------------------------------------------------------------------------
72 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
75 fSecVtxWithKF(source.fSecVtxWithKF),
76 fUseTRef(source.fUseTRef),
77 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
78 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
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)
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];
97 for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
98 for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
100 //--------------------------------------------------------------------------
101 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
104 // assignment operator
106 if(&source == this) return *this;
107 fBzkG = source.fBzkG;
108 fSecVtxWithKF = source.fSecVtxWithKF;
109 fUseTRef = source.fUseTRef;
110 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
111 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
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;
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];
127 for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
128 for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
132 //----------------------------------------------------------------------------
133 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
135 if(fV1) { delete fV1; fV1=0; }
137 //----------------------------------------------------------------------------
138 void AliAnalysisVertexingHF::FindCandidatesESDtoAOD(AliESDEvent *esd,
139 TClonesArray *aodVerticesHFTClArr,
140 TClonesArray *aodD0toKpiTClArr,
141 TClonesArray *aodJPSItoEleTClArr,
142 TClonesArray *aodCharm3ProngTClArr,
143 TClonesArray *aodCharm4ProngTClArr)
145 // Find heavy-flavour vertex candidates
147 // Output: AOD (additional branches added)
149 if(!aodVerticesHFTClArr) {
150 printf("ERROR: no aodVerticesHFTClArr");
153 if(fD0toKpi && !aodD0toKpiTClArr) {
154 printf("ERROR: no aodD0toKpiTClArr");
157 if(fJPSItoEle && !aodJPSItoEleTClArr) {
158 printf("ERROR: no aodJPSItoEleTClArr");
161 if(f3Prong && !aodCharm3ProngTClArr) {
162 printf("ERROR: no aodCharm3ProngTClArr");
165 if(f4Prong && !aodCharm4ProngTClArr) {
166 printf("ERROR: no aodCharm4ProngTClArr");
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;
176 aodD0toKpiTClArr->Delete();
177 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
180 aodJPSItoEleTClArr->Delete();
181 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
184 aodCharm3ProngTClArr->Delete();
185 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
188 aodCharm4ProngTClArr->Delete();
189 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
191 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
192 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
193 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
194 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
197 AliAODRecoDecayHF2Prong *io2Prong = 0;
198 AliAODRecoDecayHF3Prong *io3Prong = 0;
199 AliAODRecoDecayHF4Prong *io4Prong = 0;
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);
214 //------- SINGLE EVENT ANALYSIS ---------------------------------
216 Int_t ev = (Int_t)esd->GetEventNumberInFile();
217 printf("--- Finding candidates in event %d\n",ev);
221 fBzkG = (Double_t)esd->GetMagneticField();
223 trkEntries = (Int_t)esd->GetNumberOfTracks();
224 printf(" Number of tracks: %d\n",trkEntries);
226 if(trkEntries<2 || !esd->GetPrimaryVertex()) {
230 // retrieve primary vertex from the AliESDEvent
231 AliESDVertex copy(*(esd->GetPrimaryVertex()));
232 SetPrimaryVertex(©);
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);
240 printf(" Pos. tracks: %d Neg. tracks: %d\n",nTrksP,nTrksN);
242 TObjArray *twoTrackArray1 = new TObjArray(2);
243 TObjArray *twoTrackArray2 = new TObjArray(2);
244 TObjArray *threeTrackArray = new TObjArray(3);
245 TObjArray *fourTrackArray = new TObjArray(4);
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);
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; }
262 twoTrackArray1->AddAt(postrack1,0);
263 twoTrackArray1->AddAt(negtrack1,1);
264 AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1);
266 twoTrackArray1->Clear();
270 if(fD0toKpi || fJPSItoEle) {
271 io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI);
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()};
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);
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);
295 if(!okD0) v->SetParent(rd); // do something better here...
297 //printf("DCA: %f\n",rd->GetDCA());
299 if(okD0) new(aodD0toKpiRef[iD0toKpi++]) AliAODRecoDecayHF2Prong(*io2Prong);
300 if(okJPSI) new(aodJPSItoEleRef[iJPSItoEle++]) AliAODRecoDecayHF2Prong(*io2Prong);
307 twoTrackArray1->Clear();
308 if(!f3Prong && !f4Prong) {
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; }
325 twoTrackArray2->AddAt(postrack2,0);
326 twoTrackArray2->AddAt(negtrack1,1);
327 AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2);
329 twoTrackArray2->Clear();
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);
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);
356 new(aodCharm3ProngRef[i3Prong++]) AliAODRecoDecayHF3Prong(*io3Prong);
359 if(io3Prong) { /*delete io3Prong;*/ io3Prong=NULL; }
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; }
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);
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);
397 new(aodCharm4ProngRef[i4Prong++]) AliAODRecoDecayHF4Prong(*io4Prong);
400 if(io4Prong) { /*delete io4Prong;*/ io4Prong=NULL; }
401 fourTrackArray->Clear();
403 } // end loop on negative tracks
407 } // end 2nd loop on positive tracks
408 twoTrackArray2->Clear();
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; }
420 twoTrackArray2->AddAt(postrack1,0);
421 twoTrackArray2->AddAt(negtrack2,1);
422 AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2);
424 twoTrackArray2->Clear();
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);
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);
451 new(aodCharm3ProngRef[i3Prong++]) AliAODRecoDecayHF3Prong(*io3Prong);
454 if(io3Prong) { /*delete io3Prong;*/ io3Prong=NULL; }
458 } // end 2nd loop on negative tracks
459 twoTrackArray2->Clear();
463 } // end 1st loop on negative tracks
466 } // end 1st loop on positive tracks
470 printf(" D0->Kpi in event %d = %d;\n",
471 (Int_t)esd->GetEventNumberInFile(),
472 (Int_t)aodD0toKpiTClArr->GetEntriesFast());
475 printf(" JPSI->ee in event %d = %d;\n",
476 (Int_t)esd->GetEventNumberInFile(),
477 (Int_t)aodJPSItoEleTClArr->GetEntriesFast());
480 printf(" Charm->3Prong in event %d = %d;\n",
481 (Int_t)esd->GetEventNumberInFile(),
482 (Int_t)aodCharm3ProngTClArr->GetEntriesFast());
485 printf(" Charm->4Prong in event %d = %d;\n",
486 (Int_t)esd->GetEventNumberInFile(),
487 (Int_t)aodCharm4ProngTClArr->GetEntriesFast());
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;
501 //------- END SINGLE EVENT ANALYSIS --------------------------------
505 //----------------------------------------------------------------------------
506 void AliAnalysisVertexingHF::FindCandidates(AliESDEvent *esd,TTree *treeout[])
508 // Find heavy-flavour vertex candidates
510 // DEPRECATED: use FindCandidatesESDtoAOD!
512 fUseTRef=kFALSE; // cannot use TRefs outside AOD
514 AliAODRecoDecayHF2Prong *io2Prong = new AliAODRecoDecayHF2Prong();
515 AliAODRecoDecayHF3Prong *io3Prong = new AliAODRecoDecayHF3Prong();
516 AliAODRecoDecayHF4Prong *io4Prong = new AliAODRecoDecayHF4Prong();
518 Int_t itreeD0toKpi=-1,itreeJPSItoEle=-1,itree3Prong=-1,itree4Prong=-1;
519 Int_t initEntriesD0toKpi=0,initEntriesJPSItoEle=0,initEntries3Prong=0,initEntries4Prong=0;
522 treeout[itree]->SetBranchAddress("D0toKpi",&io2Prong);
524 initEntriesD0toKpi = treeout[itreeD0toKpi]->GetEntries();
527 itreeJPSItoEle=itree;
528 treeout[itree]->SetBranchAddress("JPSItoEle",&io2Prong);
530 initEntriesJPSItoEle = treeout[itreeJPSItoEle]->GetEntries();
534 treeout[itree]->SetBranchAddress("Charmto3Prong",&io3Prong);
536 initEntries3Prong = treeout[itree3Prong]->GetEntries();
540 treeout[itree]->SetBranchAddress("D0to4Prong",&io4Prong);
542 initEntries4Prong = treeout[itree4Prong]->GetEntries();
544 delete io2Prong; io2Prong = NULL;
545 delete io3Prong; io3Prong = NULL;
546 delete io4Prong; io4Prong = NULL;
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);
561 Int_t ev = (Int_t)esd->GetEventNumberInFile();
562 printf("--- Finding candidates in event %d\n",ev);
564 fBzkG = (Double_t)esd->GetMagneticField();
566 trkEntries = (Int_t)esd->GetNumberOfTracks();
567 printf(" Number of tracks: %d\n",trkEntries);
568 if(trkEntries<2) return;
570 // retrieve primary vertex from the AliESDEvent
571 if(!esd->GetPrimaryVertex()) {
572 printf(" No vertex in AliESD\n");
575 AliESDVertex copy(*(esd->GetPrimaryVertex()));
576 SetPrimaryVertex(©);
578 // call function which applies sigle-track selection and
579 // separetes positives and negatives
580 TObjArray trksP(trkEntries/2);
581 TObjArray trksN(trkEntries/2);
582 SelectTracks(esd,trksP,nTrksP,
585 printf(" Pos. tracks: %d Neg. tracks: %d\n",nTrksP,nTrksN);
587 TObjArray *twoTrackArray1 = new TObjArray(2);
588 TObjArray *twoTrackArray2 = new TObjArray(2);
589 TObjArray *threeTrackArray = new TObjArray(3);
590 TObjArray *fourTrackArray = new TObjArray(4);
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);
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
604 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
605 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
607 twoTrackArray1->AddAt(postrack1,0);
608 twoTrackArray1->AddAt(negtrack1,1);
609 AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1);
611 twoTrackArray1->Clear();
615 if(fD0toKpi || fJPSItoEle) {
616 io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI);
617 if(okD0) treeout[itreeD0toKpi]->Fill();
618 if(okJPSI) treeout[itreeJPSItoEle]->Fill();
619 delete io2Prong; io2Prong=NULL;
622 twoTrackArray1->Clear();
623 if(!f3Prong && !f4Prong) {
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);
633 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
634 if(dcap2n1>dcaMax) { postrack2=0; continue; }
635 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
636 if(dcap1p2>dcaMax) { postrack2=0; continue; }
639 twoTrackArray2->AddAt(postrack2,0);
640 twoTrackArray2->AddAt(negtrack1,1);
641 AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2);
643 twoTrackArray2->Clear();
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);
652 if(ok3Prong) treeout[itree3Prong]->Fill();
653 if(io3Prong) delete io3Prong; io3Prong=NULL;
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);
660 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
661 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
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);
668 if(ok4Prong) treeout[itree4Prong]->Fill();
669 delete io4Prong; io4Prong=NULL;
670 fourTrackArray->Clear();
672 } // end loop on negative tracks
676 } // end 2nd loop on positive tracks
677 twoTrackArray2->Clear();
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);
683 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
684 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
685 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
686 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
689 twoTrackArray2->AddAt(postrack1,0);
690 twoTrackArray2->AddAt(negtrack2,1);
691 AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2);
693 twoTrackArray2->Clear();
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);
702 if(ok3Prong) treeout[itree3Prong]->Fill();
703 if(io3Prong) delete io3Prong; io3Prong=NULL;
707 } // end 2nd loop on negative tracks
708 twoTrackArray2->Clear();
712 } // end 1st loop on negative tracks
715 } // end 1st loop on positive tracks
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;
730 // create a copy of this class to be written to output file
731 //AliAnalysisVertexingHF *copy = (AliAnalysisVertexingHF*)this->Clone("AnalysisVertexingHF");
735 printf(" D0->Kpi: event %d = %d; total = %d;\n",
736 (Int_t)esd->GetEventNumberInFile(),
737 (Int_t)treeout[itreeD0toKpi]->GetEntries()-initEntriesD0toKpi,
738 (Int_t)treeout[itreeD0toKpi]->GetEntries());
741 printf(" JPSI->ee: event %d = %d; total = %d;\n",
742 (Int_t)esd->GetEventNumberInFile(),
743 (Int_t)treeout[itreeJPSItoEle]->GetEntries()-initEntriesJPSItoEle,
744 (Int_t)treeout[itreeJPSItoEle]->GetEntries());
747 printf(" Charm->3Prong: event %d = %d; total = %d;\n",
748 (Int_t)esd->GetEventNumberInFile(),
749 (Int_t)treeout[itree3Prong]->GetEntries()-initEntries3Prong,
750 (Int_t)treeout[itree3Prong]->GetEntries());
753 printf(" Charm->4Prong: event %d = %d; total = %d;\n",
754 (Int_t)esd->GetEventNumberInFile(),
755 (Int_t)treeout[itree4Prong]->GetEntries()-initEntries4Prong,
756 (Int_t)treeout[itree4Prong]->GetEntries());
762 //----------------------------------------------------------------------------
763 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
764 TObjArray *twoTrackArray1,AliESDEvent *esd,
765 AliESDVertex *secVertexESD,Double_t dca,
766 Bool_t &okD0,Bool_t &okJPSI) const
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)
772 okD0=kFALSE; okJPSI=kFALSE;
774 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
776 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0);
777 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1);
779 // propagate tracks to secondary vertex, to compute inv. mass
780 postrack->RelateToVertex(secVertexESD,fBzkG,10.);
781 negtrack->RelateToVertex(secVertexESD,fBzkG,10.);
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];
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;
795 if(fDebug) printf(" candidate didn't pass mass cut\n");
800 AliESDVertex *primVertex = fV1;
801 AliESDVertex *ownPrimVertex=0;
803 // primary vertex from *other* tracks in the event
804 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
805 ownPrimVertex = OwnPrimaryVertex(2,twoTrackArray1,esd);
809 if(ownPrimVertex->GetNContributors()<2) {
810 delete ownPrimVertex;
813 primVertex = ownPrimVertex;
818 Float_t d0z0[2],covd0z0[3];
819 postrack->RelateToVertex(primVertex,fBzkG,10.);
820 postrack->GetImpactParameters(d0z0,covd0z0);
822 d0err[0] = TMath::Sqrt(covd0z0[0]);
823 negtrack->RelateToVertex(primVertex,fBzkG,10.);
824 negtrack->GetImpactParameters(d0z0,covd0z0);
826 d0err[1] = TMath::Sqrt(covd0z0[0]);
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());
833 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVertexAOD,px,py,pz,d0,d0err,dca);
834 the2Prong->SetOwnSecondaryVtx(secVertexAOD);
835 primVertex->GetXYZ(pos); // position
836 primVertex->GetCovMatrix(cov); //covariance matrix
837 AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
838 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
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
847 if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
848 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
852 // get PID info from ESD
854 postrack->GetESDpid(esdpid0);
856 negtrack->GetESDpid(esdpid1);
858 for(Int_t i=0;i<5;i++) {
859 esdpid[i] = esdpid0[i];
860 esdpid[5+i] = esdpid1[i];
862 the2Prong->SetPID(2,esdpid);
865 if(ownPrimVertex) delete ownPrimVertex;
869 //----------------------------------------------------------------------------
870 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
871 TObjArray *threeTrackArray,AliESDEvent *esd,
872 AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
873 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
874 Bool_t &ok3Prong) const
876 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
877 // reconstruction cuts
881 Double_t px[3],py[3],pz[3],d0[3],d0err[3];//d0z[3];
882 Float_t d0z0[2],covd0z0[3];
885 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
886 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
887 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
889 AliESDVertex *primVertex = fV1;
891 postrack1->RelateToVertex(primVertex,fBzkG,10.);
892 negtrack->RelateToVertex(primVertex,fBzkG,10.);
893 postrack2->RelateToVertex(primVertex,fBzkG,10.);
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];
903 postrack1->GetImpactParameters(d0z0,covd0z0);
905 d0err[0] = TMath::Sqrt(covd0z0[0]);
906 negtrack->GetImpactParameters(d0z0,covd0z0);
908 d0err[1] = TMath::Sqrt(covd0z0[0]);
909 postrack2->GetImpactParameters(d0z0,covd0z0);
911 d0err[2] = TMath::Sqrt(covd0z0[0]);
914 // invariant mass cut for D+, Ds, Lc
915 Bool_t okMassCut=kFALSE;
916 if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
918 if(fDebug) printf(" candidate didn't pass mass cut\n");
923 Short_t charge=(Short_t)(postrack1->GetSign()*postrack2->GetSign()*negtrack->GetSign());
925 AliESDVertex *ownPrimVertex = 0;
926 // primary vertex from *other* tracks in the event
927 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
928 ownPrimVertex = OwnPrimaryVertex(3,threeTrackArray,esd);
932 if(ownPrimVertex->GetNContributors()<2) {
933 delete ownPrimVertex;
936 primVertex = ownPrimVertex;
941 // create the object AliAODRecoDecayHF3Prong
942 AliESDVertex* secVert3Prong = ReconstructSecondaryVertex(threeTrackArray);
944 if(ownPrimVertex) delete ownPrimVertex;
947 Double_t pos[3],cov[6],sigmavert;
948 secVert3Prong->GetXYZ(pos); // position
949 secVert3Prong->GetCovMatrix(cov); //covariance matrix
950 sigmavert=secVert3Prong->GetDispersion();
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};
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]));
961 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert3PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
962 the3Prong->SetOwnSecondaryVtx(secVert3PrAOD);
963 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
966 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
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;
974 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
976 // get PID info from ESD
978 postrack1->GetESDpid(esdpid0);
980 negtrack->GetESDpid(esdpid1);
982 postrack2->GetESDpid(esdpid2);
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];
991 the3Prong->SetPID(3,esdpid);
994 if(ownPrimVertex) delete ownPrimVertex;
998 //----------------------------------------------------------------------------
999 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1000 TObjArray *fourTrackArray,AliESDEvent *esd,
1001 AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
1002 Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
1003 Bool_t &ok4Prong) const
1005 // Make 4Prong candidates and check if they pass D0toKpipipi
1006 // reconstruction cuts
1007 // G.E.Bruno, R.Romita
1011 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1012 //Float_t d0z0[2],covd0z0[3];
1014 px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
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);
1024 AliESDVertex *primVertex = fV1;
1026 postrack1->RelateToVertex(primVertex,fBzkG,10.);
1027 negtrack1->RelateToVertex(primVertex,fBzkG,10.);
1028 postrack2->RelateToVertex(primVertex,fBzkG,10.);
1029 negtrack2->RelateToVertex(primVertex,fBzkG,10.);
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];
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;
1045 // if(fDebug) printf(" candidate didn't pass mass cut\n");
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) {
1056 if(ownPrimVertex->GetNContributors()<2) {
1057 delete ownPrimVertex;
1060 primVertex = ownPrimVertex;
1065 // create the object AliAODRecoDecayHF4Prong
1066 AliESDVertex* secVert4Prong = ReconstructSecondaryVertex(fourTrackArray);
1067 if(!secVert4Prong) {
1068 if(ownPrimVertex) delete ownPrimVertex;
1071 Double_t pos[3],cov[6],sigmavert;
1072 secVert4Prong->GetXYZ(pos); // position
1073 secVert4Prong->GetCovMatrix(cov); //covariance matrix
1074 sigmavert=secVert4Prong->GetDispersion();
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
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
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);
1091 the4Prong->SetOwnSecondaryVtx(secVert4PrAOD);
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 ...
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);
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];
1118 the4Prong->SetPID(4,esdpid);
1120 if(ownPrimVertex) delete ownPrimVertex;
1124 //-----------------------------------------------------------------------------
1125 AliESDVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(Int_t ntrks,
1126 TObjArray *trkArray,
1127 AliESDEvent *esd) const
1129 // Returns primary vertex specific to this candidate
1131 AliVertexerTracks *vertexer1 = new AliVertexerTracks(esd->GetMagneticField());
1132 AliESDVertex *ownPrimVertex = 0;
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();
1149 for(Int_t i=0; i<ntrks; i++) {
1150 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1151 skipped[i] = (Int_t)t->GetID();
1153 vertexer1->SetSkipTracks(ntrks,skipped);
1154 ownPrimVertex = (AliESDVertex*)vertexer1->FindPrimaryVertex(esd);
1157 // removing the prongs tracks
1158 if(fRmTrksFromPrimVtx) {
1159 TObjArray rmArray(ntrks);
1160 UShort_t *rmId = new UShort_t[ntrks];
1161 AliESDtrack *esdTrack = 0;
1163 for(Int_t i=0; i<ntrks; i++) {
1164 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1165 esdTrack = new AliESDtrack(*t);
1166 rmArray.AddLast(esdTrack);
1167 rmId[i]=(UShort_t)esdTrack->GetID();
1169 Float_t diamondxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
1170 ownPrimVertex = vertexer1->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1171 delete [] rmId; rmId=NULL;
1175 delete vertexer1; vertexer1=NULL;
1177 return ownPrimVertex;
1179 //-----------------------------------------------------------------------------
1180 void AliAnalysisVertexingHF::PrintStatus() const {
1181 // Print parameters being used
1183 printf("Preselections:\n");
1184 printf(" fITSrefit = %d\n",(Int_t)fITSrefit);
1185 printf(" fBothSPD = %d\n",(Int_t)fBothSPD);
1186 printf(" fMinITSCls = %d\n",fMinITSCls);
1187 printf(" fMinPtCut = %f GeV/c\n",fMinPtCut);
1188 printf(" fMind0rphiCut = %f cm\n",fMind0rphiCut);
1190 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1192 printf("Secondary vertex with AliVertexerTracks\n");
1194 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1195 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1197 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1198 printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
1199 printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
1200 printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
1201 printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
1202 printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
1203 printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
1204 printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
1205 printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
1206 printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
1209 printf("Reconstruct J/psi from B candidates with cuts:\n");
1210 printf(" |M-MJPSI| [GeV] < %f\n",fBtoJPSICuts[0]);
1211 printf(" dca [cm] < %f\n",fBtoJPSICuts[1]);
1212 printf(" cosThetaStar < %f\n",fBtoJPSICuts[2]);
1213 printf(" pTP [GeV/c] > %f\n",fBtoJPSICuts[3]);
1214 printf(" pTN [GeV/c] > %f\n",fBtoJPSICuts[4]);
1215 printf(" |d0P| [cm] < %f\n",fBtoJPSICuts[5]);
1216 printf(" |d0N| [cm] < %f\n",fBtoJPSICuts[6]);
1217 printf(" d0d0 [cm^2] < %f\n",fBtoJPSICuts[7]);
1218 printf(" cosThetaPoint > %f\n",fBtoJPSICuts[8]);
1221 printf("Reconstruct 3 prong candidates.\n");
1222 printf(" D+->Kpipi cuts:\n");
1223 printf(" |M-MD+| [GeV] < %f\n",fDplusCuts[0]);
1224 printf(" pTK [GeV/c] > %f\n",fDplusCuts[1]);
1225 printf(" pTPi [GeV/c] > %f\n",fDplusCuts[2]);
1226 printf(" |d0K| [cm] > %f\n",fDplusCuts[3]);
1227 printf(" |d0Pi| [cm] > %f\n",fDplusCuts[4]);
1228 printf(" dist12 [cm] < %f\n",fDplusCuts[5]);
1229 printf(" sigmavert [cm] < %f\n",fDplusCuts[6]);
1230 printf(" dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1231 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1232 printf(" cosThetaPoint > %f\n",fDplusCuts[9]);
1233 printf(" Sum d0^2 [cm^2] > %f\n",fDplusCuts[10]);
1234 printf(" dca cut [cm] < %f\n",fDplusCuts[11]);
1235 printf(" Ds->KKpi cuts:\n");
1236 printf(" |M-MDs| [GeV] < %f\n",fDsCuts[0]);
1237 printf(" pTK [GeV/c] > %f\n",fDsCuts[1]);
1238 printf(" pTPi [GeV/c] > %f\n",fDsCuts[2]);
1239 printf(" |d0K| [cm] > %f\n",fDsCuts[3]);
1240 printf(" |d0Pi| [cm] > %f\n",fDsCuts[4]);
1241 printf(" dist12 [cm] < %f\n",fDsCuts[5]);
1242 printf(" sigmavert [cm] < %f\n",fDsCuts[6]);
1243 printf(" dist prim-sec [cm] > %f\n",fDsCuts[7]);
1244 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
1245 printf(" cosThetaPoint > %f\n",fDsCuts[9]);
1246 printf(" Sum d0^2 [cm^2] > %f\n",fDsCuts[10]);
1247 printf(" dca cut [cm] < %f\n",fDsCuts[11]);
1248 printf(" Inv. Mass phi/K0* [GeV] < %f\n",fDsCuts[12]);
1249 printf(" Lc->pKpi cuts:\n");
1250 printf(" |M-MLc| [GeV] < %f\n",fLcCuts[0]);
1251 printf(" pTP [GeV/c] > %f\n",fLcCuts[1]);
1252 printf(" pTPi and pTK [GeV/c] > %f\n",fLcCuts[2]);
1253 printf(" |d0P| [cm] > %f\n",fLcCuts[3]);
1254 printf(" |d0Pi| and |d0K| [cm] > %f\n",fLcCuts[4]);
1255 printf(" dist12 [cm] < %f\n",fLcCuts[5]);
1256 printf(" sigmavert [cm] < %f\n",fLcCuts[6]);
1257 printf(" dist prim-sec [cm] > %f\n",fLcCuts[7]);
1258 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
1259 printf(" cosThetaPoint > %f\n",fLcCuts[9]);
1260 printf(" Sum d0^2 [cm^2] > %f\n",fLcCuts[10]);
1261 printf(" dca cut [cm] < %f\n",fLcCuts[11]);
1262 printf(" Ds->KKpi cuts:\n");
1267 //-----------------------------------------------------------------------------
1268 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1272 Double_t *pz) const {
1273 // Check invariant mass cut
1275 Short_t dummycharge=0;
1276 Double_t *dummyd0 = new Double_t[nprongs];
1277 for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1278 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1281 UInt_t pdg2[2],pdg3[3];
1284 Bool_t retval=kFALSE;
1288 pdg2[0]=211; pdg2[1]=321;
1289 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1290 minv = rd->InvMass(nprongs,pdg2);
1291 if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1292 pdg2[0]=321; pdg2[1]=211;
1293 minv = rd->InvMass(nprongs,pdg2);
1294 if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1297 pdg2[0]=11; pdg2[1]=11;
1298 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1299 minv = rd->InvMass(nprongs,pdg2);
1300 if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1302 case 2: // D+->Kpipi
1303 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1304 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1305 minv = rd->InvMass(nprongs,pdg3);
1306 if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1308 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1309 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1310 minv = rd->InvMass(nprongs,pdg3);
1311 if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1312 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1313 minv = rd->InvMass(nprongs,pdg3);
1314 if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1316 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1317 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1318 minv = rd->InvMass(nprongs,pdg3);
1319 if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1320 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1321 minv = rd->InvMass(nprongs,pdg3);
1322 if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1332 //-----------------------------------------------------------------------------
1333 void AliAnalysisVertexingHF::SelectTracks(AliESDEvent *esd,
1334 TObjArray &trksP,Int_t &nTrksP,
1335 TObjArray &trksN,Int_t &nTrksN) const
1337 // Fill two TObjArrays with positive and negative tracks and
1338 // apply single-track preselection
1342 Int_t entries = (Int_t)esd->GetNumberOfTracks();
1344 // transfer ITS tracks from ESD to arrays and to a tree
1345 for(Int_t i=0; i<entries; i++) {
1347 AliESDtrack *esdtrack = esd->GetTrack(i);
1348 UInt_t status = esdtrack->GetStatus();
1350 // require refit in ITS
1351 if(fITSrefit && !(status&AliESDtrack::kITSrefit)) {
1352 if(fDebug) printf("track %d is not kITSrefit\n",i);
1356 // require minimum # of ITS points
1357 if(esdtrack->GetNcls(0)<fMinITSCls) {
1358 if(fDebug) printf("track %d has %d ITS cls\n",i,esdtrack->GetNcls(0));
1361 // require points on the 2 pixel layers
1363 if(!TESTBIT(esdtrack->GetITSClusterMap(),0) ||
1364 !TESTBIT(esdtrack->GetITSClusterMap(),1)) continue;
1366 // single track selection
1367 if(!SingleTrkCuts(*esdtrack)) continue;
1369 if(esdtrack->GetSign()<0) { // negative track
1370 trksN.AddLast(esdtrack);
1372 } else { // positive track
1373 trksP.AddLast(esdtrack);
1377 } // loop on ESD tracks
1381 //-----------------------------------------------------------------------------
1382 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1383 Double_t cut2,Double_t cut3,Double_t cut4,
1384 Double_t cut5,Double_t cut6,
1385 Double_t cut7,Double_t cut8)
1387 // Set the cuts for D0 selection
1388 fD0toKpiCuts[0] = cut0;
1389 fD0toKpiCuts[1] = cut1;
1390 fD0toKpiCuts[2] = cut2;
1391 fD0toKpiCuts[3] = cut3;
1392 fD0toKpiCuts[4] = cut4;
1393 fD0toKpiCuts[5] = cut5;
1394 fD0toKpiCuts[6] = cut6;
1395 fD0toKpiCuts[7] = cut7;
1396 fD0toKpiCuts[8] = cut8;
1400 //-----------------------------------------------------------------------------
1401 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
1403 // Set the cuts for D0 selection
1405 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1409 //-----------------------------------------------------------------------------
1410 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1411 Double_t cut2,Double_t cut3,Double_t cut4,
1412 Double_t cut5,Double_t cut6,
1413 Double_t cut7,Double_t cut8)
1415 // Set the cuts for J/psi from B selection
1416 fBtoJPSICuts[0] = cut0;
1417 fBtoJPSICuts[1] = cut1;
1418 fBtoJPSICuts[2] = cut2;
1419 fBtoJPSICuts[3] = cut3;
1420 fBtoJPSICuts[4] = cut4;
1421 fBtoJPSICuts[5] = cut5;
1422 fBtoJPSICuts[6] = cut6;
1423 fBtoJPSICuts[7] = cut7;
1424 fBtoJPSICuts[8] = cut8;
1428 //-----------------------------------------------------------------------------
1429 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1431 // Set the cuts for J/psi from B selection
1433 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1437 //-----------------------------------------------------------------------------
1438 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1439 Double_t cut2,Double_t cut3,Double_t cut4,
1440 Double_t cut5,Double_t cut6,
1441 Double_t cut7,Double_t cut8,
1442 Double_t cut9,Double_t cut10,Double_t cut11)
1444 // Set the cuts for Dplus->Kpipi selection
1445 fDplusCuts[0] = cut0;
1446 fDplusCuts[1] = cut1;
1447 fDplusCuts[2] = cut2;
1448 fDplusCuts[3] = cut3;
1449 fDplusCuts[4] = cut4;
1450 fDplusCuts[5] = cut5;
1451 fDplusCuts[6] = cut6;
1452 fDplusCuts[7] = cut7;
1453 fDplusCuts[8] = cut8;
1454 fDplusCuts[9] = cut9;
1455 fDplusCuts[10] = cut10;
1456 fDplusCuts[11] = cut11;
1460 //-----------------------------------------------------------------------------
1461 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1463 // Set the cuts for Dplus->Kpipi selection
1465 for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1469 //-----------------------------------------------------------------------------
1470 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
1471 Double_t cut2,Double_t cut3,Double_t cut4,
1472 Double_t cut5,Double_t cut6,
1473 Double_t cut7,Double_t cut8,
1474 Double_t cut9,Double_t cut10,
1475 Double_t cut11,Double_t cut12)
1477 // Set the cuts for Ds->KKpi selection
1488 fDsCuts[10] = cut10;
1489 fDsCuts[11] = cut11;
1490 fDsCuts[12] = cut12;
1494 //-----------------------------------------------------------------------------
1495 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13])
1497 // Set the cuts for Ds->KKpi selection
1499 for(Int_t i=0; i<13; i++) fDsCuts[i] = cuts[i];
1503 //-----------------------------------------------------------------------------
1504 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
1505 Double_t cut2,Double_t cut3,Double_t cut4,
1506 Double_t cut5,Double_t cut6,
1507 Double_t cut7,Double_t cut8,
1508 Double_t cut9,Double_t cut10,Double_t cut11)
1510 // Set the cuts for Lc->pKpi selection
1521 fLcCuts[10] = cut10;
1522 fLcCuts[11] = cut11;
1526 //-----------------------------------------------------------------------------
1527 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
1529 // Set the cuts for Lc->pKpi selection
1531 for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
1535 //-----------------------------------------------------------------------------
1536 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk) const
1538 // Check if track passes some kinematical cuts
1540 if(TMath::Abs(trk.Pt()) < fMinPtCut) {
1541 //printf("pt %f\n",1./trk.GetParameter()[4]);
1544 trk.RelateToVertex(fV1,fBzkG,10.);
1545 Float_t d0z0[2],covd0z0[3];
1546 trk.GetImpactParameters(d0z0,covd0z0);
1547 if(TMath::Abs(d0z0[0]) < fMind0rphiCut) {
1548 printf("d0rphi %f\n",TMath::Abs(d0z0[0]));
1554 //-----------------------------------------------------------------------------
1555 AliESDVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray) const
1557 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1559 AliESDVertex *vertex = 0;
1561 if(!fSecVtxWithKF) { // AliVertexerTracks
1563 AliVertexerTracks *vertexer2 = new AliVertexerTracks(fBzkG);
1564 vertexer2->SetDebug(0);
1565 vertexer2->SetVtxStart(fV1);
1566 vertex = (AliESDVertex*)vertexer2->VertexForSelectedESDTracks(trkArray);
1569 if(vertex->GetNContributors()!=trkArray->GetEntriesFast()) {
1570 if(fDebug) printf("vertexing failed\n");
1571 delete vertex; vertex=0;
1574 } else { // Kalman Filter vertexer (AliKFParticle)
1576 AliKFParticle::SetField(fBzkG);
1578 AliKFParticle vertexKF;
1580 Int_t nTrks = trkArray->GetEntriesFast();
1581 for(Int_t i=0; i<nTrks; i++) {
1582 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1583 AliKFParticle daughterKF(*esdTrack,211);
1584 vertexKF.AddDaughter(daughterKF);
1586 vertex = new AliESDVertex();
1587 vertexKF.CopyToESDVertex(*vertex);
1593 //-----------------------------------------------------------------------------