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