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