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