]> 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
New method to fix the references of the daughter tracks
[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; }
c8ab4e4f 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
18e7db05 1374 Bool_t okLcksp=0, okLcLpi=0;
a07ad8e0 1375 if(fCascades && fInputAOD){
cb8088a2 1376 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1377 if(okCascades==1) okLcksp=1;
1378 if(okCascades==2) okLcLpi=1;
1379 if(okCascades==3) { okLcksp=1; okLcLpi=1;}
a07ad8e0 1380 }
18e7db05 1381 else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=kTRUE; }// no cuts implemented from ESDs
a07ad8e0 1382 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1383 tmpCascade->UnsetOwnPrimaryVtx();
1384 delete tmpCascade; tmpCascade=NULL;
1385 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1386 //---
1387
1388 return theCascade;
1389}
1390
2ff20727 1391//-----------------------------------------------------------------------------
6a213b59 1392AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
2ff20727 1393 TObjArray *twoTrackArray,AliVEvent *event,
dcb444c9 1394 AliAODVertex *secVert,Double_t dca,
2ff20727 1395 Bool_t &okD0,Bool_t &okJPSI,
1396 Bool_t &okD0fromDstar) const
6a213b59 1397{
1398 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1399 // reconstruction cuts
1400 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1401
2ff20727 1402 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
6a213b59 1403
1404 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
2ff20727 1405 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1406 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
6a213b59 1407
1408 // propagate tracks to secondary vertex, to compute inv. mass
dcb444c9 1409 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1410 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
6a213b59 1411
1412 Double_t momentum[3];
1413 postrack->GetPxPyPz(momentum);
1414 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1415 negtrack->GetPxPyPz(momentum);
1416 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1417
1418
1419 // invariant mass cut (try to improve coding here..)
1420 Bool_t okMassCut=kFALSE;
dcb444c9 1421 if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
6a213b59 1422 if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
2ff20727 1423 if(!okMassCut && fDstar) if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
6a213b59 1424 if(!okMassCut) {
dcb444c9 1425 AliDebug(2," candidate didn't pass mass cut");
6a213b59 1426 return 0x0;
1427 }
dcb444c9 1428 // primary vertex to be used by this candidate
2ff20727 1429 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
dcb444c9 1430 if(!primVertexAOD) return 0x0;
6a213b59 1431
a25935a9 1432
dcb444c9 1433 Double_t d0z0[2],covd0z0[3];
1434 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
6a213b59 1435 d0[0] = d0z0[0];
1436 d0err[0] = TMath::Sqrt(covd0z0[0]);
dcb444c9 1437 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
6a213b59 1438 d0[1] = d0z0[0];
1439 d0err[1] = TMath::Sqrt(covd0z0[0]);
b6e66d86 1440
6a213b59 1441 // create the object AliAODRecoDecayHF2Prong
dcb444c9 1442 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
6a213b59 1443 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
dcb444c9 1444 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1445 the2Prong->SetProngIDs(2,id);
0a65d33f 1446 delete primVertexAOD; primVertexAOD=NULL;
6a213b59 1447
f3014dc3 1448
1449 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1450 // select D0->Kpi
a9b75906 1451 if(fD0toKpi) okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
f3014dc3 1452 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1453 // select J/psi from B
a9b75906 1454 if(fJPSItoEle) okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
f3014dc3 1455 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1456 // select D0->Kpi from Dstar
aa0fbb4a 1457 if(fDstar) okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
f3014dc3 1458 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1459 }
1460
6a213b59 1461
2ff20727 1462 // remove the primary vertex (was used only for selection)
a25935a9 1463 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
b6e66d86 1464 the2Prong->UnsetOwnPrimaryVtx();
b6e66d86 1465 }
1466
2ff20727 1467 // get PID info from ESD
5770b08b 1468 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1469 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1470 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1471 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2ff20727 1472 Double_t esdpid[10];
1473 for(Int_t i=0;i<5;i++) {
1474 esdpid[i] = esdpid0[i];
1475 esdpid[5+i] = esdpid1[i];
6a213b59 1476 }
2ff20727 1477 the2Prong->SetPID(2,esdpid);
6a213b59 1478
6a213b59 1479 return the2Prong;
1480}
1481//----------------------------------------------------------------------------
1482AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
dcb444c9 1483 TObjArray *threeTrackArray,AliVEvent *event,
1484 AliAODVertex *secVert,Double_t dispersion,
c4cdf105 1485 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
6a213b59 1486 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1487 Bool_t &ok3Prong) const
1488{
1489 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1490 // reconstruction cuts
1491 // E.Bruna, F.Prino
1492
dcb444c9 1493
6a213b59 1494 ok3Prong=kFALSE;
673ac607 1495 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
dcb444c9 1496
1497 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1498 Double_t momentum[3];
6a213b59 1499
6a213b59 1500
1501 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
dcb444c9 1502 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
6a213b59 1503 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1504
dcb444c9 1505 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1506 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1507 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
6a213b59 1508 postrack1->GetPxPyPz(momentum);
1509 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1510 negtrack->GetPxPyPz(momentum);
1511 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1512 postrack2->GetPxPyPz(momentum);
1513 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1514
b82f6d67 1515 // invariant mass cut for D+, Ds, Lc
6a213b59 1516 Bool_t okMassCut=kFALSE;
e0732246 1517 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
6a213b59 1518 if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1519 if(!okMassCut) {
dcb444c9 1520 AliDebug(2," candidate didn't pass mass cut");
6a213b59 1521 return 0x0;
1522 }
1523
dcb444c9 1524 // primary vertex to be used by this candidate
2ff20727 1525 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
dcb444c9 1526 if(!primVertexAOD) return 0x0;
1527
1528 Double_t d0z0[2],covd0z0[3];
1529 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1530 d0[0]=d0z0[0];
1531 d0err[0] = TMath::Sqrt(covd0z0[0]);
1532 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1533 d0[1]=d0z0[0];
1534 d0err[1] = TMath::Sqrt(covd0z0[0]);
1535 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1536 d0[2]=d0z0[0];
1537 d0err[2] = TMath::Sqrt(covd0z0[0]);
6a213b59 1538
6a213b59 1539
1540 // create the object AliAODRecoDecayHF3Prong
dcb444c9 1541 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
6a213b59 1542 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
dcb444c9 1543 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]));
1544 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 1545 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
dcb444c9 1546 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
6a213b59 1547 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
dcb444c9 1548 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1549 the3Prong->SetProngIDs(3,id);
6a213b59 1550
0a65d33f 1551 delete primVertexAOD; primVertexAOD=NULL;
6a213b59 1552
b82f6d67 1553 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1554 if(f3Prong) {
1555 ok3Prong = kFALSE;
a9b75906 1556
e3d40058 1557 if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1558 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1559 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
a9b75906 1560
b82f6d67 1561 }
6a213b59 1562 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
dcb444c9 1563
a25935a9 1564 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
b6e66d86 1565 the3Prong->UnsetOwnPrimaryVtx();
b6e66d86 1566 }
dcb444c9 1567
2ff20727 1568 // get PID info from ESD
5770b08b 1569 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1570 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1571 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1572 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1573 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1574 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2ff20727 1575
1576 Double_t esdpid[15];
1577 for(Int_t i=0;i<5;i++) {
1578 esdpid[i] = esdpid0[i];
1579 esdpid[5+i] = esdpid1[i];
1580 esdpid[10+i] = esdpid2[i];
6a213b59 1581 }
2ff20727 1582 the3Prong->SetPID(3,esdpid);
6a213b59 1583
6a213b59 1584 return the3Prong;
1585}
1586//----------------------------------------------------------------------------
1587AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
dcb444c9 1588 TObjArray *fourTrackArray,AliVEvent *event,
721c0b8f 1589 AliAODVertex *secVert,
c4cdf105 1590 const AliAODVertex *vertexp1n1,
1591 const AliAODVertex *vertexp1n1p2,
721c0b8f 1592 Double_t dcap1n1,Double_t dcap1n2,
1593 Double_t dcap2n1,Double_t dcap2n2,
1594 Bool_t &ok4Prong) const
6a213b59 1595{
1596 // Make 4Prong candidates and check if they pass D0toKpipipi
1597 // reconstruction cuts
1598 // G.E.Bruno, R.Romita
1599
1600 ok4Prong=kFALSE;
673ac607 1601 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
6a213b59 1602
1603 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
6a213b59 1604
6a213b59 1605 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1606 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1607 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1608 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1609
dcb444c9 1610 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1611 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1612 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1613 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
6a213b59 1614
1615 Double_t momentum[3];
1616 postrack1->GetPxPyPz(momentum);
1617 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1618 negtrack1->GetPxPyPz(momentum);
1619 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1620 postrack2->GetPxPyPz(momentum);
1621 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1622 negtrack2->GetPxPyPz(momentum);
1623 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1624
1625 // invariant mass cut for rho or D0 (try to improve coding here..)
0a65d33f 1626 Bool_t okMassCut=kFALSE;
e0732246 1627 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
cb8088a2 1628 if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
0a65d33f 1629 if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1630 }
1631 if(!okMassCut) {
1632 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1633 //printf(" candidate didn't pass mass cut\n");
1634 return 0x0;
1635 }
6a213b59 1636
dcb444c9 1637 // primary vertex to be used by this candidate
2ff20727 1638 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
dcb444c9 1639 if(!primVertexAOD) return 0x0;
1640
721c0b8f 1641 Double_t d0z0[2],covd0z0[3];
1642 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1643 d0[0]=d0z0[0];
1644 d0err[0] = TMath::Sqrt(covd0z0[0]);
1645 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1646 d0[1]=d0z0[0];
1647 d0err[1] = TMath::Sqrt(covd0z0[0]);
1648 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1649 d0[2]=d0z0[0];
1650 d0err[2] = TMath::Sqrt(covd0z0[0]);
1651 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1652 d0[3]=d0z0[0];
1653 d0err[3] = TMath::Sqrt(covd0z0[0]);
1654
6a213b59 1655
1656 // create the object AliAODRecoDecayHF4Prong
dcb444c9 1657 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
721c0b8f 1658 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
dcb444c9 1659 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 1660 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]));
1661 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 1662 Short_t charge=0;
721c0b8f 1663 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
6a213b59 1664 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
dcb444c9 1665 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1666 the4Prong->SetProngIDs(4,id);
6a213b59 1667
0a65d33f 1668 delete primVertexAOD; primVertexAOD=NULL;
6a213b59 1669
a9b75906 1670 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
e3d40058 1671
6a213b59 1672
a25935a9 1673 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
b6e66d86 1674 the4Prong->UnsetOwnPrimaryVtx();
b6e66d86 1675 }
6a213b59 1676
721c0b8f 1677
6a213b59 1678 // get PID info from ESD
721c0b8f 1679 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1680 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1681 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1682 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1683 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1684 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1685 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1686 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
6a213b59 1687
1688 Double_t esdpid[20];
1689 for(Int_t i=0;i<5;i++) {
2ff20727 1690 esdpid[i] = esdpid0[i];
1691 esdpid[5+i] = esdpid1[i];
6a213b59 1692 esdpid[10+i] = esdpid2[i];
1693 esdpid[15+i] = esdpid3[i];
1694 }
1695 the4Prong->SetPID(4,esdpid);
721c0b8f 1696
6a213b59 1697 return the4Prong;
1698}
1699//-----------------------------------------------------------------------------
c4cdf105 1700AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2ff20727 1701 AliVEvent *event) const
6a213b59 1702{
2ff20727 1703 // Returns primary vertex to be used for this candidate
6a213b59 1704
dcb444c9 1705 AliESDVertex *vertexESD = 0;
1706 AliAODVertex *vertexAOD = 0;
1707
dcb444c9 1708
1709 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1710 // primary vertex from the input event
1711
1712 vertexESD = new AliESDVertex(*fV1);
1713
1714 } else {
1715 // primary vertex specific to this candidate
1716
2ff20727 1717 Int_t nTrks = trkArray->GetEntriesFast();
dcb444c9 1718 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1719
1720 if(fRecoPrimVtxSkippingTrks) {
1721 // recalculating the vertex
1722
1723 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1724 Float_t diamondcovxy[3];
1725 event->GetDiamondCovXY(diamondcovxy);
1726 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
8bb09289 1727 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
dcb444c9 1728 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1729 vertexer->SetVtxStart(diamond);
1730 delete diamond; diamond=NULL;
1731 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1732 vertexer->SetOnlyFitter();
1733 }
0a65d33f 1734 Int_t skipped[1000];
e18fbfa7 1735 Int_t nTrksToSkip=0,id;
1736 AliExternalTrackParam *t = 0;
dcb444c9 1737 for(Int_t i=0; i<nTrks; i++) {
e18fbfa7 1738 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1739 id = (Int_t)t->GetID();
1740 if(id<0) continue;
1741 skipped[nTrksToSkip++] = id;
dcb444c9 1742 }
0a65d33f 1743 // TEMPORARY FIX
1744 // For AOD, skip also tracks without covariance matrix
1745 if(fInputAOD) {
1746 Double_t covtest[21];
1747 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1748 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1749 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
3ecb887f 1750 id = (Int_t)vtrack->GetID();
0a65d33f 1751 if(id<0) continue;
1752 skipped[nTrksToSkip++] = id;
1753 }
1754 }
1755 }
1756 //
e18fbfa7 1757 vertexer->SetSkipTracks(nTrksToSkip,skipped);
dcb444c9 1758 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1759
1760 } else if(fRmTrksFromPrimVtx) {
1761 // removing the prongs tracks
1762
1763 TObjArray rmArray(nTrks);
1764 UShort_t *rmId = new UShort_t[nTrks];
1765 AliESDtrack *esdTrack = 0;
1766 AliESDtrack *t = 0;
1767 for(Int_t i=0; i<nTrks; i++) {
1768 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1769 esdTrack = new AliESDtrack(*t);
1770 rmArray.AddLast(esdTrack);
2ff20727 1771 if(esdTrack->GetID()>=0) {
1772 rmId[i]=(UShort_t)esdTrack->GetID();
1773 } else {
1774 rmId[i]=9999;
1775 }
dcb444c9 1776 }
1777 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1778 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1779 delete [] rmId; rmId=NULL;
1780 rmArray.Delete();
1781
6a213b59 1782 }
6a213b59 1783
dcb444c9 1784 if(!vertexESD) return vertexAOD;
1785 if(vertexESD->GetNContributors()<=0) {
1786 AliDebug(2,"vertexing failed");
1787 delete vertexESD; vertexESD=NULL;
1788 return vertexAOD;
6a213b59 1789 }
dcb444c9 1790
1791 delete vertexer; vertexer=NULL;
1792
6a213b59 1793 }
1794
dcb444c9 1795 // convert to AliAODVertex
1796 Double_t pos[3],cov[6],chi2perNDF;
1797 vertexESD->GetXYZ(pos); // position
1798 vertexESD->GetCovMatrix(cov); //covariance matrix
1799 chi2perNDF = vertexESD->GetChi2toNDF();
1800 delete vertexESD; vertexESD=NULL;
1801
1802 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
6a213b59 1803
dcb444c9 1804 return vertexAOD;
6a213b59 1805}
1806//-----------------------------------------------------------------------------
1807void AliAnalysisVertexingHF::PrintStatus() const {
1808 // Print parameters being used
1809
f8fa4595 1810 //printf("Preselections:\n");
1811 // fTrackFilter->Dump();
4d5c9633 1812 if(fSecVtxWithKF) {
1813 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1814 } else {
1815 printf("Secondary vertex with AliVertexerTracks\n");
1816 }
6a213b59 1817 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1818 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1819 if(fD0toKpi) {
1820 printf("Reconstruct D0->Kpi candidates with cuts:\n");
e3d40058 1821 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
6a213b59 1822 }
2ff20727 1823 if(fDstar) {
da6fefc3 1824 printf("Reconstruct D*->D0pi candidates with cuts:\n");
eee95e18 1825 if(fFindVertexForDstar) {
1826 printf(" Reconstruct a secondary vertex for the D*\n");
1827 } else {
1828 printf(" Assume the D* comes from the primary vertex\n");
1829 }
aa0fbb4a 1830 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
2ff20727 1831 }
6a213b59 1832 if(fJPSItoEle) {
1833 printf("Reconstruct J/psi from B candidates with cuts:\n");
e3d40058 1834 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
6a213b59 1835 }
1836 if(f3Prong) {
1837 printf("Reconstruct 3 prong candidates.\n");
4d5c9633 1838 printf(" D+->Kpipi cuts:\n");
e3d40058 1839 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
4d5c9633 1840 printf(" Ds->KKpi cuts:\n");
e3d40058 1841 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
6ea608bf 1842 printf(" Lc->pKpi cuts:\n");
e3d40058 1843 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
6a213b59 1844 }
e3d40058 1845 if(f4Prong) {
1846 printf("Reconstruct 4 prong candidates.\n");
1847 printf(" D0->Kpipipi cuts:\n");
1848 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1849 }
a07ad8e0 1850 if(fCascades) {
1851 printf("Reconstruct cascades candidates formed with v0s.\n");
1852 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
1853 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1854 }
6a213b59 1855
1856 return;
1857}
1858//-----------------------------------------------------------------------------
dcb444c9 1859AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
dcfa35b3 1860 Double_t &dispersion,Bool_t useTRefArray) const
dcb444c9 1861{
1862 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1863
1864 AliESDVertex *vertexESD = 0;
1865 AliAODVertex *vertexAOD = 0;
1866
1867 if(!fSecVtxWithKF) { // AliVertexerTracks
1868
1869 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1870 vertexer->SetVtxStart(fV1);
1871 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1872 delete vertexer; vertexer=NULL;
1873
1874 if(!vertexESD) return vertexAOD;
1875
1876 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1877 AliDebug(2,"vertexing failed");
1878 delete vertexESD; vertexESD=NULL;
1879 return vertexAOD;
1880 }
1881
1882 } else { // Kalman Filter vertexer (AliKFParticle)
1883
1884 AliKFParticle::SetField(fBzkG);
1885
a353d2e4 1886 AliKFVertex vertexKF;
dcb444c9 1887
1888 Int_t nTrks = trkArray->GetEntriesFast();
1889 for(Int_t i=0; i<nTrks; i++) {
1890 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1891 AliKFParticle daughterKF(*esdTrack,211);
1892 vertexKF.AddDaughter(daughterKF);
1893 }
a353d2e4 1894 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1895 vertexKF.CovarianceMatrix(),
1896 vertexKF.GetChi2(),
1897 vertexKF.GetNContributors());
dcb444c9 1898
1899 }
1900
1901 // convert to AliAODVertex
1902 Double_t pos[3],cov[6],chi2perNDF;
1903 vertexESD->GetXYZ(pos); // position
1904 vertexESD->GetCovMatrix(cov); //covariance matrix
1905 chi2perNDF = vertexESD->GetChi2toNDF();
1906 dispersion = vertexESD->GetDispersion();
1907 delete vertexESD; vertexESD=NULL;
1908
dcfa35b3 1909 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1910 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
dcb444c9 1911
1912 return vertexAOD;
1913}
1914//-----------------------------------------------------------------------------
e0732246 1915Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
1916 // Invariant mass cut on tracks
1917
1918 Int_t nProngs=trkArray->GetEntriesFast();
1919 Int_t retval=kFALSE;
1920
1921 Double_t momentum[3];
1922 Double_t px3[3],py3[3],pz3[3];
1923 Double_t px4[4],py4[4],pz4[4];
1924 AliESDtrack *postrack1=0;
1925 AliESDtrack *postrack2=0;
1926 AliESDtrack *negtrack1=0;
1927 AliESDtrack *negtrack2=0;
1928
1929 switch(nProngs) {
1930 case 3:
1931 // invariant mass cut for D+, Ds, Lc
1932 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1933 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
1934 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1935 postrack1->GetPxPyPz(momentum);
1936 px3[0] = momentum[0]; py3[0] = momentum[1]; pz3[0] = momentum[2];
1937 negtrack1->GetPxPyPz(momentum);
1938 px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2];
1939 postrack2->GetPxPyPz(momentum);
1940 px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2];
1941 retval = SelectInvMass(2,3,px3,py3,pz3);
1942 break;
1943 case 4:
1944 // invariant mass cut for D0->4prong
1945 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1946 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
1947 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1948 negtrack2 = (AliESDtrack*)trkArray->UncheckedAt(3);
1949 postrack1->GetPxPyPz(momentum);
1950 px4[0] = momentum[0]; py4[0] = momentum[1]; pz4[0] = momentum[2];
1951 negtrack1->GetPxPyPz(momentum);
1952 px4[1] = momentum[0]; py4[1] = momentum[1]; pz4[1] = momentum[2];
1953 postrack2->GetPxPyPz(momentum);
1954 px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
1955 negtrack2->GetPxPyPz(momentum);
1956 px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
1957 retval = SelectInvMass(4,4,px4,py4,pz4);
1958 break;
1959 default:
1960 printf("SelectInvMass(): wrong decay selection\n");
1961 break;
1962 }
1963
1964 return retval;
1965}
1966//-----------------------------------------------------------------------------
6a213b59 1967Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1968 Int_t nprongs,
1969 Double_t *px,
1970 Double_t *py,
1971 Double_t *pz) const {
1972 // Check invariant mass cut
1973
721c0b8f 1974 UInt_t pdg2[2],pdg3[3],pdg4[4];
6a213b59 1975 Double_t mPDG,minv;
1976
1977 Bool_t retval=kFALSE;
1978 switch (decay)
1979 {
1980 case 0: // D0->Kpi
e0732246 1981 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
6a213b59 1982 pdg2[0]=211; pdg2[1]=321;
1983 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
e0732246 1984 minv = fMassCalc2->InvMass(nprongs,pdg2);
a9b75906 1985 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
6a213b59 1986 pdg2[0]=321; pdg2[1]=211;
e0732246 1987 minv = fMassCalc2->InvMass(nprongs,pdg2);
a9b75906 1988 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
6a213b59 1989 break;
1990 case 1: // JPSI->ee
e0732246 1991 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
6a213b59 1992 pdg2[0]=11; pdg2[1]=11;
1993 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
e0732246 1994 minv = fMassCalc2->InvMass(nprongs,pdg2);
a9b75906 1995 if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
6a213b59 1996 break;
1997 case 2: // D+->Kpipi
e0732246 1998 fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
6a213b59 1999 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2000 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
e0732246 2001 minv = fMassCalc3->InvMass(nprongs,pdg3);
a9b75906 2002 if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
b82f6d67 2003 // Ds+->KKpi
2004 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2005 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
e0732246 2006 minv = fMassCalc3->InvMass(nprongs,pdg3);
a9b75906 2007 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
b82f6d67 2008 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
e0732246 2009 minv = fMassCalc3->InvMass(nprongs,pdg3);
a9b75906 2010 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
b82f6d67 2011 // Lc->pKpi
2012 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2013 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
e0732246 2014 minv = fMassCalc3->InvMass(nprongs,pdg3);
a9b75906 2015 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
b82f6d67 2016 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
e0732246 2017 minv = fMassCalc3->InvMass(nprongs,pdg3);
a9b75906 2018 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
6a213b59 2019 break;
2ff20727 2020 case 3: // D*->D0pi
e0732246 2021 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
682639e9 2022 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2ff20727 2023 mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
e0732246 2024 minv = fMassCalc2->InvMass(nprongs,pdg2);
aa0fbb4a 2025 if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
2ff20727 2026 break;
e3d40058 2027 case 4: // D0->Kpipipi without PID
e0732246 2028 fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
721c0b8f 2029 pdg4[0]=321; pdg4[1]=211; 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]=321; pdg4[2]=211; 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]=321; pdg4[3]=211;
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 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2042 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
e0732246 2043 minv = fMassCalc4->InvMass(nprongs,pdg4);
a9b75906 2044 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
721c0b8f 2045 break;
6a213b59 2046 default:
2ff20727 2047 printf("SelectInvMass(): wrong decay selection\n");
6a213b59 2048 break;
2049 }
2050
6a213b59 2051 return retval;
2052}
2053//-----------------------------------------------------------------------------
c4cdf105 2054void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
e11ae259 2055 Int_t trkEntries,
2ff20727 2056 TObjArray &seleTrksArray,Int_t &nSeleTrks,
a25935a9 2057 UChar_t *seleFlags,Int_t *evtNumber)
6a213b59 2058{
2ff20727 2059 // Apply single-track preselection.
2060 // Fill a TObjArray with selected tracks (for displaced vertices or
2061 // soft pion from D*). Selection flag stored in seleFlags.
dcb444c9 2062 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2063 // In case of AOD input, also fill fAODMap for track index<->ID
2064
2065 const AliVVertex *vprimary = event->GetPrimaryVertex();
2066
2067 if(fV1) { delete fV1; fV1=NULL; }
c8ab4e4f 2068 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
dcb444c9 2069
2070 Int_t nindices=0;
2071 UShort_t *indices = 0;
2072 Double_t pos[3],cov[6];
2073
2074 if(!fInputAOD) { // ESD
2075 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2076 } else { // AOD
2077 vprimary->GetXYZ(pos);
2078 vprimary->GetCovarianceMatrix(cov);
2079 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2080 indices = new UShort_t[event->GetNumberOfTracks()];
c8ab4e4f 2081 fAODMapSize = 100000;
2082 fAODMap = new Int_t[fAODMapSize];
dcb444c9 2083 }
2084
2085
dcb444c9 2086 Int_t entries = (Int_t)event->GetNumberOfTracks();
2ff20727 2087 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2088 nSeleTrks=0;
6a213b59 2089
dcb444c9 2090 // transfer ITS tracks from event to arrays
6a213b59 2091 for(Int_t i=0; i<entries; i++) {
a25935a9 2092 AliVTrack *track;
2093 track = (AliVTrack*)event->GetTrack(i);
dcb444c9 2094
f7c170b4 2095 // skip pure ITS SA tracks
ac325db6 2096 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
f7c170b4 2097
7162b6d2 2098 // skip tracks without ITS
2099 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2100
65dad444 2101 // TEMPORARY: check that the cov matrix is there
0ab1bd93 2102 Double_t covtest[21];
2103 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2104 //
65dad444 2105
dcb444c9 2106 if(fInputAOD) {
2107 AliAODTrack *aodt = (AliAODTrack*)track;
2108 if(aodt->GetUsedForPrimVtxFit()) {
2109 indices[nindices]=aodt->GetID(); nindices++;
2110 }
2111 fAODMap[(Int_t)aodt->GetID()] = i;
2112 }
6a213b59 2113
dcb444c9 2114 AliESDtrack *esdt = 0;
b6e66d86 2115
dcb444c9 2116 if(!fInputAOD) {
2117 esdt = (AliESDtrack*)track;
2118 } else {
2119 esdt = new AliESDtrack(track);
2120 }
6a213b59 2121
2122 // single track selection
2ff20727 2123 okDisplaced=kFALSE; okSoftPi=kFALSE;
e11ae259 2124 if(fMixEvent && i<trkEntries){
a25935a9 2125 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2126 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2127 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2128 eventVtx->GetXYZ(vtxPos);
2129 vprimary->GetXYZ(primPos);
2130 eventVtx->GetCovarianceMatrix(primCov);
2131 for(Int_t ind=0;ind<3;ind++){
2132 trasl[ind]=vtxPos[ind]-primPos[ind];
2133 }
2134
2135 Bool_t isTransl=esdt->Translate(trasl,primCov);
2136 if(!isTransl) {
2137 delete esdt;
2138 esdt = NULL;
2139 continue;
2140 }
2141 }
2142
e11ae259 2143 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi) && nSeleTrks<trkEntries) {
2ff20727 2144 seleTrksArray.AddLast(esdt);
2145 seleFlags[nSeleTrks]=0;
2146 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2147 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2148 nSeleTrks++;
2149 } else {
dcb444c9 2150 if(fInputAOD) delete esdt;
2151 esdt = NULL;
2152 continue;
2ff20727 2153 }
6a213b59 2154
dcb444c9 2155 } // end loop on tracks
2156
2157 // primary vertex from AOD
2158 if(fInputAOD) {
2159 delete fV1; fV1=NULL;
2160 vprimary->GetXYZ(pos);
2161 vprimary->GetCovarianceMatrix(cov);
2162 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2163 Int_t ncontr=nindices;
a6e0ebfe 2164 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
dcb444c9 2165 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2166 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2167 fV1->SetTitle(vprimary->GetTitle());
2168 fV1->SetIndices(nindices,indices);
2169 }
2170 if(indices) { delete [] indices; indices=NULL; }
2171
6a213b59 2172
721c0b8f 2173 return;
2174}
2175//-----------------------------------------------------------------------------
2ff20727 2176Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2177 Bool_t &okDisplaced,Bool_t &okSoftPi) const
6a213b59 2178{
2179 // Check if track passes some kinematical cuts
6a213b59 2180
f8fa4595 2181 // this is needed to store the impact parameters
2182 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2183
2ff20727 2184 UInt_t selectInfo;
f8fa4595 2185 //
2ff20727 2186 // Track selection, displaced tracks
2187 selectInfo = 0;
f8fa4595 2188 if(fTrackFilter) {
2189 selectInfo = fTrackFilter->IsSelected(trk);
6a213b59 2190 }
2ff20727 2191 if(selectInfo) okDisplaced=kTRUE;
e11ae259 2192
2ff20727 2193 // Track selection, soft pions
2194 selectInfo = 0;
2195 if(fDstar && fTrackFilterSoftPi) {
2196 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2197 }
2198 if(selectInfo) okSoftPi=kTRUE;
6a213b59 2199
2ff20727 2200 if(okDisplaced || okSoftPi) return kTRUE;
f8fa4595 2201
2ff20727 2202 return kFALSE;
6a213b59 2203}
a07ad8e0 2204
2205
2206//-----------------------------------------------------------------------------
2207AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2208 //
2209 // Transform ESDv0 to AODv0
2210 //
2211 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2212 // and creates an AODv0 out of them
2213 //
18e7db05 2214 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
a07ad8e0 2215 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2216
2217 // create the v0 neutral track to compute the DCA to the primary vertex
2218 Double_t xyz[3], pxpypz[3];
2219 esdV0->XvYvZv(xyz);
2220 esdV0->PxPyPz(pxpypz);
2221 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2222 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2223 if(!trackesdV0) return 0;
2224 Double_t d0z0[2],covd0z0[3];
18e7db05 2225 AliAODVertex *primVertexAOD = PrimaryVertex();
2226 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
a07ad8e0 2227 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2228 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2229 Double_t dcaV0DaughterToPrimVertex[2];
2230 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2231 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
18e7db05 2232 if( !posV0track || !negV0track) {
2233 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2234 return 0;
2235 }
2236 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
a07ad8e0 2237 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2238 // else
2239 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
18e7db05 2240 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
a07ad8e0 2241 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2242 // else
2243 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
18e7db05 2244 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2245 Double_t pmom[3],nmom[3];
a07ad8e0 2246 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2247 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2248
2249 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2250 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2251
18e7db05 2252 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2253 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
a07ad8e0 2254
2255 return aodV0;
2256}
6a213b59 2257//-----------------------------------------------------------------------------