reader and analysis container pointers need to be streamed, if not plugin does not...
[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//----------------------------------------------------------------------------
b82f6d67 27#include <TFile.h>
6a213b59 28#include <TDatabasePDG.h>
29#include <TString.h>
dcb444c9 30#include "AliLog.h"
31#include "AliVEvent.h"
32#include "AliVVertex.h"
33#include "AliVTrack.h"
6a213b59 34#include "AliVertexerTracks.h"
a353d2e4 35#include "AliKFVertex.h"
dcb444c9 36#include "AliESDEvent.h"
6a213b59 37#include "AliESDVertex.h"
e18fbfa7 38#include "AliExternalTrackParam.h"
dcfa35b3 39#include "AliNeutralTrackParam.h"
dcb444c9 40#include "AliESDtrack.h"
f8fa4595 41#include "AliESDtrackCuts.h"
b056c5e3 42#include "AliAODEvent.h"
6a213b59 43#include "AliAODRecoDecay.h"
44#include "AliAODRecoDecayHF.h"
45#include "AliAODRecoDecayHF2Prong.h"
46#include "AliAODRecoDecayHF3Prong.h"
47#include "AliAODRecoDecayHF4Prong.h"
2ff20727 48#include "AliAODRecoCascadeHF.h"
f3014dc3 49#include "AliRDHFCutsD0toKpi.h"
f8fa4595 50#include "AliAnalysisFilter.h"
6a213b59 51#include "AliAnalysisVertexingHF.h"
a25935a9 52#include "AliMixedEvent.h"
6a213b59 53
54ClassImp(AliAnalysisVertexingHF)
55
56//----------------------------------------------------------------------------
57AliAnalysisVertexingHF::AliAnalysisVertexingHF():
dcb444c9 58fInputAOD(kFALSE),
c8ab4e4f 59fAODMapSize(0),
60fAODMap(0),
b82f6d67 61fBzkG(0.),
62fSecVtxWithKF(kFALSE),
6a213b59 63fRecoPrimVtxSkippingTrks(kFALSE),
64fRmTrksFromPrimVtx(kFALSE),
65fV1(0x0),
b82f6d67 66fD0toKpi(kTRUE),
67fJPSItoEle(kTRUE),
68f3Prong(kTRUE),
69f4Prong(kTRUE),
2ff20727 70fDstar(kTRUE),
dc963de9 71fLikeSign(kFALSE),
a25935a9 72fMixEvent(kFALSE),
2ff20727 73fTrackFilter(0x0),
eee95e18 74fTrackFilterSoftPi(0x0),
f3014dc3 75fCutsD0toKpi(0x0),
eee95e18 76fFindVertexForDstar(kTRUE)
6a213b59 77{
78 // Default constructor
79
80 SetD0toKpiCuts();
81 SetBtoJPSICuts();
82 SetDplusCuts();
b82f6d67 83 SetDsCuts();
84 SetLcCuts();
2ff20727 85 SetDstarCuts();
721c0b8f 86 SetD0to4ProngsCuts();
6a213b59 87}
13977a79 88//--------------------------------------------------------------------------
89AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
90TNamed(source),
dcb444c9 91fInputAOD(source.fInputAOD),
c8ab4e4f 92fAODMapSize(source.fAODMapSize),
93fAODMap(source.fAODMap),
b82f6d67 94fBzkG(source.fBzkG),
95fSecVtxWithKF(source.fSecVtxWithKF),
13977a79 96fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
97fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
98fV1(source.fV1),
13977a79 99fD0toKpi(source.fD0toKpi),
100fJPSItoEle(source.fJPSItoEle),
101f3Prong(source.f3Prong),
102f4Prong(source.f4Prong),
2ff20727 103fDstar(source.fDstar),
dc963de9 104fLikeSign(source.fLikeSign),
a25935a9 105fMixEvent(source.fMixEvent),
2ff20727 106fTrackFilter(source.fTrackFilter),
eee95e18 107fTrackFilterSoftPi(source.fTrackFilterSoftPi),
f3014dc3 108fCutsD0toKpi(source.fCutsD0toKpi),
eee95e18 109fFindVertexForDstar(source.fFindVertexForDstar)
13977a79 110{
111 //
112 // Copy constructor
113 //
114 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
115 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
116 for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
f7c170b4 117 for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
2ff20727 118 for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
119 for(Int_t i=0; i<5; i++) fDstarCuts[i]=source.fDstarCuts[i];
721c0b8f 120 for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
13977a79 121}
122//--------------------------------------------------------------------------
123AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
124{
125 //
126 // assignment operator
127 //
128 if(&source == this) return *this;
dcb444c9 129 fInputAOD = source.fInputAOD;
c8ab4e4f 130 fAODMapSize = source.fAODMapSize;
b82f6d67 131 fBzkG = source.fBzkG;
132 fSecVtxWithKF = source.fSecVtxWithKF;
13977a79 133 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
134 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
135 fV1 = source.fV1;
13977a79 136 fD0toKpi = source.fD0toKpi;
137 fJPSItoEle = source.fJPSItoEle;
138 f3Prong = source.f3Prong;
139 f4Prong = source.f4Prong;
2ff20727 140 fDstar = source.fDstar;
dc963de9 141 fLikeSign = source.fLikeSign;
a25935a9 142 fMixEvent = source.fMixEvent;
f8fa4595 143 fTrackFilter = source.fTrackFilter;
2ff20727 144 fTrackFilterSoftPi = source.fTrackFilterSoftPi;
f3014dc3 145 fCutsD0toKpi = source.fCutsD0toKpi;
eee95e18 146 fFindVertexForDstar = source.fFindVertexForDstar;
13977a79 147
148 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
149 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
150 for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
f7c170b4 151 for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
2ff20727 152 for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
153 for(Int_t i=0; i<5; i++) fDstarCuts[i]=source.fDstarCuts[i];
721c0b8f 154 for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
13977a79 155
156 return *this;
157}
6a213b59 158//----------------------------------------------------------------------------
159AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
160 // Destructor
13977a79 161 if(fV1) { delete fV1; fV1=0; }
f8fa4595 162 if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
2ff20727 163 if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
f3014dc3 164 if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
c8ab4e4f 165 if(fAODMap) { delete fAODMap; fAODMap=0; }
6a213b59 166}
167//----------------------------------------------------------------------------
dcb444c9 168void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
169 TClonesArray *aodVerticesHFTClArr,
170 TClonesArray *aodD0toKpiTClArr,
171 TClonesArray *aodJPSItoEleTClArr,
172 TClonesArray *aodCharm3ProngTClArr,
2ff20727 173 TClonesArray *aodCharm4ProngTClArr,
dc963de9 174 TClonesArray *aodDstarTClArr,
423fb9ae 175 TClonesArray *aodLikeSign2ProngTClArr,
176 TClonesArray *aodLikeSign3ProngTClArr)
b056c5e3 177{
178 // Find heavy-flavour vertex candidates
dcb444c9 179 // Input: ESD or AOD
b056c5e3 180 // Output: AOD (additional branches added)
181
a25935a9 182 if(!fMixEvent){
183 TString evtype = event->IsA()->GetName();
184 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
185 } // if we do mixing AliVEvent is a AliMixedEvent
dcb444c9 186
e15f55cb 187 if(fInputAOD) {
188 AliDebug(2,"Creating HF candidates from AOD");
189 } else {
190 AliDebug(2,"Creating HF candidates from ESD");
191 }
192
b056c5e3 193 if(!aodVerticesHFTClArr) {
194 printf("ERROR: no aodVerticesHFTClArr");
195 return;
196 }
fa925a06 197 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
b056c5e3 198 printf("ERROR: no aodD0toKpiTClArr");
199 return;
200 }
201 if(fJPSItoEle && !aodJPSItoEleTClArr) {
202 printf("ERROR: no aodJPSItoEleTClArr");
203 return;
204 }
205 if(f3Prong && !aodCharm3ProngTClArr) {
206 printf("ERROR: no aodCharm3ProngTClArr");
207 return;
208 }
209 if(f4Prong && !aodCharm4ProngTClArr) {
210 printf("ERROR: no aodCharm4ProngTClArr");
211 return;
212 }
2ff20727 213 if(fDstar && !aodDstarTClArr) {
214 printf("ERROR: no aodDstarTClArr");
215 return;
216 }
423fb9ae 217 if(fLikeSign && !aodLikeSign2ProngTClArr) {
218 printf("ERROR: no aodLikeSign2ProngTClArr");
219 return;
220 }
221 if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
222 printf("ERROR: no aodLikeSign2ProngTClArr");
dc963de9 223 return;
224 }
b056c5e3 225
226 // delete candidates from previous event and create references
423fb9ae 227 Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
b056c5e3 228 aodVerticesHFTClArr->Delete();
229 iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
230 TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
fa925a06 231 if(fD0toKpi || fDstar) {
b056c5e3 232 aodD0toKpiTClArr->Delete();
233 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
234 }
235 if(fJPSItoEle) {
236 aodJPSItoEleTClArr->Delete();
237 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
238 }
239 if(f3Prong) {
240 aodCharm3ProngTClArr->Delete();
241 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
242 }
243 if(f4Prong) {
244 aodCharm4ProngTClArr->Delete();
245 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
246 }
2ff20727 247 if(fDstar) {
248 aodDstarTClArr->Delete();
249 iDstar = aodDstarTClArr->GetEntriesFast();
250 }
423fb9ae 251 if(fLikeSign) {
252 aodLikeSign2ProngTClArr->Delete();
253 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
254 }
255 if(fLikeSign && f3Prong) {
256 aodLikeSign3ProngTClArr->Delete();
257 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
dc963de9 258 }
423fb9ae 259
260 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
261 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
262 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
263 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
264 TClonesArray &aodDstarRef = *aodDstarTClArr;
265 TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
266 TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
dc963de9 267
b056c5e3 268
2ff20727 269 AliAODRecoDecayHF2Prong *io2Prong = 0;
270 AliAODRecoDecayHF3Prong *io3Prong = 0;
271 AliAODRecoDecayHF4Prong *io4Prong = 0;
272 AliAODRecoCascadeHF *ioCascade = 0;
b056c5e3 273
2ff20727 274 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries;
721c0b8f 275 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaCasc;
b056c5e3 276 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
2ff20727 277 Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
b056c5e3 278 AliESDtrack *postrack1 = 0;
279 AliESDtrack *postrack2 = 0;
280 AliESDtrack *negtrack1 = 0;
281 AliESDtrack *negtrack2 = 0;
2ff20727 282 AliESDtrack *trackPi = 0;
3cbc09cd 283 //Double_t dcaMax = fCutsD0toKpi->GetDCACut();
b056c5e3 284 Double_t dcaMax = fD0toKpiCuts[1];
285 if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
286 if(dcaMax < fDplusCuts[11]) dcaMax=fDplusCuts[11];
721c0b8f 287 if(dcaMax < fD0to4ProngsCuts[1]) dcaMax=fD0to4ProngsCuts[1];
dcb444c9 288 AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
b056c5e3 289
b056c5e3 290
dcb444c9 291 // get Bz
292 fBzkG = (Double_t)event->GetMagneticField();
b056c5e3 293
dcb444c9 294 trkEntries = (Int_t)event->GetNumberOfTracks();
295 AliDebug(1,Form(" Number of tracks: %d",trkEntries));
357cd934 296
297 if(trkEntries<2) {
298 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
299 return;
300 }
b056c5e3 301
3cbc09cd 302 // event selection
303 //if(!fCutsD0toKpi->IsEventSelected(event)) return;
357cd934 304 const AliVVertex *primary = event->GetPrimaryVertex();
305 if(!primary) {
306 AliDebug(1," No primary vertex from tracks");
307 return;
308 }
a25935a9 309
310 TString primTitle=primary->GetTitle();
b6e66d86 311 if(!primTitle.Contains("VertexerTracks") ||
312 primary->GetNContributors()<=0) {
357cd934 313 AliDebug(1," No primary vertex from tracks");
b056c5e3 314 return;
315 }
a25935a9 316
b056c5e3 317
dcb444c9 318 // call function that applies sigle-track selection,
2ff20727 319 // for displaced tracks and soft pions (both charges) for D*,
dcb444c9 320 // and retrieves primary vertex
2ff20727 321 TObjArray seleTrksArray(trkEntries);
322 UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
323 Int_t nSeleTrks=0;
a25935a9 324 Int_t *evtNumber = new Int_t[trkEntries];
325 SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
b056c5e3 326
2ff20727 327 AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
b056c5e3 328
2ff20727 329 TObjArray *twoTrackArray1 = new TObjArray(2);
330 TObjArray *twoTrackArray2 = new TObjArray(2);
331 TObjArray *twoTrackArrayCasc = new TObjArray(2);
332 TObjArray *threeTrackArray = new TObjArray(3);
333 TObjArray *fourTrackArray = new TObjArray(4);
b056c5e3 334
dcb444c9 335 Double_t dispersion;
423fb9ae 336 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
dcb444c9 337
2ff20727 338 AliAODRecoDecayHF *rd = 0;
339 AliAODRecoCascadeHF *rc = 0;
340
b056c5e3 341 // LOOP ON POSITIVE TRACKS
2ff20727 342 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
dc963de9 343
2ff20727 344 if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
dc963de9 345
b056c5e3 346 // get track from tracks array
2ff20727 347 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
dc963de9 348
349 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
350
351 if(postrack1->Charge()<0 && !fLikeSign) continue;
352
b056c5e3 353 // LOOP ON NEGATIVE TRACKS
2ff20727 354 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
dc963de9 355
2ff20727 356 if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
dc963de9 357
358 if(iTrkN1==iTrkP1) continue;
359
b056c5e3 360 // get track from tracks array
2ff20727 361 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
dc963de9 362
363 if(negtrack1->Charge()>0 && !fLikeSign) continue;
364
365 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
366
a25935a9 367 if(fMixEvent) {
368 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
369 }
370
dc963de9 371 if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
423fb9ae 372 isLikeSign2Prong=kTRUE;
dc963de9 373 if(!fLikeSign) continue;
374 if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
375 } else { // unlike-sign
423fb9ae 376 isLikeSign2Prong=kFALSE;
dc963de9 377 if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
a25935a9 378 if(fMixEvent) {
379 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
380 }
381
dc963de9 382 }
383
dcb444c9 384 // back to primary vertex
385 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
386 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
dc963de9 387
b056c5e3 388 // DCA between the two tracks
389 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
390 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
dc963de9 391
b056c5e3 392 // Vertexing
393 twoTrackArray1->AddAt(postrack1,0);
394 twoTrackArray1->AddAt(negtrack1,1);
dcb444c9 395 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
b056c5e3 396 if(!vertexp1n1) {
397 twoTrackArray1->Clear();
398 negtrack1=0;
399 continue;
400 }
dc963de9 401
402 // 2 prong candidate
403 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
b6e66d86 404
2ff20727 405 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
b6e66d86 406
423fb9ae 407 if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
2ff20727 408 // add the vertex and the decay to the AOD
409 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
410 if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
423fb9ae 411 if(!isLikeSign2Prong) {
dc963de9 412 if(okD0) {
413 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
414 rd->SetSecondaryVtx(v2Prong);
415 v2Prong->SetParent(rd);
460cd990 416 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
dc963de9 417 }
418 if(okJPSI) {
419 rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
420 rd->SetSecondaryVtx(v2Prong);
421 if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
460cd990 422 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
dc963de9 423 }
423fb9ae 424 } else { // isLikeSign2Prong
425 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
2ff20727 426 rd->SetSecondaryVtx(v2Prong);
427 v2Prong->SetParent(rd);
460cd990 428 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
dcb444c9 429 }
b056c5e3 430 }
dc963de9 431 // D* candidates
423fb9ae 432 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
dc963de9 433 // write references in io2Prong
434 if(fInputAOD) {
2ff20727 435 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
dc963de9 436 } else {
437 vertexp1n1->AddDaughter(postrack1);
438 vertexp1n1->AddDaughter(negtrack1);
2ff20727 439 }
dc963de9 440 io2Prong->SetSecondaryVtx(vertexp1n1);
441 //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
2ff20727 442 // create a track from the D0
dcfa35b3 443 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
dc963de9 444
2ff20727 445 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
446 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
dc963de9 447
2ff20727 448 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
dc963de9 449
2ff20727 450 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
dc963de9 451
a25935a9 452 if(fMixEvent) {
453 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
454 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
455 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
456 }
457
2ff20727 458 if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
dc963de9 459
eee95e18 460 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
461 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
462
2ff20727 463 // get track from tracks array
464 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
2ff20727 465 trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
dc963de9 466
2ff20727 467 twoTrackArrayCasc->AddAt(trackPi,0);
468 twoTrackArrayCasc->AddAt(trackD0,1);
eee95e18 469
470 AliAODVertex *vertexCasc = 0;
471
472 if(fFindVertexForDstar) {
473 // DCA between the two tracks
474 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
475 // Vertexing
476 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
477 } else {
478 // assume Dstar decays at the primary vertex
479 Double_t pos[3],cov[6],chi2perNDF;
480 fV1->GetXYZ(pos);
481 fV1->GetCovMatrix(cov);
482 chi2perNDF = fV1->GetChi2toNDF();
483 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
484 dcaCasc = 0.;
485 }
2ff20727 486 if(!vertexCasc) {
487 twoTrackArrayCasc->Clear();
488 trackPi=0;
489 continue;
490 }
dc963de9 491
2ff20727 492 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
493 if(okDstar) {
dc963de9 494 // add the D0 to the AOD (if not already done)
2ff20727 495 if(!okD0) {
2ff20727 496 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
497 if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
498 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
499 rd->SetSecondaryVtx(v2Prong);
500 v2Prong->SetParent(rd);
460cd990 501 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
2ff20727 502 okD0=kTRUE; // this is done to add it only once
503 }
504 // add the vertex and the cascade to the AOD
505 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
dcfa35b3 506 if(fInputAOD) {
507 AddDaughterRefs(vCasc,event,twoTrackArrayCasc); // add the pion
508 } else {
509 vCasc->AddDaughter(rd); // just to fill ref #0
510 }
511 vCasc->AddDaughter(rd); // add the D0 (in ref #1)
2ff20727 512 rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
2ff20727 513 rc->SetSecondaryVtx(vCasc);
514 vCasc->SetParent(rc);
460cd990 515 if(fInputAOD) rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
2ff20727 516 }
517 twoTrackArrayCasc->Clear();
518 trackPi=0;
0a65d33f 519 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
b6e66d86 520 delete vertexCasc; vertexCasc=NULL;
2ff20727 521 } // end loop on soft pi tracks
dc963de9 522
2ff20727 523 if(trackD0) {delete trackD0; trackD0=NULL;}
dc963de9 524
2ff20727 525 }
0a65d33f 526 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
2ff20727 527 }
dc963de9 528
b056c5e3 529 twoTrackArray1->Clear();
423fb9ae 530 if( (!f3Prong && !f4Prong) ||
531 (isLikeSign2Prong && !f3Prong) ) {
b056c5e3 532 negtrack1=0;
533 delete vertexp1n1;
534 continue;
535 }
536
537
538 // 2nd LOOP ON POSITIVE TRACKS
2ff20727 539 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
dc963de9 540
2ff20727 541 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
dc963de9 542
2ff20727 543 if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
dc963de9 544
b056c5e3 545 // get track from tracks array
2ff20727 546 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
dc963de9 547
548 if(postrack2->Charge()<0) continue;
549
550 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
551
a25935a9 552 if(fMixEvent) {
553 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
554 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
555 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
556 }
557
423fb9ae 558 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
559 if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
560 isLikeSign3Prong=kTRUE;
561 } else { // not ok
562 continue;
563 }
564 } else { // normal triplet (+-+)
565 isLikeSign3Prong=kFALSE;
a25935a9 566 if(fMixEvent) {
567 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
568 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
569 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
570 }
423fb9ae 571 }
572
dcb444c9 573 // back to primary vertex
574 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
575 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
576 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
2ff20727 577 //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
dc963de9 578
b056c5e3 579 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
580 if(dcap2n1>dcaMax) { postrack2=0; continue; }
581 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
582 if(dcap1p2>dcaMax) { postrack2=0; continue; }
583
584 // Vertexing
585 twoTrackArray2->AddAt(postrack2,0);
586 twoTrackArray2->AddAt(negtrack1,1);
dcb444c9 587 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
b056c5e3 588 if(!vertexp2n1) {
589 twoTrackArray2->Clear();
590 postrack2=0;
423fb9ae 591 continue;
b056c5e3 592 }
dc963de9 593
594 // 3 prong candidates
b056c5e3 595 if(f3Prong) {
423fb9ae 596 if(postrack2->Charge()>0) {
597 threeTrackArray->AddAt(postrack1,0);
598 threeTrackArray->AddAt(negtrack1,1);
599 threeTrackArray->AddAt(postrack2,2);
600 } else {
601 threeTrackArray->AddAt(negtrack1,0);
602 threeTrackArray->AddAt(postrack1,1);
603 threeTrackArray->AddAt(postrack2,2);
604 }
dc963de9 605
dcb444c9 606 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
607 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
b056c5e3 608 if(ok3Prong) {
2ff20727 609 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
610 if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
423fb9ae 611 if(!isLikeSign3Prong) {
612 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
613 rd->SetSecondaryVtx(v3Prong);
614 v3Prong->SetParent(rd);
615 } else { // isLikeSign3Prong
616 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
617 rd->SetSecondaryVtx(v3Prong);
618 v3Prong->SetParent(rd);
460cd990 619 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
423fb9ae 620 }
b056c5e3 621 }
b6e66d86 622 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
623 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
b056c5e3 624 }
dc963de9 625
626 // 4 prong candidates
423fb9ae 627 if(f4Prong
628 // don't make 4 prong with like-sign pairs and triplets
629 && !isLikeSign2Prong && !isLikeSign3Prong
fae0f82d 630 // track-to-track dca cuts already now
631 && dcap1n1 < fD0to4ProngsCuts[1]
632 && dcap2n1 < fD0to4ProngsCuts[1]) {
dc963de9 633
721c0b8f 634 // back to primary vertex
635 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
636 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
637 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
638 // Vertexing for these 3 (can be taken from above?)
639 threeTrackArray->AddAt(postrack1,0);
640 threeTrackArray->AddAt(negtrack1,1);
641 threeTrackArray->AddAt(postrack2,2);
642 AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
643
b056c5e3 644 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
2ff20727 645 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
dc963de9 646
2ff20727 647 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
dc963de9 648
2ff20727 649 if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
dc963de9 650
b056c5e3 651 // get track from tracks array
2ff20727 652 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
dc963de9 653
654 if(negtrack2->Charge()>0) continue;
655
656 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
a25935a9 657 if(fMixEvent){
658 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
659 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
660 evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
661 evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
662 evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
663 evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
664 }
dc963de9 665
dcb444c9 666 // back to primary vertex
667 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
668 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
669 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
670 negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
b056c5e3 671 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
fae0f82d 672 if(dcap1n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
721c0b8f 673 dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
fae0f82d 674 if(dcap2n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
dc963de9 675
b056c5e3 676 // Vertexing
677 fourTrackArray->AddAt(postrack1,0);
678 fourTrackArray->AddAt(negtrack1,1);
679 fourTrackArray->AddAt(postrack2,2);
680 fourTrackArray->AddAt(negtrack2,3);
dc963de9 681
dcb444c9 682 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
721c0b8f 683 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
b056c5e3 684 if(ok4Prong) {
2ff20727 685 AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
686 if(fInputAOD) AddDaughterRefs(v4Prong,event,fourTrackArray);
687 rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
688 rd->SetSecondaryVtx(v4Prong);
689 v4Prong->SetParent(rd);
460cd990 690 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
b056c5e3 691 }
dc963de9 692
b6e66d86 693 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
694 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
b056c5e3 695 fourTrackArray->Clear();
696 negtrack2 = 0;
dc963de9 697
b056c5e3 698 } // end loop on negative tracks
dc963de9 699
721c0b8f 700 threeTrackArray->Clear();
701 delete vertexp1n1p2;
702
b056c5e3 703 }
dc963de9 704
b056c5e3 705 postrack2 = 0;
706 delete vertexp2n1;
dc963de9 707
b056c5e3 708 } // end 2nd loop on positive tracks
dc963de9 709
b056c5e3 710 twoTrackArray2->Clear();
711
423fb9ae 712 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
2ff20727 713 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
dc963de9 714
2ff20727 715 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
dc963de9 716
2ff20727 717 if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
dc963de9 718
b056c5e3 719 // get track from tracks array
2ff20727 720 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
dc963de9 721
722 if(negtrack2->Charge()>0) continue;
723
724 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
725
a25935a9 726 if(fMixEvent) {
727 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
728 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
729 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
730 }
731
423fb9ae 732 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
733 if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
734 isLikeSign3Prong=kTRUE;
735 } else { // not ok
736 continue;
737 }
738 } else { // normal triplet (-+-)
a25935a9 739 isLikeSign3Prong=kFALSE;
423fb9ae 740 }
741
dcb444c9 742 // back to primary vertex
743 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
744 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
745 negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
2ff20727 746 //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
dc963de9 747
b056c5e3 748 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
749 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
750 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
751 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
752
753 // Vertexing
754 twoTrackArray2->AddAt(postrack1,0);
755 twoTrackArray2->AddAt(negtrack2,1);
dcb444c9 756
757 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
b056c5e3 758 if(!vertexp1n2) {
759 twoTrackArray2->Clear();
760 negtrack2=0;
761 continue;
762 }
dc963de9 763
b056c5e3 764 if(f3Prong) {
765 threeTrackArray->AddAt(negtrack1,0);
766 threeTrackArray->AddAt(postrack1,1);
767 threeTrackArray->AddAt(negtrack2,2);
dcb444c9 768 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
769 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
b056c5e3 770 if(ok3Prong) {
2ff20727 771 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
772 if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
423fb9ae 773 if(!isLikeSign3Prong) {
774 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
775 rd->SetSecondaryVtx(v3Prong);
776 v3Prong->SetParent(rd);
777 } else { // isLikeSign3Prong
778 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
779 rd->SetSecondaryVtx(v3Prong);
780 v3Prong->SetParent(rd);
460cd990 781 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
423fb9ae 782 }
b056c5e3 783 }
b6e66d86 784 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
785 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
b056c5e3 786 }
dc963de9 787
b056c5e3 788 negtrack2 = 0;
789 delete vertexp1n2;
dc963de9 790
b056c5e3 791 } // end 2nd loop on negative tracks
423fb9ae 792
b056c5e3 793 twoTrackArray2->Clear();
794
795 negtrack1 = 0;
796 delete vertexp1n1;
797 } // end 1st loop on negative tracks
798
799 postrack1 = 0;
800 } // end 1st loop on positive tracks
801
802
2ff20727 803 AliDebug(1,Form(" Total HF vertices in event = %d;",
804 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
b056c5e3 805 if(fD0toKpi) {
dcb444c9 806 AliDebug(1,Form(" D0->Kpi in event = %d;",
807 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
b056c5e3 808 }
809 if(fJPSItoEle) {
dcb444c9 810 AliDebug(1,Form(" JPSI->ee in event = %d;",
811 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
b056c5e3 812 }
813 if(f3Prong) {
dcb444c9 814 AliDebug(1,Form(" Charm->3Prong in event = %d;",
815 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
b056c5e3 816 }
817 if(f4Prong) {
dcb444c9 818 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
819 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
b056c5e3 820 }
2ff20727 821 if(fDstar) {
822 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
823 (Int_t)aodDstarTClArr->GetEntriesFast()));
824 }
dc963de9 825 if(fLikeSign) {
423fb9ae 826 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
827 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
828 }
829 if(fLikeSign && f3Prong) {
830 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
831 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
dc963de9 832 }
b056c5e3 833
834
dcb444c9 835 twoTrackArray1->Delete(); delete twoTrackArray1;
836 twoTrackArray2->Delete(); delete twoTrackArray2;
2ff20727 837 twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
b056c5e3 838 threeTrackArray->Clear();
839 threeTrackArray->Delete(); delete threeTrackArray;
dcb444c9 840 fourTrackArray->Delete(); delete fourTrackArray;
2ff20727 841 delete [] seleFlags; seleFlags=NULL;
a25935a9 842 if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
b056c5e3 843
dcb444c9 844 if(fInputAOD) {
b6e66d86 845 seleTrksArray.Delete();
c8ab4e4f 846 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
dcb444c9 847 }
b056c5e3 848
849 return;
850}
851//----------------------------------------------------------------------------
c4cdf105 852void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
853 const AliVEvent *event,
854 const TObjArray *trkArray) const
6a213b59 855{
dcb444c9 856 // Add the AOD tracks as daughters of the vertex (TRef)
6a213b59 857
dcb444c9 858 Int_t nTrks = trkArray->GetEntriesFast();
6a213b59 859
e18fbfa7 860 AliExternalTrackParam *track = 0;
dcb444c9 861 AliAODTrack *aodTrack = 0;
862 Int_t id;
6a213b59 863
dcb444c9 864 for(Int_t i=0; i<nTrks; i++) {
e18fbfa7 865 track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
866 id = (Int_t)track->GetID();
2ff20727 867 //printf("---> %d\n",id);
e18fbfa7 868 if(id<0) continue; // this track is a AliAODRecoDecay
dcb444c9 869 aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
870 v->AddDaughter(aodTrack);
6a213b59 871 }
6a213b59 872
873 return;
dcb444c9 874}
6a213b59 875//----------------------------------------------------------------------------
2ff20727 876AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
877 TObjArray *twoTrackArray,AliVEvent *event,
878 AliAODVertex *secVert,
879 AliAODRecoDecayHF2Prong *rd2Prong,
880 Double_t dca,
881 Bool_t &okDstar) const
882{
883 // Make the cascade as a 2Prong decay and check if it passes Dstar
884 // reconstruction cuts
885
886 okDstar = kFALSE;
887
888 Bool_t dummy1,dummy2,dummy3;
889
890 // We use Make2Prong to construct the AliAODRecoCascadeHF
891 // (which inherits from AliAODRecoDecayHF2Prong)
892 AliAODRecoCascadeHF *theCascade =
893 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
894 dummy1,dummy2,dummy3);
895 if(!theCascade) return 0x0;
896
897 // charge
898 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
899 theCascade->SetCharge(trackPi->Charge());
900
901 //--- selection cuts
902 //
903 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
dcfa35b3 904 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
905 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
2ff20727 906 AliAODVertex *primVertexAOD=0;
907 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
908 // take event primary vertex
909 primVertexAOD = PrimaryVertex();
910 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
911 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
912 }
913 // select D*->D0pi
914 if(fDstar) {
915 Bool_t testD0=kTRUE;
916 okDstar = tmpCascade->SelectDstar(fDstarCuts,fD0fromDstarCuts,testD0);
917 }
dcfa35b3 918 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2ff20727 919 tmpCascade->UnsetOwnPrimaryVtx();
2ff20727 920 delete tmpCascade; tmpCascade=NULL;
a25935a9 921 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
2ff20727 922 rd2Prong->UnsetOwnPrimaryVtx();
2ff20727 923 }
b6e66d86 924 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2ff20727 925 //---
926
927 return theCascade;
928}
929//-----------------------------------------------------------------------------
6a213b59 930AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
2ff20727 931 TObjArray *twoTrackArray,AliVEvent *event,
dcb444c9 932 AliAODVertex *secVert,Double_t dca,
2ff20727 933 Bool_t &okD0,Bool_t &okJPSI,
934 Bool_t &okD0fromDstar) const
6a213b59 935{
936 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
937 // reconstruction cuts
938 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
939
2ff20727 940 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
6a213b59 941
942 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
2ff20727 943 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
944 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
6a213b59 945
946 // propagate tracks to secondary vertex, to compute inv. mass
dcb444c9 947 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
948 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
6a213b59 949
950 Double_t momentum[3];
951 postrack->GetPxPyPz(momentum);
952 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
953 negtrack->GetPxPyPz(momentum);
954 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
955
956
957 // invariant mass cut (try to improve coding here..)
958 Bool_t okMassCut=kFALSE;
dcb444c9 959 if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
6a213b59 960 if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
2ff20727 961 if(!okMassCut && fDstar) if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
6a213b59 962 if(!okMassCut) {
dcb444c9 963 AliDebug(2," candidate didn't pass mass cut");
6a213b59 964 return 0x0;
965 }
dcb444c9 966 // primary vertex to be used by this candidate
2ff20727 967 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
dcb444c9 968 if(!primVertexAOD) return 0x0;
6a213b59 969
a25935a9 970
dcb444c9 971 Double_t d0z0[2],covd0z0[3];
972 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
6a213b59 973 d0[0] = d0z0[0];
974 d0err[0] = TMath::Sqrt(covd0z0[0]);
dcb444c9 975 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
6a213b59 976 d0[1] = d0z0[0];
977 d0err[1] = TMath::Sqrt(covd0z0[0]);
b6e66d86 978
6a213b59 979 // create the object AliAODRecoDecayHF2Prong
dcb444c9 980 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
6a213b59 981 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
dcb444c9 982 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
983 the2Prong->SetProngIDs(2,id);
0a65d33f 984 delete primVertexAOD; primVertexAOD=NULL;
6a213b59 985
f3014dc3 986
987 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
988 // select D0->Kpi
989 Int_t checkD0,checkD0bar;
990 if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
991 //if(fD0toKpi) okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
992 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
993 // select J/psi from B
994 Int_t checkJPSI;
995 if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
996 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
997 // select D0->Kpi from Dstar
998 if(fDstar) okD0fromDstar = the2Prong->SelectD0(fD0fromDstarCuts,checkD0,checkD0bar);
999 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1000 }
1001
6a213b59 1002
2ff20727 1003 // remove the primary vertex (was used only for selection)
a25935a9 1004 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
b6e66d86 1005 the2Prong->UnsetOwnPrimaryVtx();
b6e66d86 1006 }
1007
2ff20727 1008 // get PID info from ESD
5770b08b 1009 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1010 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1011 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1012 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2ff20727 1013 Double_t esdpid[10];
1014 for(Int_t i=0;i<5;i++) {
1015 esdpid[i] = esdpid0[i];
1016 esdpid[5+i] = esdpid1[i];
6a213b59 1017 }
2ff20727 1018 the2Prong->SetPID(2,esdpid);
6a213b59 1019
6a213b59 1020 return the2Prong;
1021}
1022//----------------------------------------------------------------------------
1023AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
dcb444c9 1024 TObjArray *threeTrackArray,AliVEvent *event,
1025 AliAODVertex *secVert,Double_t dispersion,
c4cdf105 1026 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
6a213b59 1027 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1028 Bool_t &ok3Prong) const
1029{
1030 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1031 // reconstruction cuts
1032 // E.Bruna, F.Prino
1033
dcb444c9 1034
6a213b59 1035 ok3Prong=kFALSE;
673ac607 1036 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
dcb444c9 1037
1038 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1039 Double_t momentum[3];
6a213b59 1040
6a213b59 1041
1042 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
dcb444c9 1043 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
6a213b59 1044 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1045
dcb444c9 1046 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1047 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1048 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
6a213b59 1049 postrack1->GetPxPyPz(momentum);
1050 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1051 negtrack->GetPxPyPz(momentum);
1052 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1053 postrack2->GetPxPyPz(momentum);
1054 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1055
b82f6d67 1056 // invariant mass cut for D+, Ds, Lc
6a213b59 1057 Bool_t okMassCut=kFALSE;
1058 if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1059 if(!okMassCut) {
dcb444c9 1060 AliDebug(2," candidate didn't pass mass cut");
6a213b59 1061 return 0x0;
1062 }
1063
dcb444c9 1064 // primary vertex to be used by this candidate
2ff20727 1065 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
dcb444c9 1066 if(!primVertexAOD) return 0x0;
1067
1068 Double_t d0z0[2],covd0z0[3];
1069 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1070 d0[0]=d0z0[0];
1071 d0err[0] = TMath::Sqrt(covd0z0[0]);
1072 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1073 d0[1]=d0z0[0];
1074 d0err[1] = TMath::Sqrt(covd0z0[0]);
1075 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1076 d0[2]=d0z0[0];
1077 d0err[2] = TMath::Sqrt(covd0z0[0]);
6a213b59 1078
6a213b59 1079
1080 // create the object AliAODRecoDecayHF3Prong
dcb444c9 1081 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
6a213b59 1082 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
dcb444c9 1083 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]));
1084 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]));
1085 Short_t charge=(Short_t)(postrack1->Charge()*postrack2->Charge()*negtrack->Charge());
1086 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
6a213b59 1087 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
dcb444c9 1088 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1089 the3Prong->SetProngIDs(3,id);
6a213b59 1090
0a65d33f 1091 delete primVertexAOD; primVertexAOD=NULL;
6a213b59 1092
b82f6d67 1093 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1094 if(f3Prong) {
1095 ok3Prong = kFALSE;
1096 Int_t ok1,ok2;
f7c170b4 1097 Int_t dum1,dum2;
b82f6d67 1098 if(the3Prong->SelectDplus(fDplusCuts)) ok3Prong = kTRUE;
f7c170b4 1099 if(the3Prong->SelectDs(fDsCuts,ok1,ok2,dum1,dum2)) ok3Prong = kTRUE;
b82f6d67 1100 if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
1101 }
6a213b59 1102 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
dcb444c9 1103
a25935a9 1104 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
b6e66d86 1105 the3Prong->UnsetOwnPrimaryVtx();
b6e66d86 1106 }
dcb444c9 1107
2ff20727 1108 // get PID info from ESD
5770b08b 1109 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1110 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1111 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1112 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1113 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1114 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2ff20727 1115
1116 Double_t esdpid[15];
1117 for(Int_t i=0;i<5;i++) {
1118 esdpid[i] = esdpid0[i];
1119 esdpid[5+i] = esdpid1[i];
1120 esdpid[10+i] = esdpid2[i];
6a213b59 1121 }
2ff20727 1122 the3Prong->SetPID(3,esdpid);
6a213b59 1123
6a213b59 1124 return the3Prong;
1125}
1126//----------------------------------------------------------------------------
1127AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
dcb444c9 1128 TObjArray *fourTrackArray,AliVEvent *event,
721c0b8f 1129 AliAODVertex *secVert,
c4cdf105 1130 const AliAODVertex *vertexp1n1,
1131 const AliAODVertex *vertexp1n1p2,
721c0b8f 1132 Double_t dcap1n1,Double_t dcap1n2,
1133 Double_t dcap2n1,Double_t dcap2n2,
1134 Bool_t &ok4Prong) const
6a213b59 1135{
1136 // Make 4Prong candidates and check if they pass D0toKpipipi
1137 // reconstruction cuts
1138 // G.E.Bruno, R.Romita
1139
1140 ok4Prong=kFALSE;
673ac607 1141 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
6a213b59 1142
1143 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
6a213b59 1144
6a213b59 1145 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1146 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1147 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1148 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1149
dcb444c9 1150 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1151 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1152 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1153 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
6a213b59 1154
1155 Double_t momentum[3];
1156 postrack1->GetPxPyPz(momentum);
1157 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1158 negtrack1->GetPxPyPz(momentum);
1159 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1160 postrack2->GetPxPyPz(momentum);
1161 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1162 negtrack2->GetPxPyPz(momentum);
1163 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1164
1165 // invariant mass cut for rho or D0 (try to improve coding here..)
0a65d33f 1166 Bool_t okMassCut=kFALSE;
c4cdf105 1167 if(!okMassCut && TMath::Abs(fD0to4ProngsCuts[8])<1.e-13){ //no PID, to be implemented with PID
0a65d33f 1168 if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1169 }
1170 if(!okMassCut) {
1171 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1172 //printf(" candidate didn't pass mass cut\n");
1173 return 0x0;
1174 }
6a213b59 1175
dcb444c9 1176 // primary vertex to be used by this candidate
2ff20727 1177 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
dcb444c9 1178 if(!primVertexAOD) return 0x0;
1179
721c0b8f 1180 Double_t d0z0[2],covd0z0[3];
1181 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1182 d0[0]=d0z0[0];
1183 d0err[0] = TMath::Sqrt(covd0z0[0]);
1184 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1185 d0[1]=d0z0[0];
1186 d0err[1] = TMath::Sqrt(covd0z0[0]);
1187 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1188 d0[2]=d0z0[0];
1189 d0err[2] = TMath::Sqrt(covd0z0[0]);
1190 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1191 d0[3]=d0z0[0];
1192 d0err[3] = TMath::Sqrt(covd0z0[0]);
1193
6a213b59 1194
1195 // create the object AliAODRecoDecayHF4Prong
dcb444c9 1196 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
721c0b8f 1197 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
dcb444c9 1198 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 1199 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]));
1200 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 1201 Short_t charge=0;
721c0b8f 1202 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
6a213b59 1203 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
dcb444c9 1204 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1205 the4Prong->SetProngIDs(4,id);
6a213b59 1206
0a65d33f 1207 delete primVertexAOD; primVertexAOD=NULL;
6a213b59 1208
721c0b8f 1209 Int_t checkD0,checkD0bar;
1210 ok4Prong=the4Prong->SelectD0(fD0to4ProngsCuts,checkD0,checkD0bar);
6a213b59 1211
a25935a9 1212 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
b6e66d86 1213 the4Prong->UnsetOwnPrimaryVtx();
b6e66d86 1214 }
6a213b59 1215
721c0b8f 1216
6a213b59 1217 // get PID info from ESD
721c0b8f 1218 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1219 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1220 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1221 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1222 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1223 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1224 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1225 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
6a213b59 1226
1227 Double_t esdpid[20];
1228 for(Int_t i=0;i<5;i++) {
2ff20727 1229 esdpid[i] = esdpid0[i];
1230 esdpid[5+i] = esdpid1[i];
6a213b59 1231 esdpid[10+i] = esdpid2[i];
1232 esdpid[15+i] = esdpid3[i];
1233 }
1234 the4Prong->SetPID(4,esdpid);
721c0b8f 1235
6a213b59 1236 return the4Prong;
1237}
1238//-----------------------------------------------------------------------------
c4cdf105 1239AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2ff20727 1240 AliVEvent *event) const
6a213b59 1241{
2ff20727 1242 // Returns primary vertex to be used for this candidate
6a213b59 1243
dcb444c9 1244 AliESDVertex *vertexESD = 0;
1245 AliAODVertex *vertexAOD = 0;
1246
dcb444c9 1247
1248 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1249 // primary vertex from the input event
1250
1251 vertexESD = new AliESDVertex(*fV1);
1252
1253 } else {
1254 // primary vertex specific to this candidate
1255
2ff20727 1256 Int_t nTrks = trkArray->GetEntriesFast();
dcb444c9 1257 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1258
1259 if(fRecoPrimVtxSkippingTrks) {
1260 // recalculating the vertex
1261
1262 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1263 Float_t diamondcovxy[3];
1264 event->GetDiamondCovXY(diamondcovxy);
1265 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
8bb09289 1266 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
dcb444c9 1267 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1268 vertexer->SetVtxStart(diamond);
1269 delete diamond; diamond=NULL;
1270 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1271 vertexer->SetOnlyFitter();
1272 }
0a65d33f 1273 Int_t skipped[1000];
e18fbfa7 1274 Int_t nTrksToSkip=0,id;
1275 AliExternalTrackParam *t = 0;
dcb444c9 1276 for(Int_t i=0; i<nTrks; i++) {
e18fbfa7 1277 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1278 id = (Int_t)t->GetID();
1279 if(id<0) continue;
1280 skipped[nTrksToSkip++] = id;
dcb444c9 1281 }
0a65d33f 1282 // TEMPORARY FIX
1283 // For AOD, skip also tracks without covariance matrix
1284 if(fInputAOD) {
1285 Double_t covtest[21];
1286 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1287 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1288 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1289 id = (Int_t)t->GetID();
1290 if(id<0) continue;
1291 skipped[nTrksToSkip++] = id;
1292 }
1293 }
1294 }
1295 //
e18fbfa7 1296 vertexer->SetSkipTracks(nTrksToSkip,skipped);
dcb444c9 1297 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1298
1299 } else if(fRmTrksFromPrimVtx) {
1300 // removing the prongs tracks
1301
1302 TObjArray rmArray(nTrks);
1303 UShort_t *rmId = new UShort_t[nTrks];
1304 AliESDtrack *esdTrack = 0;
1305 AliESDtrack *t = 0;
1306 for(Int_t i=0; i<nTrks; i++) {
1307 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1308 esdTrack = new AliESDtrack(*t);
1309 rmArray.AddLast(esdTrack);
2ff20727 1310 if(esdTrack->GetID()>=0) {
1311 rmId[i]=(UShort_t)esdTrack->GetID();
1312 } else {
1313 rmId[i]=9999;
1314 }
dcb444c9 1315 }
1316 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1317 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1318 delete [] rmId; rmId=NULL;
1319 rmArray.Delete();
1320
6a213b59 1321 }
6a213b59 1322
dcb444c9 1323 if(!vertexESD) return vertexAOD;
1324 if(vertexESD->GetNContributors()<=0) {
1325 AliDebug(2,"vertexing failed");
1326 delete vertexESD; vertexESD=NULL;
1327 return vertexAOD;
6a213b59 1328 }
dcb444c9 1329
1330 delete vertexer; vertexer=NULL;
1331
6a213b59 1332 }
1333
dcb444c9 1334 // convert to AliAODVertex
1335 Double_t pos[3],cov[6],chi2perNDF;
1336 vertexESD->GetXYZ(pos); // position
1337 vertexESD->GetCovMatrix(cov); //covariance matrix
1338 chi2perNDF = vertexESD->GetChi2toNDF();
1339 delete vertexESD; vertexESD=NULL;
1340
1341 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
6a213b59 1342
dcb444c9 1343 return vertexAOD;
6a213b59 1344}
1345//-----------------------------------------------------------------------------
1346void AliAnalysisVertexingHF::PrintStatus() const {
1347 // Print parameters being used
1348
f8fa4595 1349 //printf("Preselections:\n");
1350 // fTrackFilter->Dump();
4d5c9633 1351 if(fSecVtxWithKF) {
1352 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1353 } else {
1354 printf("Secondary vertex with AliVertexerTracks\n");
1355 }
6a213b59 1356 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1357 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1358 if(fD0toKpi) {
f3014dc3 1359 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
6a213b59 1360 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1361 printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
1362 printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
1363 printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
1364 printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
1365 printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
1366 printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
1367 printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
1368 printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
1369 printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
1370 }
2ff20727 1371 if(fDstar) {
eee95e18 1372 if(fFindVertexForDstar) {
1373 printf(" Reconstruct a secondary vertex for the D*\n");
1374 } else {
1375 printf(" Assume the D* comes from the primary vertex\n");
1376 }
2ff20727 1377 printf(" |M-MD*| [GeV] < %f\n",fDstarCuts[0]);
1378 printf(" |M_Kpipi-M_Kpi-(MD*-MD0)| [GeV] < %f\n",fDstarCuts[1]);
1379 printf(" pTpisoft [GeV/c] > %f\n",fDstarCuts[2]);
1380 printf(" pTpisoft [GeV/c] < %f\n",fDstarCuts[3]);
1381 printf(" Theta(pisoft,D0plane) < %f\n",fDstarCuts[4]);
1382 printf("Reconstruct D*->D0pi candidates with cuts:\n");
1383 printf(" D0 from D* cuts:\n");
1384 printf(" |M-MD0| [GeV] < %f\n",fD0fromDstarCuts[0]);
1385 printf(" dca [cm] < %f\n",fD0fromDstarCuts[1]);
1386 printf(" cosThetaStar < %f\n",fD0fromDstarCuts[2]);
1387 printf(" pTK [GeV/c] > %f\n",fD0fromDstarCuts[3]);
1388 printf(" pTpi [GeV/c] > %f\n",fD0fromDstarCuts[4]);
1389 printf(" |d0K| [cm] < %f\n",fD0fromDstarCuts[5]);
1390 printf(" |d0pi| [cm] < %f\n",fD0fromDstarCuts[6]);
1391 printf(" d0d0 [cm^2] < %f\n",fD0fromDstarCuts[7]);
1392 printf(" cosThetaPoint > %f\n",fD0fromDstarCuts[8]);
1393 }
6a213b59 1394 if(fJPSItoEle) {
1395 printf("Reconstruct J/psi from B candidates with cuts:\n");
1396 printf(" |M-MJPSI| [GeV] < %f\n",fBtoJPSICuts[0]);
1397 printf(" dca [cm] < %f\n",fBtoJPSICuts[1]);
1398 printf(" cosThetaStar < %f\n",fBtoJPSICuts[2]);
1399 printf(" pTP [GeV/c] > %f\n",fBtoJPSICuts[3]);
1400 printf(" pTN [GeV/c] > %f\n",fBtoJPSICuts[4]);
1401 printf(" |d0P| [cm] < %f\n",fBtoJPSICuts[5]);
1402 printf(" |d0N| [cm] < %f\n",fBtoJPSICuts[6]);
1403 printf(" d0d0 [cm^2] < %f\n",fBtoJPSICuts[7]);
1404 printf(" cosThetaPoint > %f\n",fBtoJPSICuts[8]);
1405 }
1406 if(f3Prong) {
1407 printf("Reconstruct 3 prong candidates.\n");
4d5c9633 1408 printf(" D+->Kpipi cuts:\n");
6a213b59 1409 printf(" |M-MD+| [GeV] < %f\n",fDplusCuts[0]);
1410 printf(" pTK [GeV/c] > %f\n",fDplusCuts[1]);
1411 printf(" pTPi [GeV/c] > %f\n",fDplusCuts[2]);
1412 printf(" |d0K| [cm] > %f\n",fDplusCuts[3]);
1413 printf(" |d0Pi| [cm] > %f\n",fDplusCuts[4]);
1414 printf(" dist12 [cm] < %f\n",fDplusCuts[5]);
1415 printf(" sigmavert [cm] < %f\n",fDplusCuts[6]);
1416 printf(" dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1417 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1418 printf(" cosThetaPoint > %f\n",fDplusCuts[9]);
1419 printf(" Sum d0^2 [cm^2] > %f\n",fDplusCuts[10]);
1420 printf(" dca cut [cm] < %f\n",fDplusCuts[11]);
4d5c9633 1421 printf(" Ds->KKpi cuts:\n");
1422 printf(" |M-MDs| [GeV] < %f\n",fDsCuts[0]);
81679460 1423 printf(" pTK [GeV/c] > %f\n",fDsCuts[1]);
1424 printf(" pTPi [GeV/c] > %f\n",fDsCuts[2]);
1425 printf(" |d0K| [cm] > %f\n",fDsCuts[3]);
1426 printf(" |d0Pi| [cm] > %f\n",fDsCuts[4]);
1427 printf(" dist12 [cm] < %f\n",fDsCuts[5]);
1428 printf(" sigmavert [cm] < %f\n",fDsCuts[6]);
1429 printf(" dist prim-sec [cm] > %f\n",fDsCuts[7]);
1430 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
1431 printf(" cosThetaPoint > %f\n",fDsCuts[9]);
1432 printf(" Sum d0^2 [cm^2] > %f\n",fDsCuts[10]);
1433 printf(" dca cut [cm] < %f\n",fDsCuts[11]);
f7c170b4 1434 printf(" Inv. Mass phi [GeV] < %f\n",fDsCuts[12]);
1435 printf(" Inv. Mass K0* [GeV] < %f\n",fDsCuts[13]);
6ea608bf 1436 printf(" Lc->pKpi cuts:\n");
1437 printf(" |M-MLc| [GeV] < %f\n",fLcCuts[0]);
81679460 1438 printf(" pTP [GeV/c] > %f\n",fLcCuts[1]);
1439 printf(" pTPi and pTK [GeV/c] > %f\n",fLcCuts[2]);
1440 printf(" |d0P| [cm] > %f\n",fLcCuts[3]);
1441 printf(" |d0Pi| and |d0K| [cm] > %f\n",fLcCuts[4]);
1442 printf(" dist12 [cm] < %f\n",fLcCuts[5]);
1443 printf(" sigmavert [cm] < %f\n",fLcCuts[6]);
1444 printf(" dist prim-sec [cm] > %f\n",fLcCuts[7]);
1445 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
1446 printf(" cosThetaPoint > %f\n",fLcCuts[9]);
1447 printf(" Sum d0^2 [cm^2] > %f\n",fLcCuts[10]);
1448 printf(" dca cut [cm] < %f\n",fLcCuts[11]);
6a213b59 1449 }
1450
1451 return;
1452}
1453//-----------------------------------------------------------------------------
dcb444c9 1454AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
dcfa35b3 1455 Double_t &dispersion,Bool_t useTRefArray) const
dcb444c9 1456{
1457 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1458
1459 AliESDVertex *vertexESD = 0;
1460 AliAODVertex *vertexAOD = 0;
1461
1462 if(!fSecVtxWithKF) { // AliVertexerTracks
1463
1464 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1465 vertexer->SetVtxStart(fV1);
1466 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1467 delete vertexer; vertexer=NULL;
1468
1469 if(!vertexESD) return vertexAOD;
1470
1471 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1472 AliDebug(2,"vertexing failed");
1473 delete vertexESD; vertexESD=NULL;
1474 return vertexAOD;
1475 }
1476
1477 } else { // Kalman Filter vertexer (AliKFParticle)
1478
1479 AliKFParticle::SetField(fBzkG);
1480
a353d2e4 1481 AliKFVertex vertexKF;
dcb444c9 1482
1483 Int_t nTrks = trkArray->GetEntriesFast();
1484 for(Int_t i=0; i<nTrks; i++) {
1485 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1486 AliKFParticle daughterKF(*esdTrack,211);
1487 vertexKF.AddDaughter(daughterKF);
1488 }
a353d2e4 1489 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1490 vertexKF.CovarianceMatrix(),
1491 vertexKF.GetChi2(),
1492 vertexKF.GetNContributors());
dcb444c9 1493
1494 }
1495
1496 // convert to AliAODVertex
1497 Double_t pos[3],cov[6],chi2perNDF;
1498 vertexESD->GetXYZ(pos); // position
1499 vertexESD->GetCovMatrix(cov); //covariance matrix
1500 chi2perNDF = vertexESD->GetChi2toNDF();
1501 dispersion = vertexESD->GetDispersion();
1502 delete vertexESD; vertexESD=NULL;
1503
dcfa35b3 1504 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1505 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
dcb444c9 1506
1507 return vertexAOD;
1508}
1509//-----------------------------------------------------------------------------
6a213b59 1510Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1511 Int_t nprongs,
1512 Double_t *px,
1513 Double_t *py,
1514 Double_t *pz) const {
1515 // Check invariant mass cut
1516
1517 Short_t dummycharge=0;
1518 Double_t *dummyd0 = new Double_t[nprongs];
b056c5e3 1519 for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
6a213b59 1520 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
b6e66d86 1521 delete [] dummyd0; dummyd0=NULL;
6a213b59 1522
721c0b8f 1523 UInt_t pdg2[2],pdg3[3],pdg4[4];
6a213b59 1524 Double_t mPDG,minv;
1525
1526 Bool_t retval=kFALSE;
1527 switch (decay)
1528 {
1529 case 0: // D0->Kpi
1530 pdg2[0]=211; pdg2[1]=321;
1531 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1532 minv = rd->InvMass(nprongs,pdg2);
1533 if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
f3014dc3 1534 //if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
6a213b59 1535 pdg2[0]=321; pdg2[1]=211;
1536 minv = rd->InvMass(nprongs,pdg2);
1537 if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
f3014dc3 1538 //if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
6a213b59 1539 break;
1540 case 1: // JPSI->ee
1541 pdg2[0]=11; pdg2[1]=11;
1542 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1543 minv = rd->InvMass(nprongs,pdg2);
1544 if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1545 break;
1546 case 2: // D+->Kpipi
1547 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1548 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1549 minv = rd->InvMass(nprongs,pdg3);
1550 if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
b82f6d67 1551 // Ds+->KKpi
1552 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1553 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1554 minv = rd->InvMass(nprongs,pdg3);
1555 if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1556 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1557 minv = rd->InvMass(nprongs,pdg3);
1558 if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1559 // Lc->pKpi
1560 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1561 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1562 minv = rd->InvMass(nprongs,pdg3);
1563 if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1564 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1565 minv = rd->InvMass(nprongs,pdg3);
1566 if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
6a213b59 1567 break;
2ff20727 1568 case 3: // D*->D0pi
682639e9 1569 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2ff20727 1570 mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1571 minv = rd->InvMass(nprongs,pdg2);
1572 if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
1573 break;
721c0b8f 1574 case 4: // D0->K3pi without PID
1575 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1576 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1577 minv = rd->InvMass(nprongs,pdg4);
1578 if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0])
1579 retval=kTRUE;
1580 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1581 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1582 minv = rd->InvMass(nprongs,pdg4);
1583 if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1584 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1585 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1586 minv = rd->InvMass(nprongs,pdg4);
1587 if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1588 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1589 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1590 minv = rd->InvMass(nprongs,pdg4);
1591 if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1592 break;
6a213b59 1593 default:
2ff20727 1594 printf("SelectInvMass(): wrong decay selection\n");
6a213b59 1595 break;
1596 }
1597
b6e66d86 1598 delete rd; rd=NULL;
6a213b59 1599
1600 return retval;
1601}
1602//-----------------------------------------------------------------------------
c4cdf105 1603void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2ff20727 1604 TObjArray &seleTrksArray,Int_t &nSeleTrks,
a25935a9 1605 UChar_t *seleFlags,Int_t *evtNumber)
6a213b59 1606{
2ff20727 1607 // Apply single-track preselection.
1608 // Fill a TObjArray with selected tracks (for displaced vertices or
1609 // soft pion from D*). Selection flag stored in seleFlags.
dcb444c9 1610 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1611 // In case of AOD input, also fill fAODMap for track index<->ID
1612
1613 const AliVVertex *vprimary = event->GetPrimaryVertex();
1614
1615 if(fV1) { delete fV1; fV1=NULL; }
c8ab4e4f 1616 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
dcb444c9 1617
1618 Int_t nindices=0;
1619 UShort_t *indices = 0;
1620 Double_t pos[3],cov[6];
1621
1622 if(!fInputAOD) { // ESD
1623 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
1624 } else { // AOD
1625 vprimary->GetXYZ(pos);
1626 vprimary->GetCovarianceMatrix(cov);
1627 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1628 indices = new UShort_t[event->GetNumberOfTracks()];
c8ab4e4f 1629 fAODMapSize = 100000;
1630 fAODMap = new Int_t[fAODMapSize];
dcb444c9 1631 }
1632
1633
dcb444c9 1634 Int_t entries = (Int_t)event->GetNumberOfTracks();
2ff20727 1635 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
1636 nSeleTrks=0;
6a213b59 1637
dcb444c9 1638 // transfer ITS tracks from event to arrays
6a213b59 1639 for(Int_t i=0; i<entries; i++) {
a25935a9 1640 AliVTrack *track;
1641 track = (AliVTrack*)event->GetTrack(i);
dcb444c9 1642
f7c170b4 1643 // skip pure ITS SA tracks
ac325db6 1644 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
f7c170b4 1645
65dad444 1646 // TEMPORARY: check that the cov matrix is there
0ab1bd93 1647 Double_t covtest[21];
1648 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1649 //
65dad444 1650
dcb444c9 1651 if(fInputAOD) {
1652 AliAODTrack *aodt = (AliAODTrack*)track;
1653 if(aodt->GetUsedForPrimVtxFit()) {
1654 indices[nindices]=aodt->GetID(); nindices++;
1655 }
1656 fAODMap[(Int_t)aodt->GetID()] = i;
1657 }
6a213b59 1658
dcb444c9 1659 AliESDtrack *esdt = 0;
b6e66d86 1660
dcb444c9 1661 if(!fInputAOD) {
1662 esdt = (AliESDtrack*)track;
1663 } else {
1664 esdt = new AliESDtrack(track);
1665 }
6a213b59 1666
1667 // single track selection
2ff20727 1668 okDisplaced=kFALSE; okSoftPi=kFALSE;
a25935a9 1669 if(fMixEvent){
1670 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
1671 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
1672 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
1673 eventVtx->GetXYZ(vtxPos);
1674 vprimary->GetXYZ(primPos);
1675 eventVtx->GetCovarianceMatrix(primCov);
1676 for(Int_t ind=0;ind<3;ind++){
1677 trasl[ind]=vtxPos[ind]-primPos[ind];
1678 }
1679
1680 Bool_t isTransl=esdt->Translate(trasl,primCov);
1681 if(!isTransl) {
1682 delete esdt;
1683 esdt = NULL;
1684 continue;
1685 }
1686 }
1687
2ff20727 1688 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
1689 seleTrksArray.AddLast(esdt);
1690 seleFlags[nSeleTrks]=0;
1691 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
1692 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
1693 nSeleTrks++;
1694 } else {
dcb444c9 1695 if(fInputAOD) delete esdt;
1696 esdt = NULL;
1697 continue;
2ff20727 1698 }
6a213b59 1699
dcb444c9 1700 } // end loop on tracks
1701
1702 // primary vertex from AOD
1703 if(fInputAOD) {
1704 delete fV1; fV1=NULL;
1705 vprimary->GetXYZ(pos);
1706 vprimary->GetCovarianceMatrix(cov);
1707 Double_t chi2toNDF = vprimary->GetChi2perNDF();
1708 Int_t ncontr=nindices;
a6e0ebfe 1709 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
dcb444c9 1710 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
1711 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
1712 fV1->SetTitle(vprimary->GetTitle());
1713 fV1->SetIndices(nindices,indices);
1714 }
1715 if(indices) { delete [] indices; indices=NULL; }
1716
6a213b59 1717
1718 return;
1719}
1720//-----------------------------------------------------------------------------
1721void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1722 Double_t cut2,Double_t cut3,Double_t cut4,
1723 Double_t cut5,Double_t cut6,
1724 Double_t cut7,Double_t cut8)
1725{
1726 // Set the cuts for D0 selection
1727 fD0toKpiCuts[0] = cut0;
1728 fD0toKpiCuts[1] = cut1;
1729 fD0toKpiCuts[2] = cut2;
1730 fD0toKpiCuts[3] = cut3;
1731 fD0toKpiCuts[4] = cut4;
1732 fD0toKpiCuts[5] = cut5;
1733 fD0toKpiCuts[6] = cut6;
1734 fD0toKpiCuts[7] = cut7;
1735 fD0toKpiCuts[8] = cut8;
1736
1737 return;
1738}
1739//-----------------------------------------------------------------------------
1740void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
1741{
1742 // Set the cuts for D0 selection
1743
1744 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1745
1746 return;
1747}
1748//-----------------------------------------------------------------------------
2ff20727 1749void AliAnalysisVertexingHF::SetD0fromDstarCuts(Double_t cut0,Double_t cut1,
1750 Double_t cut2,Double_t cut3,Double_t cut4,
1751 Double_t cut5,Double_t cut6,
1752 Double_t cut7,Double_t cut8)
1753{
1754 // Set the cuts for D0 from D* selection
1755 fD0fromDstarCuts[0] = cut0;
1756 fD0fromDstarCuts[1] = cut1;
1757 fD0fromDstarCuts[2] = cut2;
1758 fD0fromDstarCuts[3] = cut3;
1759 fD0fromDstarCuts[4] = cut4;
1760 fD0fromDstarCuts[5] = cut5;
1761 fD0fromDstarCuts[6] = cut6;
1762 fD0fromDstarCuts[7] = cut7;
1763 fD0fromDstarCuts[8] = cut8;
1764
1765 return;
1766}
1767//-----------------------------------------------------------------------------
1768void AliAnalysisVertexingHF::SetD0fromDstarCuts(const Double_t cuts[9])
1769{
1770 // Set the cuts for D0 from D* selection
1771
1772 for(Int_t i=0; i<9; i++) fD0fromDstarCuts[i] = cuts[i];
1773
1774 return;
1775}
1776//-----------------------------------------------------------------------------
1777void AliAnalysisVertexingHF::SetDstarCuts(Double_t cut0,Double_t cut1,
1778 Double_t cut2,Double_t cut3,
1779 Double_t cut4)
1780{
1781 // Set the cuts for D* selection
1782 fDstarCuts[0] = cut0;
1783 fDstarCuts[1] = cut1;
1784 fDstarCuts[2] = cut2;
1785 fDstarCuts[3] = cut3;
1786 fDstarCuts[4] = cut4;
1787
1788 return;
1789}
1790//-----------------------------------------------------------------------------
1791void AliAnalysisVertexingHF::SetDstarCuts(const Double_t cuts[5])
1792{
1793 // Set the cuts for D* selection
1794
1795 for(Int_t i=0; i<5; i++) fDstarCuts[i] = cuts[i];
1796
1797 return;
1798}
1799//-----------------------------------------------------------------------------
6a213b59 1800void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1801 Double_t cut2,Double_t cut3,Double_t cut4,
1802 Double_t cut5,Double_t cut6,
1803 Double_t cut7,Double_t cut8)
1804{
1805 // Set the cuts for J/psi from B selection
1806 fBtoJPSICuts[0] = cut0;
1807 fBtoJPSICuts[1] = cut1;
1808 fBtoJPSICuts[2] = cut2;
1809 fBtoJPSICuts[3] = cut3;
1810 fBtoJPSICuts[4] = cut4;
1811 fBtoJPSICuts[5] = cut5;
1812 fBtoJPSICuts[6] = cut6;
1813 fBtoJPSICuts[7] = cut7;
1814 fBtoJPSICuts[8] = cut8;
1815
1816 return;
1817}
1818//-----------------------------------------------------------------------------
1819void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1820{
1821 // Set the cuts for J/psi from B selection
1822
1823 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1824
1825 return;
1826}
1827//-----------------------------------------------------------------------------
1828void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1829 Double_t cut2,Double_t cut3,Double_t cut4,
1830 Double_t cut5,Double_t cut6,
1831 Double_t cut7,Double_t cut8,
1832 Double_t cut9,Double_t cut10,Double_t cut11)
1833{
1834 // Set the cuts for Dplus->Kpipi selection
1835 fDplusCuts[0] = cut0;
1836 fDplusCuts[1] = cut1;
1837 fDplusCuts[2] = cut2;
1838 fDplusCuts[3] = cut3;
1839 fDplusCuts[4] = cut4;
1840 fDplusCuts[5] = cut5;
1841 fDplusCuts[6] = cut6;
1842 fDplusCuts[7] = cut7;
1843 fDplusCuts[8] = cut8;
1844 fDplusCuts[9] = cut9;
1845 fDplusCuts[10] = cut10;
1846 fDplusCuts[11] = cut11;
1847
1848 return;
1849}
1850//-----------------------------------------------------------------------------
1851void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1852{
b82f6d67 1853 // Set the cuts for Dplus->Kpipi selection
6a213b59 1854
1855 for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1856
1857 return;
1858}
1859//-----------------------------------------------------------------------------
6ea608bf 1860void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
f7c170b4 1861 Double_t cut2,Double_t cut3,
1862 Double_t cut4,Double_t cut5,
1863 Double_t cut6,Double_t cut7,
1864 Double_t cut8,Double_t cut9,
1865 Double_t cut10,Double_t cut11,
1866 Double_t cut12,Double_t cut13)
b82f6d67 1867{
1868 // Set the cuts for Ds->KKpi selection
1869 fDsCuts[0] = cut0;
6ea608bf 1870 fDsCuts[1] = cut1;
1871 fDsCuts[2] = cut2;
1872 fDsCuts[3] = cut3;
1873 fDsCuts[4] = cut4;
1874 fDsCuts[5] = cut5;
1875 fDsCuts[6] = cut6;
1876 fDsCuts[7] = cut7;
1877 fDsCuts[8] = cut8;
1878 fDsCuts[9] = cut9;
1879 fDsCuts[10] = cut10;
1880 fDsCuts[11] = cut11;
81679460 1881 fDsCuts[12] = cut12;
f7c170b4 1882 fDsCuts[13] = cut13;
b82f6d67 1883
1884 return;
1885}
1886//-----------------------------------------------------------------------------
81679460 1887void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13])
b82f6d67 1888{
1889 // Set the cuts for Ds->KKpi selection
1890
f7c170b4 1891 for(Int_t i=0; i<14; i++) fDsCuts[i] = cuts[i];
b82f6d67 1892
1893 return;
1894}
1895//-----------------------------------------------------------------------------
6ea608bf 1896void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
1897 Double_t cut2,Double_t cut3,Double_t cut4,
1898 Double_t cut5,Double_t cut6,
1899 Double_t cut7,Double_t cut8,
1900 Double_t cut9,Double_t cut10,Double_t cut11)
b82f6d67 1901{
1902 // Set the cuts for Lc->pKpi selection
1903 fLcCuts[0] = cut0;
6ea608bf 1904 fLcCuts[1] = cut1;
1905 fLcCuts[2] = cut2;
1906 fLcCuts[3] = cut3;
1907 fLcCuts[4] = cut4;
1908 fLcCuts[5] = cut5;
1909 fLcCuts[6] = cut6;
1910 fLcCuts[7] = cut7;
1911 fLcCuts[8] = cut8;
1912 fLcCuts[9] = cut9;
1913 fLcCuts[10] = cut10;
1914 fLcCuts[11] = cut11;
b82f6d67 1915
1916 return;
1917}
1918//-----------------------------------------------------------------------------
6ea608bf 1919void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
b82f6d67 1920{
1921 // Set the cuts for Lc->pKpi selection
1922
6ea608bf 1923 for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
b82f6d67 1924
1925 return;
1926}
1927//-----------------------------------------------------------------------------
721c0b8f 1928void AliAnalysisVertexingHF::SetD0to4ProngsCuts(Double_t cut0,Double_t cut1,
1929 Double_t cut2,Double_t cut3,Double_t cut4,
1930 Double_t cut5,Double_t cut6,
1931 Double_t cut7,Double_t cut8)
1932{
1933 // Set the cuts for D0->Kpipipi selection
1934
1935 fD0to4ProngsCuts[0] = cut0;
1936 fD0to4ProngsCuts[1] = cut1;
1937 fD0to4ProngsCuts[2] = cut2;
1938 fD0to4ProngsCuts[3] = cut3;
1939 fD0to4ProngsCuts[4] = cut4;
1940 fD0to4ProngsCuts[5] = cut5;
1941 fD0to4ProngsCuts[6] = cut6;
1942 fD0to4ProngsCuts[7] = cut7;
1943 fD0to4ProngsCuts[8] = cut8;
1944
1945 return;
1946}
1947//-----------------------------------------------------------------------------
1948void AliAnalysisVertexingHF::SetD0to4ProngsCuts(const Double_t cuts[9])
1949{
1950 // Set the cuts for D0->Kpipipi selection
1951
1952 for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i] = cuts[i];
1953
1954 return;
1955}
1956//-----------------------------------------------------------------------------
2ff20727 1957Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
1958 Bool_t &okDisplaced,Bool_t &okSoftPi) const
6a213b59 1959{
1960 // Check if track passes some kinematical cuts
6a213b59 1961
f8fa4595 1962 // this is needed to store the impact parameters
1963 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
1964
2ff20727 1965 UInt_t selectInfo;
f8fa4595 1966 //
2ff20727 1967 // Track selection, displaced tracks
1968 selectInfo = 0;
f8fa4595 1969 if(fTrackFilter) {
1970 selectInfo = fTrackFilter->IsSelected(trk);
6a213b59 1971 }
2ff20727 1972 if(selectInfo) okDisplaced=kTRUE;
1973 // Track selection, soft pions
1974 selectInfo = 0;
1975 if(fDstar && fTrackFilterSoftPi) {
1976 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
1977 }
1978 if(selectInfo) okSoftPi=kTRUE;
6a213b59 1979
2ff20727 1980 if(okDisplaced || okSoftPi) return kTRUE;
f8fa4595 1981
2ff20727 1982 return kFALSE;
6a213b59 1983}
1984//-----------------------------------------------------------------------------