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