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 = dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[id]));
1503 if(!aodTrack) AliFatal("Not a standard AOD");
1504 v->AddDaughter(aodTrack);
1509 //---------------------------------------------------------------------------
1510 void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)
1512 // Checks that the references to the daughter tracks are properly
1513 // assigned and reassigns them if needed
1515 //AliCodeTimerAuto("",0);
1518 TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1519 if(!inputArray) return;
1521 AliAODTrack *track = 0;
1522 AliAODVertex *vertex = 0;
1524 Bool_t needtofix=kFALSE;
1525 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1526 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1527 for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1528 track = (AliAODTrack*)vertex->GetDaughter(id);
1529 if(!track->GetStatus()) needtofix=kTRUE;
1531 if(needtofix) break;
1534 if(!needtofix) return;
1537 printf("Fixing references\n");
1539 fAODMapSize = 100000;
1540 fAODMap = new Int_t[fAODMapSize];
1541 memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
1543 for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1544 track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
1545 if(!track) AliFatal("Not a standard AOD");
1547 // skip pure ITS SA tracks
1548 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1550 // skip tracks without ITS
1551 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1553 // TEMPORARY: check that the cov matrix is there
1554 Double_t covtest[21];
1555 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1558 Int_t ind = (Int_t)track->GetID();
1559 if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
1563 Int_t ids[4]={-1,-1,-1,-1};
1564 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1565 Bool_t cascade=kFALSE;
1566 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1568 Int_t nDgs = vertex->GetNDaughters();
1569 for(id=0; id<nDgs; id++) {
1570 track = (AliAODTrack*)vertex->GetDaughter(id);
1571 if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1572 ids[id]=(Int_t)track->GetID();
1573 vertex->RemoveDaughter(track);
1575 if(cascade) continue;
1576 for(id=0; id<nDgs; id++) {
1577 if (ids[id]>-1 && ids[id] < fAODMapSize) {
1578 track = dynamic_cast<AliAODTrack*>(aod->GetTrack(fAODMap[ids[id]]));
1579 if(!track) AliFatal("Not a standard AOD");
1580 vertex->AddDaughter(track);
1588 //----------------------------------------------------------------------------
1589 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1590 TObjArray *twoTrackArray,AliVEvent *event,
1591 AliAODVertex *secVert,
1592 AliAODRecoDecayHF2Prong *rd2Prong,
1596 // Make the cascade as a 2Prong decay and check if it passes Dstar
1597 // reconstruction cuts
1598 //AliCodeTimerAuto("",0);
1602 Bool_t dummy1,dummy2,dummy3;
1604 // We use Make2Prong to construct the AliAODRecoCascadeHF
1605 // (which inherits from AliAODRecoDecayHF2Prong)
1606 AliAODRecoCascadeHF *theCascade =
1607 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1608 dummy1,dummy2,dummy3);
1609 if(!theCascade) return 0x0;
1612 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1613 theCascade->SetCharge(trackPi->Charge());
1615 //--- selection cuts
1617 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1619 Int_t idSoftPi=(Int_t)trackPi->GetID();
1620 if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
1621 AliAODTrack* trackPiAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idSoftPi]));
1622 if(!trackPiAOD) AliFatal("Not a standard AOD");
1623 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
1626 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1628 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1630 AliAODVertex *primVertexAOD=0;
1631 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1632 // take event primary vertex
1633 primVertexAOD = PrimaryVertex();
1634 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1635 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1639 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1640 if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1642 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1643 tmpCascade->UnsetOwnPrimaryVtx();
1644 delete tmpCascade; tmpCascade=NULL;
1645 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1646 rd2Prong->UnsetOwnPrimaryVtx();
1648 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1656 //----------------------------------------------------------------------------
1657 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1658 TObjArray *twoTrackArray,AliVEvent *event,
1659 AliAODVertex *secVert,
1665 // Make the cascade as a 2Prong decay and check if it passes
1666 // cascades reconstruction cuts
1667 //AliCodeTimerAuto("",0);
1669 // AliDebug(2,Form(" building the cascade"));
1671 Bool_t dummy1,dummy2,dummy3;
1673 // We use Make2Prong to construct the AliAODRecoCascadeHF
1674 // (which inherits from AliAODRecoDecayHF2Prong)
1675 AliAODRecoCascadeHF *theCascade =
1676 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1677 dummy1,dummy2,dummy3);
1678 if(!theCascade) return 0x0;
1680 // bachelor track and charge
1681 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1682 theCascade->SetCharge(trackBachelor->Charge());
1684 //--- selection cuts
1687 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1689 Int_t idBachelor=(Int_t)trackBachelor->GetID();
1690 if (idBachelor > -1 && idBachelor < fAODMapSize) {
1691 AliAODTrack* trackBachelorAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idBachelor]));
1692 if(!trackBachelorAOD) AliFatal("Not a standard AOD");
1693 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
1696 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1698 tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1700 AliAODVertex *primVertexAOD=0;
1701 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1702 // take event primary vertex
1703 primVertexAOD = PrimaryVertex();
1704 if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1705 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1709 if(fCascades && fInputAOD){
1710 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1713 //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
1715 }// no cuts implemented from ESDs
1716 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1717 tmpCascade->UnsetOwnPrimaryVtx();
1718 delete tmpCascade; tmpCascade=NULL;
1719 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1725 //-----------------------------------------------------------------------------
1726 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1727 TObjArray *twoTrackArray,AliVEvent *event,
1728 AliAODVertex *secVert,Double_t dca,
1729 Bool_t &okD0,Bool_t &okJPSI,
1730 Bool_t &okD0fromDstar)
1732 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1733 // reconstruction cuts
1734 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1735 //AliCodeTimerAuto("",0);
1737 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1739 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1740 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1741 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1743 // propagate tracks to secondary vertex, to compute inv. mass
1744 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1745 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1747 Double_t momentum[3];
1748 postrack->GetPxPyPz(momentum);
1749 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1750 negtrack->GetPxPyPz(momentum);
1751 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1754 // invariant mass cut (try to improve coding here..)
1755 Bool_t okMassCut=kFALSE;
1756 if(!okMassCut && fD0toKpi) if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
1757 if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
1758 if(!okMassCut && fDstar) if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
1759 if(!okMassCut && fCascades) if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
1761 //AliDebug(2," candidate didn't pass mass cut");
1764 // primary vertex to be used by this candidate
1765 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1766 if(!primVertexAOD) return 0x0;
1769 Double_t d0z0[2],covd0z0[3];
1770 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1772 d0err[0] = TMath::Sqrt(covd0z0[0]);
1773 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1775 d0err[1] = TMath::Sqrt(covd0z0[0]);
1777 // create the object AliAODRecoDecayHF2Prong
1778 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1779 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1780 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1781 the2Prong->SetProngIDs(2,id);
1782 delete primVertexAOD; primVertexAOD=NULL;
1785 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1787 // Add daughter references already here
1788 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
1792 okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
1793 if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1795 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1796 // select J/psi from B
1798 okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1800 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1801 // select D0->Kpi from Dstar
1803 okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1804 if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1806 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1810 // remove the primary vertex (was used only for selection)
1811 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1812 the2Prong->UnsetOwnPrimaryVtx();
1815 // get PID info from ESD
1816 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1817 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1818 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1819 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1820 Double_t esdpid[10];
1821 for(Int_t i=0;i<5;i++) {
1822 esdpid[i] = esdpid0[i];
1823 esdpid[5+i] = esdpid1[i];
1825 the2Prong->SetPID(2,esdpid);
1829 //----------------------------------------------------------------------------
1830 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1831 TObjArray *threeTrackArray,AliVEvent *event,
1832 AliAODVertex *secVert,Double_t dispersion,
1833 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1834 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1835 Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
1837 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1838 // reconstruction cuts
1841 //AliCodeTimerAuto("",0);
1844 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1846 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1847 Double_t momentum[3];
1850 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1851 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1852 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1854 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1855 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1856 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1857 postrack1->GetPxPyPz(momentum);
1858 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1859 negtrack->GetPxPyPz(momentum);
1860 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1861 postrack2->GetPxPyPz(momentum);
1862 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1864 // invariant mass cut for D+, Ds, Lc
1865 Bool_t okMassCut=kFALSE;
1866 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1867 if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
1869 //AliDebug(2," candidate didn't pass mass cut");
1873 // primary vertex to be used by this candidate
1874 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1875 if(!primVertexAOD) return 0x0;
1877 Double_t d0z0[2],covd0z0[3];
1878 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1880 d0err[0] = TMath::Sqrt(covd0z0[0]);
1881 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1883 d0err[1] = TMath::Sqrt(covd0z0[0]);
1884 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1886 d0err[2] = TMath::Sqrt(covd0z0[0]);
1889 // create the object AliAODRecoDecayHF3Prong
1890 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1891 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1892 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]));
1893 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]));
1894 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1895 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1896 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1897 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1898 the3Prong->SetProngIDs(3,id);
1900 delete primVertexAOD; primVertexAOD=NULL;
1902 // Add daughter references already here
1903 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
1905 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1909 if(fOKInvMassDplus && fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1911 the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1913 if(useForDs && fOKInvMassDs){
1914 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1916 the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1919 if(useForLc && fOKInvMassLc){
1920 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1922 the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1926 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1928 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1929 the3Prong->UnsetOwnPrimaryVtx();
1932 // get PID info from ESD
1933 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1934 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1935 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1936 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1937 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1938 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1940 Double_t esdpid[15];
1941 for(Int_t i=0;i<5;i++) {
1942 esdpid[i] = esdpid0[i];
1943 esdpid[5+i] = esdpid1[i];
1944 esdpid[10+i] = esdpid2[i];
1946 the3Prong->SetPID(3,esdpid);
1950 //----------------------------------------------------------------------------
1951 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1952 TObjArray *fourTrackArray,AliVEvent *event,
1953 AliAODVertex *secVert,
1954 const AliAODVertex *vertexp1n1,
1955 const AliAODVertex *vertexp1n1p2,
1956 Double_t dcap1n1,Double_t dcap1n2,
1957 Double_t dcap2n1,Double_t dcap2n2,
1960 // Make 4Prong candidates and check if they pass D0toKpipipi
1961 // reconstruction cuts
1962 // G.E.Bruno, R.Romita
1963 //AliCodeTimerAuto("",0);
1966 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1968 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1970 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1971 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1972 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1973 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1975 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1976 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1977 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1978 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1980 Double_t momentum[3];
1981 postrack1->GetPxPyPz(momentum);
1982 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1983 negtrack1->GetPxPyPz(momentum);
1984 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1985 postrack2->GetPxPyPz(momentum);
1986 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1987 negtrack2->GetPxPyPz(momentum);
1988 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1990 // invariant mass cut for rho or D0 (try to improve coding here..)
1991 Bool_t okMassCut=kFALSE;
1992 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1993 if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
1994 if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
1997 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1998 //printf(" candidate didn't pass mass cut\n");
2002 // primary vertex to be used by this candidate
2003 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
2004 if(!primVertexAOD) return 0x0;
2006 Double_t d0z0[2],covd0z0[3];
2007 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2009 d0err[0] = TMath::Sqrt(covd0z0[0]);
2010 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2012 d0err[1] = TMath::Sqrt(covd0z0[0]);
2013 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2015 d0err[2] = TMath::Sqrt(covd0z0[0]);
2016 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2018 d0err[3] = TMath::Sqrt(covd0z0[0]);
2021 // create the object AliAODRecoDecayHF4Prong
2022 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2023 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2024 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]));
2025 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]));
2026 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]));
2028 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2029 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2030 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2031 the4Prong->SetProngIDs(4,id);
2033 delete primVertexAOD; primVertexAOD=NULL;
2035 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
2038 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
2039 the4Prong->UnsetOwnPrimaryVtx();
2043 // get PID info from ESD
2044 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2045 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2046 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2047 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2048 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2049 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2050 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2051 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2053 Double_t esdpid[20];
2054 for(Int_t i=0;i<5;i++) {
2055 esdpid[i] = esdpid0[i];
2056 esdpid[5+i] = esdpid1[i];
2057 esdpid[10+i] = esdpid2[i];
2058 esdpid[15+i] = esdpid3[i];
2060 the4Prong->SetPID(4,esdpid);
2064 //-----------------------------------------------------------------------------
2065 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2066 AliVEvent *event) const
2068 // Returns primary vertex to be used for this candidate
2069 //AliCodeTimerAuto("",0);
2071 AliESDVertex *vertexESD = 0;
2072 AliAODVertex *vertexAOD = 0;
2075 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
2076 // primary vertex from the input event
2078 vertexESD = new AliESDVertex(*fV1);
2081 // primary vertex specific to this candidate
2083 Int_t nTrks = trkArray->GetEntriesFast();
2084 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2086 if(fRecoPrimVtxSkippingTrks) {
2087 // recalculating the vertex
2089 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2090 Float_t diamondcovxy[3];
2091 event->GetDiamondCovXY(diamondcovxy);
2092 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2093 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2094 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2095 vertexer->SetVtxStart(diamond);
2096 delete diamond; diamond=NULL;
2097 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2098 vertexer->SetOnlyFitter();
2100 Int_t skipped[1000];
2101 Int_t nTrksToSkip=0,id;
2102 AliExternalTrackParam *t = 0;
2103 for(Int_t i=0; i<nTrks; i++) {
2104 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2105 id = (Int_t)t->GetID();
2107 skipped[nTrksToSkip++] = id;
2110 // For AOD, skip also tracks without covariance matrix
2112 Double_t covtest[21];
2113 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2114 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2115 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2116 id = (Int_t)vtrack->GetID();
2118 skipped[nTrksToSkip++] = id;
2122 for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2124 vertexer->SetSkipTracks(nTrksToSkip,skipped);
2125 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2127 } else if(fRmTrksFromPrimVtx && nTrks>0) {
2128 // removing the prongs tracks
2130 TObjArray rmArray(nTrks);
2131 UShort_t *rmId = new UShort_t[nTrks];
2132 AliESDtrack *esdTrack = 0;
2134 for(Int_t i=0; i<nTrks; i++) {
2135 t = (AliESDtrack*)trkArray->UncheckedAt(i);
2136 esdTrack = new AliESDtrack(*t);
2137 rmArray.AddLast(esdTrack);
2138 if(esdTrack->GetID()>=0) {
2139 rmId[i]=(UShort_t)esdTrack->GetID();
2144 Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2145 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2146 delete [] rmId; rmId=NULL;
2151 if(!vertexESD) return vertexAOD;
2152 if(vertexESD->GetNContributors()<=0) {
2153 //AliDebug(2,"vertexing failed");
2154 delete vertexESD; vertexESD=NULL;
2158 delete vertexer; vertexer=NULL;
2162 // convert to AliAODVertex
2163 Double_t pos[3],cov[6],chi2perNDF;
2164 vertexESD->GetXYZ(pos); // position
2165 vertexESD->GetCovMatrix(cov); //covariance matrix
2166 chi2perNDF = vertexESD->GetChi2toNDF();
2167 delete vertexESD; vertexESD=NULL;
2169 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2173 //-----------------------------------------------------------------------------
2174 void AliAnalysisVertexingHF::PrintStatus() const {
2175 // Print parameters being used
2177 //printf("Preselections:\n");
2178 // fTrackFilter->Dump();
2180 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2182 printf("Secondary vertex with AliVertexerTracks\n");
2184 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2185 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2187 printf("Reconstruct D0->Kpi candidates with cuts:\n");
2188 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
2191 printf("Reconstruct D*->D0pi candidates with cuts:\n");
2192 if(fFindVertexForDstar) {
2193 printf(" Reconstruct a secondary vertex for the D*\n");
2195 printf(" Assume the D* comes from the primary vertex\n");
2197 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
2200 printf("Reconstruct J/psi from B candidates with cuts:\n");
2201 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
2204 printf("Reconstruct 3 prong candidates.\n");
2205 printf(" D+->Kpipi cuts:\n");
2206 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
2207 printf(" Ds->KKpi cuts:\n");
2208 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
2209 printf(" Lc->pKpi cuts:\n");
2210 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
2213 printf("Reconstruct 4 prong candidates.\n");
2214 printf(" D0->Kpipipi cuts:\n");
2215 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
2218 printf("Reconstruct cascades candidates formed with v0s.\n");
2219 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
2220 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
2225 //-----------------------------------------------------------------------------
2226 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
2227 Double_t &dispersion,Bool_t useTRefArray) const
2229 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
2230 //AliCodeTimerAuto("",0);
2232 AliESDVertex *vertexESD = 0;
2233 AliAODVertex *vertexAOD = 0;
2235 if(!fSecVtxWithKF) { // AliVertexerTracks
2237 fVertexerTracks->SetVtxStart(fV1);
2238 vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2240 if(!vertexESD) return vertexAOD;
2242 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2243 //AliDebug(2,"vertexing failed");
2244 delete vertexESD; vertexESD=NULL;
2248 Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
2250 // vertex outside beam pipe, reject candidate to avoid propagation through material
2251 delete vertexESD; vertexESD=NULL;
2255 } else { // Kalman Filter vertexer (AliKFParticle)
2257 AliKFParticle::SetField(fBzkG);
2259 AliKFVertex vertexKF;
2261 Int_t nTrks = trkArray->GetEntriesFast();
2262 for(Int_t i=0; i<nTrks; i++) {
2263 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2264 AliKFParticle daughterKF(*esdTrack,211);
2265 vertexKF.AddDaughter(daughterKF);
2267 vertexESD = new AliESDVertex(vertexKF.Parameters(),
2268 vertexKF.CovarianceMatrix(),
2270 vertexKF.GetNContributors());
2274 // convert to AliAODVertex
2275 Double_t pos[3],cov[6],chi2perNDF;
2276 vertexESD->GetXYZ(pos); // position
2277 vertexESD->GetCovMatrix(cov); //covariance matrix
2278 chi2perNDF = vertexESD->GetChi2toNDF();
2279 dispersion = vertexESD->GetDispersion();
2280 delete vertexESD; vertexESD=NULL;
2282 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2283 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2287 //-----------------------------------------------------------------------------
2288 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(TObjArray *trkArray){
2289 // Invariant mass cut on tracks
2290 //AliCodeTimerAuto("",0);
2292 Int_t retval=kFALSE;
2293 Double_t momentum[3];
2294 Double_t px[3],py[3],pz[3];
2295 for(Int_t iTrack=0; iTrack<3; iTrack++){
2296 AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2297 track->GetPxPyPz(momentum);
2298 px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2300 retval = SelectInvMassAndPt3prong(px,py,pz);
2305 //-----------------------------------------------------------------------------
2306 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(TObjArray *trkArray){
2307 // Invariant mass cut on tracks
2308 //AliCodeTimerAuto("",0);
2310 Int_t retval=kFALSE;
2311 Double_t momentum[3];
2312 Double_t px[4],py[4],pz[4];
2314 for(Int_t iTrack=0; iTrack<4; iTrack++){
2315 AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2316 track->GetPxPyPz(momentum);
2317 px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2320 retval = SelectInvMassAndPt4prong(px,py,pz);
2324 //-----------------------------------------------------------------------------
2325 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(TObjArray *trkArray){
2326 // Invariant mass cut on tracks
2327 //AliCodeTimerAuto("",0);
2329 Int_t retval=kFALSE;
2330 Double_t momentum[3];
2331 Double_t px[2],py[2],pz[2];
2333 for(Int_t iTrack=0; iTrack<2; iTrack++){
2334 AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2335 track->GetPxPyPz(momentum);
2336 px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2338 retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
2342 //-----------------------------------------------------------------------------
2343 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtD0Kpi(Double_t *px,
2346 // Check invariant mass cut and pt candidate cut
2347 //AliCodeTimerAuto("",0);
2351 Double_t minv2,mrange;
2352 Double_t lolim,hilim;
2354 Bool_t retval=kFALSE;
2356 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2357 fOKInvMassD0=kFALSE;
2359 minPt=fCutsD0toKpi->GetMinPtCandidate();
2361 if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2363 mrange=fCutsD0toKpi->GetMassCut();
2364 lolim=fMassDzero-mrange;
2365 hilim=fMassDzero+mrange;
2366 pdg2[0]=211; pdg2[1]=321;
2367 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2368 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2372 pdg2[0]=321; pdg2[1]=211;
2373 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2374 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2381 //-----------------------------------------------------------------------------
2382 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtJpsiee(Double_t *px,
2385 // Check invariant mass cut and pt candidate cut
2386 //AliCodeTimerAuto("",0);
2390 Double_t minv2,mrange;
2391 Double_t lolim,hilim;
2393 Bool_t retval=kFALSE;
2395 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2396 fOKInvMassJpsi=kFALSE;
2398 minPt=fCutsJpsitoee->GetMinPtCandidate();
2400 if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2402 mrange=fCutsJpsitoee->GetMassCut();
2403 lolim=fMassJpsi-mrange;
2404 hilim=fMassJpsi+mrange;
2406 pdg2[0]=11; pdg2[1]=11;
2407 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2408 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2410 fOKInvMassJpsi=kTRUE;
2415 //-----------------------------------------------------------------------------
2416 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(Double_t *px,
2420 // Check invariant mass cut and pt candidate cut
2421 //AliCodeTimerAuto("",0);
2425 Double_t minv2,mrange;
2426 Double_t lolim,hilim;
2428 Bool_t retval=kFALSE;
2431 fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2432 fOKInvMassDplus=kFALSE;
2433 fOKInvMassDs=kFALSE;
2434 fOKInvMassLc=kFALSE;
2436 minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
2437 minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
2439 if(fMassCalc3->Pt2() < minPt*minPt) return retval;
2441 mrange=fCutsDplustoKpipi->GetMassCut();
2442 lolim=fMassDplus-mrange;
2443 hilim=fMassDplus+mrange;
2444 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2445 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2446 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2448 fOKInvMassDplus=kTRUE;
2451 mrange=fCutsDstoKKpi->GetMassCut();
2452 lolim=fMassDs-mrange;
2453 hilim=fMassDs+mrange;
2454 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2455 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2456 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2460 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2461 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2462 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2467 mrange=fCutsLctopKpi->GetMassCut();
2468 lolim=fMassLambdaC-mrange;
2469 hilim=fMassLambdaC+mrange;
2471 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2472 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2473 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2479 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2480 minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2481 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2490 //-----------------------------------------------------------------------------
2491 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(Double_t *px,
2494 // Check invariant mass cut and pt candidate cut
2495 //AliCodeTimerAuto("",0);
2499 Double_t minv2,mrange;
2500 Double_t lolim,hilim;
2502 Bool_t retval=kFALSE;
2504 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2505 fOKInvMassDstar=kFALSE;
2507 minPt=fCutsDStartoKpipi->GetMinPtCandidate();
2509 if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2511 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2512 mrange=fCutsDStartoKpipi->GetMassCut();
2513 lolim=fMassDstar-mrange;
2514 hilim=fMassDstar+mrange;
2515 minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2516 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2518 fOKInvMassDstar=kTRUE;
2524 //-----------------------------------------------------------------------------
2525 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(Double_t *px,
2528 // Check invariant mass cut and pt candidate cut
2529 //AliCodeTimerAuto("",0);
2533 Double_t minv2,mrange;
2534 Double_t lolim,hilim;
2536 Bool_t retval=kFALSE;
2538 // D0->Kpipipi without PID
2539 fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2540 fOKInvMassD0to4p=kFALSE;
2542 minPt=fCutsD0toKpipipi->GetMinPtCandidate();
2544 if(fMassCalc4->Pt2() < minPt*minPt) return retval;
2546 mrange=fCutsD0toKpipipi->GetMassCut();
2547 lolim=fMassDzero-mrange;
2548 hilim=fMassDzero+mrange;
2550 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2551 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2552 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2554 fOKInvMassD0to4p=kTRUE;
2557 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2558 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2559 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2561 fOKInvMassD0to4p=kTRUE;
2564 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2565 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2566 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2568 fOKInvMassD0to4p=kTRUE;
2571 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2572 minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2573 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2575 fOKInvMassD0to4p=kTRUE;
2580 //-----------------------------------------------------------------------------
2581 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtCascade(Double_t *px,
2584 // Check invariant mass cut and pt candidate cut
2585 //AliCodeTimerAuto("",0);
2589 Double_t minv2,mrange;
2590 Double_t lolim,hilim;
2591 // Double_t minPt=0;
2592 Bool_t retval=kFALSE;
2594 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2595 // minPt=fCutsLctoV0->GetMinPtCandidate();
2596 fOKInvMassLctoV0=kFALSE;
2597 mrange=fCutsLctoV0->GetMassCut();
2598 lolim=fMassLambdaC-mrange;
2599 hilim=fMassLambdaC+mrange;
2600 pdg2[0]=2212;pdg2[1]=310;
2601 minv2=fMassCalc2->InvMass2(2,pdg2);
2602 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2604 fOKInvMassLctoV0=kTRUE;
2606 pdg2[0]=211;pdg2[1]=3122;
2607 minv2=fMassCalc2->InvMass2(2,pdg2);
2608 if(minv2>lolim*lolim && minv2<hilim*hilim ){
2610 fOKInvMassLctoV0=kTRUE;
2615 //-----------------------------------------------------------------------------
2616 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2618 TObjArray &seleTrksArray,
2619 TObjArray &tracksAtVertex,
2621 UChar_t *seleFlags,Int_t *evtNumber)
2623 // Apply single-track preselection.
2624 // Fill a TObjArray with selected tracks (for displaced vertices or
2625 // soft pion from D*). Selection flag stored in seleFlags.
2626 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2627 // In case of AOD input, also fill fAODMap for track index<->ID
2628 //AliCodeTimerAuto("",0);
2630 const AliVVertex *vprimary = event->GetPrimaryVertex();
2632 if(fV1) { delete fV1; fV1=NULL; }
2633 if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
2636 UShort_t *indices = 0;
2637 Double_t pos[3],cov[6];
2638 const Int_t entries = event->GetNumberOfTracks();
2639 AliCentrality* cent;
2641 if(!fInputAOD) { // ESD
2642 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2643 cent=((AliESDEvent*)event)->GetCentrality();
2645 vprimary->GetXYZ(pos);
2646 vprimary->GetCovarianceMatrix(cov);
2647 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2648 if(entries<=0) return;
2649 indices = new UShort_t[entries];
2650 memset(indices,0,sizeof(UShort_t)*entries);
2651 fAODMapSize = 100000;
2652 fAODMap = new Int_t[fAODMapSize];
2653 memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2654 cent=((AliAODEvent*)event)->GetCentrality();
2656 Float_t centperc=cent->GetCentralityPercentile("V0M");
2658 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE;
2661 // transfer ITS tracks from event to arrays
2662 for(Int_t i=0; i<entries; i++) {
2664 track = (AliVTrack*)event->GetTrack(i);
2666 // skip pure ITS SA tracks
2667 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2669 // skip tracks without ITS
2670 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2672 // skip tracks with negative ID
2673 // (these are duplicated TPC-only AOD tracks, for jet analysis...)
2674 if(track->GetID()<0) continue;
2676 // TEMPORARY: check that the cov matrix is there
2677 Double_t covtest[21];
2678 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2682 AliAODTrack *aodt = (AliAODTrack*)track;
2683 if(aodt->GetUsedForPrimVtxFit()) {
2684 indices[nindices]=aodt->GetID(); nindices++;
2686 Int_t ind = (Int_t)aodt->GetID();
2687 if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2690 AliESDtrack *esdt = 0;
2693 esdt = (AliESDtrack*)track;
2695 esdt = new AliESDtrack(track);
2698 // single track selection
2699 okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE;
2700 if(fMixEvent && i<trkEntries){
2701 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2702 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2703 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2704 eventVtx->GetXYZ(vtxPos);
2705 vprimary->GetXYZ(primPos);
2706 eventVtx->GetCovarianceMatrix(primCov);
2707 for(Int_t ind=0;ind<3;ind++){
2708 trasl[ind]=vtxPos[ind]-primPos[ind];
2711 Bool_t isTransl=esdt->Translate(trasl,primCov);
2719 if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong) && nSeleTrks<trkEntries) {
2720 esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
2721 seleTrksArray.AddLast(esdt);
2722 tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
2723 seleFlags[nSeleTrks]=0;
2724 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2725 if(okFor3Prong) SETBIT(seleFlags[nSeleTrks],kBit3Prong);
2726 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2729 SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
2730 SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2731 SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2732 Bool_t useTPC=kTRUE;
2734 Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
2735 if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
2736 CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
2738 Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
2739 if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
2740 CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2742 Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
2743 if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
2744 CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2746 if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
2748 if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
2749 Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
2750 if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
2751 CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
2753 Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
2754 if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
2755 CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2757 Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
2758 if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
2759 CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2764 if(fInputAOD) delete esdt;
2769 } // end loop on tracks
2771 // primary vertex from AOD
2773 delete fV1; fV1=NULL;
2774 vprimary->GetXYZ(pos);
2775 vprimary->GetCovarianceMatrix(cov);
2776 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2777 Int_t ncontr=nindices;
2778 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2779 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2780 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2781 fV1->SetTitle(vprimary->GetTitle());
2782 fV1->SetIndices(nindices,indices);
2784 if(indices) { delete [] indices; indices=NULL; }
2789 //-----------------------------------------------------------------------------
2790 void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
2792 // Set the selection bit for PID
2794 //AliCodeTimerAuto("",0);
2795 if(cuts->GetPidHF()) {
2796 Bool_t usepid=cuts->GetIsUsePID();
2797 cuts->SetUsePID(kTRUE);
2798 if(cuts->IsSelectedPID(rd))
2799 rd->SetSelectionBit(bit);
2800 cuts->SetUsePID(usepid);
2804 //-----------------------------------------------------------------------------
2805 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2806 Float_t centralityperc,
2807 Bool_t &okDisplaced,
2809 Bool_t &okFor3Prong) const
2811 // Check if track passes some kinematical cuts
2813 // this is needed to store the impact parameters
2814 //AliCodeTimerAuto("",0);
2816 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2820 // Track selection, displaced tracks -- 2 prong
2822 if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
2823 && fTrackFilter2prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2824 // central PbPb events, tighter cuts
2825 selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
2829 selectInfo = fTrackFilter->IsSelected(trk);
2832 if(selectInfo) okDisplaced=kTRUE;
2834 // Track selection, displaced tracks -- 3 prong
2836 if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
2837 && fTrackFilter3prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2838 // central PbPb events, tighter cuts
2839 selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
2843 selectInfo = fTrackFilter->IsSelected(trk);
2846 if(selectInfo) okFor3Prong=kTRUE;
2848 // Track selection, soft pions
2850 if(fDstar && fTrackFilterSoftPi) {
2851 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2853 if(selectInfo) okSoftPi=kTRUE;
2855 if(okDisplaced || okSoftPi || okFor3Prong) return kTRUE;
2861 //-----------------------------------------------------------------------------
2862 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2864 // Transform ESDv0 to AODv0
2866 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2867 // and creates an AODv0 out of them
2869 //AliCodeTimerAuto("",0);
2870 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2871 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2873 // create the v0 neutral track to compute the DCA to the primary vertex
2874 Double_t xyz[3], pxpypz[3];
2876 esdV0->PxPyPz(pxpypz);
2877 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2878 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2883 Double_t d0z0[2],covd0z0[3];
2884 AliAODVertex *primVertexAOD = PrimaryVertex();
2885 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2886 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2887 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2888 Double_t dcaV0DaughterToPrimVertex[2];
2889 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2890 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2891 if( !posV0track || !negV0track) {
2892 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2894 delete primVertexAOD;
2897 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2898 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2900 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2901 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2902 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2904 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2905 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2906 Double_t pmom[3],nmom[3];
2907 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2908 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2910 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2911 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2914 delete primVertexAOD;
2918 //-----------------------------------------------------------------------------
2919 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
2920 // Set the stored track parameters at primary vertex into AliESDtrack
2921 //AliCodeTimerAuto("",0);
2923 const Double_t *par=extpar->GetParameter();
2924 const Double_t *cov=extpar->GetCovariance();
2925 Double_t alpha=extpar->GetAlpha();
2926 Double_t x=extpar->GetX();
2927 esdt->Set(x,alpha,par,cov);
2930 //-----------------------------------------------------------------------------
2931 void AliAnalysisVertexingHF::SetMasses(){
2932 // Set the hadron mass values from TDatabasePDG
2934 fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2935 fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2936 fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2937 fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2938 fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2939 fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2941 //-----------------------------------------------------------------------------
2942 Bool_t AliAnalysisVertexingHF::CheckCutsConsistency(){
2944 // Check the Vertexer and the analysts task consitstecny
2948 //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
2951 if ( fCutsLctoV0 && fV0TypeForCascadeVertex != fCutsLctoV0->GetV0Type())
2953 printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
2958 //-----------------------------------------------------------------------------