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