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