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