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