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