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