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