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