1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //----------------------------------------------------------------------------
19 // Implementation of the heavy-flavour vertexing analysis class
20 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
21 // To be used as a task of AliAnalysisManager by means of the interface
22 // class AliAnalysisTaskSEVertexingHF.
23 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
25 // Contact: andrea.dainese@pd.infn.it
26 // Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
27 // F.Prino, R.Romita, X.M.Zhang
28 //----------------------------------------------------------------------------
29 #include <Riostream.h>
31 #include <TDatabasePDG.h>
35 #include "AliVEvent.h"
36 #include "AliVVertex.h"
37 #include "AliVTrack.h"
38 #include "AliVertexerTracks.h"
39 #include "AliKFVertex.h"
40 #include "AliESDEvent.h"
41 #include "AliESDVertex.h"
42 #include "AliExternalTrackParam.h"
43 #include "AliNeutralTrackParam.h"
44 #include "AliESDtrack.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliAODEvent.h"
47 #include "AliPIDResponse.h"
48 #include "AliAODRecoDecay.h"
49 #include "AliAODRecoDecayHF.h"
50 #include "AliAODRecoDecayHF2Prong.h"
51 #include "AliAODRecoDecayHF3Prong.h"
52 #include "AliAODRecoDecayHF4Prong.h"
53 #include "AliAODRecoCascadeHF.h"
54 #include "AliRDHFCutsD0toKpi.h"
55 #include "AliRDHFCutsJpsitoee.h"
56 #include "AliRDHFCutsDplustoKpipi.h"
57 #include "AliRDHFCutsDstoKKpi.h"
58 #include "AliRDHFCutsLctopKpi.h"
59 #include "AliRDHFCutsLctoV0.h"
60 #include "AliRDHFCutsD0toKpipipi.h"
61 #include "AliRDHFCutsDStartoKpipi.h"
62 #include "AliAnalysisFilter.h"
63 #include "AliAnalysisVertexingHF.h"
64 #include "AliMixedEvent.h"
67 #include "AliCodeTimer.h"
70 ClassImp(AliAnalysisVertexingHF)
72 //----------------------------------------------------------------------------
73 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
79 fSecVtxWithKF(kFALSE),
80 fRecoPrimVtxSkippingTrks(kFALSE),
81 fRmTrksFromPrimVtx(kFALSE),
90 fLikeSign3prong(kFALSE),
93 fUseKaonPIDfor3Prong(kFALSE),
96 fUseKaonPIDforDs(kFALSE),
99 fUseTPCPIDOnlyIfNoTOF(kFALSE),
100 fMaxMomForTPCPid(1.),
101 fnSigmaTPCPionLow(5.),
102 fnSigmaTPCPionHi(5.),
103 fnSigmaTOFPionLow(5.),
104 fnSigmaTOFPionHi(5.),
105 fnSigmaTPCKaonLow(5.),
106 fnSigmaTPCKaonHi(5.),
107 fnSigmaTOFKaonLow(5.),
108 fnSigmaTOFKaonHi(5.),
109 fnSigmaTPCProtonLow(5.),
110 fnSigmaTPCProtonHi(5.),
111 fnSigmaTOFProtonLow(5.),
112 fnSigmaTOFProtonHi(5.),
113 fMaxCentPercentileForTightCuts(-9999),
115 fTrackFilter2prongCentral(0x0),
116 fTrackFilter3prongCentral(0x0),
117 fTrackFilterSoftPi(0x0),
120 fCutsDplustoKpipi(0x0),
124 fCutsD0toKpipipi(0x0),
125 fCutsDStartoKpipi(0x0),
127 fFindVertexForDstar(kTRUE),
128 fFindVertexForCascades(kTRUE),
129 fV0TypeForCascadeVertex(0),
130 fMassCutBeforeVertexing(kFALSE),
134 fOKInvMassD0(kFALSE),
135 fOKInvMassJpsi(kFALSE),
136 fOKInvMassDplus(kFALSE),
137 fOKInvMassDs(kFALSE),
138 fOKInvMassLc(kFALSE),
139 fOKInvMassDstar(kFALSE),
140 fOKInvMassD0to4p(kFALSE),
141 fOKInvMassLctoV0(kFALSE),
151 // Default constructor
153 Double_t d02[2]={0.,0.};
154 Double_t d03[3]={0.,0.,0.};
155 Double_t d04[4]={0.,0.,0.,0.};
156 fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
157 fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
158 fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
161 //--------------------------------------------------------------------------
162 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
164 fInputAOD(source.fInputAOD),
165 fAODMapSize(source.fAODMapSize),
166 fAODMap(source.fAODMap),
167 fVertexerTracks(source.fVertexerTracks),
169 fSecVtxWithKF(source.fSecVtxWithKF),
170 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
171 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
173 fD0toKpi(source.fD0toKpi),
174 fJPSItoEle(source.fJPSItoEle),
175 f3Prong(source.f3Prong),
176 f4Prong(source.f4Prong),
177 fDstar(source.fDstar),
178 fCascades(source.fCascades),
179 fLikeSign(source.fLikeSign),
180 fLikeSign3prong(source.fLikeSign3prong),
181 fMixEvent(source.fMixEvent),
182 fPidResponse(source.fPidResponse),
183 fUseKaonPIDfor3Prong(source.fUseKaonPIDfor3Prong),
184 fUsePIDforLc(source.fUsePIDforLc),
185 fUsePIDforLc2V0(source.fUsePIDforLc2V0),
186 fUseKaonPIDforDs(source.fUseKaonPIDforDs),
187 fUseTPCPID(source.fUseTPCPID),
188 fUseTOFPID(source.fUseTOFPID),
189 fUseTPCPIDOnlyIfNoTOF(source.fUseTPCPIDOnlyIfNoTOF),
190 fMaxMomForTPCPid(source.fMaxMomForTPCPid),
191 fnSigmaTPCPionLow(source.fnSigmaTPCPionLow),
192 fnSigmaTPCPionHi(source.fnSigmaTPCPionHi),
193 fnSigmaTOFPionLow(source.fnSigmaTOFPionLow),
194 fnSigmaTOFPionHi(source.fnSigmaTOFPionHi),
195 fnSigmaTPCKaonLow(source.fnSigmaTPCKaonLow),
196 fnSigmaTPCKaonHi(source.fnSigmaTPCKaonHi),
197 fnSigmaTOFKaonLow(source.fnSigmaTOFKaonLow),
198 fnSigmaTOFKaonHi(source.fnSigmaTOFKaonHi),
199 fnSigmaTPCProtonLow(source.fnSigmaTPCProtonLow),
200 fnSigmaTPCProtonHi(source.fnSigmaTPCProtonHi),
201 fnSigmaTOFProtonLow(source.fnSigmaTOFProtonLow),
202 fnSigmaTOFProtonHi(source.fnSigmaTOFProtonHi),
203 fMaxCentPercentileForTightCuts(source.fMaxCentPercentileForTightCuts),
204 fTrackFilter(source.fTrackFilter),
205 fTrackFilter2prongCentral(source.fTrackFilter2prongCentral),
206 fTrackFilter3prongCentral(source.fTrackFilter3prongCentral),
207 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
208 fCutsD0toKpi(source.fCutsD0toKpi),
209 fCutsJpsitoee(source.fCutsJpsitoee),
210 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
211 fCutsDstoKKpi(source.fCutsDstoKKpi),
212 fCutsLctopKpi(source.fCutsLctopKpi),
213 fCutsLctoV0(source.fCutsLctoV0),
214 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
215 fCutsDStartoKpipi(source.fCutsDStartoKpipi),
216 fListOfCuts(source.fListOfCuts),
217 fFindVertexForDstar(source.fFindVertexForDstar),
218 fFindVertexForCascades(source.fFindVertexForCascades),
219 fV0TypeForCascadeVertex(source.fV0TypeForCascadeVertex),
220 fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
221 fMassCalc2(source.fMassCalc2),
222 fMassCalc3(source.fMassCalc3),
223 fMassCalc4(source.fMassCalc4),
224 fOKInvMassD0(source.fOKInvMassD0),
225 fOKInvMassJpsi(source.fOKInvMassJpsi),
226 fOKInvMassDplus(source.fOKInvMassDplus),
227 fOKInvMassDs(source.fOKInvMassDs),
228 fOKInvMassLc(source.fOKInvMassLc),
229 fOKInvMassDstar(source.fOKInvMassDstar),
230 fOKInvMassD0to4p(source.fOKInvMassD0to4p),
231 fOKInvMassLctoV0(source.fOKInvMassLctoV0),
234 fMassDzero(source.fMassDzero),
235 fMassDplus(source.fMassDplus),
236 fMassDs(source.fMassDs),
237 fMassLambdaC(source.fMassLambdaC),
238 fMassDstar(source.fMassDstar),
239 fMassJpsi(source.fMassJpsi)
245 //--------------------------------------------------------------------------
246 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
249 // assignment operator
251 if(&source == this) return *this;
252 fInputAOD = source.fInputAOD;
253 fAODMapSize = source.fAODMapSize;
254 fVertexerTracks = source.fVertexerTracks;
255 fBzkG = source.fBzkG;
256 fSecVtxWithKF = source.fSecVtxWithKF;
257 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
258 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
260 fD0toKpi = source.fD0toKpi;
261 fJPSItoEle = source.fJPSItoEle;
262 f3Prong = source.f3Prong;
263 f4Prong = source.f4Prong;
264 fDstar = source.fDstar;
265 fCascades = source.fCascades;
266 fLikeSign = source.fLikeSign;
267 fLikeSign3prong = source.fLikeSign3prong;
268 fMixEvent = source.fMixEvent;
269 fPidResponse = source.fPidResponse;
270 fUseKaonPIDfor3Prong = source.fUseKaonPIDfor3Prong;
271 fUsePIDforLc = source.fUsePIDforLc;
272 fUsePIDforLc2V0 = source.fUsePIDforLc2V0;
273 fUseKaonPIDforDs = source.fUseKaonPIDforDs;
274 fUseTPCPID = source.fUseTPCPID;
275 fUseTOFPID = source.fUseTOFPID;
276 fUseTPCPIDOnlyIfNoTOF = source.fUseTPCPIDOnlyIfNoTOF;
277 fMaxMomForTPCPid = source.fMaxMomForTPCPid;
278 fnSigmaTPCPionLow = source.fnSigmaTPCPionLow;
279 fnSigmaTPCPionHi = source.fnSigmaTPCPionHi;
280 fnSigmaTOFPionLow = source.fnSigmaTOFPionLow;
281 fnSigmaTOFPionHi = source.fnSigmaTOFPionHi;
282 fnSigmaTPCKaonLow = source.fnSigmaTPCKaonLow;
283 fnSigmaTPCKaonHi = source.fnSigmaTPCKaonHi;
284 fnSigmaTOFKaonLow = source.fnSigmaTOFKaonLow;
285 fnSigmaTOFKaonHi = source.fnSigmaTOFKaonHi;
286 fnSigmaTPCProtonLow = source.fnSigmaTPCProtonLow;
287 fnSigmaTPCProtonHi = source.fnSigmaTPCProtonHi;
288 fnSigmaTOFProtonLow = source.fnSigmaTOFProtonLow;
289 fnSigmaTOFProtonHi = source.fnSigmaTOFProtonHi;
290 fnSigmaTOFKaonLow = source.fnSigmaTOFKaonLow;
291 fnSigmaTOFKaonHi = source.fnSigmaTOFKaonHi;
292 fMaxCentPercentileForTightCuts = source.fMaxCentPercentileForTightCuts;
293 fTrackFilter = source.fTrackFilter;
294 fTrackFilter2prongCentral = source.fTrackFilter2prongCentral;
295 fTrackFilter3prongCentral = source.fTrackFilter3prongCentral;
296 fTrackFilterSoftPi = source.fTrackFilterSoftPi;
297 fCutsD0toKpi = source.fCutsD0toKpi;
298 fCutsJpsitoee = source.fCutsJpsitoee;
299 fCutsDplustoKpipi = source.fCutsDplustoKpipi;
300 fCutsDstoKKpi = source.fCutsDstoKKpi;
301 fCutsLctopKpi = source.fCutsLctopKpi;
302 fCutsLctoV0 = source.fCutsLctoV0;
303 fCutsD0toKpipipi = source.fCutsD0toKpipipi;
304 fCutsDStartoKpipi = source.fCutsDStartoKpipi;
305 fListOfCuts = source.fListOfCuts;
306 fFindVertexForDstar = source.fFindVertexForDstar;
307 fFindVertexForCascades = source.fFindVertexForCascades;
308 fV0TypeForCascadeVertex = source.fV0TypeForCascadeVertex;
309 fMassCutBeforeVertexing = source.fMassCutBeforeVertexing;
310 fMassCalc2 = source.fMassCalc2;
311 fMassCalc3 = source.fMassCalc3;
312 fMassCalc4 = source.fMassCalc4;
313 fOKInvMassD0 = source.fOKInvMassD0;
314 fOKInvMassJpsi = source.fOKInvMassJpsi;
315 fOKInvMassDplus = source.fOKInvMassDplus;
316 fOKInvMassDs = source.fOKInvMassDs;
317 fOKInvMassLc = source.fOKInvMassLc;
318 fOKInvMassDstar = source.fOKInvMassDstar;
319 fOKInvMassD0to4p = source.fOKInvMassD0to4p;
320 fOKInvMassLctoV0 = source.fOKInvMassLctoV0;
321 fMassDzero = source.fMassDzero;
322 fMassDplus = source.fMassDplus;
323 fMassDs = source.fMassDs;
324 fMassLambdaC = source.fMassLambdaC;
325 fMassDstar = source.fMassDstar;
326 fMassJpsi = source.fMassJpsi;
330 //----------------------------------------------------------------------------
331 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
333 if(fV1) { delete fV1; fV1=0; }
334 delete fVertexerTracks;
335 if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
336 if(fTrackFilter2prongCentral) { delete fTrackFilter2prongCentral; fTrackFilter2prongCentral=0; }
337 if(fTrackFilter3prongCentral) { delete fTrackFilter3prongCentral; fTrackFilter3prongCentral=0; }
338 if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
339 if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
340 if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
341 if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
342 if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
343 if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
344 if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
345 if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
346 if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
347 if(fAODMap) { delete [] fAODMap; fAODMap=0; }
348 if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
349 if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
350 if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
352 //----------------------------------------------------------------------------
353 TList *AliAnalysisVertexingHF::FillListOfCuts() {
354 // Fill list of analysis cuts
356 TList *list = new TList();
358 list->SetName("ListOfCuts");
361 AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
362 list->Add(cutsD0toKpi);
365 AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
366 list->Add(cutsJpsitoee);
368 if(fCutsDplustoKpipi) {
369 AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
370 list->Add(cutsDplustoKpipi);
373 AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
374 list->Add(cutsDstoKKpi);
377 AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
378 list->Add(cutsLctopKpi);
381 AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
382 list->Add(cutsLctoV0);
384 if(fCutsD0toKpipipi) {
385 AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
386 list->Add(cutsD0toKpipipi);
388 if(fCutsDStartoKpipi) {
389 AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
390 list->Add(cutsDStartoKpipi);
393 //___ Check consitstency of cuts between vertexer and analysis tasks
394 Bool_t bCutsOk = CheckCutsConsistency();
395 if (bCutsOk == kFALSE) {AliFatal("AliAnalysisVertexingHF::FillListOfCuts vertexing and the analysis task cuts are not consistent!");}
397 // keep a pointer to the list
402 //----------------------------------------------------------------------------
403 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
404 TClonesArray *aodVerticesHFTClArr,
405 TClonesArray *aodD0toKpiTClArr,
406 TClonesArray *aodJPSItoEleTClArr,
407 TClonesArray *aodCharm3ProngTClArr,
408 TClonesArray *aodCharm4ProngTClArr,
409 TClonesArray *aodDstarTClArr,
410 TClonesArray *aodCascadesTClArr,
411 TClonesArray *aodLikeSign2ProngTClArr,
412 TClonesArray *aodLikeSign3ProngTClArr)
414 // Find heavy-flavour vertex candidates
416 // Output: AOD (additional branches added)
417 //AliCodeTimerAuto("",0);
420 TString evtype = event->IsA()->GetName();
421 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
422 } // if we do mixing AliVEvent is a AliMixedEvent
425 AliDebug(2,"Creating HF candidates from AOD");
427 AliDebug(2,"Creating HF candidates from ESD");
430 if(!aodVerticesHFTClArr) {
431 printf("ERROR: no aodVerticesHFTClArr");
434 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
435 printf("ERROR: no aodD0toKpiTClArr");
438 if(fJPSItoEle && !aodJPSItoEleTClArr) {
439 printf("ERROR: no aodJPSItoEleTClArr");
442 if(f3Prong && !aodCharm3ProngTClArr) {
443 printf("ERROR: no aodCharm3ProngTClArr");
446 if(f4Prong && !aodCharm4ProngTClArr) {
447 printf("ERROR: no aodCharm4ProngTClArr");
450 if(fDstar && !aodDstarTClArr) {
451 printf("ERROR: no aodDstarTClArr");
454 if(fCascades && !aodCascadesTClArr){
455 printf("ERROR: no aodCascadesTClArr ");
458 if(fLikeSign && !aodLikeSign2ProngTClArr) {
459 printf("ERROR: no aodLikeSign2ProngTClArr");
462 if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
463 printf("ERROR: no aodLikeSign3ProngTClArr");
467 // delete candidates from previous event and create references
468 Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
469 aodVerticesHFTClArr->Delete();
470 iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
471 TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
472 if(fD0toKpi || fDstar) {
473 aodD0toKpiTClArr->Delete();
474 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
477 aodJPSItoEleTClArr->Delete();
478 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
481 aodCharm3ProngTClArr->Delete();
482 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
485 aodCharm4ProngTClArr->Delete();
486 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
489 aodDstarTClArr->Delete();
490 iDstar = aodDstarTClArr->GetEntriesFast();
493 aodCascadesTClArr->Delete();
494 iCascades = aodCascadesTClArr->GetEntriesFast();
497 aodLikeSign2ProngTClArr->Delete();
498 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
500 if(fLikeSign3prong && f3Prong) {
501 aodLikeSign3ProngTClArr->Delete();
502 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
505 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
506 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
507 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
508 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
509 TClonesArray &aodDstarRef = *aodDstarTClArr;
510 TClonesArray &aodCascadesRef = *aodCascadesTClArr;
511 TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
512 TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
515 AliAODRecoDecayHF2Prong *io2Prong = 0;
516 AliAODRecoDecayHF3Prong *io3Prong = 0;
517 AliAODRecoDecayHF4Prong *io4Prong = 0;
518 AliAODRecoCascadeHF *ioCascade = 0;
520 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
521 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
522 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
523 Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
524 Bool_t okCascades=kFALSE;
525 AliESDtrack *postrack1 = 0;
526 AliESDtrack *postrack2 = 0;
527 AliESDtrack *negtrack1 = 0;
528 AliESDtrack *negtrack2 = 0;
529 AliESDtrack *trackPi = 0;
530 Double_t mompos1[3],mompos2[3],momneg1[3],momneg2[3];
531 // AliESDtrack *posV0track = 0;
532 // AliESDtrack *negV0track = 0;
533 Float_t dcaMax = fCutsD0toKpi->GetDCACut();
534 if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
535 if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
536 if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
537 if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
538 if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
539 if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
541 AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
545 fBzkG = (Double_t)event->GetMagneticField();
546 if(!fVertexerTracks){
547 fVertexerTracks=new AliVertexerTracks(fBzkG);
549 Double_t oldField=fVertexerTracks->GetFieldkG();
550 if(oldField!=fBzkG) fVertexerTracks->SetFieldkG(fBzkG);
553 trkEntries = (Int_t)event->GetNumberOfTracks();
554 AliDebug(1,Form(" Number of tracks: %d",trkEntries));
555 fnTrksTotal += trkEntries;
557 nv0 = (Int_t)event->GetNumberOfV0s();
558 AliDebug(1,Form(" Number of V0s: %d",nv0));
560 if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
561 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
565 // event selection + PID configuration
566 if(!fCutsD0toKpi->IsEventSelected(event)) return;
567 if(fCutsJpsitoee) fCutsJpsitoee->SetupPID(event);
568 if(fCutsDplustoKpipi) fCutsDplustoKpipi->SetupPID(event);
569 if(fCutsDstoKKpi) fCutsDstoKKpi->SetupPID(event);
570 if(fCutsLctopKpi) fCutsLctopKpi->SetupPID(event);
571 if(fCutsLctoV0) fCutsLctoV0->SetupPID(event);
572 if(fCutsD0toKpipipi) fCutsD0toKpipipi->SetupPID(event);
573 if(fCutsDStartoKpipi) fCutsDStartoKpipi->SetupPID(event);
575 // call function that applies sigle-track selection,
576 // for displaced tracks and soft pions (both charges) for D*,
577 // and retrieves primary vertex
578 TObjArray seleTrksArray(trkEntries);
579 TObjArray tracksAtVertex(trkEntries);
580 UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi, bit 2: 3 prong, bits 3-4-5: for PID
582 Int_t *evtNumber = new Int_t[trkEntries];
583 SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
585 AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
586 fnSeleTrksTotal += nSeleTrks;
589 TObjArray *twoTrackArray1 = new TObjArray(2);
590 TObjArray *twoTrackArray2 = new TObjArray(2);
591 TObjArray *twoTrackArrayV0 = new TObjArray(2);
592 TObjArray *twoTrackArrayCasc = new TObjArray(2);
593 TObjArray *threeTrackArray = new TObjArray(3);
594 TObjArray *fourTrackArray = new TObjArray(4);
597 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
599 AliAODRecoDecayHF *rd = 0;
600 AliAODRecoCascadeHF *rc = 0;
604 Bool_t massCutOK=kTRUE;
606 // LOOP ON POSITIVE TRACKS
607 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
609 //if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
610 //if(iTrkP1%1==0) printf(" 1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);
612 // get track from tracks array
613 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
614 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
615 postrack1->GetPxPyPz(mompos1);
617 // Make cascades with V0+track
621 for(iv0=0; iv0<nv0; iv0++){
623 //AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
625 if ( fUsePIDforLc2V0 && !TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) ) continue; //clm
629 v0 = ((AliAODEvent*)event)->GetV0(iv0);
631 esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
633 if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
634 (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
636 if ( v0 && ((v0->GetOnFlyStatus() == kTRUE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOfflineV0s) ||
637 (v0->GetOnFlyStatus() == kFALSE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOnTheFlyV0s)) ) continue;
639 if ( esdV0 && ((esdV0->GetOnFlyStatus() == kTRUE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOfflineV0s) ||
640 ( esdV0->GetOnFlyStatus() == kFALSE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOnTheFlyV0s)) ) continue;
642 // Get the tracks that form the V0
643 // ( parameters at primary vertex )
644 // and define an AliExternalTrackParam out of them
645 AliExternalTrackParam * posV0track;
646 AliExternalTrackParam * negV0track;
649 AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
650 AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
651 if( !posVV0track || !negVV0track ) continue;
653 // Apply some basic V0 daughter criteria
655 // bachelor must not be a v0-track
656 if (posVV0track->GetID() == postrack1->GetID() ||
657 negVV0track->GetID() == postrack1->GetID()) continue;
658 // reject like-sign v0
659 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
660 // avoid ghost TPC tracks
661 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
662 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
663 // Get AliExternalTrackParam out of the AliAODTracks
664 Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
665 posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
666 posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
667 posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
668 negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
669 negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
670 negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
672 AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
673 AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
674 if( !posVV0track || !negVV0track ) continue;
676 // Apply some basic V0 daughter criteria
678 // bachelor must not be a v0-track
679 if (posVV0track->GetID() == postrack1->GetID() ||
680 negVV0track->GetID() == postrack1->GetID()) continue;
681 // reject like-sign v0
682 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
683 // avoid ghost TPC tracks
684 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
685 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
686 // reject kinks (only necessary on AliESDtracks)
687 if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
688 // Get AliExternalTrackParam out of the AliESDtracks
689 posV0track = new AliExternalTrackParam(*posVV0track);
690 negV0track = new AliExternalTrackParam(*negVV0track);
692 // Define the AODv0 from ESDv0 if reading ESDs
693 v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
695 if( !posV0track || !negV0track ){
696 AliDebug(1,Form(" Couldn't get the V0 daughters"));
700 // fill in the v0 two-external-track-param array
701 twoTrackArrayV0->AddAt(posV0track,0);
702 twoTrackArrayV0->AddAt(negV0track,1);
705 dcaV0 = v0->DcaV0Daughters();
707 // Define the V0 (neutral) track
708 AliNeutralTrackParam *trackV0=NULL;
710 const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
711 if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
713 Double_t xyz[3], pxpypz[3];
715 esdV0->PxPyPz(pxpypz);
716 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
717 trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
721 // Fill in the object array to create the cascade
722 twoTrackArrayCasc->AddAt(postrack1,0);
723 twoTrackArrayCasc->AddAt(trackV0,1);
724 // Compute the cascade vertex
725 AliAODVertex *vertexCasc = 0;
726 if(fFindVertexForCascades) {
727 // DCA between the two tracks
728 dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
730 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
732 // assume Cascade decays at the primary vertex
733 Double_t pos[3],cov[6],chi2perNDF;
735 fV1->GetCovMatrix(cov);
736 chi2perNDF = fV1->GetChi2toNDF();
737 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
741 delete posV0track; posV0track=NULL;
742 delete negV0track; negV0track=NULL;
743 delete trackV0; trackV0=NULL;
744 if(!fInputAOD) {delete v0; v0=NULL;}
745 twoTrackArrayV0->Clear();
746 twoTrackArrayCasc->Clear();
750 // Create and store the Cascade if passed the cuts
751 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
752 if(okCascades && ioCascade) {
753 //AliDebug(1,Form("Storing a cascade object... "));
754 // add the vertex and the cascade to the AOD
755 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
756 rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
757 rc->SetSecondaryVtx(vCasc);
758 vCasc->SetParent(rc);
759 rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
760 if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
761 AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
762 vCasc->AddDaughter(v0); // fill the 2prong V0
766 delete posV0track; posV0track=NULL;
767 delete negV0track; negV0track=NULL;
768 delete trackV0; trackV0=NULL;
769 twoTrackArrayV0->Clear();
770 twoTrackArrayCasc->Clear();
771 if(ioCascade) { delete ioCascade; ioCascade=NULL; }
772 if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
773 if(!fInputAOD) {delete v0; v0=NULL;}
775 } // end loop on V0's
778 // If there is less than 2 particles continue
780 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
784 if(postrack1->Charge()<0 && !fLikeSign) continue;
786 // LOOP ON NEGATIVE TRACKS
787 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
789 //if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
790 //if(iTrkN1%1==0) printf(" 1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);
792 if(iTrkN1==iTrkP1) continue;
794 // get track from tracks array
795 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
797 if(negtrack1->Charge()>0 && !fLikeSign) continue;
799 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
802 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
805 if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
806 isLikeSign2Prong=kTRUE;
807 if(!fLikeSign) continue;
808 if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
809 } else { // unlike-sign
810 isLikeSign2Prong=kFALSE;
811 if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
813 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
818 // back to primary vertex
819 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
820 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
821 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
822 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
823 negtrack1->GetPxPyPz(momneg1);
825 // DCA between the two tracks
826 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
827 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
830 twoTrackArray1->AddAt(postrack1,0);
831 twoTrackArray1->AddAt(negtrack1,1);
832 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
834 twoTrackArray1->Clear();
840 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
842 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
844 if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
845 // add the vertex and the decay to the AOD
846 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
847 if(!isLikeSign2Prong) {
849 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
850 rd->SetSecondaryVtx(v2Prong);
851 v2Prong->SetParent(rd);
852 AddRefs(v2Prong,rd,event,twoTrackArray1);
855 rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
856 rd->SetSecondaryVtx(v2Prong);
857 if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
858 AddRefs(v2Prong,rd,event,twoTrackArray1);
860 } else { // isLikeSign2Prong
861 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
862 rd->SetSecondaryVtx(v2Prong);
863 v2Prong->SetParent(rd);
864 AddRefs(v2Prong,rd,event,twoTrackArray1);
866 // Set selection bit for PID
867 if(okD0) SetSelectionBitForPID(fCutsD0toKpi,rd,AliRDHFCuts::kD0toKpiPID);
870 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
871 // write references in io2Prong
873 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
875 vertexp1n1->AddDaughter(postrack1);
876 vertexp1n1->AddDaughter(negtrack1);
878 io2Prong->SetSecondaryVtx(vertexp1n1);
879 //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
880 // create a track from the D0
881 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
883 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
884 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
886 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
888 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
891 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
892 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
893 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
896 //if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
898 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
899 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
901 // get track from tracks array
902 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
903 // trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
904 SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
905 twoTrackArrayCasc->AddAt(trackPi,0);
906 twoTrackArrayCasc->AddAt(trackD0,1);
907 if(!SelectInvMassAndPtDstarD0pi(twoTrackArrayCasc)){
908 twoTrackArrayCasc->Clear();
913 AliAODVertex *vertexCasc = 0;
915 if(fFindVertexForDstar) {
916 // DCA between the two tracks
917 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
919 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
921 // assume Dstar decays at the primary vertex
922 Double_t pos[3],cov[6],chi2perNDF;
924 fV1->GetCovMatrix(cov);
925 chi2perNDF = fV1->GetChi2toNDF();
926 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
930 twoTrackArrayCasc->Clear();
935 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
937 // add the D0 to the AOD (if not already done)
939 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
940 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
941 rd->SetSecondaryVtx(v2Prong);
942 v2Prong->SetParent(rd);
943 AddRefs(v2Prong,rd,event,twoTrackArray1);
944 okD0=kTRUE; // this is done to add it only once
946 // add the vertex and the cascade to the AOD
947 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
948 rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
949 rc->SetSecondaryVtx(vCasc);
950 vCasc->SetParent(rc);
951 if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
952 AddRefs(vCasc,rc,event,twoTrackArrayCasc);
953 vCasc->AddDaughter(rd); // add the D0 (in ref #1)
954 // Set selection bit for PID
955 SetSelectionBitForPID(fCutsDStartoKpipi,rc,AliRDHFCuts::kDstarPID);
957 twoTrackArrayCasc->Clear();
959 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
960 delete vertexCasc; vertexCasc=NULL;
961 } // end loop on soft pi tracks
963 if(trackD0) {delete trackD0; trackD0=NULL;}
966 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
969 twoTrackArray1->Clear();
970 if( (!f3Prong && !f4Prong) ||
971 (isLikeSign2Prong && !f3Prong) ) {
978 // 2nd LOOP ON POSITIVE TRACKS
979 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
981 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
983 //if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
985 // get track from tracks array
986 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
988 if(postrack2->Charge()<0) continue;
990 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
992 // Check single tracks cuts specific for 3 prongs
993 if(!TESTBIT(seleFlags[iTrkP2],kBit3Prong)) continue;
994 if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
995 if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
998 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
999 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
1000 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1003 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1004 if(!fLikeSign3prong) continue;
1005 if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
1006 isLikeSign3Prong=kTRUE;
1010 } else { // normal triplet (+-+)
1011 isLikeSign3Prong=kFALSE;
1013 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1014 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
1015 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1019 if(fUseKaonPIDfor3Prong){
1020 if(!TESTBIT(seleFlags[iTrkN1],kBitKaonCompat)) continue;
1022 Bool_t okForLcTopKpi=kTRUE;
1023 Int_t pidLcStatus=3; // 3= OK as pKpi and Kpipi
1025 if(!TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) &&
1026 !TESTBIT(seleFlags[iTrkP2],kBitProtonCompat) ){
1027 okForLcTopKpi=kFALSE;
1030 if(okForLcTopKpi && fUsePIDforLc>1){
1031 okForLcTopKpi=kFALSE;
1033 if(TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) &&
1034 TESTBIT(seleFlags[iTrkP2],kBitPionCompat) ){
1035 okForLcTopKpi=kTRUE;
1036 pidLcStatus+=1; // 1= OK as pKpi
1038 if(TESTBIT(seleFlags[iTrkP2],kBitProtonCompat) &&
1039 TESTBIT(seleFlags[iTrkP1],kBitPionCompat) ){
1040 okForLcTopKpi=kTRUE;
1041 pidLcStatus+=2; // 2= OK as piKp
1045 Bool_t okForDsToKKpi=kTRUE;
1046 if(fUseKaonPIDforDs){
1047 if(!TESTBIT(seleFlags[iTrkP1],kBitKaonCompat) &&
1048 !TESTBIT(seleFlags[iTrkP2],kBitKaonCompat) ) okForDsToKKpi=kFALSE;
1050 // back to primary vertex
1051 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1052 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1053 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1054 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1055 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1056 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1058 //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
1060 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1061 if(dcap2n1>dcaMax) { postrack2=0; continue; }
1062 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1063 if(dcap1p2>dcaMax) { postrack2=0; continue; }
1065 // check invariant mass cuts for D+,Ds,Lc
1068 if(postrack2->Charge()>0) {
1069 threeTrackArray->AddAt(postrack1,0);
1070 threeTrackArray->AddAt(negtrack1,1);
1071 threeTrackArray->AddAt(postrack2,2);
1073 threeTrackArray->AddAt(negtrack1,0);
1074 threeTrackArray->AddAt(postrack1,1);
1075 threeTrackArray->AddAt(postrack2,2);
1077 if(fMassCutBeforeVertexing){
1078 postrack2->GetPxPyPz(mompos2);
1079 Double_t pxDau[3]={mompos1[0],momneg1[0],mompos2[0]};
1080 Double_t pyDau[3]={mompos1[1],momneg1[1],mompos2[1]};
1081 Double_t pzDau[3]={mompos1[2],momneg1[2],mompos2[2]};
1082 // massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1083 massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau,pidLcStatus);
1087 if(f3Prong && !massCutOK) {
1088 threeTrackArray->Clear();
1096 twoTrackArray2->AddAt(postrack2,0);
1097 twoTrackArray2->AddAt(negtrack1,1);
1098 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1100 twoTrackArray2->Clear();
1105 // 3 prong candidates
1106 if(f3Prong && massCutOK) {
1108 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1109 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,okForLcTopKpi,okForDsToKKpi,ok3Prong);
1111 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1112 if(!isLikeSign3Prong) {
1113 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1114 rd->SetSecondaryVtx(v3Prong);
1115 v3Prong->SetParent(rd);
1116 AddRefs(v3Prong,rd,event,threeTrackArray);
1117 } else { // isLikeSign3Prong
1118 if(fLikeSign3prong){
1119 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1120 rd->SetSecondaryVtx(v3Prong);
1121 v3Prong->SetParent(rd);
1122 AddRefs(v3Prong,rd,event,threeTrackArray);
1125 // Set selection bit for PID
1126 SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
1127 SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
1128 SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
1130 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1131 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1134 // 4 prong candidates
1136 // don't make 4 prong with like-sign pairs and triplets
1137 && !isLikeSign2Prong && !isLikeSign3Prong
1138 // track-to-track dca cuts already now
1139 && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
1140 && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
1142 // back to primary vertex
1143 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1144 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1145 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1146 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1147 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1148 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1150 // Vertexing for these 3 (can be taken from above?)
1151 threeTrackArray->AddAt(postrack1,0);
1152 threeTrackArray->AddAt(negtrack1,1);
1153 threeTrackArray->AddAt(postrack2,2);
1154 AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1156 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
1157 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1159 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1161 //if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1163 // get track from tracks array
1164 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1166 if(negtrack2->Charge()>0) continue;
1168 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1170 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1171 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1172 evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
1173 evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
1174 evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1175 evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
1178 // back to primary vertex
1179 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1180 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1181 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1182 // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1183 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1184 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1185 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1186 SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1188 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1189 if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1190 dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1191 if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1194 fourTrackArray->AddAt(postrack1,0);
1195 fourTrackArray->AddAt(negtrack1,1);
1196 fourTrackArray->AddAt(postrack2,2);
1197 fourTrackArray->AddAt(negtrack2,3);
1199 // check invariant mass cuts for D0
1201 if(fMassCutBeforeVertexing)
1202 massCutOK = SelectInvMassAndPt4prong(fourTrackArray);
1205 fourTrackArray->Clear();
1211 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
1212 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
1214 AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
1215 rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
1216 rd->SetSecondaryVtx(v4Prong);
1217 v4Prong->SetParent(rd);
1218 AddRefs(v4Prong,rd,event,fourTrackArray);
1221 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
1222 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
1223 fourTrackArray->Clear();
1226 } // end loop on negative tracks
1228 threeTrackArray->Clear();
1229 delete vertexp1n1p2;
1236 } // end 2nd loop on positive tracks
1238 twoTrackArray2->Clear();
1240 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
1241 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1243 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1245 //if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1247 // get track from tracks array
1248 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1250 if(negtrack2->Charge()>0) continue;
1252 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1254 // Check single tracks cuts specific for 3 prongs
1255 if(!TESTBIT(seleFlags[iTrkN2],kBit3Prong)) continue;
1256 if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
1257 if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
1260 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1261 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1262 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1265 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1266 if(!fLikeSign3prong) continue;
1267 if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1268 isLikeSign3Prong=kTRUE;
1272 } else { // normal triplet (-+-)
1273 isLikeSign3Prong=kFALSE;
1276 if(fUseKaonPIDfor3Prong){
1277 if(!TESTBIT(seleFlags[iTrkP1],kBitKaonCompat)) continue;
1279 Bool_t okForLcTopKpi=kTRUE;
1280 Int_t pidLcStatus=3; // 3= OK as pKpi and Kpipi
1282 if(!TESTBIT(seleFlags[iTrkN1],kBitProtonCompat) &&
1283 !TESTBIT(seleFlags[iTrkN2],kBitProtonCompat) ){
1284 okForLcTopKpi=kFALSE;
1287 if(okForLcTopKpi && fUsePIDforLc>1){
1288 okForLcTopKpi=kFALSE;
1290 if(TESTBIT(seleFlags[iTrkN1],kBitProtonCompat) &&
1291 TESTBIT(seleFlags[iTrkN2],kBitPionCompat) ){
1292 okForLcTopKpi=kTRUE;
1293 pidLcStatus+=1; // 1= OK as pKpi
1295 if(TESTBIT(seleFlags[iTrkN2],kBitProtonCompat) &&
1296 TESTBIT(seleFlags[iTrkN1],kBitPionCompat) ){
1297 okForLcTopKpi=kTRUE;
1298 pidLcStatus+=2; // 2= OK as piKp
1302 Bool_t okForDsToKKpi=kTRUE;
1303 if(fUseKaonPIDforDs){
1304 if(!TESTBIT(seleFlags[iTrkN1],kBitKaonCompat) &&
1305 !TESTBIT(seleFlags[iTrkN2],kBitKaonCompat) ) okForDsToKKpi=kFALSE;
1308 // back to primary vertex
1309 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1310 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1311 // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1312 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1313 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1314 SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1315 //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1317 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1318 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1319 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1320 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1322 threeTrackArray->AddAt(negtrack1,0);
1323 threeTrackArray->AddAt(postrack1,1);
1324 threeTrackArray->AddAt(negtrack2,2);
1326 // check invariant mass cuts for D+,Ds,Lc
1328 if(fMassCutBeforeVertexing && f3Prong){
1329 negtrack2->GetPxPyPz(momneg2);
1330 Double_t pxDau[3]={momneg1[0],mompos1[0],momneg2[0]};
1331 Double_t pyDau[3]={momneg1[1],mompos1[1],momneg2[1]};
1332 Double_t pzDau[3]={momneg1[2],mompos1[2],momneg2[2]};
1333 // massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1334 massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau,pidLcStatus);
1337 threeTrackArray->Clear();
1343 twoTrackArray2->AddAt(postrack1,0);
1344 twoTrackArray2->AddAt(negtrack2,1);
1346 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1348 twoTrackArray2->Clear();
1354 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1355 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,okForLcTopKpi,okForDsToKKpi,ok3Prong);
1357 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1358 if(!isLikeSign3Prong) {
1359 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1360 rd->SetSecondaryVtx(v3Prong);
1361 v3Prong->SetParent(rd);
1362 AddRefs(v3Prong,rd,event,threeTrackArray);
1363 } else { // isLikeSign3Prong
1364 if(fLikeSign3prong){
1365 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1366 rd->SetSecondaryVtx(v3Prong);
1367 v3Prong->SetParent(rd);
1368 AddRefs(v3Prong,rd,event,threeTrackArray);
1371 // Set selection bit for PID
1372 SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
1373 SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
1374 SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
1376 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1377 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1379 threeTrackArray->Clear();
1383 } // end 2nd loop on negative tracks
1385 twoTrackArray2->Clear();
1389 } // end 1st loop on negative tracks
1392 } // end 1st loop on positive tracks
1395 AliDebug(1,Form(" Total HF vertices in event = %d;",
1396 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1398 AliDebug(1,Form(" D0->Kpi in event = %d;",
1399 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1402 AliDebug(1,Form(" JPSI->ee in event = %d;",
1403 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1406 AliDebug(1,Form(" Charm->3Prong in event = %d;",
1407 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1410 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1411 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1414 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1415 (Int_t)aodDstarTClArr->GetEntriesFast()));
1418 AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1419 (Int_t)aodCascadesTClArr->GetEntriesFast()));
1422 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1423 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1425 if(fLikeSign3prong && f3Prong) {
1426 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1427 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1431 twoTrackArray1->Delete(); delete twoTrackArray1;
1432 twoTrackArray2->Delete(); delete twoTrackArray2;
1433 twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1434 twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1435 threeTrackArray->Clear();
1436 threeTrackArray->Delete(); delete threeTrackArray;
1437 fourTrackArray->Delete(); delete fourTrackArray;
1438 delete [] seleFlags; seleFlags=NULL;
1439 if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1440 tracksAtVertex.Delete();
1443 seleTrksArray.Delete();
1444 if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
1448 //printf("Trks: total %d sele %d\n",fnTrksTotal,fnSeleTrksTotal);
1452 //----------------------------------------------------------------------------
1453 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1454 const AliVEvent *event,
1455 const TObjArray *trkArray) const
1457 // Add the AOD tracks as daughters of the vertex (TRef)
1458 // Also add the references to the primary vertex and to the cuts
1459 //AliCodeTimerAuto("",0);
1462 AddDaughterRefs(v,event,trkArray);
1463 rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1467 rd->SetListOfCutsRef((TList*)fListOfCuts);
1468 //fListOfCuts->Print();
1469 cout<<fListOfCuts<<endl;
1470 TList *l=(TList*)rd->GetListOfCuts();
1472 if(l) {l->Print(); }else{printf("error\n");}
1477 //----------------------------------------------------------------------------
1478 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1479 const AliVEvent *event,
1480 const TObjArray *trkArray) const
1482 // Add the AOD tracks as daughters of the vertex (TRef)
1483 //AliCodeTimerAuto("",0);
1485 Int_t nDg = v->GetNDaughters();
1487 if(nDg) dg = v->GetDaughter(0);
1489 if(dg) return; // daughters already added
1491 Int_t nTrks = trkArray->GetEntriesFast();
1493 AliExternalTrackParam *track = 0;
1494 AliAODTrack *aodTrack = 0;
1497 for(Int_t i=0; i<nTrks; i++) {
1498 track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1499 id = (Int_t)track->GetID();
1500 //printf("---> %d\n",id);
1501 if(id<0) continue; // this track is a AliAODRecoDecay
1502 aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1503 v->AddDaughter(aodTrack);
1508 //---------------------------------------------------------------------------
1509 void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)
1511 // Checks that the references to the daughter tracks are properly
1512 // assigned and reassigns them if needed
1514 //AliCodeTimerAuto("",0);
1517 TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1518 if(!inputArray) return;
1520 AliAODTrack *track = 0;
1521 AliAODVertex *vertex = 0;
1523 Bool_t needtofix=kFALSE;
1524 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1525 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1526 for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1527 track = (AliAODTrack*)vertex->GetDaughter(id);
1528 if(!track->GetStatus()) needtofix=kTRUE;
1530 if(needtofix) break;
1533 if(!needtofix) return;
1536 printf("Fixing references\n");
1538 fAODMapSize = 100000;
1539 fAODMap = new Int_t[fAODMapSize];
1540 memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
1542 for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1543 track = aod->GetTrack(i);
1545 // skip pure ITS SA tracks
1546 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1548 // skip tracks without ITS
1549 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1551 // TEMPORARY: check that the cov matrix is there
1552 Double_t covtest[21];
1553 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1556 Int_t ind = (Int_t)track->GetID();
1557 if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
1561 Int_t ids[4]={-1,-1,-1,-1};
1562 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1563 Bool_t cascade=kFALSE;
1564 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1566 Int_t nDgs = vertex->GetNDaughters();
1567 for(id=0; id<nDgs; id++) {
1568 track = (AliAODTrack*)vertex->GetDaughter(id);
1569 if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1570 ids[id]=(Int_t)track->GetID();
1571 vertex->RemoveDaughter(track);
1573 if(cascade) continue;
1574 for(id=0; id<nDgs; id++) {
1575 if (ids[id]>-1 && ids[id] < fAODMapSize) {
1576 track = aod->GetTrack(fAODMap[ids[id]]);
1577 vertex->AddDaughter(track);
1585 //----------------------------------------------------------------------------
1586 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1587 TObjArray *twoTrackArray,AliVEvent *event,
1588 AliAODVertex *secVert,
1589 AliAODRecoDecayHF2Prong *rd2Prong,
1593 // Make the cascade as a 2Prong decay and check if it passes Dstar
1594 // reconstruction cuts
1595 //AliCodeTimerAuto("",0);
1599 Bool_t dummy1,dummy2,dummy3;
1601 // We use Make2Prong to construct the AliAODRecoCascadeHF
1602 // (which inherits from AliAODRecoDecayHF2Prong)
1603 AliAODRecoCascadeHF *theCascade =
1604 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1605 dummy1,dummy2,dummy3);
1606 if(!theCascade) return 0x0;
1609 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1610 theCascade->SetCharge(trackPi->Charge());
1612 //--- selection cuts
1614 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1616 Int_t idSoftPi=(Int_t)trackPi->GetID();
1617 if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
1618 AliAODTrack* trackPiAOD=(AliAODTrack*)event->GetTrack(fAODMap[idSoftPi]);
1619 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
1622 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1624 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1626 AliAODVertex *primVertexAOD=0;
1627 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1628 // take event primary vertex
1629 primVertexAOD = PrimaryVertex();
1630 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1631 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1635 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1636 if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1638 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1639 tmpCascade->UnsetOwnPrimaryVtx();
1640 delete tmpCascade; tmpCascade=NULL;
1641 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1642 rd2Prong->UnsetOwnPrimaryVtx();
1644 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1652 //----------------------------------------------------------------------------
1653 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1654 TObjArray *twoTrackArray,AliVEvent *event,
1655 AliAODVertex *secVert,
1661 // Make the cascade as a 2Prong decay and check if it passes
1662 // cascades reconstruction cuts
1663 //AliCodeTimerAuto("",0);
1665 // AliDebug(2,Form(" building the cascade"));
1667 Bool_t dummy1,dummy2,dummy3;
1669 // We use Make2Prong to construct the AliAODRecoCascadeHF
1670 // (which inherits from AliAODRecoDecayHF2Prong)
1671 AliAODRecoCascadeHF *theCascade =
1672 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1673 dummy1,dummy2,dummy3);
1674 if(!theCascade) return 0x0;
1676 // bachelor track and charge
1677 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1678 theCascade->SetCharge(trackBachelor->Charge());
1680 //--- selection cuts
1683 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1685 Int_t idBachelor=(Int_t)trackBachelor->GetID();
1686 if (idBachelor > -1 && idBachelor < fAODMapSize) {
1687 AliAODTrack* trackBachelorAOD=(AliAODTrack*)event->GetTrack(fAODMap[idBachelor]);
1688 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
1691 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1693 tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1695 AliAODVertex *primVertexAOD=0;
1696 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1697 // take event primary vertex
1698 primVertexAOD = PrimaryVertex();
1699 if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1700 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1704 if(fCascades && fInputAOD){
1705 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1708 //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
1710 }// no cuts implemented from ESDs
1711 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1712 tmpCascade->UnsetOwnPrimaryVtx();
1713 delete tmpCascade; tmpCascade=NULL;
1714 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1720 //-----------------------------------------------------------------------------
1721 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1722 TObjArray *twoTrackArray,AliVEvent *event,
1723 AliAODVertex *secVert,Double_t dca,
1724 Bool_t &okD0,Bool_t &okJPSI,
1725 Bool_t &okD0fromDstar)
1727 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1728 // reconstruction cuts
1729 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1730 //AliCodeTimerAuto("",0);
1732 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1734 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1735 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1736 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1738 // propagate tracks to secondary vertex, to compute inv. mass
1739 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1740 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1742 Double_t momentum[3];
1743 postrack->GetPxPyPz(momentum);
1744 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1745 negtrack->GetPxPyPz(momentum);
1746 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1749 // invariant mass cut (try to improve coding here..)
1750 Bool_t okMassCut=kFALSE;
1751 if(!okMassCut && fD0toKpi) if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
1752 if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
1753 if(!okMassCut && fDstar) if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
1754 if(!okMassCut && fCascades) if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
1756 //AliDebug(2," candidate didn't pass mass cut");
1759 // primary vertex to be used by this candidate
1760 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1761 if(!primVertexAOD) return 0x0;
1764 Double_t d0z0[2],covd0z0[3];
1765 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1767 d0err[0] = TMath::Sqrt(covd0z0[0]);
1768 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1770 d0err[1] = TMath::Sqrt(covd0z0[0]);
1772 // create the object AliAODRecoDecayHF2Prong
1773 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1774 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1775 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1776 the2Prong->SetProngIDs(2,id);
1777 delete primVertexAOD; primVertexAOD=NULL;
1780 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1782 // Add daughter references already here
1783 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
1787 okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
1788 if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1790 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1791 // select J/psi from B
1793 okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1795 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1796 // select D0->Kpi from Dstar
1798 okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1799 if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1801 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1805 // remove the primary vertex (was used only for selection)
1806 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1807 the2Prong->UnsetOwnPrimaryVtx();
1810 // get PID info from ESD
1811 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1812 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1813 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1814 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1815 Double_t esdpid[10];
1816 for(Int_t i=0;i<5;i++) {
1817 esdpid[i] = esdpid0[i];
1818 esdpid[5+i] = esdpid1[i];
1820 the2Prong->SetPID(2,esdpid);
1824 //----------------------------------------------------------------------------
1825 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1826 TObjArray *threeTrackArray,AliVEvent *event,
1827 AliAODVertex *secVert,Double_t dispersion,
1828 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1829 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1830 Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
1832 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1833 // reconstruction cuts
1836 //AliCodeTimerAuto("",0);
1839 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1841 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1842 Double_t momentum[3];
1845 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1846 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1847 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1849 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1850 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1851 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1852 postrack1->GetPxPyPz(momentum);
1853 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1854 negtrack->GetPxPyPz(momentum);
1855 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1856 postrack2->GetPxPyPz(momentum);
1857 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1859 // invariant mass cut for D+, Ds, Lc
1860 Bool_t okMassCut=kFALSE;
1861 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1862 if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
1864 //AliDebug(2," candidate didn't pass mass cut");
1868 // primary vertex to be used by this candidate
1869 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1870 if(!primVertexAOD) return 0x0;
1872 Double_t d0z0[2],covd0z0[3];
1873 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1875 d0err[0] = TMath::Sqrt(covd0z0[0]);
1876 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1878 d0err[1] = TMath::Sqrt(covd0z0[0]);
1879 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1881 d0err[2] = TMath::Sqrt(covd0z0[0]);
1884 // create the object AliAODRecoDecayHF3Prong
1885 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1886 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1887 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]));
1888 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]));
1889 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1890 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1891 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1892 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1893 the3Prong->SetProngIDs(3,id);
1895 delete primVertexAOD; primVertexAOD=NULL;
1897 // Add daughter references already here
1898 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
1900 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1904 if(fOKInvMassDplus && fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1906 the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1908 if(useForDs && fOKInvMassDs){
1909 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1911 the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1914 if(useForLc && fOKInvMassLc){
1915 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1917 the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1921 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1923 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1924 the3Prong->UnsetOwnPrimaryVtx();
1927 // get PID info from ESD
1928 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1929 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1930 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1931 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1932 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1933 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1935 Double_t esdpid[15];
1936 for(Int_t i=0;i<5;i++) {
1937 esdpid[i] = esdpid0[i];
1938 esdpid[5+i] = esdpid1[i];
1939 esdpid[10+i] = esdpid2[i];
1941 the3Prong->SetPID(3,esdpid);
1945 //----------------------------------------------------------------------------
1946 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1947 TObjArray *fourTrackArray,AliVEvent *event,
1948 AliAODVertex *secVert,
1949 const AliAODVertex *vertexp1n1,
1950 const AliAODVertex *vertexp1n1p2,
1951 Double_t dcap1n1,Double_t dcap1n2,
1952 Double_t dcap2n1,Double_t dcap2n2,
1955 // Make 4Prong candidates and check if they pass D0toKpipipi
1956 // reconstruction cuts
1957 // G.E.Bruno, R.Romita
1958 //AliCodeTimerAuto("",0);
1961 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1963 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1965 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1966 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1967 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1968 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1970 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1971 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1972 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1973 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1975 Double_t momentum[3];
1976 postrack1->GetPxPyPz(momentum);
1977 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1978 negtrack1->GetPxPyPz(momentum);
1979 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1980 postrack2->GetPxPyPz(momentum);
1981 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1982 negtrack2->GetPxPyPz(momentum);
1983 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1985 // invariant mass cut for rho or D0 (try to improve coding here..)
1986 Bool_t okMassCut=kFALSE;
1987 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1988 if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
1989 if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
1992 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1993 //printf(" candidate didn't pass mass cut\n");
1997 // primary vertex to be used by this candidate
1998 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
1999 if(!primVertexAOD) return 0x0;
2001 Double_t d0z0[2],covd0z0[3];
2002 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2004 d0err[0] = TMath::Sqrt(covd0z0[0]);
2005 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2007 d0err[1] = TMath::Sqrt(covd0z0[0]);
2008 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2010 d0err[2] = TMath::Sqrt(covd0z0[0]);
2011 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2013 d0err[3] = TMath::Sqrt(covd0z0[0]);
2016 // create the object AliAODRecoDecayHF4Prong
2017 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2018 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2019 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]));
2020 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]));
2021 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]));
2023 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2024 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2025 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2026 the4Prong->SetProngIDs(4,id);
2028 delete primVertexAOD; primVertexAOD=NULL;
2030 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
2033 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
2034 the4Prong->UnsetOwnPrimaryVtx();
2038 // get PID info from ESD
2039 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2040 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2041 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2042 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2043 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2044 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2045 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2046 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2048 Double_t esdpid[20];
2049 for(Int_t i=0;i<5;i++) {
2050 esdpid[i] = esdpid0[i];
2051 esdpid[5+i] = esdpid1[i];
2052 esdpid[10+i] = esdpid2[i];
2053 esdpid[15+i] = esdpid3[i];
2055 the4Prong->SetPID(4,esdpid);
2059 //-----------------------------------------------------------------------------
2060 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2061 AliVEvent *event) const
2063 // Returns primary vertex to be used for this candidate
2064 //AliCodeTimerAuto("",0);
2066 AliESDVertex *vertexESD = 0;
2067 AliAODVertex *vertexAOD = 0;
2070 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
2071 // primary vertex from the input event
2073 vertexESD = new AliESDVertex(*fV1);
2076 // primary vertex specific to this candidate
2078 Int_t nTrks = trkArray->GetEntriesFast();
2079 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2081 if(fRecoPrimVtxSkippingTrks) {
2082 // recalculating the vertex
2084 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2085 Float_t diamondcovxy[3];
2086 event->GetDiamondCovXY(diamondcovxy);
2087 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2088 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2089 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2090 vertexer->SetVtxStart(diamond);
2091 delete diamond; diamond=NULL;
2092 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2093 vertexer->SetOnlyFitter();
2095 Int_t skipped[1000];
2096 Int_t nTrksToSkip=0,id;
2097 AliExternalTrackParam *t = 0;
2098 for(Int_t i=0; i<nTrks; i++) {
2099 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2100 id = (Int_t)t->GetID();
2102 skipped[nTrksToSkip++] = id;
2105 // For AOD, skip also tracks without covariance matrix
2107 Double_t covtest[21];
2108 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2109 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2110 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2111 id = (Int_t)vtrack->GetID();
2113 skipped[nTrksToSkip++] = id;
2117 for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2119 vertexer->SetSkipTracks(nTrksToSkip,skipped);
2120 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2122 } else if(fRmTrksFromPrimVtx && nTrks>0) {
2123 // removing the prongs tracks
2125 TObjArray rmArray(nTrks);
2126 UShort_t *rmId = new UShort_t[nTrks];
2127 AliESDtrack *esdTrack = 0;
2129 for(Int_t i=0; i<nTrks; i++) {
2130 t = (AliESDtrack*)trkArray->UncheckedAt(i);
2131 esdTrack = new AliESDtrack(*t);
2132 rmArray.AddLast(esdTrack);
2133 if(esdTrack->GetID()>=0) {
2134 rmId[i]=(UShort_t)esdTrack->GetID();
2139 Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2140 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2141 delete [] rmId; rmId=NULL;
2146 if(!vertexESD) return vertexAOD;
2147 if(vertexESD->GetNContributors()<=0) {
2148 //AliDebug(2,"vertexing failed");
2149 delete vertexESD; vertexESD=NULL;
2153 delete vertexer; vertexer=NULL;
2157 // convert to AliAODVertex
2158 Double_t pos[3],cov[6],chi2perNDF;
2159 vertexESD->GetXYZ(pos); // position
2160 vertexESD->GetCovMatrix(cov); //covariance matrix
2161 chi2perNDF = vertexESD->GetChi2toNDF();
2162 delete vertexESD; vertexESD=NULL;
2164 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2168 //-----------------------------------------------------------------------------
2169 void AliAnalysisVertexingHF::PrintStatus() const {
2170 // Print parameters being used
2172 //printf("Preselections:\n");
2173 // fTrackFilter->Dump();
2175 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2177 printf("Secondary vertex with AliVertexerTracks\n");
2179 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2180 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2182 printf("Reconstruct D0->Kpi candidates with cuts:\n");
2183 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
2186 printf("Reconstruct D*->D0pi candidates with cuts:\n");
2187 if(fFindVertexForDstar) {
2188 printf(" Reconstruct a secondary vertex for the D*\n");
2190 printf(" Assume the D* comes from the primary vertex\n");
2192 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
2195 printf("Reconstruct J/psi from B candidates with cuts:\n");
2196 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
2199 printf("Reconstruct 3 prong candidates.\n");
2200 printf(" D+->Kpipi cuts:\n");
2201 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
2202 printf(" Ds->KKpi cuts:\n");
2203 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
2204 printf(" Lc->pKpi cuts:\n");
2205 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
2208 printf("Reconstruct 4 prong candidates.\n");
2209 printf(" D0->Kpipipi cuts:\n");
2210 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
2213 printf("Reconstruct cascades candidates formed with v0s.\n");
2214 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
2215 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
2220 //-----------------------------------------------------------------------------
2221 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
2222 Double_t &dispersion,Bool_t useTRefArray) const
2224 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
2225 //AliCodeTimerAuto("",0);
2227 AliESDVertex *vertexESD = 0;
2228 AliAODVertex *vertexAOD = 0;
2230 if(!fSecVtxWithKF) { // AliVertexerTracks
2232 fVertexerTracks->SetVtxStart(fV1);
2233 vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2235 if(!vertexESD) return vertexAOD;
2237 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2238 //AliDebug(2,"vertexing failed");
2239 delete vertexESD; vertexESD=NULL;
2243 Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
2245 // vertex outside beam pipe, reject candidate to avoid propagation through material
2246 delete vertexESD; vertexESD=NULL;
2250 } else { // Kalman Filter vertexer (AliKFParticle)
2252 AliKFParticle::SetField(fBzkG);
2254 AliKFVertex vertexKF;
2256 Int_t nTrks = trkArray->GetEntriesFast();
2257 for(Int_t i=0; i<nTrks; i++) {
2258 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2259 AliKFParticle daughterKF(*esdTrack,211);
2260 vertexKF.AddDaughter(daughterKF);
2262 vertexESD = new AliESDVertex(vertexKF.Parameters(),
2263 vertexKF.CovarianceMatrix(),
2265 vertexKF.GetNContributors());
2269 // convert to AliAODVertex
2270 Double_t pos[3],cov[6],chi2perNDF;
2271 vertexESD->GetXYZ(pos); // position
2272 vertexESD->GetCovMatrix(cov); //covariance matrix
2273 chi2perNDF = vertexESD->GetChi2toNDF();
2274 dispersion = vertexESD->GetDispersion();
2275 delete vertexESD; vertexESD=NULL;
2277 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2278 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2282 //-----------------------------------------------------------------------------
2283 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(TObjArray *trkArray){
2284 // Invariant mass cut on tracks
2285 //AliCodeTimerAuto("",0);
2287 Int_t retval=kFALSE;
2288 Double_t momentum[3];
2289 Double_t px[3],py[3],pz[3];
2290 for(Int_t iTrack=0; iTrack<3; iTrack++){
2291 AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2292 track->GetPxPyPz(momentum);
2293 px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2295 retval = SelectInvMassAndPt3prong(px,py,pz);
2300 //-----------------------------------------------------------------------------
2301 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(TObjArray *trkArray){
2302 // Invariant mass cut on tracks
2303 //AliCodeTimerAuto("",0);
2305 Int_t retval=kFALSE;
2306 Double_t momentum[3];
2307 Double_t px[4],py[4],pz[4];
2309 for(Int_t iTrack=0; iTrack<4; iTrack++){
2310 AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2311 track->GetPxPyPz(momentum);
2312 px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2315 retval = SelectInvMassAndPt4prong(px,py,pz);
2319 //-----------------------------------------------------------------------------
2320 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(TObjArray *trkArray){
2321 // Invariant mass cut on tracks
2322 //AliCodeTimerAuto("",0);
2324 Int_t retval=kFALSE;
2325 Double_t momentum[3];
2326 Double_t px[2],py[2],pz[2];
2328 for(Int_t iTrack=0; iTrack<2; iTrack++){
2329 AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2330 track->GetPxPyPz(momentum);
2331 px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2333 retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
2337 //-----------------------------------------------------------------------------
2338 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtD0Kpi(Double_t *px,
2341 // Check invariant mass cut and pt candidate cut
2342 //AliCodeTimerAuto("",0);
2346 Double_t minv2,mrange;
2347 Double_t lolim,hilim;
2349 Bool_t retval=kFALSE;
2351 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2352 fOKInvMassD0=kFALSE;
2354 minPt=fCutsD0toKpi->GetMinPtCandidate();
2356 if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2358 mrange=fCutsD0toKpi->GetMassCut();
2359 lolim=fMassDzero-mrange;
2360 hilim=fMassDzero+mrange;
2361 pdg2[0]=211; pdg2[1]=321;
2362 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2363 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2367 pdg2[0]=321; pdg2[1]=211;
2368 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2369 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2376 //-----------------------------------------------------------------------------
2377 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtJpsiee(Double_t *px,
2380 // Check invariant mass cut and pt candidate cut
2381 //AliCodeTimerAuto("",0);
2385 Double_t minv2,mrange;
2386 Double_t lolim,hilim;
2388 Bool_t retval=kFALSE;
2390 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2391 fOKInvMassJpsi=kFALSE;
2393 minPt=fCutsJpsitoee->GetMinPtCandidate();
2395 if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2397 mrange=fCutsJpsitoee->GetMassCut();
2398 lolim=fMassJpsi-mrange;
2399 hilim=fMassJpsi+mrange;
2401 pdg2[0]=11; pdg2[1]=11;
2402 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2403 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2405 fOKInvMassJpsi=kTRUE;
2410 //-----------------------------------------------------------------------------
2411 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(Double_t *px,
2415 // Check invariant mass cut and pt candidate cut
2416 //AliCodeTimerAuto("",0);
2420 Double_t minv2,mrange;
2421 Double_t lolim,hilim;
2423 Bool_t retval=kFALSE;
2426 fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2427 fOKInvMassDplus=kFALSE;
2428 fOKInvMassDs=kFALSE;
2429 fOKInvMassLc=kFALSE;
2431 minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
2432 minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
2434 if(fMassCalc3->Pt2() < minPt*minPt) return retval;
2436 mrange=fCutsDplustoKpipi->GetMassCut();
2437 lolim=fMassDplus-mrange;
2438 hilim=fMassDplus+mrange;
2439 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2440 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2441 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2443 fOKInvMassDplus=kTRUE;
2446 mrange=fCutsDstoKKpi->GetMassCut();
2447 lolim=fMassDs-mrange;
2448 hilim=fMassDs+mrange;
2449 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2450 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2451 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2455 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2456 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2457 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2462 mrange=fCutsLctopKpi->GetMassCut();
2463 lolim=fMassLambdaC-mrange;
2464 hilim=fMassLambdaC+mrange;
2466 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2467 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2468 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2474 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2475 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2476 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2485 //-----------------------------------------------------------------------------
2486 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(Double_t *px,
2489 // Check invariant mass cut and pt candidate cut
2490 //AliCodeTimerAuto("",0);
2494 Double_t minv2,mrange;
2495 Double_t lolim,hilim;
2497 Bool_t retval=kFALSE;
2499 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2500 fOKInvMassDstar=kFALSE;
2502 minPt=fCutsDStartoKpipi->GetMinPtCandidate();
2504 if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2506 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2507 mrange=fCutsDStartoKpipi->GetMassCut();
2508 lolim=fMassDstar-mrange;
2509 hilim=fMassDstar+mrange;
2510 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2511 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2513 fOKInvMassDstar=kTRUE;
2519 //-----------------------------------------------------------------------------
2520 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(Double_t *px,
2523 // Check invariant mass cut and pt candidate cut
2524 //AliCodeTimerAuto("",0);
2528 Double_t minv2,mrange;
2529 Double_t lolim,hilim;
2531 Bool_t retval=kFALSE;
2533 // D0->Kpipipi without PID
2534 fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2535 fOKInvMassD0to4p=kFALSE;
2537 minPt=fCutsD0toKpipipi->GetMinPtCandidate();
2539 if(fMassCalc4->Pt2() < minPt*minPt) return retval;
2541 mrange=fCutsD0toKpipipi->GetMassCut();
2542 lolim=fMassDzero-mrange;
2543 hilim=fMassDzero+mrange;
2545 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2546 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2547 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2549 fOKInvMassD0to4p=kTRUE;
2552 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2553 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2554 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2556 fOKInvMassD0to4p=kTRUE;
2559 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2560 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2561 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2563 fOKInvMassD0to4p=kTRUE;
2566 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2567 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2568 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2570 fOKInvMassD0to4p=kTRUE;
2575 //-----------------------------------------------------------------------------
2576 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtCascade(Double_t *px,
2579 // Check invariant mass cut and pt candidate cut
2580 //AliCodeTimerAuto("",0);
2584 Double_t minv2,mrange;
2585 Double_t lolim,hilim;
2586 // Double_t minPt=0;
2587 Bool_t retval=kFALSE;
2589 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2590 // minPt=fCutsLctoV0->GetMinPtCandidate();
2591 fOKInvMassLctoV0=kFALSE;
2592 mrange=fCutsLctoV0->GetMassCut();
2593 lolim=fMassLambdaC-mrange;
2594 hilim=fMassLambdaC+mrange;
2595 pdg2[0]=2212;pdg2[1]=310;
2596 minv2=fMassCalc2->InvMass2(2,pdg2);
2597 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2599 fOKInvMassLctoV0=kTRUE;
2601 pdg2[0]=211;pdg2[1]=3122;
2602 minv2=fMassCalc2->InvMass2(2,pdg2);
2603 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2605 fOKInvMassLctoV0=kTRUE;
2610 //-----------------------------------------------------------------------------
2611 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2613 TObjArray &seleTrksArray,
2614 TObjArray &tracksAtVertex,
2616 UChar_t *seleFlags,Int_t *evtNumber)
2618 // Apply single-track preselection.
2619 // Fill a TObjArray with selected tracks (for displaced vertices or
2620 // soft pion from D*). Selection flag stored in seleFlags.
2621 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2622 // In case of AOD input, also fill fAODMap for track index<->ID
2623 //AliCodeTimerAuto("",0);
2625 const AliVVertex *vprimary = event->GetPrimaryVertex();
2627 if(fV1) { delete fV1; fV1=NULL; }
2628 if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
2631 UShort_t *indices = 0;
2632 Double_t pos[3],cov[6];
2633 const Int_t entries = event->GetNumberOfTracks();
2634 AliCentrality* cent;
2636 if(!fInputAOD) { // ESD
2637 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2638 cent=((AliESDEvent*)event)->GetCentrality();
2640 vprimary->GetXYZ(pos);
2641 vprimary->GetCovarianceMatrix(cov);
2642 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2643 if(entries<=0) return;
2644 indices = new UShort_t[entries];
2645 memset(indices,0,sizeof(UShort_t)*entries);
2646 fAODMapSize = 100000;
2647 fAODMap = new Int_t[fAODMapSize];
2648 memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2649 cent=((AliAODEvent*)event)->GetCentrality();
2651 Float_t centperc=cent->GetCentralityPercentile("V0M");
2653 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE;
2656 // transfer ITS tracks from event to arrays
2657 for(Int_t i=0; i<entries; i++) {
2659 track = (AliVTrack*)event->GetTrack(i);
2661 // skip pure ITS SA tracks
2662 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2664 // skip tracks without ITS
2665 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2667 // skip tracks with negative ID
2668 // (these are duplicated TPC-only AOD tracks, for jet analysis...)
2669 if(track->GetID()<0) continue;
2671 // TEMPORARY: check that the cov matrix is there
2672 Double_t covtest[21];
2673 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2677 AliAODTrack *aodt = (AliAODTrack*)track;
2678 if(aodt->GetUsedForPrimVtxFit()) {
2679 indices[nindices]=aodt->GetID(); nindices++;
2681 Int_t ind = (Int_t)aodt->GetID();
2682 if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2685 AliESDtrack *esdt = 0;
2688 esdt = (AliESDtrack*)track;
2690 esdt = new AliESDtrack(track);
2693 // single track selection
2694 okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE;
2695 if(fMixEvent && i<trkEntries){
2696 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2697 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2698 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2699 eventVtx->GetXYZ(vtxPos);
2700 vprimary->GetXYZ(primPos);
2701 eventVtx->GetCovarianceMatrix(primCov);
2702 for(Int_t ind=0;ind<3;ind++){
2703 trasl[ind]=vtxPos[ind]-primPos[ind];
2706 Bool_t isTransl=esdt->Translate(trasl,primCov);
2714 if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong) && nSeleTrks<trkEntries) {
2715 esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
2716 seleTrksArray.AddLast(esdt);
2717 tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
2718 seleFlags[nSeleTrks]=0;
2719 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2720 if(okFor3Prong) SETBIT(seleFlags[nSeleTrks],kBit3Prong);
2721 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2724 SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
2725 SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2726 SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2727 Bool_t useTPC=kTRUE;
2729 Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
2730 if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
2731 CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
2733 Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
2734 if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
2735 CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2737 Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
2738 if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
2739 CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2741 if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
2743 if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
2744 Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
2745 if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
2746 CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
2748 Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
2749 if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
2750 CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2752 Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
2753 if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
2754 CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2759 if(fInputAOD) delete esdt;
2764 } // end loop on tracks
2766 // primary vertex from AOD
2768 delete fV1; fV1=NULL;
2769 vprimary->GetXYZ(pos);
2770 vprimary->GetCovarianceMatrix(cov);
2771 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2772 Int_t ncontr=nindices;
2773 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2774 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2775 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2776 fV1->SetTitle(vprimary->GetTitle());
2777 fV1->SetIndices(nindices,indices);
2779 if(indices) { delete [] indices; indices=NULL; }
2784 //-----------------------------------------------------------------------------
2785 void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
2787 // Set the selection bit for PID
2789 //AliCodeTimerAuto("",0);
2790 if(cuts->GetPidHF()) {
2791 Bool_t usepid=cuts->GetIsUsePID();
2792 cuts->SetUsePID(kTRUE);
2793 if(cuts->IsSelectedPID(rd))
2794 rd->SetSelectionBit(bit);
2795 cuts->SetUsePID(usepid);
2799 //-----------------------------------------------------------------------------
2800 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2801 Float_t centralityperc,
2802 Bool_t &okDisplaced,
2804 Bool_t &okFor3Prong) const
2806 // Check if track passes some kinematical cuts
2808 // this is needed to store the impact parameters
2809 //AliCodeTimerAuto("",0);
2811 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2815 // Track selection, displaced tracks -- 2 prong
2817 if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
2818 && fTrackFilter2prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2819 // central PbPb events, tighter cuts
2820 selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
2824 selectInfo = fTrackFilter->IsSelected(trk);
2827 if(selectInfo) okDisplaced=kTRUE;
2829 // Track selection, displaced tracks -- 3 prong
2831 if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
2832 && fTrackFilter3prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2833 // central PbPb events, tighter cuts
2834 selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
2838 selectInfo = fTrackFilter->IsSelected(trk);
2841 if(selectInfo) okFor3Prong=kTRUE;
2843 // Track selection, soft pions
2845 if(fDstar && fTrackFilterSoftPi) {
2846 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2848 if(selectInfo) okSoftPi=kTRUE;
2850 if(okDisplaced || okSoftPi || okFor3Prong) return kTRUE;
2856 //-----------------------------------------------------------------------------
2857 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2859 // Transform ESDv0 to AODv0
2861 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2862 // and creates an AODv0 out of them
2864 //AliCodeTimerAuto("",0);
2865 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2866 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2868 // create the v0 neutral track to compute the DCA to the primary vertex
2869 Double_t xyz[3], pxpypz[3];
2871 esdV0->PxPyPz(pxpypz);
2872 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2873 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2878 Double_t d0z0[2],covd0z0[3];
2879 AliAODVertex *primVertexAOD = PrimaryVertex();
2880 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2881 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2882 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2883 Double_t dcaV0DaughterToPrimVertex[2];
2884 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2885 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2886 if( !posV0track || !negV0track) {
2887 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2889 delete primVertexAOD;
2892 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2893 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2895 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2896 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2897 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2899 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2900 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2901 Double_t pmom[3],nmom[3];
2902 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2903 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2905 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2906 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2909 delete primVertexAOD;
2913 //-----------------------------------------------------------------------------
2914 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
2915 // Set the stored track parameters at primary vertex into AliESDtrack
2916 //AliCodeTimerAuto("",0);
2918 const Double_t *par=extpar->GetParameter();
2919 const Double_t *cov=extpar->GetCovariance();
2920 Double_t alpha=extpar->GetAlpha();
2921 Double_t x=extpar->GetX();
2922 esdt->Set(x,alpha,par,cov);
2925 //-----------------------------------------------------------------------------
2926 void AliAnalysisVertexingHF::SetMasses(){
2927 // Set the hadron mass values from TDatabasePDG
2929 fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2930 fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2931 fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2932 fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2933 fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2934 fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2936 //-----------------------------------------------------------------------------
2937 Bool_t AliAnalysisVertexingHF::CheckCutsConsistency(){
2939 // Check the Vertexer and the analysts task consitstecny
2943 //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
2946 if ( fCutsLctoV0 && fV0TypeForCascadeVertex != fCutsLctoV0->GetV0Type())
2948 printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
2953 //-----------------------------------------------------------------------------