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