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