]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG3/vertexingHF/AliAnalysisVertexingHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisVertexingHF.cxx
... / ...
CommitLineData
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/* $Id$ */
17
18//----------------------------------------------------------------------------
19// Implementation of the heavy-flavour vertexing analysis class
20// Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
21// To be used as a task of AliAnalysisManager by means of the interface
22// class AliAnalysisTaskSEVertexingHF.
23// An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
24//
25// Contact: andrea.dainese@pd.infn.it
26// Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
27// F.Prino, R.Romita, X.M.Zhang
28//----------------------------------------------------------------------------
29#include <Riostream.h>
30#include <TFile.h>
31#include <TDatabasePDG.h>
32#include <TString.h>
33#include <TList.h>
34#include "AliLog.h"
35#include "AliVEvent.h"
36#include "AliVVertex.h"
37#include "AliVTrack.h"
38#include "AliVertexerTracks.h"
39#include "AliKFVertex.h"
40#include "AliESDEvent.h"
41#include "AliESDVertex.h"
42#include "AliExternalTrackParam.h"
43#include "AliNeutralTrackParam.h"
44#include "AliESDtrack.h"
45#include "AliESDtrackCuts.h"
46#include "AliAODEvent.h"
47#include "AliAODRecoDecay.h"
48#include "AliAODRecoDecayHF.h"
49#include "AliAODRecoDecayHF2Prong.h"
50#include "AliAODRecoDecayHF3Prong.h"
51#include "AliAODRecoDecayHF4Prong.h"
52#include "AliAODRecoCascadeHF.h"
53#include "AliRDHFCutsD0toKpi.h"
54#include "AliRDHFCutsJpsitoee.h"
55#include "AliRDHFCutsDplustoKpipi.h"
56#include "AliRDHFCutsDstoKKpi.h"
57#include "AliRDHFCutsLctopKpi.h"
58#include "AliRDHFCutsLctoV0.h"
59#include "AliRDHFCutsD0toKpipipi.h"
60#include "AliRDHFCutsDStartoKpipi.h"
61#include "AliAnalysisFilter.h"
62#include "AliAnalysisVertexingHF.h"
63#include "AliMixedEvent.h"
64#include "AliESDv0.h"
65#include "AliAODv0.h"
66
67ClassImp(AliAnalysisVertexingHF)
68
69//----------------------------------------------------------------------------
70AliAnalysisVertexingHF::AliAnalysisVertexingHF():
71fInputAOD(kFALSE),
72fAODMapSize(0),
73fAODMap(0),
74fBzkG(0.),
75fSecVtxWithKF(kFALSE),
76fRecoPrimVtxSkippingTrks(kFALSE),
77fRmTrksFromPrimVtx(kFALSE),
78fV1(0x0),
79fD0toKpi(kTRUE),
80fJPSItoEle(kTRUE),
81f3Prong(kTRUE),
82f4Prong(kTRUE),
83fDstar(kTRUE),
84fCascades(kTRUE),
85fLikeSign(kFALSE),
86fLikeSign3prong(kFALSE),
87fMixEvent(kFALSE),
88fTrackFilter(0x0),
89fTrackFilterSoftPi(0x0),
90fCutsD0toKpi(0x0),
91fCutsJpsitoee(0x0),
92fCutsDplustoKpipi(0x0),
93fCutsDstoKKpi(0x0),
94fCutsLctopKpi(0x0),
95fCutsLctoV0(0x0),
96fCutsD0toKpipipi(0x0),
97fCutsDStartoKpipi(0x0),
98fListOfCuts(0x0),
99fFindVertexForDstar(kTRUE),
100fFindVertexForCascades(kTRUE),
101fMassCutBeforeVertexing(kFALSE),
102fMassCalc2(0),
103fMassCalc3(0),
104fMassCalc4(0),
105fnTrksTotal(0),
106fnSeleTrksTotal(0)
107{
108 // Default constructor
109
110 Double_t d02[2]={0.,0.};
111 Double_t d03[3]={0.,0.,0.};
112 Double_t d04[4]={0.,0.,0.,0.};
113 fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
114 fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
115 fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
116
117}
118//--------------------------------------------------------------------------
119AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
120TNamed(source),
121fInputAOD(source.fInputAOD),
122fAODMapSize(source.fAODMapSize),
123fAODMap(source.fAODMap),
124fBzkG(source.fBzkG),
125fSecVtxWithKF(source.fSecVtxWithKF),
126fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
127fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
128fV1(source.fV1),
129fD0toKpi(source.fD0toKpi),
130fJPSItoEle(source.fJPSItoEle),
131f3Prong(source.f3Prong),
132f4Prong(source.f4Prong),
133fDstar(source.fDstar),
134fCascades(source.fCascades),
135fLikeSign(source.fLikeSign),
136fLikeSign3prong(source.fLikeSign3prong),
137fMixEvent(source.fMixEvent),
138fTrackFilter(source.fTrackFilter),
139fTrackFilterSoftPi(source.fTrackFilterSoftPi),
140fCutsD0toKpi(source.fCutsD0toKpi),
141fCutsJpsitoee(source.fCutsJpsitoee),
142fCutsDplustoKpipi(source.fCutsDplustoKpipi),
143fCutsDstoKKpi(source.fCutsDstoKKpi),
144fCutsLctopKpi(source.fCutsLctopKpi),
145fCutsLctoV0(source.fCutsLctoV0),
146fCutsD0toKpipipi(source.fCutsD0toKpipipi),
147fCutsDStartoKpipi(source.fCutsDStartoKpipi),
148fListOfCuts(source.fListOfCuts),
149fFindVertexForDstar(source.fFindVertexForDstar),
150fFindVertexForCascades(source.fFindVertexForCascades),
151fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
152fMassCalc2(source.fMassCalc2),
153fMassCalc3(source.fMassCalc3),
154fMassCalc4(source.fMassCalc4),
155fnTrksTotal(0),
156fnSeleTrksTotal(0)
157{
158 //
159 // Copy constructor
160 //
161}
162//--------------------------------------------------------------------------
163AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
164{
165 //
166 // assignment operator
167 //
168 if(&source == this) return *this;
169 fInputAOD = source.fInputAOD;
170 fAODMapSize = source.fAODMapSize;
171 fBzkG = source.fBzkG;
172 fSecVtxWithKF = source.fSecVtxWithKF;
173 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
174 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
175 fV1 = source.fV1;
176 fD0toKpi = source.fD0toKpi;
177 fJPSItoEle = source.fJPSItoEle;
178 f3Prong = source.f3Prong;
179 f4Prong = source.f4Prong;
180 fDstar = source.fDstar;
181 fCascades = source.fCascades;
182 fLikeSign = source.fLikeSign;
183 fLikeSign3prong = source.fLikeSign3prong;
184 fMixEvent = source.fMixEvent;
185 fTrackFilter = source.fTrackFilter;
186 fTrackFilterSoftPi = source.fTrackFilterSoftPi;
187 fCutsD0toKpi = source.fCutsD0toKpi;
188 fCutsJpsitoee = source.fCutsJpsitoee;
189 fCutsDplustoKpipi = source.fCutsDplustoKpipi;
190 fCutsDstoKKpi = source.fCutsDstoKKpi;
191 fCutsLctopKpi = source.fCutsLctopKpi;
192 fCutsLctoV0 = source.fCutsLctoV0;
193 fCutsD0toKpipipi = source.fCutsD0toKpipipi;
194 fCutsDStartoKpipi = source.fCutsDStartoKpipi;
195 fListOfCuts = source.fListOfCuts;
196 fFindVertexForDstar = source.fFindVertexForDstar;
197 fFindVertexForCascades = source.fFindVertexForCascades;
198 fMassCutBeforeVertexing = source.fMassCutBeforeVertexing;
199 fMassCalc2 = source.fMassCalc2;
200 fMassCalc3 = source.fMassCalc3;
201 fMassCalc4 = source.fMassCalc4;
202
203 return *this;
204}
205//----------------------------------------------------------------------------
206AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
207 // Destructor
208 if(fV1) { delete fV1; fV1=0; }
209 if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
210 if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
211 if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
212 if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
213 if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
214 if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
215 if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
216 if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
217 if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
218 if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
219 if(fAODMap) { delete [] fAODMap; fAODMap=0; }
220 if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
221 if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
222 if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
223}
224//----------------------------------------------------------------------------
225TList *AliAnalysisVertexingHF::FillListOfCuts() {
226 // Fill list of analysis cuts
227
228 TList *list = new TList();
229 list->SetOwner();
230 list->SetName("ListOfCuts");
231
232 if(fCutsD0toKpi) {
233 AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
234 list->Add(cutsD0toKpi);
235 }
236 if(fCutsJpsitoee) {
237 AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
238 list->Add(cutsJpsitoee);
239 }
240 if(fCutsDplustoKpipi) {
241 AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
242 list->Add(cutsDplustoKpipi);
243 }
244 if(fCutsDstoKKpi) {
245 AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
246 list->Add(cutsDstoKKpi);
247 }
248 if(fCutsLctopKpi) {
249 AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
250 list->Add(cutsLctopKpi);
251 }
252 if(fCutsLctoV0){
253 AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
254 list->Add(cutsLctoV0);
255 }
256 if(fCutsD0toKpipipi) {
257 AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
258 list->Add(cutsD0toKpipipi);
259 }
260 if(fCutsDStartoKpipi) {
261 AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
262 list->Add(cutsDStartoKpipi);
263 }
264
265 // keep a pointer to the list
266 fListOfCuts = list;
267
268 return list;
269}
270//----------------------------------------------------------------------------
271void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
272 TClonesArray *aodVerticesHFTClArr,
273 TClonesArray *aodD0toKpiTClArr,
274 TClonesArray *aodJPSItoEleTClArr,
275 TClonesArray *aodCharm3ProngTClArr,
276 TClonesArray *aodCharm4ProngTClArr,
277 TClonesArray *aodDstarTClArr,
278 TClonesArray *aodCascadesTClArr,
279 TClonesArray *aodLikeSign2ProngTClArr,
280 TClonesArray *aodLikeSign3ProngTClArr)
281{
282 // Find heavy-flavour vertex candidates
283 // Input: ESD or AOD
284 // Output: AOD (additional branches added)
285
286 if(!fMixEvent){
287 TString evtype = event->IsA()->GetName();
288 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
289 } // if we do mixing AliVEvent is a AliMixedEvent
290
291 if(fInputAOD) {
292 AliDebug(2,"Creating HF candidates from AOD");
293 } else {
294 AliDebug(2,"Creating HF candidates from ESD");
295 }
296
297 if(!aodVerticesHFTClArr) {
298 printf("ERROR: no aodVerticesHFTClArr");
299 return;
300 }
301 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
302 printf("ERROR: no aodD0toKpiTClArr");
303 return;
304 }
305 if(fJPSItoEle && !aodJPSItoEleTClArr) {
306 printf("ERROR: no aodJPSItoEleTClArr");
307 return;
308 }
309 if(f3Prong && !aodCharm3ProngTClArr) {
310 printf("ERROR: no aodCharm3ProngTClArr");
311 return;
312 }
313 if(f4Prong && !aodCharm4ProngTClArr) {
314 printf("ERROR: no aodCharm4ProngTClArr");
315 return;
316 }
317 if(fDstar && !aodDstarTClArr) {
318 printf("ERROR: no aodDstarTClArr");
319 return;
320 }
321 if(fCascades && !aodCascadesTClArr){
322 printf("ERROR: no aodCascadesTClArr ");
323 return;
324 }
325 if(fLikeSign && !aodLikeSign2ProngTClArr) {
326 printf("ERROR: no aodLikeSign2ProngTClArr");
327 return;
328 }
329 if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
330 printf("ERROR: no aodLikeSign3ProngTClArr");
331 return;
332 }
333
334 // delete candidates from previous event and create references
335 Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
336 aodVerticesHFTClArr->Delete();
337 iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
338 TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
339 if(fD0toKpi || fDstar) {
340 aodD0toKpiTClArr->Delete();
341 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
342 }
343 if(fJPSItoEle) {
344 aodJPSItoEleTClArr->Delete();
345 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
346 }
347 if(f3Prong) {
348 aodCharm3ProngTClArr->Delete();
349 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
350 }
351 if(f4Prong) {
352 aodCharm4ProngTClArr->Delete();
353 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
354 }
355 if(fDstar) {
356 aodDstarTClArr->Delete();
357 iDstar = aodDstarTClArr->GetEntriesFast();
358 }
359 if(fCascades) {
360 aodCascadesTClArr->Delete();
361 iCascades = aodCascadesTClArr->GetEntriesFast();
362 }
363 if(fLikeSign) {
364 aodLikeSign2ProngTClArr->Delete();
365 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
366 }
367 if(fLikeSign3prong && f3Prong) {
368 aodLikeSign3ProngTClArr->Delete();
369 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
370 }
371
372 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
373 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
374 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
375 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
376 TClonesArray &aodDstarRef = *aodDstarTClArr;
377 TClonesArray &aodCascadesRef = *aodCascadesTClArr;
378 TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
379 TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
380
381
382 AliAODRecoDecayHF2Prong *io2Prong = 0;
383 AliAODRecoDecayHF3Prong *io3Prong = 0;
384 AliAODRecoDecayHF4Prong *io4Prong = 0;
385 AliAODRecoCascadeHF *ioCascade = 0;
386
387 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
388 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
389 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
390 Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
391 Bool_t okCascades=kFALSE;
392 AliESDtrack *postrack1 = 0;
393 AliESDtrack *postrack2 = 0;
394 AliESDtrack *negtrack1 = 0;
395 AliESDtrack *negtrack2 = 0;
396 AliESDtrack *trackPi = 0;
397 // AliESDtrack *posV0track = 0;
398 // AliESDtrack *negV0track = 0;
399 Float_t dcaMax = fCutsD0toKpi->GetDCACut();
400 if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
401 if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
402 if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
403 if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
404 if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
405 if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
406
407 AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
408
409
410 // get Bz
411 fBzkG = (Double_t)event->GetMagneticField();
412
413 trkEntries = (Int_t)event->GetNumberOfTracks();
414 AliDebug(1,Form(" Number of tracks: %d",trkEntries));
415 fnTrksTotal += trkEntries;
416
417 nv0 = (Int_t)event->GetNumberOfV0s();
418 AliDebug(1,Form(" Number of V0s: %d",nv0));
419
420 if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
421 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
422 return;
423 }
424
425 // event selection
426 if(!fCutsD0toKpi->IsEventSelected(event)) return;
427
428 // call function that applies sigle-track selection,
429 // for displaced tracks and soft pions (both charges) for D*,
430 // and retrieves primary vertex
431 TObjArray seleTrksArray(trkEntries);
432 TObjArray tracksAtVertex(trkEntries);
433 UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
434 Int_t nSeleTrks=0;
435 Int_t *evtNumber = new Int_t[trkEntries];
436 SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
437
438 AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
439 fnSeleTrksTotal += nSeleTrks;
440
441
442 TObjArray *twoTrackArray1 = new TObjArray(2);
443 TObjArray *twoTrackArray2 = new TObjArray(2);
444 TObjArray *twoTrackArrayV0 = new TObjArray(2);
445 TObjArray *twoTrackArrayCasc = new TObjArray(2);
446 TObjArray *threeTrackArray = new TObjArray(3);
447 TObjArray *fourTrackArray = new TObjArray(4);
448
449 Double_t dispersion;
450 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
451
452 AliAODRecoDecayHF *rd = 0;
453 AliAODRecoCascadeHF *rc = 0;
454 AliAODv0 *v0 = 0;
455 AliESDv0 *esdV0 = 0;
456
457 Bool_t massCutOK=kTRUE;
458
459 // LOOP ON POSITIVE TRACKS
460 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
461
462 //if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
463 //if(iTrkP1%1==0) printf(" 1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);
464
465 // get track from tracks array
466 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
467
468 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
469
470 // Make cascades with V0+track
471 //
472 if(fCascades) {
473 // loop on V0's
474 for(iv0=0; iv0<nv0; iv0++){
475
476 //AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
477
478 // Get the V0
479 if(fInputAOD) {
480 v0 = ((AliAODEvent*)event)->GetV0(iv0);
481 } else {
482 esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
483 }
484 if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
485 (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
486
487
488 // Get the tracks that form the V0
489 // ( parameters at primary vertex )
490 // and define an AliExternalTrackParam out of them
491 AliExternalTrackParam * posV0track;
492 AliExternalTrackParam * negV0track;
493
494 if(fInputAOD){
495 AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
496 AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
497 if( !posVV0track || !negVV0track ) continue;
498 //
499 // Apply some basic V0 daughter criteria
500 //
501 // bachelor must not be a v0-track
502 if (posVV0track->GetID() == postrack1->GetID() ||
503 negVV0track->GetID() == postrack1->GetID()) continue;
504 // reject like-sign v0
505 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
506 // avoid ghost TPC tracks
507 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
508 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
509 // Get AliExternalTrackParam out of the AliAODTracks
510 Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
511 posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
512 posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
513 posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
514 negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
515 negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
516 negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
517 } else {
518 AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
519 AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
520 if( !posVV0track || !negVV0track ) continue;
521 //
522 // Apply some basic V0 daughter criteria
523 //
524 // bachelor must not be a v0-track
525 if (posVV0track->GetID() == postrack1->GetID() ||
526 negVV0track->GetID() == postrack1->GetID()) continue;
527 // reject like-sign v0
528 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
529 // avoid ghost TPC tracks
530 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
531 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
532 // reject kinks (only necessary on AliESDtracks)
533 if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
534 // Get AliExternalTrackParam out of the AliESDtracks
535 posV0track = new AliExternalTrackParam(*posVV0track);
536 negV0track = new AliExternalTrackParam(*negVV0track);
537
538 // Define the AODv0 from ESDv0 if reading ESDs
539 v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
540 }
541 if( !posV0track || !negV0track ){
542 AliDebug(1,Form(" Couldn't get the V0 daughters"));
543 continue;
544 }
545
546 // fill in the v0 two-external-track-param array
547 twoTrackArrayV0->AddAt(posV0track,0);
548 twoTrackArrayV0->AddAt(negV0track,1);
549
550 // Get the V0 dca
551 dcaV0 = v0->DcaV0Daughters();
552
553 // Define the V0 (neutral) track
554 AliNeutralTrackParam *trackV0;
555 if(fInputAOD) {
556 const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
557 if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
558 } else {
559 Double_t xyz[3], pxpypz[3];
560 esdV0->XvYvZv(xyz);
561 esdV0->PxPyPz(pxpypz);
562 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
563 trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
564 }
565 // Fill in the object array to create the cascade
566 twoTrackArrayCasc->AddAt(postrack1,0);
567 twoTrackArrayCasc->AddAt(trackV0,1);
568 // Compute the cascade vertex
569 AliAODVertex *vertexCasc = 0;
570 if(fFindVertexForCascades) {
571 // DCA between the two tracks
572 dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
573 // Vertexing+
574 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
575 } else {
576 // assume Cascade decays at the primary vertex
577 Double_t pos[3],cov[6],chi2perNDF;
578 fV1->GetXYZ(pos);
579 fV1->GetCovMatrix(cov);
580 chi2perNDF = fV1->GetChi2toNDF();
581 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
582 dcaCasc = 0.;
583 }
584 if(!vertexCasc) {
585 delete posV0track; posV0track=NULL;
586 delete negV0track; negV0track=NULL;
587 delete trackV0; trackV0=NULL;
588 if(!fInputAOD) {delete v0; v0=NULL;}
589 twoTrackArrayV0->Clear();
590 twoTrackArrayCasc->Clear();
591 continue;
592 }
593
594 // Create and store the Cascade if passed the cuts
595 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
596 if(okCascades && ioCascade) {
597 //AliDebug(1,Form("Storing a cascade object... "));
598 // add the vertex and the cascade to the AOD
599 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
600 rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
601 rc->SetSecondaryVtx(vCasc);
602 vCasc->SetParent(rc);
603 rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
604 if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
605 AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
606 vCasc->AddDaughter(v0); // fill the 2prong V0
607 }
608
609 // Clean up
610 delete posV0track; posV0track=NULL;
611 delete negV0track; negV0track=NULL;
612 delete trackV0; trackV0=NULL;
613 twoTrackArrayV0->Clear();
614 twoTrackArrayCasc->Clear();
615 if(ioCascade) { delete ioCascade; ioCascade=NULL; }
616 if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
617 if(!fInputAOD) {delete v0; v0=NULL;}
618
619 } // end loop on V0's
620 }
621
622 // If there is less than 2 particles continue
623 if(trkEntries<2) {
624 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
625 continue;
626 }
627
628 if(postrack1->Charge()<0 && !fLikeSign) continue;
629
630 // LOOP ON NEGATIVE TRACKS
631 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
632
633 //if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
634 //if(iTrkN1%1==0) printf(" 1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);
635
636 if(iTrkN1==iTrkP1) continue;
637
638 // get track from tracks array
639 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
640
641 if(negtrack1->Charge()>0 && !fLikeSign) continue;
642
643 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
644
645 if(fMixEvent) {
646 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
647 }
648
649 if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
650 isLikeSign2Prong=kTRUE;
651 if(!fLikeSign) continue;
652 if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
653 } else { // unlike-sign
654 isLikeSign2Prong=kFALSE;
655 if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
656 if(fMixEvent) {
657 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
658 }
659
660 }
661
662 // back to primary vertex
663 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
664 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
665 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
666 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
667
668 // DCA between the two tracks
669 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
670 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
671
672 // Vertexing
673 twoTrackArray1->AddAt(postrack1,0);
674 twoTrackArray1->AddAt(negtrack1,1);
675 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
676 if(!vertexp1n1) {
677 twoTrackArray1->Clear();
678 negtrack1=0;
679 continue;
680 }
681
682 // 2 prong candidate
683 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
684
685 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
686
687 if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
688 // add the vertex and the decay to the AOD
689 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
690 if(!isLikeSign2Prong) {
691 if(okD0) {
692 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
693 rd->SetSecondaryVtx(v2Prong);
694 v2Prong->SetParent(rd);
695 AddRefs(v2Prong,rd,event,twoTrackArray1);
696 }
697 if(okJPSI) {
698 rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
699 rd->SetSecondaryVtx(v2Prong);
700 if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
701 AddRefs(v2Prong,rd,event,twoTrackArray1);
702 }
703 } else { // isLikeSign2Prong
704 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
705 rd->SetSecondaryVtx(v2Prong);
706 v2Prong->SetParent(rd);
707 AddRefs(v2Prong,rd,event,twoTrackArray1);
708 }
709 // Set selection bit for PID
710 if(okD0) SetSelectionBitForPID(fCutsD0toKpi,rd,AliRDHFCuts::kD0toKpiPID);
711 }
712 // D* candidates
713 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
714 // write references in io2Prong
715 if(fInputAOD) {
716 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
717 } else {
718 vertexp1n1->AddDaughter(postrack1);
719 vertexp1n1->AddDaughter(negtrack1);
720 }
721 io2Prong->SetSecondaryVtx(vertexp1n1);
722 //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
723 // create a track from the D0
724 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
725
726 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
727 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
728
729 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
730
731 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
732
733 if(fMixEvent) {
734 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
735 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
736 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
737 }
738
739 //if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
740
741 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
742 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
743
744 // get track from tracks array
745 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
746 // trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
747 SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
748 twoTrackArrayCasc->AddAt(trackPi,0);
749 twoTrackArrayCasc->AddAt(trackD0,1);
750
751 AliAODVertex *vertexCasc = 0;
752
753 if(fFindVertexForDstar) {
754 // DCA between the two tracks
755 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
756 // Vertexing
757 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
758 } else {
759 // assume Dstar decays at the primary vertex
760 Double_t pos[3],cov[6],chi2perNDF;
761 fV1->GetXYZ(pos);
762 fV1->GetCovMatrix(cov);
763 chi2perNDF = fV1->GetChi2toNDF();
764 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
765 dcaCasc = 0.;
766 }
767 if(!vertexCasc) {
768 twoTrackArrayCasc->Clear();
769 trackPi=0;
770 continue;
771 }
772
773 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
774 if(okDstar) {
775 // add the D0 to the AOD (if not already done)
776 if(!okD0) {
777 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
778 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
779 rd->SetSecondaryVtx(v2Prong);
780 v2Prong->SetParent(rd);
781 AddRefs(v2Prong,rd,event,twoTrackArray1);
782 okD0=kTRUE; // this is done to add it only once
783 }
784 // add the vertex and the cascade to the AOD
785 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
786 rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
787 rc->SetSecondaryVtx(vCasc);
788 vCasc->SetParent(rc);
789 if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
790 AddRefs(vCasc,rc,event,twoTrackArrayCasc);
791 vCasc->AddDaughter(rd); // add the D0 (in ref #1)
792 // Set selection bit for PID
793 SetSelectionBitForPID(fCutsDStartoKpipi,rc,AliRDHFCuts::kDstarPID);
794 }
795 twoTrackArrayCasc->Clear();
796 trackPi=0;
797 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
798 delete vertexCasc; vertexCasc=NULL;
799 } // end loop on soft pi tracks
800
801 if(trackD0) {delete trackD0; trackD0=NULL;}
802
803 }
804 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
805 }
806
807 twoTrackArray1->Clear();
808 if( (!f3Prong && !f4Prong) ||
809 (isLikeSign2Prong && !f3Prong) ) {
810 negtrack1=0;
811 delete vertexp1n1;
812 continue;
813 }
814
815
816 // 2nd LOOP ON POSITIVE TRACKS
817 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
818
819 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
820
821 //if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
822
823 // get track from tracks array
824 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
825
826 if(postrack2->Charge()<0) continue;
827
828 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
829
830 if(fMixEvent) {
831 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
832 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
833 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
834 }
835
836 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
837 if(!fLikeSign3prong) continue;
838 if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
839 isLikeSign3Prong=kTRUE;
840 } else { // not ok
841 continue;
842 }
843 } else { // normal triplet (+-+)
844 isLikeSign3Prong=kFALSE;
845 if(fMixEvent) {
846 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
847 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
848 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
849 }
850 }
851
852 // back to primary vertex
853 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
854 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
855 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
856 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
857 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
858 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
859
860 //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
861
862 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
863 if(dcap2n1>dcaMax) { postrack2=0; continue; }
864 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
865 if(dcap1p2>dcaMax) { postrack2=0; continue; }
866
867 // check invariant mass cuts for D+,Ds,Lc
868 massCutOK=kTRUE;
869 if(f3Prong) {
870 if(postrack2->Charge()>0) {
871 threeTrackArray->AddAt(postrack1,0);
872 threeTrackArray->AddAt(negtrack1,1);
873 threeTrackArray->AddAt(postrack2,2);
874 } else {
875 threeTrackArray->AddAt(negtrack1,0);
876 threeTrackArray->AddAt(postrack1,1);
877 threeTrackArray->AddAt(postrack2,2);
878 }
879 if(fMassCutBeforeVertexing)
880 massCutOK = SelectInvMassAndPt(threeTrackArray);
881 }
882
883 if(f3Prong && !massCutOK) {
884 threeTrackArray->Clear();
885 if(!f4Prong) {
886 postrack2=0;
887 continue;
888 }
889 }
890
891 // Vertexing
892 twoTrackArray2->AddAt(postrack2,0);
893 twoTrackArray2->AddAt(negtrack1,1);
894 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
895 if(!vertexp2n1) {
896 twoTrackArray2->Clear();
897 postrack2=0;
898 continue;
899 }
900
901 // 3 prong candidates
902 if(f3Prong && massCutOK) {
903
904 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
905 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
906 if(ok3Prong) {
907 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
908 if(!isLikeSign3Prong) {
909 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
910 rd->SetSecondaryVtx(v3Prong);
911 v3Prong->SetParent(rd);
912 AddRefs(v3Prong,rd,event,threeTrackArray);
913 } else { // isLikeSign3Prong
914 if(fLikeSign3prong){
915 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
916 rd->SetSecondaryVtx(v3Prong);
917 v3Prong->SetParent(rd);
918 AddRefs(v3Prong,rd,event,threeTrackArray);
919 }
920 }
921 // Set selection bit for PID
922 SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
923 SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
924 SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
925 }
926 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
927 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
928 }
929
930 // 4 prong candidates
931 if(f4Prong
932 // don't make 4 prong with like-sign pairs and triplets
933 && !isLikeSign2Prong && !isLikeSign3Prong
934 // track-to-track dca cuts already now
935 && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
936 && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
937
938 // back to primary vertex
939 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
940 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
941 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
942 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
943 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
944 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
945
946 // Vertexing for these 3 (can be taken from above?)
947 threeTrackArray->AddAt(postrack1,0);
948 threeTrackArray->AddAt(negtrack1,1);
949 threeTrackArray->AddAt(postrack2,2);
950 AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
951
952 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
953 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
954
955 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
956
957 //if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
958
959 // get track from tracks array
960 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
961
962 if(negtrack2->Charge()>0) continue;
963
964 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
965 if(fMixEvent){
966 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
967 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
968 evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
969 evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
970 evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
971 evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
972 }
973
974 // back to primary vertex
975 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
976 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
977 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
978 // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
979 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
980 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
981 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
982 SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
983
984 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
985 if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
986 dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
987 if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
988
989
990 fourTrackArray->AddAt(postrack1,0);
991 fourTrackArray->AddAt(negtrack1,1);
992 fourTrackArray->AddAt(postrack2,2);
993 fourTrackArray->AddAt(negtrack2,3);
994
995 // check invariant mass cuts for D0
996 massCutOK=kTRUE;
997 if(fMassCutBeforeVertexing)
998 massCutOK = SelectInvMassAndPt(fourTrackArray);
999
1000 if(!massCutOK) {
1001 fourTrackArray->Clear();
1002 negtrack2=0;
1003 continue;
1004 }
1005
1006 // Vertexing
1007 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
1008 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
1009 if(ok4Prong) {
1010 AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
1011 rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
1012 rd->SetSecondaryVtx(v4Prong);
1013 v4Prong->SetParent(rd);
1014 AddRefs(v4Prong,rd,event,fourTrackArray);
1015 }
1016
1017 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
1018 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
1019 fourTrackArray->Clear();
1020 negtrack2 = 0;
1021
1022 } // end loop on negative tracks
1023
1024 threeTrackArray->Clear();
1025 delete vertexp1n1p2;
1026
1027 }
1028
1029 postrack2 = 0;
1030 delete vertexp2n1;
1031
1032 } // end 2nd loop on positive tracks
1033
1034 twoTrackArray2->Clear();
1035
1036 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
1037 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1038
1039 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1040
1041 //if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1042
1043 // get track from tracks array
1044 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1045
1046 if(negtrack2->Charge()>0) continue;
1047
1048 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1049
1050 if(fMixEvent) {
1051 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1052 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1053 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1054 }
1055
1056 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1057 if(!fLikeSign3prong) continue;
1058 if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1059 isLikeSign3Prong=kTRUE;
1060 } else { // not ok
1061 continue;
1062 }
1063 } else { // normal triplet (-+-)
1064 isLikeSign3Prong=kFALSE;
1065 }
1066
1067 // back to primary vertex
1068 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1069 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1070 // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1071 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1072 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1073 SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1074 //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1075
1076 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1077 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1078 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1079 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1080
1081 threeTrackArray->AddAt(negtrack1,0);
1082 threeTrackArray->AddAt(postrack1,1);
1083 threeTrackArray->AddAt(negtrack2,2);
1084
1085 // check invariant mass cuts for D+,Ds,Lc
1086 massCutOK=kTRUE;
1087 if(fMassCutBeforeVertexing && f3Prong)
1088 massCutOK = SelectInvMassAndPt(threeTrackArray);
1089
1090 if(!massCutOK) {
1091 threeTrackArray->Clear();
1092 negtrack2=0;
1093 continue;
1094 }
1095
1096 // Vertexing
1097 twoTrackArray2->AddAt(postrack1,0);
1098 twoTrackArray2->AddAt(negtrack2,1);
1099
1100 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1101 if(!vertexp1n2) {
1102 twoTrackArray2->Clear();
1103 negtrack2=0;
1104 continue;
1105 }
1106
1107 if(f3Prong) {
1108 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1109 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1110 if(ok3Prong) {
1111 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1112 if(!isLikeSign3Prong) {
1113 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1114 rd->SetSecondaryVtx(v3Prong);
1115 v3Prong->SetParent(rd);
1116 AddRefs(v3Prong,rd,event,threeTrackArray);
1117 } else { // isLikeSign3Prong
1118 if(fLikeSign3prong){
1119 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1120 rd->SetSecondaryVtx(v3Prong);
1121 v3Prong->SetParent(rd);
1122 AddRefs(v3Prong,rd,event,threeTrackArray);
1123 }
1124 }
1125 // Set selection bit for PID
1126 SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
1127 SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
1128 SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
1129 }
1130 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1131 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1132 }
1133 threeTrackArray->Clear();
1134 negtrack2 = 0;
1135 delete vertexp1n2;
1136
1137 } // end 2nd loop on negative tracks
1138
1139 twoTrackArray2->Clear();
1140
1141 negtrack1 = 0;
1142 delete vertexp1n1;
1143 } // end 1st loop on negative tracks
1144
1145 postrack1 = 0;
1146 } // end 1st loop on positive tracks
1147
1148
1149 AliDebug(1,Form(" Total HF vertices in event = %d;",
1150 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1151 if(fD0toKpi) {
1152 AliDebug(1,Form(" D0->Kpi in event = %d;",
1153 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1154 }
1155 if(fJPSItoEle) {
1156 AliDebug(1,Form(" JPSI->ee in event = %d;",
1157 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1158 }
1159 if(f3Prong) {
1160 AliDebug(1,Form(" Charm->3Prong in event = %d;",
1161 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1162 }
1163 if(f4Prong) {
1164 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1165 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1166 }
1167 if(fDstar) {
1168 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1169 (Int_t)aodDstarTClArr->GetEntriesFast()));
1170 }
1171 if(fCascades){
1172 AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1173 (Int_t)aodCascadesTClArr->GetEntriesFast()));
1174 }
1175 if(fLikeSign) {
1176 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1177 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1178 }
1179 if(fLikeSign3prong && f3Prong) {
1180 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1181 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1182 }
1183
1184
1185 twoTrackArray1->Delete(); delete twoTrackArray1;
1186 twoTrackArray2->Delete(); delete twoTrackArray2;
1187 twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1188 twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1189 threeTrackArray->Clear();
1190 threeTrackArray->Delete(); delete threeTrackArray;
1191 fourTrackArray->Delete(); delete fourTrackArray;
1192 delete [] seleFlags; seleFlags=NULL;
1193 if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1194 tracksAtVertex.Delete();
1195
1196 if(fInputAOD) {
1197 seleTrksArray.Delete();
1198 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1199 }
1200
1201
1202 //printf("Trks: total %d sele %d\n",fnTrksTotal,fnSeleTrksTotal);
1203
1204 return;
1205}
1206//----------------------------------------------------------------------------
1207void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1208 const AliVEvent *event,
1209 const TObjArray *trkArray) const
1210{
1211 // Add the AOD tracks as daughters of the vertex (TRef)
1212 // Also add the references to the primary vertex and to the cuts
1213
1214 if(fInputAOD) {
1215 AddDaughterRefs(v,event,trkArray);
1216 rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1217 }
1218
1219 /*
1220 rd->SetListOfCutsRef((TList*)fListOfCuts);
1221 //fListOfCuts->Print();
1222 cout<<fListOfCuts<<endl;
1223 TList *l=(TList*)rd->GetListOfCuts();
1224 cout<<l<<endl;
1225 if(l) {l->Print(); }else{printf("error\n");}
1226 */
1227
1228 return;
1229}
1230//----------------------------------------------------------------------------
1231void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1232 const AliVEvent *event,
1233 const TObjArray *trkArray) const
1234{
1235 // Add the AOD tracks as daughters of the vertex (TRef)
1236
1237 Int_t nDg = v->GetNDaughters();
1238 TObject *dg = 0;
1239 if(nDg) dg = v->GetDaughter(0);
1240
1241 if(dg) return; // daughters already added
1242
1243 Int_t nTrks = trkArray->GetEntriesFast();
1244
1245 AliExternalTrackParam *track = 0;
1246 AliAODTrack *aodTrack = 0;
1247 Int_t id;
1248
1249 for(Int_t i=0; i<nTrks; i++) {
1250 track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1251 id = (Int_t)track->GetID();
1252 //printf("---> %d\n",id);
1253 if(id<0) continue; // this track is a AliAODRecoDecay
1254 aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1255 v->AddDaughter(aodTrack);
1256 }
1257
1258 return;
1259}
1260//---------------------------------------------------------------------------
1261void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)
1262{
1263 // Checks that the references to the daughter tracks are properly
1264 // assigned and reassigns them if needed
1265 //
1266
1267
1268 TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1269 if(!inputArray) return;
1270
1271 AliAODTrack *track = 0;
1272 AliAODVertex *vertex = 0;
1273
1274 Bool_t needtofix=kFALSE;
1275 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1276 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1277 for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1278 track = (AliAODTrack*)vertex->GetDaughter(id);
1279 if(!track->GetStatus()) needtofix=kTRUE;
1280 }
1281 if(needtofix) break;
1282 }
1283
1284 if(!needtofix) return;
1285
1286
1287 printf("Fixing references\n");
1288
1289 fAODMapSize = 100000;
1290 fAODMap = new Int_t[fAODMapSize];
1291
1292 for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1293 track = aod->GetTrack(i);
1294
1295 // skip pure ITS SA tracks
1296 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1297
1298 // skip tracks without ITS
1299 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1300
1301 // TEMPORARY: check that the cov matrix is there
1302 Double_t covtest[21];
1303 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1304 //
1305
1306 fAODMap[(Int_t)track->GetID()] = i;
1307 }
1308
1309
1310 Int_t ids[4]={-1,-1,-1,-1};
1311 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1312 Bool_t cascade=kFALSE;
1313 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1314 Int_t id=0;
1315 Int_t nDgs = vertex->GetNDaughters();
1316 for(id=0; id<nDgs; id++) {
1317 track = (AliAODTrack*)vertex->GetDaughter(id);
1318 if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1319 ids[id]=(Int_t)track->GetID();
1320 vertex->RemoveDaughter(track);
1321 }
1322 if(cascade) continue;
1323 for(id=0; id<nDgs; id++) {
1324 track = aod->GetTrack(fAODMap[ids[id]]);
1325 vertex->AddDaughter(track);
1326 }
1327
1328 }
1329
1330 return;
1331}
1332//----------------------------------------------------------------------------
1333AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1334 TObjArray *twoTrackArray,AliVEvent *event,
1335 AliAODVertex *secVert,
1336 AliAODRecoDecayHF2Prong *rd2Prong,
1337 Double_t dca,
1338 Bool_t &okDstar) const
1339{
1340 // Make the cascade as a 2Prong decay and check if it passes Dstar
1341 // reconstruction cuts
1342
1343 okDstar = kFALSE;
1344
1345 Bool_t dummy1,dummy2,dummy3;
1346
1347 // We use Make2Prong to construct the AliAODRecoCascadeHF
1348 // (which inherits from AliAODRecoDecayHF2Prong)
1349 AliAODRecoCascadeHF *theCascade =
1350 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1351 dummy1,dummy2,dummy3);
1352 if(!theCascade) return 0x0;
1353
1354 // charge
1355 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1356 theCascade->SetCharge(trackPi->Charge());
1357
1358 //--- selection cuts
1359 //
1360 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1361 if(fInputAOD){
1362 Int_t idSoftPi=(Int_t)trackPi->GetID();
1363 AliAODTrack* trackPiAOD=(AliAODTrack*)event->GetTrack(fAODMap[idSoftPi]);
1364 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
1365 }else{
1366 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1367 }
1368 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1369
1370 AliAODVertex *primVertexAOD=0;
1371 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1372 // take event primary vertex
1373 primVertexAOD = PrimaryVertex();
1374 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1375 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1376 }
1377 // select D*->D0pi
1378 if(fDstar) {
1379 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1380 if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1381 }
1382 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1383 tmpCascade->UnsetOwnPrimaryVtx();
1384 delete tmpCascade; tmpCascade=NULL;
1385 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1386 rd2Prong->UnsetOwnPrimaryVtx();
1387 }
1388 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1389 //---
1390
1391
1392 return theCascade;
1393}
1394
1395
1396//----------------------------------------------------------------------------
1397AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1398 TObjArray *twoTrackArray,AliVEvent *event,
1399 AliAODVertex *secVert,
1400 AliAODv0 *v0,
1401 Double_t dca,
1402 Bool_t &okCascades) const
1403{
1404 //
1405 // Make the cascade as a 2Prong decay and check if it passes
1406 // cascades reconstruction cuts
1407
1408 // AliDebug(2,Form(" building the cascade"));
1409 okCascades= kFALSE;
1410 Bool_t dummy1,dummy2,dummy3;
1411
1412 // We use Make2Prong to construct the AliAODRecoCascadeHF
1413 // (which inherits from AliAODRecoDecayHF2Prong)
1414 AliAODRecoCascadeHF *theCascade =
1415 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1416 dummy1,dummy2,dummy3);
1417 if(!theCascade) return 0x0;
1418
1419 // bachelor track and charge
1420 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1421 theCascade->SetCharge(trackBachelor->Charge());
1422
1423 //--- selection cuts
1424 //
1425
1426 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1427 if(fInputAOD){
1428 Int_t idBachelor=(Int_t)trackBachelor->GetID();
1429 AliAODTrack* trackBachelorAOD=(AliAODTrack*)event->GetTrack(fAODMap[idBachelor]);
1430 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
1431 }else{
1432 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1433 }
1434 tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1435
1436 AliAODVertex *primVertexAOD=0;
1437 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1438 // take event primary vertex
1439 primVertexAOD = PrimaryVertex();
1440 if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1441 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1442 }
1443
1444 // select Cascades
1445 if(fCascades && fInputAOD){
1446 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1447 }
1448 else {
1449 //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
1450 okCascades=kTRUE;
1451 }// no cuts implemented from ESDs
1452 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1453 tmpCascade->UnsetOwnPrimaryVtx();
1454 delete tmpCascade; tmpCascade=NULL;
1455 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1456 //---
1457
1458 return theCascade;
1459}
1460
1461//-----------------------------------------------------------------------------
1462AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1463 TObjArray *twoTrackArray,AliVEvent *event,
1464 AliAODVertex *secVert,Double_t dca,
1465 Bool_t &okD0,Bool_t &okJPSI,
1466 Bool_t &okD0fromDstar) const
1467{
1468 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1469 // reconstruction cuts
1470 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1471
1472 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1473
1474 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1475 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1476 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1477
1478 // propagate tracks to secondary vertex, to compute inv. mass
1479 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1480 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1481
1482 Double_t momentum[3];
1483 postrack->GetPxPyPz(momentum);
1484 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1485 negtrack->GetPxPyPz(momentum);
1486 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1487
1488
1489 // invariant mass cut (try to improve coding here..)
1490 Bool_t okMassCut=kFALSE;
1491 if(!okMassCut && fD0toKpi) if(SelectInvMassAndPt(0,2,px,py,pz)) okMassCut=kTRUE;
1492 if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPt(1,2,px,py,pz)) okMassCut=kTRUE;
1493 if(!okMassCut && fDstar) if(SelectInvMassAndPt(3,2,px,py,pz)) okMassCut=kTRUE;
1494 if(!okMassCut && fCascades) if(SelectInvMassAndPt(5,2,px,py,pz)) okMassCut=kTRUE;
1495 if(!okMassCut) {
1496 //AliDebug(2," candidate didn't pass mass cut");
1497 return 0x0;
1498 }
1499 // primary vertex to be used by this candidate
1500 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1501 if(!primVertexAOD) return 0x0;
1502
1503
1504 Double_t d0z0[2],covd0z0[3];
1505 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1506 d0[0] = d0z0[0];
1507 d0err[0] = TMath::Sqrt(covd0z0[0]);
1508 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1509 d0[1] = d0z0[0];
1510 d0err[1] = TMath::Sqrt(covd0z0[0]);
1511
1512 // create the object AliAODRecoDecayHF2Prong
1513 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1514 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1515 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1516 the2Prong->SetProngIDs(2,id);
1517 delete primVertexAOD; primVertexAOD=NULL;
1518
1519
1520 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1521
1522 // Add daughter references already here
1523 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
1524
1525 // select D0->Kpi
1526 if(fD0toKpi) {
1527 okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
1528 if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1529 }
1530 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1531 // select J/psi from B
1532 if(fJPSItoEle) {
1533 okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1534 }
1535 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1536 // select D0->Kpi from Dstar
1537 if(fDstar) {
1538 okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1539 if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1540 }
1541 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1542 }
1543
1544
1545 // remove the primary vertex (was used only for selection)
1546 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1547 the2Prong->UnsetOwnPrimaryVtx();
1548 }
1549
1550 // get PID info from ESD
1551 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1552 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1553 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1554 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1555 Double_t esdpid[10];
1556 for(Int_t i=0;i<5;i++) {
1557 esdpid[i] = esdpid0[i];
1558 esdpid[5+i] = esdpid1[i];
1559 }
1560 the2Prong->SetPID(2,esdpid);
1561
1562 return the2Prong;
1563}
1564//----------------------------------------------------------------------------
1565AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1566 TObjArray *threeTrackArray,AliVEvent *event,
1567 AliAODVertex *secVert,Double_t dispersion,
1568 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1569 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1570 Bool_t &ok3Prong) const
1571{
1572 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1573 // reconstruction cuts
1574 // E.Bruna, F.Prino
1575
1576
1577 ok3Prong=kFALSE;
1578 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1579
1580 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1581 Double_t momentum[3];
1582
1583
1584 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1585 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1586 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1587
1588 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1589 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1590 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1591 postrack1->GetPxPyPz(momentum);
1592 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1593 negtrack->GetPxPyPz(momentum);
1594 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1595 postrack2->GetPxPyPz(momentum);
1596 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1597
1598 // invariant mass cut for D+, Ds, Lc
1599 Bool_t okMassCut=kFALSE;
1600 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1601 if(!okMassCut && f3Prong) if(SelectInvMassAndPt(2,3,px,py,pz)) okMassCut=kTRUE;
1602 if(!okMassCut) {
1603 //AliDebug(2," candidate didn't pass mass cut");
1604 return 0x0;
1605 }
1606
1607 // primary vertex to be used by this candidate
1608 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1609 if(!primVertexAOD) return 0x0;
1610
1611 Double_t d0z0[2],covd0z0[3];
1612 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1613 d0[0]=d0z0[0];
1614 d0err[0] = TMath::Sqrt(covd0z0[0]);
1615 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1616 d0[1]=d0z0[0];
1617 d0err[1] = TMath::Sqrt(covd0z0[0]);
1618 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1619 d0[2]=d0z0[0];
1620 d0err[2] = TMath::Sqrt(covd0z0[0]);
1621
1622
1623 // create the object AliAODRecoDecayHF3Prong
1624 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1625 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1626 Double_t dist12=TMath::Sqrt((vertexp1n1->GetX()-pos[0])*(vertexp1n1->GetX()-pos[0])+(vertexp1n1->GetY()-pos[1])*(vertexp1n1->GetY()-pos[1])+(vertexp1n1->GetZ()-pos[2])*(vertexp1n1->GetZ()-pos[2]));
1627 Double_t dist23=TMath::Sqrt((vertexp2n1->GetX()-pos[0])*(vertexp2n1->GetX()-pos[0])+(vertexp2n1->GetY()-pos[1])*(vertexp2n1->GetY()-pos[1])+(vertexp2n1->GetZ()-pos[2])*(vertexp2n1->GetZ()-pos[2]));
1628 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1629 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1630 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1631 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1632 the3Prong->SetProngIDs(3,id);
1633
1634 delete primVertexAOD; primVertexAOD=NULL;
1635
1636 // Add daughter references already here
1637 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
1638
1639 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1640 if(f3Prong) {
1641 ok3Prong = kFALSE;
1642
1643 if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1644 ok3Prong = kTRUE;
1645 the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1646 }
1647 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1648 ok3Prong = kTRUE;
1649 the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1650 }
1651 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1652 ok3Prong = kTRUE;
1653 the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1654 }
1655 }
1656 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1657
1658 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1659 the3Prong->UnsetOwnPrimaryVtx();
1660 }
1661
1662 // get PID info from ESD
1663 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1664 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1665 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1666 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1667 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1668 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1669
1670 Double_t esdpid[15];
1671 for(Int_t i=0;i<5;i++) {
1672 esdpid[i] = esdpid0[i];
1673 esdpid[5+i] = esdpid1[i];
1674 esdpid[10+i] = esdpid2[i];
1675 }
1676 the3Prong->SetPID(3,esdpid);
1677
1678 return the3Prong;
1679}
1680//----------------------------------------------------------------------------
1681AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1682 TObjArray *fourTrackArray,AliVEvent *event,
1683 AliAODVertex *secVert,
1684 const AliAODVertex *vertexp1n1,
1685 const AliAODVertex *vertexp1n1p2,
1686 Double_t dcap1n1,Double_t dcap1n2,
1687 Double_t dcap2n1,Double_t dcap2n2,
1688 Bool_t &ok4Prong) const
1689{
1690 // Make 4Prong candidates and check if they pass D0toKpipipi
1691 // reconstruction cuts
1692 // G.E.Bruno, R.Romita
1693
1694 ok4Prong=kFALSE;
1695 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1696
1697 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1698
1699 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1700 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1701 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1702 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1703
1704 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1705 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1706 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1707 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1708
1709 Double_t momentum[3];
1710 postrack1->GetPxPyPz(momentum);
1711 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1712 negtrack1->GetPxPyPz(momentum);
1713 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1714 postrack2->GetPxPyPz(momentum);
1715 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1716 negtrack2->GetPxPyPz(momentum);
1717 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1718
1719 // invariant mass cut for rho or D0 (try to improve coding here..)
1720 Bool_t okMassCut=kFALSE;
1721 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1722 if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
1723 if(SelectInvMassAndPt(4,4,px,py,pz)) okMassCut=kTRUE;
1724 }
1725 if(!okMassCut) {
1726 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1727 //printf(" candidate didn't pass mass cut\n");
1728 return 0x0;
1729 }
1730
1731 // primary vertex to be used by this candidate
1732 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
1733 if(!primVertexAOD) return 0x0;
1734
1735 Double_t d0z0[2],covd0z0[3];
1736 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1737 d0[0]=d0z0[0];
1738 d0err[0] = TMath::Sqrt(covd0z0[0]);
1739 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1740 d0[1]=d0z0[0];
1741 d0err[1] = TMath::Sqrt(covd0z0[0]);
1742 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1743 d0[2]=d0z0[0];
1744 d0err[2] = TMath::Sqrt(covd0z0[0]);
1745 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1746 d0[3]=d0z0[0];
1747 d0err[3] = TMath::Sqrt(covd0z0[0]);
1748
1749
1750 // create the object AliAODRecoDecayHF4Prong
1751 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1752 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1753 Double_t dist12=TMath::Sqrt((vertexp1n1->GetX()-pos[0])*(vertexp1n1->GetX()-pos[0])+(vertexp1n1->GetY()-pos[1])*(vertexp1n1->GetY()-pos[1])+(vertexp1n1->GetZ()-pos[2])*(vertexp1n1->GetZ()-pos[2]));
1754 Double_t dist3=TMath::Sqrt((vertexp1n1p2->GetX()-pos[0])*(vertexp1n1p2->GetX()-pos[0])+(vertexp1n1p2->GetY()-pos[1])*(vertexp1n1p2->GetY()-pos[1])+(vertexp1n1p2->GetZ()-pos[2])*(vertexp1n1p2->GetZ()-pos[2]));
1755 Double_t dist4=TMath::Sqrt((secVert->GetX()-pos[0])*(secVert->GetX()-pos[0])+(secVert->GetY()-pos[1])*(secVert->GetY()-pos[1])+(secVert->GetZ()-pos[2])*(secVert->GetZ()-pos[2]));
1756 Short_t charge=0;
1757 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1758 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1759 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1760 the4Prong->SetProngIDs(4,id);
1761
1762 delete primVertexAOD; primVertexAOD=NULL;
1763
1764 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1765
1766
1767 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1768 the4Prong->UnsetOwnPrimaryVtx();
1769 }
1770
1771
1772 // get PID info from ESD
1773 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1774 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1775 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1776 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1777 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1778 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1779 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1780 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1781
1782 Double_t esdpid[20];
1783 for(Int_t i=0;i<5;i++) {
1784 esdpid[i] = esdpid0[i];
1785 esdpid[5+i] = esdpid1[i];
1786 esdpid[10+i] = esdpid2[i];
1787 esdpid[15+i] = esdpid3[i];
1788 }
1789 the4Prong->SetPID(4,esdpid);
1790
1791 return the4Prong;
1792}
1793//-----------------------------------------------------------------------------
1794AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1795 AliVEvent *event) const
1796{
1797 // Returns primary vertex to be used for this candidate
1798
1799 AliESDVertex *vertexESD = 0;
1800 AliAODVertex *vertexAOD = 0;
1801
1802
1803 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1804 // primary vertex from the input event
1805
1806 vertexESD = new AliESDVertex(*fV1);
1807
1808 } else {
1809 // primary vertex specific to this candidate
1810
1811 Int_t nTrks = trkArray->GetEntriesFast();
1812 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1813
1814 if(fRecoPrimVtxSkippingTrks) {
1815 // recalculating the vertex
1816
1817 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1818 Float_t diamondcovxy[3];
1819 event->GetDiamondCovXY(diamondcovxy);
1820 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1821 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1822 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1823 vertexer->SetVtxStart(diamond);
1824 delete diamond; diamond=NULL;
1825 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1826 vertexer->SetOnlyFitter();
1827 }
1828 Int_t skipped[1000];
1829 Int_t nTrksToSkip=0,id;
1830 AliExternalTrackParam *t = 0;
1831 for(Int_t i=0; i<nTrks; i++) {
1832 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1833 id = (Int_t)t->GetID();
1834 if(id<0) continue;
1835 skipped[nTrksToSkip++] = id;
1836 }
1837 // TEMPORARY FIX
1838 // For AOD, skip also tracks without covariance matrix
1839 if(fInputAOD) {
1840 Double_t covtest[21];
1841 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1842 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1843 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1844 id = (Int_t)vtrack->GetID();
1845 if(id<0) continue;
1846 skipped[nTrksToSkip++] = id;
1847 }
1848 }
1849 }
1850 for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
1851 //
1852 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1853 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1854
1855 } else if(fRmTrksFromPrimVtx && nTrks>0) {
1856 // removing the prongs tracks
1857
1858 TObjArray rmArray(nTrks);
1859 UShort_t *rmId = new UShort_t[nTrks];
1860 AliESDtrack *esdTrack = 0;
1861 AliESDtrack *t = 0;
1862 for(Int_t i=0; i<nTrks; i++) {
1863 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1864 esdTrack = new AliESDtrack(*t);
1865 rmArray.AddLast(esdTrack);
1866 if(esdTrack->GetID()>=0) {
1867 rmId[i]=(UShort_t)esdTrack->GetID();
1868 } else {
1869 rmId[i]=9999;
1870 }
1871 }
1872 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1873 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1874 delete [] rmId; rmId=NULL;
1875 rmArray.Delete();
1876
1877 }
1878
1879 if(!vertexESD) return vertexAOD;
1880 if(vertexESD->GetNContributors()<=0) {
1881 //AliDebug(2,"vertexing failed");
1882 delete vertexESD; vertexESD=NULL;
1883 return vertexAOD;
1884 }
1885
1886 delete vertexer; vertexer=NULL;
1887
1888 }
1889
1890 // convert to AliAODVertex
1891 Double_t pos[3],cov[6],chi2perNDF;
1892 vertexESD->GetXYZ(pos); // position
1893 vertexESD->GetCovMatrix(cov); //covariance matrix
1894 chi2perNDF = vertexESD->GetChi2toNDF();
1895 delete vertexESD; vertexESD=NULL;
1896
1897 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1898
1899 return vertexAOD;
1900}
1901//-----------------------------------------------------------------------------
1902void AliAnalysisVertexingHF::PrintStatus() const {
1903 // Print parameters being used
1904
1905 //printf("Preselections:\n");
1906 // fTrackFilter->Dump();
1907 if(fSecVtxWithKF) {
1908 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1909 } else {
1910 printf("Secondary vertex with AliVertexerTracks\n");
1911 }
1912 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1913 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1914 if(fD0toKpi) {
1915 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1916 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1917 }
1918 if(fDstar) {
1919 printf("Reconstruct D*->D0pi candidates with cuts:\n");
1920 if(fFindVertexForDstar) {
1921 printf(" Reconstruct a secondary vertex for the D*\n");
1922 } else {
1923 printf(" Assume the D* comes from the primary vertex\n");
1924 }
1925 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1926 }
1927 if(fJPSItoEle) {
1928 printf("Reconstruct J/psi from B candidates with cuts:\n");
1929 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1930 }
1931 if(f3Prong) {
1932 printf("Reconstruct 3 prong candidates.\n");
1933 printf(" D+->Kpipi cuts:\n");
1934 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1935 printf(" Ds->KKpi cuts:\n");
1936 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1937 printf(" Lc->pKpi cuts:\n");
1938 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1939 }
1940 if(f4Prong) {
1941 printf("Reconstruct 4 prong candidates.\n");
1942 printf(" D0->Kpipipi cuts:\n");
1943 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1944 }
1945 if(fCascades) {
1946 printf("Reconstruct cascades candidates formed with v0s.\n");
1947 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
1948 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1949 }
1950
1951 return;
1952}
1953//-----------------------------------------------------------------------------
1954AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1955 Double_t &dispersion,Bool_t useTRefArray) const
1956{
1957 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1958
1959 AliESDVertex *vertexESD = 0;
1960 AliAODVertex *vertexAOD = 0;
1961
1962 if(!fSecVtxWithKF) { // AliVertexerTracks
1963
1964 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1965 vertexer->SetVtxStart(fV1);
1966 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1967 delete vertexer; vertexer=NULL;
1968
1969 if(!vertexESD) return vertexAOD;
1970
1971 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1972 //AliDebug(2,"vertexing failed");
1973 delete vertexESD; vertexESD=NULL;
1974 return vertexAOD;
1975 }
1976
1977 Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
1978 if(vertRadius2>8.){
1979 // vertex outside beam pipe, reject candidate to avoid propagation through material
1980 delete vertexESD; vertexESD=NULL;
1981 return vertexAOD;
1982 }
1983
1984 } else { // Kalman Filter vertexer (AliKFParticle)
1985
1986 AliKFParticle::SetField(fBzkG);
1987
1988 AliKFVertex vertexKF;
1989
1990 Int_t nTrks = trkArray->GetEntriesFast();
1991 for(Int_t i=0; i<nTrks; i++) {
1992 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1993 AliKFParticle daughterKF(*esdTrack,211);
1994 vertexKF.AddDaughter(daughterKF);
1995 }
1996 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1997 vertexKF.CovarianceMatrix(),
1998 vertexKF.GetChi2(),
1999 vertexKF.GetNContributors());
2000
2001 }
2002
2003 // convert to AliAODVertex
2004 Double_t pos[3],cov[6],chi2perNDF;
2005 vertexESD->GetXYZ(pos); // position
2006 vertexESD->GetCovMatrix(cov); //covariance matrix
2007 chi2perNDF = vertexESD->GetChi2toNDF();
2008 dispersion = vertexESD->GetDispersion();
2009 delete vertexESD; vertexESD=NULL;
2010
2011 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2012 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2013
2014 return vertexAOD;
2015}
2016//-----------------------------------------------------------------------------
2017Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt(TObjArray *trkArray) const {
2018 // Invariant mass cut on tracks
2019
2020 Int_t nProngs=trkArray->GetEntriesFast();
2021 Int_t retval=kFALSE;
2022
2023 Double_t momentum[3];
2024 Double_t px3[3],py3[3],pz3[3];
2025 Double_t px4[4],py4[4],pz4[4];
2026 AliESDtrack *postrack1=0;
2027 AliESDtrack *postrack2=0;
2028 AliESDtrack *negtrack1=0;
2029 AliESDtrack *negtrack2=0;
2030
2031 switch(nProngs) {
2032 case 3:
2033 // invariant mass cut for D+, Ds, Lc
2034 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
2035 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
2036 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
2037 postrack1->GetPxPyPz(momentum);
2038 px3[0] = momentum[0]; py3[0] = momentum[1]; pz3[0] = momentum[2];
2039 negtrack1->GetPxPyPz(momentum);
2040 px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2];
2041 postrack2->GetPxPyPz(momentum);
2042 px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2];
2043 retval = SelectInvMassAndPt(2,3,px3,py3,pz3);
2044 break;
2045 case 4:
2046 // invariant mass cut for D0->4prong
2047 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
2048 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
2049 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
2050 negtrack2 = (AliESDtrack*)trkArray->UncheckedAt(3);
2051 postrack1->GetPxPyPz(momentum);
2052 px4[0] = momentum[0]; py4[0] = momentum[1]; pz4[0] = momentum[2];
2053 negtrack1->GetPxPyPz(momentum);
2054 px4[1] = momentum[0]; py4[1] = momentum[1]; pz4[1] = momentum[2];
2055 postrack2->GetPxPyPz(momentum);
2056 px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
2057 negtrack2->GetPxPyPz(momentum);
2058 px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
2059 retval = SelectInvMassAndPt(4,4,px4,py4,pz4);
2060 break;
2061 default:
2062 printf("SelectInvMassAndPt(): wrong decay selection\n");
2063 break;
2064 }
2065
2066 return retval;
2067}
2068//-----------------------------------------------------------------------------
2069Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt(Int_t decay,
2070 Int_t nprongs,
2071 Double_t *px,
2072 Double_t *py,
2073 Double_t *pz) const {
2074 // Check invariant mass cut and pt candidate cut
2075
2076 UInt_t pdg2[2],pdg3[3],pdg4[4];
2077 Double_t mPDG,minv;
2078 Double_t minPt=0;
2079
2080 Bool_t retval=kFALSE;
2081 switch (decay)
2082 {
2083 case 0: // D0->Kpi
2084 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2085 // pt cut
2086 minPt=fCutsD0toKpi->GetMinPtCandidate();
2087 if(minPt>0.1)
2088 if(fMassCalc2->Pt2() < minPt*minPt) break;
2089 // mass cut
2090 pdg2[0]=211; pdg2[1]=321;
2091 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2092 minv = fMassCalc2->InvMass(nprongs,pdg2);
2093 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
2094 pdg2[0]=321; pdg2[1]=211;
2095 minv = fMassCalc2->InvMass(nprongs,pdg2);
2096 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
2097 break;
2098 case 1: // JPSI->ee
2099 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2100 // pt cut
2101 minPt=fCutsJpsitoee->GetMinPtCandidate();
2102 if(minPt>0.1)
2103 if(fMassCalc2->Pt2() < minPt*minPt) break;
2104 // mass cut
2105 pdg2[0]=11; pdg2[1]=11;
2106 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2107 minv = fMassCalc2->InvMass(nprongs,pdg2);
2108 if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
2109 break;
2110 case 2: // D+->Kpipi
2111 fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2112 // pt cut
2113 minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
2114 minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
2115 if(minPt>0.1)
2116 if(fMassCalc3->Pt2() < minPt*minPt) break;
2117 // mass cut
2118 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2119 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2120 minv = fMassCalc3->InvMass(nprongs,pdg3);
2121 if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
2122 // Ds+->KKpi
2123 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2124 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2125 minv = fMassCalc3->InvMass(nprongs,pdg3);
2126 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
2127 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2128 minv = fMassCalc3->InvMass(nprongs,pdg3);
2129 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
2130 // Lc->pKpi
2131 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2132 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2133 minv = fMassCalc3->InvMass(nprongs,pdg3);
2134 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
2135 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2136 minv = fMassCalc3->InvMass(nprongs,pdg3);
2137 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
2138 break;
2139 case 3: // D*->D0pi
2140 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2141 // pt cut
2142 minPt=fCutsDStartoKpipi->GetMinPtCandidate();
2143 if(minPt>0.1)
2144 if(fMassCalc2->Pt2() < minPt*minPt) break;
2145 // mass cut
2146 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2147 mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2148 minv = fMassCalc2->InvMass(nprongs,pdg2);
2149 if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
2150 break;
2151 case 4: // D0->Kpipipi without PID
2152 fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2153 // pt cut
2154 minPt=fCutsD0toKpipipi->GetMinPtCandidate();
2155 if(minPt>0.1)
2156 if(fMassCalc4->Pt2() < minPt*minPt) break;
2157 // mass cut
2158 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2159 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2160 minv = fMassCalc4->InvMass(nprongs,pdg4);
2161 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2162 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2163 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2164 minv = fMassCalc4->InvMass(nprongs,pdg4);
2165 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2166 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2167 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2168 minv = fMassCalc4->InvMass(nprongs,pdg4);
2169 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2170 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2171 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2172 minv = fMassCalc4->InvMass(nprongs,pdg4);
2173 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2174 break;
2175 case 5:
2176 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2177 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2178 pdg2[0]=2212;pdg2[1]=310;
2179 minv=fMassCalc2->InvMass(2,pdg2);
2180 if(TMath::Abs(minv-mPDG)<fCutsLctoV0->GetMassCut()) retval=kTRUE;
2181 pdg2[0]=211;pdg2[1]=3122;
2182 minv=fMassCalc2->InvMass(2,pdg2);
2183 if(TMath::Abs(minv-mPDG)<fCutsLctoV0->GetMassCut()) retval=kTRUE;
2184 break;
2185 default:
2186 printf("SelectInvMassAndPt(): wrong decay selection\n");
2187 break;
2188 }
2189
2190 return retval;
2191}
2192//-----------------------------------------------------------------------------
2193void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2194 Int_t trkEntries,
2195 TObjArray &seleTrksArray,TObjArray &tracksAtVertex,
2196 Int_t &nSeleTrks,
2197 UChar_t *seleFlags,Int_t *evtNumber)
2198{
2199 // Apply single-track preselection.
2200 // Fill a TObjArray with selected tracks (for displaced vertices or
2201 // soft pion from D*). Selection flag stored in seleFlags.
2202 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2203 // In case of AOD input, also fill fAODMap for track index<->ID
2204
2205 const AliVVertex *vprimary = event->GetPrimaryVertex();
2206
2207 if(fV1) { delete fV1; fV1=NULL; }
2208 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
2209
2210 Int_t nindices=0;
2211 UShort_t *indices = 0;
2212 Double_t pos[3],cov[6];
2213
2214 if(!fInputAOD) { // ESD
2215 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2216 } else { // AOD
2217 vprimary->GetXYZ(pos);
2218 vprimary->GetCovarianceMatrix(cov);
2219 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2220 indices = new UShort_t[event->GetNumberOfTracks()];
2221 for(Int_t ijk=0; ijk<event->GetNumberOfTracks(); ijk++) indices[ijk]=0;
2222 fAODMapSize = 100000;
2223 fAODMap = new Int_t[fAODMapSize];
2224 }
2225
2226
2227 Int_t entries = (Int_t)event->GetNumberOfTracks();
2228 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2229 nSeleTrks=0;
2230
2231 // transfer ITS tracks from event to arrays
2232 for(Int_t i=0; i<entries; i++) {
2233 AliVTrack *track;
2234 track = (AliVTrack*)event->GetTrack(i);
2235
2236 // skip pure ITS SA tracks
2237 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2238
2239 // skip tracks without ITS
2240 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2241
2242 // TEMPORARY: check that the cov matrix is there
2243 Double_t covtest[21];
2244 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2245 //
2246
2247 if(fInputAOD) {
2248 AliAODTrack *aodt = (AliAODTrack*)track;
2249 if(aodt->GetUsedForPrimVtxFit()) {
2250 indices[nindices]=aodt->GetID(); nindices++;
2251 }
2252 fAODMap[(Int_t)aodt->GetID()] = i;
2253 }
2254
2255 AliESDtrack *esdt = 0;
2256
2257 if(!fInputAOD) {
2258 esdt = (AliESDtrack*)track;
2259 } else {
2260 esdt = new AliESDtrack(track);
2261 }
2262
2263 // single track selection
2264 okDisplaced=kFALSE; okSoftPi=kFALSE;
2265 if(fMixEvent && i<trkEntries){
2266 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2267 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2268 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2269 eventVtx->GetXYZ(vtxPos);
2270 vprimary->GetXYZ(primPos);
2271 eventVtx->GetCovarianceMatrix(primCov);
2272 for(Int_t ind=0;ind<3;ind++){
2273 trasl[ind]=vtxPos[ind]-primPos[ind];
2274 }
2275
2276 Bool_t isTransl=esdt->Translate(trasl,primCov);
2277 if(!isTransl) {
2278 delete esdt;
2279 esdt = NULL;
2280 continue;
2281 }
2282 }
2283
2284 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi) && nSeleTrks<trkEntries) {
2285 esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
2286 seleTrksArray.AddLast(esdt);
2287 tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
2288
2289 seleFlags[nSeleTrks]=0;
2290 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2291 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2292 nSeleTrks++;
2293 } else {
2294 if(fInputAOD) delete esdt;
2295 esdt = NULL;
2296 continue;
2297 }
2298
2299 } // end loop on tracks
2300
2301 // primary vertex from AOD
2302 if(fInputAOD) {
2303 delete fV1; fV1=NULL;
2304 vprimary->GetXYZ(pos);
2305 vprimary->GetCovarianceMatrix(cov);
2306 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2307 Int_t ncontr=nindices;
2308 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2309 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2310 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2311 fV1->SetTitle(vprimary->GetTitle());
2312 fV1->SetIndices(nindices,indices);
2313 }
2314 if(indices) { delete [] indices; indices=NULL; }
2315
2316
2317 return;
2318}
2319//-----------------------------------------------------------------------------
2320void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
2321 //
2322 // Set the selection bit for PID
2323 //
2324 if(cuts->GetPidHF()) {
2325 Bool_t usepid=cuts->GetIsUsePID();
2326 cuts->SetUsePID(kTRUE);
2327 if(cuts->IsSelectedPID(rd))
2328 rd->SetSelectionBit(bit);
2329 cuts->SetUsePID(usepid);
2330 }
2331 return;
2332}
2333//-----------------------------------------------------------------------------
2334Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2335 Bool_t &okDisplaced,Bool_t &okSoftPi) const
2336{
2337 // Check if track passes some kinematical cuts
2338
2339 // this is needed to store the impact parameters
2340 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2341
2342 UInt_t selectInfo;
2343 //
2344 // Track selection, displaced tracks
2345 selectInfo = 0;
2346 if(fTrackFilter) {
2347 selectInfo = fTrackFilter->IsSelected(trk);
2348 }
2349 if(selectInfo) okDisplaced=kTRUE;
2350
2351 // Track selection, soft pions
2352 selectInfo = 0;
2353 if(fDstar && fTrackFilterSoftPi) {
2354 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2355 }
2356 if(selectInfo) okSoftPi=kTRUE;
2357
2358 if(okDisplaced || okSoftPi) return kTRUE;
2359
2360 return kFALSE;
2361}
2362
2363
2364//-----------------------------------------------------------------------------
2365AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2366 //
2367 // Transform ESDv0 to AODv0
2368 //
2369 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2370 // and creates an AODv0 out of them
2371 //
2372 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2373 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2374
2375 // create the v0 neutral track to compute the DCA to the primary vertex
2376 Double_t xyz[3], pxpypz[3];
2377 esdV0->XvYvZv(xyz);
2378 esdV0->PxPyPz(pxpypz);
2379 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2380 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2381 if(!trackesdV0) {
2382 delete vertexV0;
2383 return 0;
2384 }
2385 Double_t d0z0[2],covd0z0[3];
2386 AliAODVertex *primVertexAOD = PrimaryVertex();
2387 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2388 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2389 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2390 Double_t dcaV0DaughterToPrimVertex[2];
2391 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2392 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2393 if( !posV0track || !negV0track) {
2394 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2395 delete vertexV0;
2396 delete primVertexAOD;
2397 return 0;
2398 }
2399 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2400 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2401 // else
2402 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2403 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2404 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2405 // else
2406 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2407 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2408 Double_t pmom[3],nmom[3];
2409 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2410 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2411
2412 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2413 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2414
2415 delete trackesdV0;
2416 delete primVertexAOD;
2417
2418 return aodV0;
2419}
2420//-----------------------------------------------------------------------------
2421void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, AliExternalTrackParam* extpar) const{
2422 // Set the stored track parameters at primary vertex into AliESDtrack
2423 Double_t xyz[3], pxpypz[3], cv[21];
2424 extpar->PxPyPz(pxpypz);
2425 extpar->XvYvZv(xyz);
2426 extpar->GetCovarianceXYZPxPyPz(cv);
2427 Short_t sign=extpar->Charge();
2428 esdt->Set(xyz,pxpypz,cv,sign);
2429 return;
2430}