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