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