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