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