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