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