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