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