AliPhysics  d497afb (d497afb)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisVertexingHF.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 /* $Id$ */
17 
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.
24 //
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>
30 #include <TFile.h>
31 #include <TDatabasePDG.h>
32 #include <TString.h>
33 #include <TList.h>
34 #include <TProcessID.h>
35 #include "AliLog.h"
36 #include "AliVEvent.h"
37 #include "AliVVertex.h"
38 #include "AliVTrack.h"
39 #include "AliVertexerTracks.h"
40 #include "AliKFVertex.h"
41 #include "AliESDEvent.h"
42 #include "AliESDVertex.h"
43 #include "AliExternalTrackParam.h"
44 #include "AliNeutralTrackParam.h"
45 #include "AliESDtrack.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliAODEvent.h"
48 #include "AliPIDResponse.h"
49 #include "AliAODRecoDecay.h"
50 #include "AliAODRecoDecayHF.h"
54 #include "AliAODRecoCascadeHF.h"
55 #include "AliRDHFCutsD0toKpi.h"
56 #include "AliRDHFCutsJpsitoee.h"
59 #include "AliRDHFCutsDstoK0sK.h"
60 #include "AliRDHFCutsDstoKKpi.h"
61 #include "AliRDHFCutsLctopKpi.h"
62 #include "AliRDHFCutsLctoV0.h"
63 #include "AliRDHFCutsD0toKpipipi.h"
65 #include "AliAnalysisFilter.h"
66 #include "AliAnalysisVertexingHF.h"
67 #include "AliMixedEvent.h"
68 #include "AliESDv0.h"
69 #include "AliAODv0.h"
70 #include "AliCodeTimer.h"
71 #include "AliMultSelection.h"
72 #include <cstring>
73 
75 ClassImp(AliAnalysisVertexingHF);
77 
78 //----------------------------------------------------------------------------
80 fInputAOD(kFALSE),
81 fAODMapSize(0),
82 fAODMap(0),
83 fVertexerTracks(0x0),
84 fBzkG(0.),
85 fSecVtxWithKF(kFALSE),
86 fRecoPrimVtxSkippingTrks(kFALSE),
87 fRmTrksFromPrimVtx(kFALSE),
88 fV1(0x0),
89 fD0toKpi(kTRUE),
90 fJPSItoEle(kTRUE),
91 f3Prong(kTRUE),
92 f4Prong(kTRUE),
93 fDstar(kTRUE),
94 fCascades(kTRUE),
95 fLikeSign(kFALSE),
96 fLikeSign3prong(kFALSE),
97 fMixEvent(kFALSE),
98 fPidResponse(0x0),
99 fUseKaonPIDfor3Prong(kFALSE),
100 fUsePIDforLc(0),
101 fUsePIDforLc2V0(0),
102 fUseKaonPIDforDs(kFALSE),
103 fUseTPCPID(kFALSE),
104 fUseTOFPID(kFALSE),
105 fUseTPCPIDOnlyIfNoTOF(kFALSE),
106 fMaxMomForTPCPid(1.),
107 fUsePidTag(kFALSE),
108 fnSigmaTPCPionLow(5.),
109 fnSigmaTPCPionHi(5.),
110 fnSigmaTOFPionLow(5.),
111 fnSigmaTOFPionHi(5.),
112 fnSigmaTPCKaonLow(5.),
113 fnSigmaTPCKaonHi(5.),
114 fnSigmaTOFKaonLow(5.),
115 fnSigmaTOFKaonHi(5.),
116 fnSigmaTPCProtonLow(5.),
117 fnSigmaTPCProtonHi(5.),
118 fnSigmaTOFProtonLow(5.),
119 fnSigmaTOFProtonHi(5.),
120 fMaxCentPercentileForTightCuts(-9999),
121 fTrackFilter(0x0),
122 fTrackFilter2prongCentral(0x0),
123 fTrackFilter3prongCentral(0x0),
124 fTrackFilterSoftPi(0x0),
125 fTrackFilterBachelor(0x0),
126 fCutsD0toKpi(0x0),
127 fCutsJpsitoee(0x0),
128 fCutsDplustoK0spi(0x0),
129 fCutsDplustoKpipi(0x0),
130 fCutsDstoK0sK(0x0),
131 fCutsDstoKKpi(0x0),
132 fCutsLctopKpi(0x0),
133 fCutsLctoV0(0x0),
134 fCutsD0toKpipipi(0x0),
135 fCutsDStartoKpipi(0x0),
136 fListOfCuts(0x0),
137 fFindVertexForDstar(kTRUE),
138 fFindVertexForCascades(kTRUE),
139 fV0TypeForCascadeVertex(0),
140 fMassCutBeforeVertexing(kFALSE),
141 fMassCalc2(0),
142 fMassCalc3(0),
143 fMassCalc4(0),
144 fOKInvMassD0(kFALSE),
145 fOKInvMassJpsi(kFALSE),
146 fOKInvMassDplus(kFALSE),
147 fOKInvMassDs(kFALSE),
148 fOKInvMassLc(kFALSE),
149 fOKInvMassDstar(kFALSE),
150 fOKInvMassD0to4p(kFALSE),
151 fOKInvMassLctoV0(kFALSE),
152 fnTrksTotal(0),
153 fnSeleTrksTotal(0),
154 fMakeReducedRHF(kFALSE),
155 fMassDzero(0.),
156 fMassDplus(0.),
157 fMassDs(0.),
158 fMassLambdaC(0.),
159 fMassDstar(0.),
160 fMassJpsi(0.)
161 {
163 
164  Double_t d02[2]={0.,0.};
165  Double_t d03[3]={0.,0.,0.};
166  Double_t d04[4]={0.,0.,0.,0.};
167  fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
168  fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
169  fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
170  SetMasses();
171 }
172 //--------------------------------------------------------------------------
174 TNamed(source),
175 fInputAOD(source.fInputAOD),
176 fAODMapSize(source.fAODMapSize),
177 fAODMap(source.fAODMap),
178 fVertexerTracks(source.fVertexerTracks),
179 fBzkG(source.fBzkG),
180 fSecVtxWithKF(source.fSecVtxWithKF),
181 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
182 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
183 fV1(source.fV1),
184 fD0toKpi(source.fD0toKpi),
185 fJPSItoEle(source.fJPSItoEle),
186 f3Prong(source.f3Prong),
187 f4Prong(source.f4Prong),
188 fDstar(source.fDstar),
189 fCascades(source.fCascades),
190 fLikeSign(source.fLikeSign),
191 fLikeSign3prong(source.fLikeSign3prong),
192 fMixEvent(source.fMixEvent),
193 fPidResponse(source.fPidResponse),
194 fUseKaonPIDfor3Prong(source.fUseKaonPIDfor3Prong),
195 fUsePIDforLc(source.fUsePIDforLc),
196 fUsePIDforLc2V0(source.fUsePIDforLc2V0),
197 fUseKaonPIDforDs(source.fUseKaonPIDforDs),
198 fUseTPCPID(source.fUseTPCPID),
199 fUseTOFPID(source.fUseTOFPID),
200 fUseTPCPIDOnlyIfNoTOF(source.fUseTPCPIDOnlyIfNoTOF),
201 fMaxMomForTPCPid(source.fMaxMomForTPCPid),
202 fUsePidTag(source.fUsePidTag),
203 fnSigmaTPCPionLow(source.fnSigmaTPCPionLow),
204 fnSigmaTPCPionHi(source.fnSigmaTPCPionHi),
205 fnSigmaTOFPionLow(source.fnSigmaTOFPionLow),
206 fnSigmaTOFPionHi(source.fnSigmaTOFPionHi),
207 fnSigmaTPCKaonLow(source.fnSigmaTPCKaonLow),
208 fnSigmaTPCKaonHi(source.fnSigmaTPCKaonHi),
209 fnSigmaTOFKaonLow(source.fnSigmaTOFKaonLow),
210 fnSigmaTOFKaonHi(source.fnSigmaTOFKaonHi),
211 fnSigmaTPCProtonLow(source.fnSigmaTPCProtonLow),
212 fnSigmaTPCProtonHi(source.fnSigmaTPCProtonHi),
213 fnSigmaTOFProtonLow(source.fnSigmaTOFProtonLow),
214 fnSigmaTOFProtonHi(source.fnSigmaTOFProtonHi),
215 fMaxCentPercentileForTightCuts(source.fMaxCentPercentileForTightCuts),
216 fTrackFilter(source.fTrackFilter),
217 fTrackFilter2prongCentral(source.fTrackFilter2prongCentral),
218 fTrackFilter3prongCentral(source.fTrackFilter3prongCentral),
219 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
220 fTrackFilterBachelor(source.fTrackFilterBachelor),
221 fCutsD0toKpi(source.fCutsD0toKpi),
222 fCutsJpsitoee(source.fCutsJpsitoee),
223 fCutsDplustoK0spi(source.fCutsDplustoK0spi),
224 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
225 fCutsDstoK0sK(source.fCutsDstoK0sK),
226 fCutsDstoKKpi(source.fCutsDstoKKpi),
227 fCutsLctopKpi(source.fCutsLctopKpi),
228 fCutsLctoV0(source.fCutsLctoV0),
229 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
230 fCutsDStartoKpipi(source.fCutsDStartoKpipi),
231 fListOfCuts(source.fListOfCuts),
232 fFindVertexForDstar(source.fFindVertexForDstar),
233 fFindVertexForCascades(source.fFindVertexForCascades),
234 fV0TypeForCascadeVertex(source.fV0TypeForCascadeVertex),
235 fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
236 fMassCalc2(source.fMassCalc2),
237 fMassCalc3(source.fMassCalc3),
238 fMassCalc4(source.fMassCalc4),
239 fOKInvMassD0(source.fOKInvMassD0),
240 fOKInvMassJpsi(source.fOKInvMassJpsi),
241 fOKInvMassDplus(source.fOKInvMassDplus),
242 fOKInvMassDs(source.fOKInvMassDs),
243 fOKInvMassLc(source.fOKInvMassLc),
244 fOKInvMassDstar(source.fOKInvMassDstar),
245 fOKInvMassD0to4p(source.fOKInvMassD0to4p),
246 fOKInvMassLctoV0(source.fOKInvMassLctoV0),
247 fnTrksTotal(0),
248 fnSeleTrksTotal(0),
249 fMakeReducedRHF(kFALSE),
250 fMassDzero(source.fMassDzero),
251 fMassDplus(source.fMassDplus),
252 fMassDs(source.fMassDs),
253 fMassLambdaC(source.fMassLambdaC),
254 fMassDstar(source.fMassDstar),
255 fMassJpsi(source.fMassJpsi)
256 {
260 }
261 //--------------------------------------------------------------------------
263 {
267  if(&source == this) return *this;
268  fInputAOD = source.fInputAOD;
269  fAODMapSize = source.fAODMapSize;
271  fBzkG = source.fBzkG;
272  fSecVtxWithKF = source.fSecVtxWithKF;
275  fV1 = source.fV1;
276  fD0toKpi = source.fD0toKpi;
277  fJPSItoEle = source.fJPSItoEle;
278  f3Prong = source.f3Prong;
279  f4Prong = source.f4Prong;
280  fDstar = source.fDstar;
281  fCascades = source.fCascades;
282  fLikeSign = source.fLikeSign;
284  fMixEvent = source.fMixEvent;
285  fPidResponse = source.fPidResponse;
287  fUsePIDforLc = source.fUsePIDforLc;
290  fUseTPCPID = source.fUseTPCPID;
291  fUseTOFPID = source.fUseTOFPID;
294  fUsePidTag = source.fUsePidTag;
310  fTrackFilter = source.fTrackFilter;
315  fCutsD0toKpi = source.fCutsD0toKpi;
316  fCutsJpsitoee = source.fCutsJpsitoee;
319  fCutsDstoK0sK = source.fCutsDstoK0sK;
320  fCutsDstoKKpi = source.fCutsDstoKKpi;
321  fCutsLctopKpi = source.fCutsLctopKpi;
322  fCutsLctoV0 = source.fCutsLctoV0;
325  fListOfCuts = source.fListOfCuts;
330  fMassCalc2 = source.fMassCalc2;
331  fMassCalc3 = source.fMassCalc3;
332  fMassCalc4 = source.fMassCalc4;
333  fOKInvMassD0 = source.fOKInvMassD0;
336  fOKInvMassDs = source.fOKInvMassDs;
337  fOKInvMassLc = source.fOKInvMassLc;
341  fMassDzero = source.fMassDzero;
342  fMassDplus = source.fMassDplus;
343  fMassDs = source.fMassDs;
344  fMassLambdaC = source.fMassLambdaC;
345  fMassDstar = source.fMassDstar;
346  fMassJpsi = source.fMassJpsi;
347 
348  return *this;
349 }
350 //----------------------------------------------------------------------------
353  if(fV1) { delete fV1; fV1=0; }
354  delete fVertexerTracks;
355  if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
360  if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
361  if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
364  if(fCutsDstoK0sK) { delete fCutsDstoK0sK; fCutsDstoK0sK=0; }
365  if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
366  if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
367  if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
370  if(fAODMap) { delete [] fAODMap; fAODMap=0; }
371  if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
372  if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
373  if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
374 }
375 //----------------------------------------------------------------------------
378 
379  TList *list = new TList();
380  list->SetOwner();
381  list->SetName("ListOfCuts");
382 
383  if(fCutsD0toKpi) {
385  list->Add(cutsD0toKpi);
386  }
387  if(fCutsJpsitoee) {
389  list->Add(cutsJpsitoee);
390  }
391  if (fCutsDplustoK0spi) {
393  list->Add(cutsDplustoK0spi);
394  }
395  if(fCutsDplustoKpipi) {
397  list->Add(cutsDplustoKpipi);
398  }
399  if (fCutsDstoK0sK) {
401  list->Add(cutsDstoK0sK);
402  }
403  if(fCutsDstoKKpi) {
405  list->Add(cutsDstoKKpi);
406  }
407  if(fCutsLctopKpi) {
409  list->Add(cutsLctopKpi);
410  }
411  if(fCutsLctoV0){
412  AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
413  list->Add(cutsLctoV0);
414  }
415  if(fCutsD0toKpipipi) {
417  list->Add(cutsD0toKpipipi);
418  }
419  if(fCutsDStartoKpipi) {
421  list->Add(cutsDStartoKpipi);
422  }
423 
424  //___ Check consitstency of cuts between vertexer and analysis tasks
425  Bool_t bCutsOk = CheckCutsConsistency();
426  if (bCutsOk == kFALSE) {AliFatal("AliAnalysisVertexingHF::FillListOfCuts vertexing and the analysis task cuts are not consistent!");}
427 
428  // keep a pointer to the list
429  fListOfCuts = list;
430 
431  return list;
432 }
433 //----------------------------------------------------------------------------
435  TClonesArray *aodVerticesHFTClArr,
436  TClonesArray *aodD0toKpiTClArr,
437  TClonesArray *aodJPSItoEleTClArr,
438  TClonesArray *aodCharm3ProngTClArr,
439  TClonesArray *aodCharm4ProngTClArr,
440  TClonesArray *aodDstarTClArr,
441  TClonesArray *aodCascadesTClArr,
442  TClonesArray *aodLikeSign2ProngTClArr,
443  TClonesArray *aodLikeSign3ProngTClArr)
444 {
448  //AliCodeTimerAuto("",0);
449 
450  if(!fMixEvent){
451  TString evtype = event->IsA()->GetName();
452  fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
453  } // if we do mixing AliVEvent is a AliMixedEvent
454 
455  if(fInputAOD) {
456  AliDebug(2,"Creating HF candidates from AOD");
457  } else {
458  AliDebug(2,"Creating HF candidates from ESD");
459  }
460 
461  if(!aodVerticesHFTClArr) {
462  printf("ERROR: no aodVerticesHFTClArr");
463  return;
464  }
465  if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
466  printf("ERROR: no aodD0toKpiTClArr");
467  return;
468  }
469  if(fJPSItoEle && !aodJPSItoEleTClArr) {
470  printf("ERROR: no aodJPSItoEleTClArr");
471  return;
472  }
473  if(f3Prong && !aodCharm3ProngTClArr) {
474  printf("ERROR: no aodCharm3ProngTClArr");
475  return;
476  }
477  if(f4Prong && !aodCharm4ProngTClArr) {
478  printf("ERROR: no aodCharm4ProngTClArr");
479  return;
480  }
481  if(fDstar && !aodDstarTClArr) {
482  printf("ERROR: no aodDstarTClArr");
483  return;
484  }
485  if(fCascades && !aodCascadesTClArr){
486  printf("ERROR: no aodCascadesTClArr ");
487  return;
488  }
489  if(fLikeSign && !aodLikeSign2ProngTClArr) {
490  printf("ERROR: no aodLikeSign2ProngTClArr");
491  return;
492  }
493  if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
494  printf("ERROR: no aodLikeSign3ProngTClArr");
495  return;
496  }
497 
498  // delete candidates from previous event and create references
499  Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
500  aodVerticesHFTClArr->Delete();
501  iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
502  TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
503  if(fD0toKpi || fDstar) {
504  aodD0toKpiTClArr->Delete();
505  iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
506  }
507  if(fJPSItoEle) {
508  aodJPSItoEleTClArr->Delete();
509  iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
510  }
511  if(f3Prong) {
512  aodCharm3ProngTClArr->Delete();
513  i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
514  }
515  if(f4Prong) {
516  aodCharm4ProngTClArr->Delete();
517  i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
518  }
519  if(fDstar) {
520  aodDstarTClArr->Delete();
521  iDstar = aodDstarTClArr->GetEntriesFast();
522  }
523  if(fCascades) {
524  aodCascadesTClArr->Delete();
525  iCascades = aodCascadesTClArr->GetEntriesFast();
526  }
527  if(fLikeSign) {
528  aodLikeSign2ProngTClArr->Delete();
529  iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
530  }
531  if(fLikeSign3prong && f3Prong) {
532  aodLikeSign3ProngTClArr->Delete();
533  iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
534  }
535 
536  TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
537  TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
538  TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
539  TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
540  TClonesArray &aodDstarRef = *aodDstarTClArr;
541  TClonesArray &aodCascadesRef = *aodCascadesTClArr;
542  TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
543  TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
544 
545 
546  AliAODRecoDecayHF2Prong *io2Prong = 0;
547  AliAODRecoDecayHF3Prong *io3Prong = 0;
548  AliAODRecoDecayHF4Prong *io4Prong = 0;
549  AliAODRecoCascadeHF *ioCascade = 0;
550 
551  Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
552  Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
553  Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
554  Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
555  Bool_t okCascades=kFALSE;
556  AliESDtrack *postrack1 = 0;
557  AliESDtrack *postrack2 = 0;
558  AliESDtrack *negtrack1 = 0;
559  AliESDtrack *negtrack2 = 0;
560  AliESDtrack *trackPi = 0;
561  Double_t mompos1[3],mompos2[3],momneg1[3],momneg2[3];
562  // AliESDtrack *posV0track = 0;
563  // AliESDtrack *negV0track = 0;
564  Float_t dcaMax = fCutsD0toKpi->GetDCACut();
565  if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
566  if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
567  if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
568  if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
569  if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
570  if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
571 
572  AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
573 
574 
575  // get Bz
576  fBzkG = (Double_t)event->GetMagneticField();
577  if(!fVertexerTracks){
578  fVertexerTracks=new AliVertexerTracks(fBzkG);
579  }else{
580  Double_t oldField=fVertexerTracks->GetFieldkG();
581  if(oldField!=fBzkG) fVertexerTracks->SetFieldkG(fBzkG);
582  }
583 
584  trkEntries = (Int_t)event->GetNumberOfTracks();
585  AliDebug(1,Form(" Number of tracks: %d",trkEntries));
586  fnTrksTotal += trkEntries;
587 
588  nv0 = (Int_t)event->GetNumberOfV0s();
589  AliDebug(1,Form(" Number of V0s: %d",nv0));
590 
591  if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
592  AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
593  return;
594  }
595 
596  // event selection + PID configuration
597  if(!fCutsD0toKpi->IsEventSelected(event)) return;
604  if(fCutsLctoV0) fCutsLctoV0->SetupPID(event);
607 
608  // call function that applies sigle-track selection,
609  // for displaced tracks and soft pions (both charges) for D*,
610  // and retrieves primary vertex
611  TObjArray seleTrksArray(trkEntries);
612  TObjArray tracksAtVertex(trkEntries);
613  UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi, bit 2: 3 prong, bits 3-4-5: for PID
614  Int_t nSeleTrks=0;
615  Int_t *evtNumber = new Int_t[trkEntries];
616  SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
617 
618  AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
619  fnSeleTrksTotal += nSeleTrks;
620 
621 
622  TObjArray *twoTrackArray1 = new TObjArray(2);
623  TObjArray *twoTrackArray2 = new TObjArray(2);
624  TObjArray *twoTrackArrayV0 = new TObjArray(2);
625  TObjArray *twoTrackArrayCasc = new TObjArray(2);
626  TObjArray *threeTrackArray = new TObjArray(3);
627  TObjArray *fourTrackArray = new TObjArray(4);
628 
629  Double_t dispersion;
630  Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
631 
632  AliAODRecoDecayHF *rd = 0;
633  AliAODRecoCascadeHF *rc = 0;
634  AliAODv0 *v0 = 0;
635  AliESDv0 *esdV0 = 0;
636 
637  Bool_t massCutOK=kTRUE;
638 
639  // LOOP ON POSITIVE TRACKS
640  for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
641 
642  //if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
643  //if(iTrkP1%1==0) printf(" 1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);
644 
645  // get track from tracks array
646  postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
647  postrack1->GetPxPyPz(mompos1);
648 
649  // Make cascades with V0+track
650  //
651  if(fCascades) {
652  // loop on V0's
653  for(iv0=0; iv0<nv0; iv0++){
654 
655  //AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
656  if ( !TESTBIT(seleFlags[iTrkP1],kBitBachelor) ) continue;
657 
658  if ( fUsePIDforLc2V0 && !TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) ) continue; //clm
659 
660  // Get the V0
661  if(fInputAOD) {
662  v0 = ((AliAODEvent*)event)->GetV0(iv0);
663  } else {
664  esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
665  }
666  if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
667  (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
668 
669  if ( v0 && ((v0->GetOnFlyStatus() == kTRUE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOfflineV0s) ||
670  (v0->GetOnFlyStatus() == kFALSE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOnTheFlyV0s)) ) continue;
671 
672  if ( esdV0 && ((esdV0->GetOnFlyStatus() == kTRUE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOfflineV0s) ||
673  ( esdV0->GetOnFlyStatus() == kFALSE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOnTheFlyV0s)) ) continue;
674 
675  // Get the tracks that form the V0
676  // ( parameters at primary vertex )
677  // and define an AliExternalTrackParam out of them
678  AliExternalTrackParam * posV0track;
679  AliExternalTrackParam * negV0track;
680 
681  if(fInputAOD){
682  AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
683  AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
684  if( !posVV0track || !negVV0track ) continue;
685  //
686  // Apply some basic V0 daughter criteria
687  //
688  // bachelor must not be a v0-track
689  if (posVV0track->GetID() == postrack1->GetID() ||
690  negVV0track->GetID() == postrack1->GetID()) continue;
691  // reject like-sign v0
692  if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
693  // avoid ghost TPC tracks
694  if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
695  !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
696  // Get AliExternalTrackParam out of the AliAODTracks
697  Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
698  posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
699  posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
700  posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
701  negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
702  negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
703  negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
704  } else {
705  AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
706  AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
707  if( !posVV0track || !negVV0track ) continue;
708  //
709  // Apply some basic V0 daughter criteria
710  //
711  // bachelor must not be a v0-track
712  if (posVV0track->GetID() == postrack1->GetID() ||
713  negVV0track->GetID() == postrack1->GetID()) continue;
714  // reject like-sign v0
715  if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
716  // avoid ghost TPC tracks
717  if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
718  !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
719  // reject kinks (only necessary on AliESDtracks)
720  if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
721  // Get AliExternalTrackParam out of the AliESDtracks
722  posV0track = new AliExternalTrackParam(*posVV0track);
723  negV0track = new AliExternalTrackParam(*negVV0track);
724 
725  // Define the AODv0 from ESDv0 if reading ESDs
726  v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
727  }
728  if( !posV0track || !negV0track ){
729  AliDebug(1,Form(" Couldn't get the V0 daughters"));
730  continue;
731  }
732 
733  // fill in the v0 two-external-track-param array
734  twoTrackArrayV0->AddAt(posV0track,0);
735  twoTrackArrayV0->AddAt(negV0track,1);
736 
737  // Get the V0 dca
738  dcaV0 = v0->DcaV0Daughters();
739 
740  // Define the V0 (neutral) track
741  AliNeutralTrackParam *trackV0=NULL;
742  if(fInputAOD) {
743  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
744  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
745  } else {
746  Double_t xyz[3], pxpypz[3];
747  esdV0->XvYvZv(xyz);
748  esdV0->PxPyPz(pxpypz);
749  Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
750  trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
751  }
752 
753 
754  // Fill in the object array to create the cascade
755  twoTrackArrayCasc->AddAt(postrack1,0);
756  twoTrackArrayCasc->AddAt(trackV0,1);
757  // Compute the cascade vertex
758  AliAODVertex *vertexCasc = 0;
760  // DCA between the two tracks
761  dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
762  // Vertexing+
763  vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
764  } else {
765  // assume Cascade decays at the primary vertex
766  Double_t pos[3],cov[6],chi2perNDF;
767  fV1->GetXYZ(pos);
768  fV1->GetCovMatrix(cov);
769  chi2perNDF = fV1->GetChi2toNDF();
770  vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
771  dcaCasc = 0.;
772  }
773  if(!vertexCasc) {
774  delete posV0track; posV0track=NULL;
775  delete negV0track; negV0track=NULL;
776  delete trackV0; trackV0=NULL;
777  if(!fInputAOD) {delete v0; v0=NULL;}
778  twoTrackArrayV0->Clear();
779  twoTrackArrayCasc->Clear();
780  continue;
781  }
782 
783  // Create and store the Cascade if passed the cuts
784  ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
785  if(okCascades && ioCascade) {
786  //AliDebug(1,Form("Storing a cascade object... "));
787  // add the vertex and the cascade to the AOD
788  rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
789  if(fMakeReducedRHF){
790  UShort_t id[2]={(UShort_t)postrack1->GetID(),(UShort_t)iv0};
791  rc->SetProngIDs(2,id);
792  rc->DeleteRecoD();
793  }else{
794  AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
795  rc->SetSecondaryVtx(vCasc);
796  vCasc->SetParent(rc);
797  if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
798  AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
799  vCasc->AddDaughter(v0); // fill the 2prong V0
800  }
801  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
802  }
803 
804 
805  // Clean up
806  delete posV0track; posV0track=NULL;
807  delete negV0track; negV0track=NULL;
808  delete trackV0; trackV0=NULL;
809  twoTrackArrayV0->Clear();
810  twoTrackArrayCasc->Clear();
811  if(ioCascade) { delete ioCascade; ioCascade=NULL; }
812  if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
813  if(!fInputAOD) {delete v0; v0=NULL;}
814 
815  } // end loop on V0's
816 
817  // re-set parameters at vertex
818  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
819  } // end fCascades
820 
821  // If there is less than 2 particles continue
822  if(trkEntries<2) {
823  AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
824  continue;
825  }
826 
827  if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
828  if(postrack1->Charge()<0 && !fLikeSign) continue;
829 
830  // LOOP ON NEGATIVE TRACKS
831  for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
832 
833  //if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
834  //if(iTrkN1%1==0) printf(" 1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);
835 
836  if(iTrkN1==iTrkP1) continue;
837 
838  // get track from tracks array
839  negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
840 
841  if(negtrack1->Charge()>0 && !fLikeSign) continue;
842 
843  if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
844 
845  if(fMixEvent) {
846  if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
847  }
848 
849  if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
850  isLikeSign2Prong=kTRUE;
851  if(!fLikeSign) continue;
852  if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
853  } else { // unlike-sign
854  isLikeSign2Prong=kFALSE;
855  if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
856  if(fMixEvent) {
857  if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
858  }
859 
860  }
861 
862  // back to primary vertex
863  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
864  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
865  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
866  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
867  negtrack1->GetPxPyPz(momneg1);
868 
869  // DCA between the two tracks
870  dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
871  if(dcap1n1>dcaMax) { negtrack1=0; continue; }
872 
873  // Vertexing
874  twoTrackArray1->AddAt(postrack1,0);
875  twoTrackArray1->AddAt(negtrack1,1);
876  AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
877  if(!vertexp1n1) {
878  twoTrackArray1->Clear();
879  negtrack1=0;
880  continue;
881  }
882  // 2 prong candidate
883  if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
884 
885  io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
886 
887  if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
888  // add the vertex and the decay to the AOD
889  AliAODVertex *v2Prong =0x0;
890  if(!fMakeReducedRHF)v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
891  if(!isLikeSign2Prong) {
892  if(okD0) {
893  rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
895 
896  if(fMakeReducedRHF){
897  rd->DeleteRecoD();
898  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
899  }else{
900  rd->SetSecondaryVtx(v2Prong);
901  v2Prong->SetParent(rd);
902  AddRefs(v2Prong,rd,event,twoTrackArray1);
903  }
904  }
905  if(okJPSI) {
906  rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
907  if(fMakeReducedRHF){
908  rd->DeleteRecoD();
909  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
910  }else{
911  if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
912  AddRefs(v2Prong,rd,event,twoTrackArray1);
913  }
914  }
915  } else { // isLikeSign2Prong
916  rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
917  //Set selection bit for PID
919  if(fMakeReducedRHF){
920  rd->DeleteRecoD();
921  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
922  }else{
923  rd->SetSecondaryVtx(v2Prong);
924  v2Prong->SetParent(rd);
925  AddRefs(v2Prong,rd,event,twoTrackArray1);
926  }
927  }
928  }
929  // D* candidates
930  if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
931  // write references in io2Prong
932  if(fInputAOD) {
933  AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
934  } else {
935  vertexp1n1->AddDaughter(postrack1);
936  vertexp1n1->AddDaughter(negtrack1);
937  }
938  io2Prong->SetSecondaryVtx(vertexp1n1);
939  //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
940  // create a track from the D0
941  AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
942 
943  // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
944  for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
945 
946  if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
947 
948  if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
949 
950  if(fMixEvent) {
951  if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
952  evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
953  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
954  }
955 
956  //if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
957 
958  trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
959  if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
960 
961  // get track from tracks array
962  trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
963  // trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
964  SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
965  twoTrackArrayCasc->AddAt(trackPi,0);
966  twoTrackArrayCasc->AddAt(trackD0,1);
967  if(!SelectInvMassAndPtDstarD0pi(twoTrackArrayCasc)){
968  twoTrackArrayCasc->Clear();
969  trackPi=0;
970  continue;
971  }
972 
973  AliAODVertex *vertexCasc = 0;
974 
975  if(fFindVertexForDstar) {
976  // DCA between the two tracks
977  dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
978  // Vertexing
979  vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
980  } else {
981  // assume Dstar decays at the primary vertex
982  Double_t pos[3],cov[6],chi2perNDF;
983  fV1->GetXYZ(pos);
984  fV1->GetCovMatrix(cov);
985  chi2perNDF = fV1->GetChi2toNDF();
986  vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
987  dcaCasc = 0.;
988  }
989  if(!vertexCasc) {
990  twoTrackArrayCasc->Clear();
991  trackPi=0;
992  continue;
993  }
994 
995  ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
996  if(okDstar) {
997  // add the D0 to the AOD (if not already done)
998  if(!okD0) {
999  rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
1000  if(fMakeReducedRHF){
1001  rd->DeleteRecoD();
1002  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1003  }else{
1004  AliAODVertex *v2Prong = new (verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
1005  rd->SetSecondaryVtx(v2Prong);
1006  v2Prong->SetParent(rd);
1007  AddRefs(v2Prong,rd,event,twoTrackArray1);
1008  }
1009  okD0=kTRUE; // this is done to add it only once
1010  }
1011  // add the vertex and the cascade to the AOD
1012  rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
1013  // Set selection bit for PID
1015  AliAODVertex *vCasc = 0x0;
1016  if(fMakeReducedRHF){
1017  //assign a ID to the D0 candidate, daughter of the Cascade. ID = position in the D0toKpi array
1018  UShort_t idCasc[2]={(UShort_t)trackPi->GetID(),(UShort_t)(iD0toKpi-1)};
1019  rc->SetProngIDs(2,idCasc);
1020  rc->DeleteRecoD();
1021  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1022  }else{
1023  AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
1024  rc->SetSecondaryVtx(vCasc);
1025  vCasc->SetParent(rc);
1026  if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
1027  AddRefs(vCasc,rc,event,twoTrackArrayCasc);
1028  vCasc->AddDaughter(rd); // add the D0 (in ref #1)
1029  }
1030  }
1031  twoTrackArrayCasc->Clear();
1032  trackPi=0;
1033  if(ioCascade) {delete ioCascade; ioCascade=NULL;}
1034  delete vertexCasc; vertexCasc=NULL;
1035  } // end loop on soft pi tracks
1036 
1037  if(trackD0) {delete trackD0; trackD0=NULL;}
1038 
1039  }
1040  if(io2Prong) {delete io2Prong; io2Prong=NULL;}
1041  }
1042 
1043  twoTrackArray1->Clear();
1044  if( (!f3Prong && !f4Prong) ||
1045  (isLikeSign2Prong && !f3Prong) ) {
1046  negtrack1=0;
1047  delete vertexp1n1;
1048  continue;
1049  }
1050 
1051 
1052  // 2nd LOOP ON POSITIVE TRACKS
1053  for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
1054 
1055  if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
1056 
1057  //if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
1058 
1059  // get track from tracks array
1060  postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
1061 
1062  if(postrack2->Charge()<0) continue;
1063 
1064  if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
1065 
1066  // Check single tracks cuts specific for 3 prongs
1067  if(!TESTBIT(seleFlags[iTrkP2],kBit3Prong)) continue;
1068  if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
1069  if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
1070 
1071  if(fMixEvent) {
1072  if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1073  evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
1074  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1075  }
1076 
1077  if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1078  if(!fLikeSign3prong) continue;
1079  if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
1080  isLikeSign3Prong=kTRUE;
1081  } else { // not ok
1082  continue;
1083  }
1084  } else { // normal triplet (+-+)
1085  isLikeSign3Prong=kFALSE;
1086  if(fMixEvent) {
1087  if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1088  evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
1089  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1090  }
1091  }
1092 
1094  if(!TESTBIT(seleFlags[iTrkN1],kBitKaonCompat)) continue;
1095  }
1096  Bool_t okForLcTopKpi=kTRUE;
1097  Int_t pidLcStatus=3; // 3= OK as pKpi and Kpipi
1098  if(fUsePIDforLc>0){
1099  if(!TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) &&
1100  !TESTBIT(seleFlags[iTrkP2],kBitProtonCompat) ){
1101  okForLcTopKpi=kFALSE;
1102  pidLcStatus=0;
1103  }
1104  if(okForLcTopKpi && fUsePIDforLc>1){
1105  okForLcTopKpi=kFALSE;
1106  pidLcStatus=0;
1107  if(TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) &&
1108  TESTBIT(seleFlags[iTrkP2],kBitPionCompat) ){
1109  okForLcTopKpi=kTRUE;
1110  pidLcStatus+=1; // 1= OK as pKpi
1111  }
1112  if(TESTBIT(seleFlags[iTrkP2],kBitProtonCompat) &&
1113  TESTBIT(seleFlags[iTrkP1],kBitPionCompat) ){
1114  okForLcTopKpi=kTRUE;
1115  pidLcStatus+=2; // 2= OK as piKp
1116  }
1117  }
1118  }
1119  Bool_t okForDsToKKpi=kTRUE;
1120  if(fUseKaonPIDforDs){
1121  if(!TESTBIT(seleFlags[iTrkP1],kBitKaonCompat) &&
1122  !TESTBIT(seleFlags[iTrkP2],kBitKaonCompat) ) okForDsToKKpi=kFALSE;
1123  }
1124  // back to primary vertex
1125  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1126  // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1127  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1128  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1129  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1130  SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1131 
1132  //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
1133 
1134  dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1135  if(dcap2n1>dcaMax) { postrack2=0; continue; }
1136  dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1137  if(dcap1p2>dcaMax) { postrack2=0; continue; }
1138 
1139  // check invariant mass cuts for D+,Ds,Lc
1140  massCutOK=kTRUE;
1141  if(f3Prong) {
1142  if(postrack2->Charge()>0) {
1143  threeTrackArray->AddAt(postrack1,0);
1144  threeTrackArray->AddAt(negtrack1,1);
1145  threeTrackArray->AddAt(postrack2,2);
1146  } else {
1147  threeTrackArray->AddAt(negtrack1,0);
1148  threeTrackArray->AddAt(postrack1,1);
1149  threeTrackArray->AddAt(postrack2,2);
1150  }
1152  postrack2->GetPxPyPz(mompos2);
1153  Double_t pxDau[3]={mompos1[0],momneg1[0],mompos2[0]};
1154  Double_t pyDau[3]={mompos1[1],momneg1[1],mompos2[1]};
1155  Double_t pzDau[3]={mompos1[2],momneg1[2],mompos2[2]};
1156  // massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1157  massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau,pidLcStatus);
1158  }
1159  }
1160 
1161  if(f3Prong && !massCutOK) {
1162  threeTrackArray->Clear();
1163  if(!f4Prong) {
1164  postrack2=0;
1165  continue;
1166  }
1167  }
1168 
1169  // Vertexing
1170  twoTrackArray2->AddAt(postrack2,0);
1171  twoTrackArray2->AddAt(negtrack1,1);
1172  AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1173  if(!vertexp2n1) {
1174  twoTrackArray2->Clear();
1175  postrack2=0;
1176  continue;
1177  }
1178 
1179  // 3 prong candidates
1180  if(f3Prong && massCutOK) {
1181 
1182  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1183  io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,okForLcTopKpi,okForDsToKKpi,ok3Prong);
1184  if(ok3Prong) {
1185  AliAODVertex *v3Prong=0x0;
1186  if(!fMakeReducedRHF)v3Prong = new (verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1187  if(!isLikeSign3Prong) {
1188  rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1189  // Set selection bit for PID
1193  if(fMakeReducedRHF){
1194  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1195  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1196  }else{
1197  v3Prong = new (verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1198  rd->SetSecondaryVtx(v3Prong);
1199  v3Prong->SetParent(rd);
1200  AddRefs(v3Prong,rd,event,threeTrackArray);
1201  }
1202  } else { // isLikeSign3Prong
1203  if(fLikeSign3prong){
1204  rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1205  // Set selection bit for PID
1209  if(fMakeReducedRHF){
1210  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1211  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1212  }else{
1213  rd->SetSecondaryVtx(v3Prong);
1214  v3Prong->SetParent(rd);
1215  AddRefs(v3Prong,rd,event,threeTrackArray);
1216  }
1217  }
1218  }
1219 
1220  }
1221  if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1222  if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1223  }
1224 
1225  // 4 prong candidates
1226  if(f4Prong
1227  // don't make 4 prong with like-sign pairs and triplets
1228  && !isLikeSign2Prong && !isLikeSign3Prong
1229  // track-to-track dca cuts already now
1230  && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
1231  && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
1232  // back to primary vertex
1233  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1234  // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1235  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1236  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1237  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1238  SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1239 
1240  // Vertexing for these 3 (can be taken from above?)
1241  threeTrackArray->AddAt(postrack1,0);
1242  threeTrackArray->AddAt(negtrack1,1);
1243  threeTrackArray->AddAt(postrack2,2);
1244  AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1245 
1246  // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
1247  for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1248 
1249  if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1250 
1251  //if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1252 
1253  // get track from tracks array
1254  negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1255 
1256  if(negtrack2->Charge()>0) continue;
1257 
1258  if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1259  if(fMixEvent){
1260  if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1261  evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1262  evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
1263  evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
1264  evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1265  evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
1266  }
1267 
1268  // back to primary vertex
1269  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1270  // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1271  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1272  // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1273  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1274  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1275  SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1276  SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1277 
1278  dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1279  if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1280  dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1281  if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1282 
1283 
1284  fourTrackArray->AddAt(postrack1,0);
1285  fourTrackArray->AddAt(negtrack1,1);
1286  fourTrackArray->AddAt(postrack2,2);
1287  fourTrackArray->AddAt(negtrack2,3);
1288 
1289  // check invariant mass cuts for D0
1290  massCutOK=kTRUE;
1292  massCutOK = SelectInvMassAndPt4prong(fourTrackArray);
1293 
1294  if(!massCutOK) {
1295  fourTrackArray->Clear();
1296  negtrack2=0;
1297  continue;
1298  }
1299 
1300  // Vertexing
1301  AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
1302  io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
1303  if(ok4Prong) {
1304  rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
1305  if(fMakeReducedRHF){
1306  rd->DeleteRecoD();
1307  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1308  }else{
1309  AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
1310  rd->SetSecondaryVtx(v4Prong);
1311  v4Prong->SetParent(rd);
1312  AddRefs(v4Prong,rd,event,fourTrackArray);
1313  }
1314  }
1315 
1316  if(io4Prong) {delete io4Prong; io4Prong=NULL;}
1317  if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
1318  fourTrackArray->Clear();
1319  negtrack2 = 0;
1320 
1321  } // end loop on negative tracks
1322 
1323  threeTrackArray->Clear();
1324  delete vertexp1n1p2;
1325 
1326  }
1327 
1328  postrack2 = 0;
1329  delete vertexp2n1;
1330 
1331  } // end 2nd loop on positive tracks
1332 
1333  twoTrackArray2->Clear();
1334 
1335  // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
1336  for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1337 
1338  if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1339 
1340  //if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1341 
1342  // get track from tracks array
1343  negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1344 
1345  if(negtrack2->Charge()>0) continue;
1346 
1347  if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1348 
1349  // Check single tracks cuts specific for 3 prongs
1350  if(!TESTBIT(seleFlags[iTrkN2],kBit3Prong)) continue;
1351  if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
1352  if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
1353 
1354  if(fMixEvent) {
1355  if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1356  evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1357  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1358  }
1359 
1360  if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1361  if(!fLikeSign3prong) continue;
1362  if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1363  isLikeSign3Prong=kTRUE;
1364  } else { // not ok
1365  continue;
1366  }
1367  } else { // normal triplet (-+-)
1368  isLikeSign3Prong=kFALSE;
1369  }
1370 
1372  if(!TESTBIT(seleFlags[iTrkP1],kBitKaonCompat)) continue;
1373  }
1374  Bool_t okForLcTopKpi=kTRUE;
1375  Int_t pidLcStatus=3; // 3= OK as pKpi and Kpipi
1376  if(fUsePIDforLc>0){
1377  if(!TESTBIT(seleFlags[iTrkN1],kBitProtonCompat) &&
1378  !TESTBIT(seleFlags[iTrkN2],kBitProtonCompat) ){
1379  okForLcTopKpi=kFALSE;
1380  pidLcStatus=0;
1381  }
1382  if(okForLcTopKpi && fUsePIDforLc>1){
1383  okForLcTopKpi=kFALSE;
1384  pidLcStatus=0;
1385  if(TESTBIT(seleFlags[iTrkN1],kBitProtonCompat) &&
1386  TESTBIT(seleFlags[iTrkN2],kBitPionCompat) ){
1387  okForLcTopKpi=kTRUE;
1388  pidLcStatus+=1; // 1= OK as pKpi
1389  }
1390  if(TESTBIT(seleFlags[iTrkN2],kBitProtonCompat) &&
1391  TESTBIT(seleFlags[iTrkN1],kBitPionCompat) ){
1392  okForLcTopKpi=kTRUE;
1393  pidLcStatus+=2; // 2= OK as piKp
1394  }
1395  }
1396  }
1397  Bool_t okForDsToKKpi=kTRUE;
1398  if(fUseKaonPIDforDs){
1399  if(!TESTBIT(seleFlags[iTrkN1],kBitKaonCompat) &&
1400  !TESTBIT(seleFlags[iTrkN2],kBitKaonCompat) ) okForDsToKKpi=kFALSE;
1401  }
1402 
1403  // back to primary vertex
1404  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1405  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1406  // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1407  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1408  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1409  SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1410  //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1411 
1412  dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1413  if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1414  dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1415  if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1416 
1417  threeTrackArray->AddAt(negtrack1,0);
1418  threeTrackArray->AddAt(postrack1,1);
1419  threeTrackArray->AddAt(negtrack2,2);
1420 
1421  // check invariant mass cuts for D+,Ds,Lc
1422  massCutOK=kTRUE;
1424  negtrack2->GetPxPyPz(momneg2);
1425  Double_t pxDau[3]={momneg1[0],mompos1[0],momneg2[0]};
1426  Double_t pyDau[3]={momneg1[1],mompos1[1],momneg2[1]};
1427  Double_t pzDau[3]={momneg1[2],mompos1[2],momneg2[2]};
1428  // massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1429  massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau,pidLcStatus);
1430  }
1431  if(!massCutOK) {
1432  threeTrackArray->Clear();
1433  negtrack2=0;
1434  continue;
1435  }
1436 
1437  // Vertexing
1438  twoTrackArray2->AddAt(postrack1,0);
1439  twoTrackArray2->AddAt(negtrack2,1);
1440 
1441  AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1442  if(!vertexp1n2) {
1443  twoTrackArray2->Clear();
1444  negtrack2=0;
1445  continue;
1446  }
1447 
1448  if(f3Prong) {
1449  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1450  io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,okForLcTopKpi,okForDsToKKpi,ok3Prong);
1451  if(ok3Prong) {
1452  AliAODVertex *v3Prong = 0x0;
1453  if(!fMakeReducedRHF) v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1454  if(!isLikeSign3Prong) {
1455  rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1459  if(fMakeReducedRHF){
1460  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1461  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1462  }else{
1463  rd->SetSecondaryVtx(v3Prong);
1464  v3Prong->SetParent(rd);
1465  AddRefs(v3Prong,rd,event,threeTrackArray);
1466  }
1467  } else { // isLikeSign3Prong
1468  if(fLikeSign3prong){
1469  rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1473  if(fMakeReducedRHF){
1474  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1475  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1476  }else{
1477  rd->SetSecondaryVtx(v3Prong);
1478  v3Prong->SetParent(rd);
1479  AddRefs(v3Prong,rd,event,threeTrackArray);
1480  }
1481  }
1482 
1483  }
1484  }
1485  if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1486  if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1487  }
1488  threeTrackArray->Clear();
1489  negtrack2 = 0;
1490  delete vertexp1n2;
1491 
1492  } // end 2nd loop on negative tracks
1493 
1494  twoTrackArray2->Clear();
1495 
1496  negtrack1 = 0;
1497  delete vertexp1n1;
1498  } // end 1st loop on negative tracks
1499 
1500  postrack1 = 0;
1501  } // end 1st loop on positive tracks
1502 
1503 
1504  // AliDebug(1,Form(" Total HF vertices in event = %d;",
1505  // (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1506  if(fD0toKpi) {
1507  AliDebug(1,Form(" D0->Kpi in event = %d;",
1508  (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1509  }
1510  if(fJPSItoEle) {
1511  AliDebug(1,Form(" JPSI->ee in event = %d;",
1512  (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1513  }
1514  if(f3Prong) {
1515  AliDebug(1,Form(" Charm->3Prong in event = %d;",
1516  (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1517  }
1518  if(f4Prong) {
1519  AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1520  (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1521  }
1522  if(fDstar) {
1523  AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1524  (Int_t)aodDstarTClArr->GetEntriesFast()));
1525  }
1526  if(fCascades){
1527  AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1528  (Int_t)aodCascadesTClArr->GetEntriesFast()));
1529  }
1530  if(fLikeSign) {
1531  AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1532  (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1533  }
1534  if(fLikeSign3prong && f3Prong) {
1535  AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1536  (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1537  }
1538 
1539 
1540  twoTrackArray1->Delete(); delete twoTrackArray1;
1541  twoTrackArray2->Delete(); delete twoTrackArray2;
1542  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1543  twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1544  threeTrackArray->Clear();
1545  threeTrackArray->Delete(); delete threeTrackArray;
1546  fourTrackArray->Delete(); delete fourTrackArray;
1547  delete [] seleFlags; seleFlags=NULL;
1548  if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1549  tracksAtVertex.Delete();
1550 
1551  if(fInputAOD) {
1552  seleTrksArray.Delete();
1553  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
1554  }
1555 
1556 
1557  //printf("Trks: total %d sele %d\n",fnTrksTotal,fnSeleTrksTotal);
1558 
1559  return;
1560 }
1561 //----------------------------------------------------------------------------
1563  const AliVEvent *event,
1564  const TObjArray *trkArray) const
1565 {
1568  //AliCodeTimerAuto("",0);
1569 
1570  if(fInputAOD) {
1571  AddDaughterRefs(v,event,trkArray);
1572  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1573  }
1574 
1575  /*
1576  rd->SetListOfCutsRef((TList*)fListOfCuts);
1577  //fListOfCuts->Print();
1578  cout<<fListOfCuts<<endl;
1579  TList *l=(TList*)rd->GetListOfCuts();
1580  cout<<l<<endl;
1581  if(l) {l->Print(); }else{printf("error\n");}
1582  */
1583 
1584  return;
1585 }
1586 //---------------------------------------------------------------------------
1588  const AliVEvent *event,
1589  const TObjArray *trkArray) const
1590 {
1592  //AliCodeTimerAuto("",0);
1593 
1594  Int_t nDg = v->GetNDaughters();
1595  TObject *dg = 0;
1596  if(nDg) dg = v->GetDaughter(0);
1597 
1598  if(dg) return; // daughters already added
1599 
1600  Int_t nTrks = trkArray->GetEntriesFast();
1601 
1602  AliExternalTrackParam *track = 0;
1603  AliAODTrack *aodTrack = 0;
1604  Int_t id;
1605 
1606  for(Int_t i=0; i<nTrks; i++) {
1607  track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1608  id = (Int_t)track->GetID();
1609  //printf("---> %d\n",id);
1610  if(id<0) continue; // this track is a AliAODRecoDecay
1611  aodTrack = dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[id]));
1612  if(!aodTrack) AliFatal("Not a standard AOD");
1613  v->AddDaughter(aodTrack);
1614  }
1615 
1616  return;
1617 }
1618 //----------------------------------------------------------------------------
1620 {
1623  //
1624  //AliCodeTimerAuto("",0);
1625 
1626 
1627  TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1628  if(!inputArray) return;
1629 
1630  AliAODTrack *track = 0;
1631  AliAODVertex *vertex = 0;
1632 
1633  Bool_t needtofix=kFALSE;
1634  for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1635  vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1636  for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1637  track = (AliAODTrack*)vertex->GetDaughter(id);
1638  if(!track->GetStatus()) needtofix=kTRUE;
1639  }
1640  if(needtofix) break;
1641  }
1642 
1643  if(!needtofix) return;
1644 
1645 
1646  printf("Fixing references\n");
1647 
1648  fAODMapSize = 100000;
1649  fAODMap = new Int_t[fAODMapSize];
1650  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
1651 
1652  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1653  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
1654  if(!track) AliFatal("Not a standard AOD");
1655 
1656  // skip pure ITS SA tracks
1657  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1658 
1659  // skip tracks without ITS
1660  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1661 
1662  // TEMPORARY: check that the cov matrix is there
1663  Double_t covtest[21];
1664  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1665  //
1666 
1667  Int_t ind = (Int_t)track->GetID();
1668  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
1669  }
1670 
1671 
1672  Int_t ids[4]={-1,-1,-1,-1};
1673  for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1674  Bool_t cascade=kFALSE;
1675  vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1676  Int_t id=0;
1677  Int_t nDgs = vertex->GetNDaughters();
1678  for(id=0; id<nDgs; id++) {
1679  track = (AliAODTrack*)vertex->GetDaughter(id);
1680  if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1681  ids[id]=(Int_t)track->GetID();
1682  vertex->RemoveDaughter(track);
1683  }
1684  if(cascade) continue;
1685  for(id=0; id<nDgs; id++) {
1686  if (ids[id]>-1 && ids[id] < fAODMapSize) {
1687  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(fAODMap[ids[id]]));
1688  if(!track) AliFatal("Not a standard AOD");
1689  vertex->AddDaughter(track);
1690  }
1691  }
1692 
1693  }
1694 
1695  return;
1696 }
1697 //----------------------------------------------------------------------------
1699  // method to retrieve daughters from trackID and reconstruct secondary vertex
1700  // save the TRefs to the candidate AliAODRecoDecayHF3Prong rd
1701  // and fill on-the-fly the data member of rd
1702  if(rd->GetIsFilled()!=0)return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2 already refilled)
1703  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1704  TObjArray *threeTrackArray = new TObjArray(3);
1705 
1706  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1707  if(!track1)return kFALSE;
1708  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1709  if(!track2)return kFALSE;
1710  AliESDtrack *postrack1 = 0;
1711  AliESDtrack *negtrack1 = 0;
1712  postrack1 = new AliESDtrack(track1);
1713  negtrack1 = new AliESDtrack(track2);
1714 
1715  // DCA between the two tracks
1716  Double_t xdummy, ydummy;
1717  fBzkG = (Double_t)event->GetMagneticField();
1718  Double_t dca12 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1719 
1720  const AliVVertex *vprimary = event->GetPrimaryVertex();
1721  Double_t pos[3];
1722  Double_t cov[6];
1723  vprimary->GetXYZ(pos);
1724  vprimary->GetCovarianceMatrix(cov);
1725  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1726  fV1->GetCovMatrix(cov);
1727  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1728 
1729  AliAODTrack *track3 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(2)]);
1730  if(!track3)return kFALSE;
1731  AliESDtrack *esdt3 = new AliESDtrack(track3);
1732 
1733  Double_t dca2;
1734  Double_t dca3;
1735  threeTrackArray->AddAt(postrack1,0);
1736  threeTrackArray->AddAt(negtrack1,1);
1737  threeTrackArray->AddAt(esdt3,2);
1738  dca2 = esdt3->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1739  dca3 = esdt3->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1740  Double_t dispersion;
1741 
1742  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray, dispersion);
1743  if (!secVert3PrAOD) {
1744  threeTrackArray->Clear();
1745  threeTrackArray->Delete(); delete threeTrackArray;
1746  delete fV1; fV1=0;
1747  delete postrack1; postrack1=NULL;
1748  delete negtrack1; negtrack1=NULL;
1749  delete esdt3; esdt3=NULL;
1750  return kFALSE;
1751  }
1752 
1753  rd->SetNProngs();
1754  Double_t vtxRec=rd->GetDist12toPrim();
1755  Double_t vertexp2n1=rd->GetDist23toPrim();
1756  rd= Make3Prong(threeTrackArray, event, secVert3PrAOD,dispersion, vtxRec, vertexp2n1, dca12, dca2, dca3, rd);
1757  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1758  rd->SetIsFilled(2);
1759  threeTrackArray->Clear();
1760  threeTrackArray->Delete(); delete threeTrackArray;
1761  delete fV1; fV1=0;
1762  delete postrack1; postrack1=NULL;
1763  delete negtrack1; negtrack1=NULL;
1764  delete esdt3; esdt3=NULL;
1765  return kTRUE;
1766 }
1767 //___________________________
1769  // method to retrieve daughters from trackID and reconstruct secondary vertex
1770  // save the TRefs to the candidate AliAODRecoDecayHF2Prong rd
1771  // and fill on-the-fly the data member of rd
1772  if(rd->GetIsFilled()!=0)return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1:standard dAOD, 2 already refilled)
1773  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1774 
1775  Double_t dispersion;
1776  TObjArray *twoTrackArray1 = new TObjArray(2);
1777 
1778  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1779  if(!track1)return kFALSE;
1780  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1781  if(!track2)return kFALSE;
1782 
1783  AliESDtrack *esdt1 = 0;
1784  AliESDtrack *esdt2 = 0;
1785  esdt1 = new AliESDtrack(track1);
1786  esdt2 = new AliESDtrack(track2);
1787 
1788  twoTrackArray1->AddAt(esdt1,0);
1789  twoTrackArray1->AddAt(esdt2,1);
1790  // DCA between the two tracks
1791  Double_t xdummy, ydummy;
1792  fBzkG = (Double_t)event->GetMagneticField();
1793  Double_t dca12 = esdt1->GetDCA(esdt2,fBzkG,xdummy,ydummy);
1794  const AliVVertex *vprimary = event->GetPrimaryVertex();
1795  Double_t pos[3];
1796  Double_t cov[6];
1797  vprimary->GetXYZ(pos);
1798  vprimary->GetCovarianceMatrix(cov);
1799  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1800  fV1->GetCovMatrix(cov);
1801  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1802 
1803 
1804  AliAODVertex *vtxRec = ReconstructSecondaryVertex(twoTrackArray1, dispersion);
1805  if(!vtxRec) {
1806  twoTrackArray1->Clear();
1807  twoTrackArray1->Delete(); delete twoTrackArray1;
1808  delete fV1; fV1=0;
1809  delete esdt1; esdt1=NULL;
1810  delete esdt2; esdt2=NULL;
1811  return kFALSE; }
1812  Bool_t okD0=kFALSE;
1813  Bool_t okJPSI=kFALSE;
1814  Bool_t okD0FromDstar=kFALSE;
1815  Bool_t refill =kTRUE;
1816  rd->SetNProngs();
1817  rd= Make2Prong(twoTrackArray1, event, vtxRec, dca12, okD0, okJPSI, okD0FromDstar,refill,rd);
1818  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1819  rd->SetIsFilled(2);
1820  delete fV1; fV1=0;
1821  twoTrackArray1->Clear();
1822  twoTrackArray1->Delete(); delete twoTrackArray1;
1823  delete esdt1; esdt1=NULL;
1824  delete esdt2; esdt2=NULL;
1825  return kTRUE;
1826 }
1827 //----------------------------------------------------------------------------
1829  // method to retrieve daughters from trackID
1830  // and fill on-the-fly the data member of rCasc and their AliAODRecoDecayHF2Prong daughters
1831  if(rCasc->GetIsFilled()!=0) return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2: already refilled)
1832  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1833  TObjArray *twoTrackArrayCasc = new TObjArray(2);
1834 
1835  AliAODTrack *trackB =(AliAODTrack*)event->GetTrack(fAODMap[rCasc->GetProngID(0)]);//retrieve bachelor from the trackID through the AOD index map
1836  if(!trackB)return kFALSE;
1837 
1838  AliNeutralTrackParam *trackV0=NULL;
1839  AliAODv0 *v0 =NULL;
1840  AliAODRecoDecayHF2Prong *trackD0=NULL;
1841  rCasc->SetNProngs();
1842 
1843  if(DStar){
1844  TClonesArray *inputArrayD0=(TClonesArray*)event->GetList()->FindObject("D0toKpi");
1845  if(!inputArrayD0) return kFALSE;
1846  trackD0=(AliAODRecoDecayHF2Prong*)inputArrayD0->At(rCasc->GetProngID(1));
1847  if(!trackD0)return kFALSE;
1848  if(!FillRecoCand(event,trackD0)) return kFALSE; //fill missing info of the corresponding D0 daughter
1849 
1850  trackV0 = new AliNeutralTrackParam(trackD0);
1851 
1852  }else{//is a V0 candidate
1853  v0 = ((AliAODEvent*)event)->GetV0(rCasc->GetProngID(1));
1854  if(!v0) return kFALSE;
1855  // Define the V0 (neutral) track
1856  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1857  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
1858  }
1859 
1860  AliESDtrack *esdB = new AliESDtrack(trackB);
1861 
1862  twoTrackArrayCasc->AddAt(esdB,0);
1863  twoTrackArrayCasc->AddAt(trackV0,1);
1864 
1865  fBzkG = (Double_t)event->GetMagneticField();
1866  const AliVVertex *vprimary = event->GetPrimaryVertex();
1867 
1868  Double_t pos[3];
1869  Double_t cov[6];
1870  vprimary->GetXYZ(pos);
1871  vprimary->GetCovarianceMatrix(cov);
1872  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1873  fV1->GetCovMatrix(cov);
1874 
1875  Double_t dca = 0.;
1876  AliAODVertex *vtxCasc = 0x0;
1877  Double_t chi2perNDF = fV1->GetChi2toNDF();
1878  if (recoSecVtx) {
1879  Double_t dispersion, xdummy, ydummy;
1880  dca = esdB->GetDCA(trackV0,fBzkG,xdummy,ydummy);
1881  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
1882  vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
1883  } else {
1884  vtxCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
1885  }
1886  if(!vtxCasc) {
1887  twoTrackArrayCasc->Clear();
1888  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1889  delete fV1; fV1=0;
1890  delete esdB; esdB=NULL;
1891  delete vtxCasc;vtxCasc=NULL;
1892  trackB=NULL;
1893  delete trackV0; trackV0=NULL;
1894  if(!DStar){
1895  v0=NULL;
1896  }
1897  return kFALSE;
1898  }
1899  vtxCasc->SetParent(rCasc);
1900  rCasc->SetSecondaryVtx(vtxCasc);
1901  AddDaughterRefs(vtxCasc,(AliAODEvent*)event,twoTrackArrayCasc);
1902  if(DStar)vtxCasc->AddDaughter(trackD0);
1903  else vtxCasc->AddDaughter(v0);
1904  rCasc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1905 
1906  Bool_t refill =kTRUE;
1907 
1908  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1909  // propagate tracks to secondary vertex, to compute inv. mass
1910  esdB->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1911  trackV0->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1912  Double_t momentum[3];
1913  esdB->GetPxPyPz(momentum);
1914  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1915  trackV0->GetPxPyPz(momentum);
1916  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1917 
1918  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc,event);
1919  if(!primVertexAOD){
1920  delete fV1; fV1=0;
1921  delete vtxCasc; vtxCasc=NULL;
1922  twoTrackArrayCasc->Clear();
1923  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1924  delete esdB; esdB=NULL;
1925  delete trackV0; trackV0=NULL;
1926  if(!DStar)v0=NULL;
1927  return kFALSE;
1928  }
1929  Double_t d0z0[2],covd0z0[3];
1930  esdB->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1931  d0[0] = d0z0[0];
1932  d0err[0] = TMath::Sqrt(covd0z0[0]);
1933  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1934  d0[1] = d0z0[0];
1935  d0err[1] = TMath::Sqrt(covd0z0[0]);
1936  rCasc->SetPxPyPzProngs(2,px,py,pz);
1937  rCasc->SetDCA(dca);
1938  rCasc->Setd0Prongs(2,d0);
1939  rCasc->Setd0errProngs(2,d0err);
1940  rCasc->SetOwnPrimaryVtx(primVertexAOD);
1941  rCasc->SetCharge(esdB->Charge());
1942  // get PID info from ESD
1943  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1944  if(esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
1945  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1946  Double_t esdpid[10];
1947  for(Int_t i=0;i<5;i++) {
1948  esdpid[i] = esdpid0[i];
1949  esdpid[5+i] = esdpid1[i];
1950  }
1951  rCasc->SetPID(2,esdpid);
1952  rCasc->SetIsFilled(2);
1953 
1954 
1955  delete fV1; fV1=0;
1956  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1957  twoTrackArrayCasc->Clear();
1958  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1959  delete esdB; esdB=NULL;
1960  delete trackV0; trackV0=NULL;
1961  if(!DStar)v0=NULL;
1962  return kTRUE;
1963 }
1964 //---------------------------------------------------------------------------
1966 {
1971 
1972 
1973  if (rc->GetIsFilled()==0) return kFALSE; // Reduced dAOD with candidates not yet refilled
1974  if (!fAODMap) MapAODtracks(event); // Fill the AOD index map if it is not yet done
1975 
1976 
1977  // - Get bachelor and V0 from the cascade
1978  AliAODTrack *trackB = dynamic_cast<AliAODTrack*> (rc->GetBachelor());
1979  if (!trackB) return kFALSE;
1980 
1981  AliAODv0 *v0 = dynamic_cast<AliAODv0*> (rc->Getv0());
1982  if (!v0) return kFALSE;
1983 
1984 
1985  // - Cast bachelor (AOD->ESD) and V0 (AliAODv0->AliNeutralTrackParam)
1986  AliESDtrack *esdB = new AliESDtrack(trackB);
1987  if (!esdB) return kFALSE;
1988 
1989  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1990  if (!trackVV0) return kFALSE;
1991 
1992  AliNeutralTrackParam *trackV0 = new AliNeutralTrackParam(trackVV0);
1993  if (!trackV0) return kFALSE;
1994 
1995  TObjArray *twoTrackArrayCasc = new TObjArray(2);
1996  twoTrackArrayCasc->AddAt(esdB, 0);
1997  twoTrackArrayCasc->AddAt(trackV0, 1);
1998 
1999 
2000  Double_t dispersion, xdummy, ydummy, pos[3], cov[6];
2001 
2002 
2003  // - Some ingredients will be needed:
2004  // magnetic field
2005  fBzkG = (Double_t) event->GetMagneticField();
2006 
2007  // primary vertex
2008  const AliVVertex *vprimary = event->GetPrimaryVertex();
2009  vprimary->GetXYZ(pos);
2010  vprimary->GetCovarianceMatrix(cov);
2011  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2012 
2013  // vertexer
2014  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
2015 
2016 
2017  // - Compute the DCA and the decay vertex
2018  Double_t dca = esdB->GetDCA(trackV0, fBzkG, xdummy, ydummy);
2019  AliAODVertex *vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc, dispersion, kFALSE);
2020 
2021  if (!vtxCasc) {
2022  twoTrackArrayCasc->Clear();
2023  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2024  delete fV1; fV1=0;
2025  delete esdB; esdB=0;
2026  delete vtxCasc; vtxCasc=0;
2027  delete trackV0; trackV0=0;
2028  trackB=0; v0=0;
2029  return kFALSE;
2030  }
2031 
2032 
2033  // - Linking the cascade with the new secondary vertex
2034  vtxCasc->SetParent(rc);
2035  rc->SetSecondaryVtx(vtxCasc);
2036  AddDaughterRefs(vtxCasc, (AliAODEvent*)event, twoTrackArrayCasc);
2037  vtxCasc->AddDaughter(v0);
2038  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
2039 
2040 
2041  // - Propagate tracks to secondary vertex, to get momenta
2042  Double_t momentum[3], px[2], py[2], pz[2], d0[2], d0err[2], d0z0[2], covd0z0[3];
2043  esdB->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2044  trackV0->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2045  esdB->GetPxPyPz(momentum);
2046  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2047  trackV0->GetPxPyPz(momentum);
2048  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2049 
2050 
2051  // - Propagate tracks to primary vertex, to get impact parameters
2052  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc, event);
2053  if (!primVertexAOD) {
2054  twoTrackArrayCasc->Clear();
2055  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2056  delete fV1; fV1=0;
2057  delete esdB; esdB=0;
2058  delete vtxCasc; vtxCasc=0;
2059  delete trackV0; trackV0=0;
2060  trackB=0; v0=0;
2061  return kFALSE;
2062  }
2063 
2064  esdB->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2065  d0[0] = d0z0[0];
2066  d0err[0] = TMath::Sqrt(covd0z0[0]);
2067  trackV0->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2068  d0[1] = d0z0[0];
2069  d0err[1] = TMath::Sqrt(covd0z0[0]);
2070 
2071 
2072  // - Get PID info from ESD
2073  Double_t esdpid[10];
2074  Double_t esdpid0[5] = {0., 0., 0., 0., 0.};
2075  Double_t esdpid1[5] = {0., 0., 0., 0., 0.};
2076  if (esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
2077  for (Int_t ipid=0; ipid<5; ipid++) {
2078  esdpid[ipid] = esdpid0[ipid];
2079  esdpid[5+ipid] = esdpid1[ipid];
2080  }
2081 
2082 
2083  // - Set data members
2084  rc->SetOwnPrimaryVtx(primVertexAOD);
2085  rc->SetDCA(dca);
2086  rc->SetPxPyPzProngs(2, px, py, pz);
2087  rc->Setd0Prongs(2, d0);
2088  rc->Setd0errProngs(2, d0err);
2089  rc->SetCharge(esdB->Charge());
2090  rc->SetPID(2, esdpid);
2091 
2092  rc->SetIsFilled(2); // To clean the newly reconstructed secondary vertex with the CleanUp task
2093 
2094 
2095  twoTrackArrayCasc->Clear();
2096  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2097  delete fV1; fV1=0;
2098  delete primVertexAOD; primVertexAOD=0;
2099  delete esdB; esdB=0;
2100  delete trackV0; trackV0=0;
2101  trackB=0; v0=0;
2102 
2103  return kTRUE;
2104 }
2105 
2106 //---------------------------------------------------------------------------
2108  TObjArray *twoTrackArray,AliVEvent *event,
2109  AliAODVertex *secVert,
2110  AliAODRecoDecayHF2Prong *rd2Prong,
2111  Double_t dca,
2112  Bool_t &okDstar)
2113 {
2116  //AliCodeTimerAuto("",0);
2117  UInt_t ntref=TProcessID::GetObjectCount();
2118  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2119  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2120  }
2121 
2122  okDstar = kFALSE;
2123 
2124  Bool_t dummy1,dummy2,dummy3;
2125 
2126  // We use Make2Prong to construct the AliAODRecoCascadeHF
2127  // (which inherits from AliAODRecoDecayHF2Prong)
2128  AliAODRecoCascadeHF *theCascade =
2129  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2130  dummy1,dummy2,dummy3);
2131  if(!theCascade) return 0x0;
2132 
2133  // charge
2134  AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2135  theCascade->SetCharge(trackPi->Charge());
2136 
2137  //--- selection cuts
2138  //
2139  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2140  if(fInputAOD){
2141  Int_t idSoftPi=(Int_t)trackPi->GetID();
2142  if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
2143  AliAODTrack* trackPiAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idSoftPi]));
2144  if(!trackPiAOD) AliFatal("Not a standard AOD");
2145  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
2146  }
2147  }else{
2148  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
2149  }
2150  tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
2151 
2152  AliAODVertex *primVertexAOD=0;
2154  // take event primary vertex
2155  primVertexAOD = PrimaryVertex();
2156  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2157  rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
2158  }
2159  // select D*->D0pi
2160  if(fDstar) {
2162  if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
2163  }
2164  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2165  tmpCascade->UnsetOwnPrimaryVtx();
2166  delete tmpCascade; tmpCascade=NULL;
2168  rd2Prong->UnsetOwnPrimaryVtx();
2169  }
2170  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2171  //---
2172 
2173 
2174  return theCascade;
2175 }
2176 
2177 
2178 //----------------------------------------------------------------------------
2180  TObjArray *twoTrackArray,AliVEvent *event,
2181  AliAODVertex *secVert,
2182  AliAODv0 *v0,
2183  Double_t dca,
2184  Bool_t &okCascades)
2185 {
2188  //AliCodeTimerAuto("",0);
2189  UInt_t ntref=TProcessID::GetObjectCount();
2190  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2191  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2192  }
2193 
2194  // preselection to reduce the number of TRefs
2195  Double_t px[2],py[2],pz[2];
2196  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2197  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2198  // propagate tracks to secondary vertex, to compute inv. mass
2199  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2200  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2201  Double_t momentum[3];
2202  postrack->GetPxPyPz(momentum);
2203  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2204  negtrack->GetPxPyPz(momentum);
2205  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2206  if(!SelectInvMassAndPtCascade(px,py,pz)) return 0x0;
2207  Double_t dummyd0[2]={0,0};
2208  Double_t dummyd0err[2]={0,0};
2209  AliAODRecoCascadeHF tmpCasc(0x0,postrack->Charge(),px, py, pz, dummyd0, dummyd0err,0.);
2210  Bool_t presel=fCutsLctoV0->PreSelect(&tmpCasc,v0,postrack);
2211  if(!presel) return 0x0;
2212 
2213  // AliDebug(2,Form(" building the cascade"));
2214  okCascades = kFALSE;
2215  Bool_t dummy1,dummy2,dummy3;
2216 
2217  // We use Make2Prong to construct the AliAODRecoCascadeHF
2218  // (which inherits from AliAODRecoDecayHF2Prong)
2219  AliAODRecoCascadeHF *theCascade =
2220  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2221  dummy1,dummy2,dummy3);
2222  if(!theCascade) return 0x0;
2223 
2224  // bachelor track and charge
2225  AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2226  theCascade->SetCharge(trackBachelor->Charge());
2227 
2228  //--- selection cuts
2229  //
2230 
2231  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2232  if(fInputAOD){
2233  Int_t idBachelor=(Int_t)trackBachelor->GetID();
2234  if (idBachelor > -1 && idBachelor < fAODMapSize) {
2235  AliAODTrack* trackBachelorAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idBachelor]));
2236  if(!trackBachelorAOD) AliFatal("Not a standard AOD");
2237  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
2238  }
2239  }else{
2240  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
2241  }
2242  tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
2243 
2244  AliAODVertex *primVertexAOD=0;
2246  // take event primary vertex
2247  primVertexAOD = PrimaryVertex();
2248  if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
2249  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2250  }
2251 
2252  // select Cascades
2253  if (fCascades && fInputAOD) {
2254  if (fCutsLctoV0->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2255  okCascades = kTRUE;
2257  }
2258  if (fCutsDplustoK0spi) {
2260  okCascades = kTRUE;
2262  }
2263  }
2264  if (fCutsDstoK0sK){
2265  if (fCutsDstoK0sK->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2266  okCascades = kTRUE;
2268  }
2269  }
2270  }
2271  else {
2272  //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
2273  okCascades=kTRUE;
2274  } // no cuts implemented from ESDs
2275  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2276  tmpCascade->UnsetOwnPrimaryVtx();
2277  delete tmpCascade; tmpCascade=NULL;
2278  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2279  //---
2280 
2281  return theCascade;
2282 }
2283 
2284 //-----------------------------------------------------------------------------
2286  TObjArray *twoTrackArray,AliVEvent *event,
2287  AliAODVertex *secVert,Double_t dca,
2288  Bool_t &okD0,Bool_t &okJPSI,
2289  Bool_t &okD0fromDstar, Bool_t refill, AliAODRecoDecayHF2Prong *rd)
2290 {
2293  // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
2294  // AliCodeTimerAuto("",0);
2295  okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
2296 
2297  UInt_t ntref=TProcessID::GetObjectCount();
2298  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2299  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2300  }
2301 
2302  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
2303  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2304  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2305 
2306  // propagate tracks to secondary vertex, to compute inv. mass
2307  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2308  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2309 
2310  Double_t momentum[3];
2311  postrack->GetPxPyPz(momentum);
2312  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2313  negtrack->GetPxPyPz(momentum);
2314  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2315 
2316  if(!refill){//skip if it is called in refill step because already checked
2317  // invariant mass cut (try to improve coding here..)
2318  Bool_t okMassCut=kFALSE;
2319  if(!okMassCut && fD0toKpi) if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
2320  if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
2321  if(!okMassCut && fDstar) if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
2322  if(!okMassCut && fCascades) if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
2323  if(!okMassCut) {
2324  //AliDebug(2," candidate didn't pass mass cut");
2325  return 0x0;
2326  }
2327  }
2328  // primary vertex to be used by this candidate
2329  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
2330  if(!primVertexAOD) return 0x0;
2331 
2332  Double_t d0z0[2],covd0z0[3];
2333  postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2334  d0[0] = d0z0[0];
2335  d0err[0] = TMath::Sqrt(covd0z0[0]);
2336  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2337  d0[1] = d0z0[0];
2338  d0err[1] = TMath::Sqrt(covd0z0[0]);
2339  AliAODRecoDecayHF2Prong *the2Prong;
2340  // create the object AliAODRecoDecayHF2Prong
2341  if(!refill){
2342  the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
2343  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2344  UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
2345  the2Prong->SetProngIDs(2,id);
2346  if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
2347  // Add daughter references already here
2348  if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2349 
2350  // select D0->Kpi
2351  if(fD0toKpi) {
2353  if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
2354  }
2355  //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
2356  // select J/psi from B
2357  if(fJPSItoEle) {
2359  }
2360  //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
2361  // select D0->Kpi from Dstar
2362  if(fDstar) {
2363  okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
2364  if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
2365  }
2366  //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
2367  }
2368  }else{
2369  the2Prong =rd;
2370  the2Prong->SetSecondaryVtx(secVert);
2371  secVert->SetParent(the2Prong);
2372  AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2373  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2374  the2Prong->SetPxPyPzProngs(2,px,py,pz);
2375  the2Prong->SetDCA(dca);
2376  the2Prong->Setd0Prongs(2,d0);
2377  the2Prong->Setd0errProngs(2,d0err);
2378  the2Prong->SetCharge(0);
2379  }
2380  delete primVertexAOD; primVertexAOD=NULL;
2381 
2382  // remove the primary vertex (was used only for selection)
2384  the2Prong->UnsetOwnPrimaryVtx();
2385  }
2386 
2387  // get PID info from ESD
2388  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2389  if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
2390  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2391  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2392  Double_t esdpid[10];
2393  for(Int_t i=0;i<5;i++) {
2394  esdpid[i] = esdpid0[i];
2395  esdpid[5+i] = esdpid1[i];
2396  }
2397  the2Prong->SetPID(2,esdpid);
2398  return the2Prong;
2399 }
2400 //----------------------------------------------------------------------------
2402  TObjArray *threeTrackArray,AliVEvent *event,
2403  AliAODVertex *secVert,Double_t dispersion,
2404  const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
2405  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2406  Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
2407 {
2410  // E.Bruna, F.Prino
2411  // AliCodeTimerAuto("",0);
2412 
2413 
2414  UInt_t ntref=TProcessID::GetObjectCount();
2415  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2416  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2417  }
2418 
2419  ok3Prong=kFALSE;
2420  if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
2421 
2422  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2423  Double_t momentum[3];
2424 
2425 
2426  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2427  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2428  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2429 
2430  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2431  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2432  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2433  postrack1->GetPxPyPz(momentum);
2434  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2435  negtrack->GetPxPyPz(momentum);
2436  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2437  postrack2->GetPxPyPz(momentum);
2438  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2439 
2440  // invariant mass cut for D+, Ds, Lc
2441  Bool_t okMassCut=kFALSE;
2442  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2443  if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
2444  if(!okMassCut) {
2445  //AliDebug(2," candidate didn't pass mass cut");
2446  return 0x0;
2447  }
2448 
2449  // primary vertex to be used by this candidate
2450  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2451  if(!primVertexAOD) return 0x0;
2452 
2453  Double_t d0z0[2],covd0z0[3];
2454  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2455  d0[0]=d0z0[0];
2456  d0err[0] = TMath::Sqrt(covd0z0[0]);
2457  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2458  d0[1]=d0z0[0];
2459  d0err[1] = TMath::Sqrt(covd0z0[0]);
2460  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2461  d0[2]=d0z0[0];
2462  d0err[2] = TMath::Sqrt(covd0z0[0]);
2463 
2464 
2465  // create the object AliAODRecoDecayHF3Prong
2466  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2467  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2468  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]));
2469  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]));
2470  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2471 
2472 
2473  // construct the candidate passing a NULL pointer for the secondary vertex to avoid creation of TRef
2474  AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(0x0,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
2475  the3Prong->SetOwnPrimaryVtx(primVertexAOD);
2476  // add a pointer to the secondary vertex via SetOwnSecondaryVtx (no TRef created)
2477  AliAODVertex* ownsecv=secVert->CloneWithoutRefs();
2478  the3Prong->SetOwnSecondaryVtx(ownsecv);
2479  UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
2480  the3Prong->SetProngIDs(3,id);
2481 
2482  delete primVertexAOD; primVertexAOD=NULL;
2483 
2484  // disable PID, which requires the TRefs to the daughter tracks
2485  fCutsDplustoKpipi->SetUsePID(kFALSE);
2486  fCutsDstoKKpi->SetUsePID(kFALSE);
2487  fCutsLctopKpi->SetUsePID(kFALSE);
2488 
2489  // select D+->Kpipi, Ds->KKpi, Lc->pKpi
2490  if(f3Prong) {
2491  ok3Prong = kFALSE;
2492 
2494  ok3Prong = kTRUE;
2496  }
2497  if(useForDs && fOKInvMassDs){
2499  ok3Prong = kTRUE;
2501  }
2502  }
2503  if(useForLc && fOKInvMassLc){
2505  ok3Prong = kTRUE;
2507  }
2508  }
2509  }
2510  //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
2511  the3Prong->UnsetOwnSecondaryVtx();
2513  the3Prong->UnsetOwnPrimaryVtx();
2514  }
2515 
2516  // Add TRefs to secondary vertex and daughter tracks only for candidates passing the filtering cuts
2517  if(ok3Prong && fInputAOD){
2518  the3Prong->SetSecondaryVtx(secVert);
2519  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2520  }
2521 
2522  // get PID info from ESD
2523  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2524  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2525  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2526  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2527  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2528  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2529 
2530  Double_t esdpid[15];
2531  for(Int_t i=0;i<5;i++) {
2532  esdpid[i] = esdpid0[i];
2533  esdpid[5+i] = esdpid1[i];
2534  esdpid[10+i] = esdpid2[i];
2535  }
2536  the3Prong->SetPID(3,esdpid);
2537 
2538  return the3Prong;
2539 }
2540 //----------------------------------------------------------------------------
2542  TObjArray *threeTrackArray,AliVEvent *event,
2543  AliAODVertex *secVert,Double_t dispersion,
2544  Double32_t dist12, Double32_t dist23,
2545  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2547 {
2548  // Fill on-the-fly the 3prong data member missing info
2549  // set TRef of vertex and daughters
2550  // do not recalculate the two-track secondary vertex dist12 and dist23
2551  // because they are stored in the AliAODRecoDecayHF3Prong candidate
2552  // reconstructed in the FindCandidate step
2553  // do not check if it is a Lambdac, Ds or Dplus because it is already check in the FindCandidate step
2554  // and the info is stored
2555  // AliCodeTimerAuto("",0);
2556 
2557  UInt_t ntref=TProcessID::GetObjectCount();
2558  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2559  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2560  }
2561 
2562  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2563  Double_t momentum[3];
2564 
2565 
2566  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2567  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2568  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2569 
2570  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2571  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2572  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2573  postrack1->GetPxPyPz(momentum);
2574  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2575  negtrack->GetPxPyPz(momentum);
2576  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2577  postrack2->GetPxPyPz(momentum);
2578  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2579  // primary vertex to be used by this candidate
2580  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2581  if(!primVertexAOD) return 0x0;
2582  Double_t d0z0[2],covd0z0[3];
2583  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2584  d0[0]=d0z0[0];
2585  d0err[0] = TMath::Sqrt(covd0z0[0]);
2586  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2587  d0[1]=d0z0[0];
2588  d0err[1] = TMath::Sqrt(covd0z0[0]);
2589  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2590  d0[2]=d0z0[0];
2591  d0err[2] = TMath::Sqrt(covd0z0[0]);
2592 
2593  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2594  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2595  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2596 
2597  rd->SetSecondaryVtx(secVert);
2598  secVert->SetParent(rd);
2599  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2600 
2601  rd->SetPxPyPzProngs(3,px,py,pz);
2602  rd->SetDCAs(3,dca);
2603  rd->Setd0Prongs(3,d0);
2604  rd->Setd0errProngs(3,d0err);
2605  rd->SetCharge(charge);
2606  rd->SetOwnPrimaryVtx(primVertexAOD);
2607  rd->SetSigmaVert(dispersion);
2608  delete primVertexAOD; primVertexAOD=NULL;
2609 
2611  rd->UnsetOwnPrimaryVtx();
2612  }
2613 
2614  // get PID info from ESD
2615  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2616  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2617  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2618  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2619  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2620  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2621 
2622  Double_t esdpid[15];
2623  for(Int_t i=0;i<5;i++) {
2624  esdpid[i] = esdpid0[i];
2625  esdpid[5+i] = esdpid1[i];
2626  esdpid[10+i] = esdpid2[i];
2627  }
2628  rd->SetPID(3,esdpid);
2629  return rd;
2630 }
2631 //----------------------------------------------------------------------------
2633  TObjArray *fourTrackArray,AliVEvent *event,
2634  AliAODVertex *secVert,
2635  const AliAODVertex *vertexp1n1,
2636  const AliAODVertex *vertexp1n1p2,
2637  Double_t dcap1n1,Double_t dcap1n2,
2638  Double_t dcap2n1,Double_t dcap2n2,
2639  Bool_t &ok4Prong)
2640 {
2643  // G.E.Bruno, R.Romita
2644  // AliCodeTimerAuto("",0);
2645 
2646  UInt_t ntref=TProcessID::GetObjectCount();
2647  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2648  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2649  }
2650 
2651  ok4Prong=kFALSE;
2652  if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
2653 
2654  Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
2655 
2656  AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
2657  AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
2658  AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
2659  AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
2660 
2661  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2662  negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2663  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2664  negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2665 
2666  Double_t momentum[3];
2667  postrack1->GetPxPyPz(momentum);
2668  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2669  negtrack1->GetPxPyPz(momentum);
2670  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2671  postrack2->GetPxPyPz(momentum);
2672  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2673  negtrack2->GetPxPyPz(momentum);
2674  px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
2675 
2676  // invariant mass cut for rho or D0 (try to improve coding here..)
2677  Bool_t okMassCut=kFALSE;
2678  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2679  if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
2680  if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
2681  }
2682  if(!okMassCut) {
2683  //if(fDebug) printf(" candidate didn't pass mass cut\n");
2684  //printf(" candidate didn't pass mass cut\n");
2685  return 0x0;
2686  }
2687 
2688  // primary vertex to be used by this candidate
2689  AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
2690  if(!primVertexAOD) return 0x0;
2691 
2692  Double_t d0z0[2],covd0z0[3];
2693  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2694  d0[0]=d0z0[0];
2695  d0err[0] = TMath::Sqrt(covd0z0[0]);
2696  negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2697  d0[1]=d0z0[0];
2698  d0err[1] = TMath::Sqrt(covd0z0[0]);
2699  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2700  d0[2]=d0z0[0];
2701  d0err[2] = TMath::Sqrt(covd0z0[0]);
2702  negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2703  d0[3]=d0z0[0];
2704  d0err[3] = TMath::Sqrt(covd0z0[0]);
2705 
2706 
2707  // create the object AliAODRecoDecayHF4Prong
2708  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2709  Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2710  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]));
2711  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]));
2712  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]));
2713  Short_t charge=0;
2714  AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2715  the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2716  UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2717  the4Prong->SetProngIDs(4,id);
2718 
2719  delete primVertexAOD; primVertexAOD=NULL;
2720 
2722 
2723 
2725  the4Prong->UnsetOwnPrimaryVtx();
2726  }
2727 
2728 
2729  // get PID info from ESD
2730  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2731  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2732  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2733  if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2734  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2735  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2736  Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2737  if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2738 
2739  Double_t esdpid[20];
2740  for(Int_t i=0;i<5;i++) {
2741  esdpid[i] = esdpid0[i];
2742  esdpid[5+i] = esdpid1[i];
2743  esdpid[10+i] = esdpid2[i];
2744  esdpid[15+i] = esdpid3[i];
2745  }
2746  the4Prong->SetPID(4,esdpid);
2747 
2748  return the4Prong;
2749 }
2750 //----------------------------------------------------------------------------------
2752  //assign and save in fAODMap the index of the AliAODTrack track
2753  //ordering them on the basis of selected criteria
2754 
2755  fAODMapSize = 100000;
2756  fAODMap = new Int_t[fAODMapSize];
2757  AliAODTrack *track=0;
2758  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2759  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
2760  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
2761  if(!track) AliFatal("Not a standard AOD");
2762  // skip pure ITS SA tracks
2763  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2764 
2765  // skip tracks without ITS
2766  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2767 
2768  // TEMPORARY: check that the cov matrix is there
2769  Double_t covtest[21];
2770  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2771  //
2772 
2773  Int_t ind = (Int_t)track->GetID();
2774  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2775  }
2776  return;
2777 }
2778 //-----------------------------------------------------------------------------
2779 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2780  AliVEvent *event) const
2781 {
2783  //AliCodeTimerAuto("",0);
2784 
2785  AliESDVertex *vertexESD = 0;
2786  AliAODVertex *vertexAOD = 0;
2787 
2788 
2790  // primary vertex from the input event
2791 
2792  vertexESD = new AliESDVertex(*fV1);
2793 
2794  } else {
2795  // primary vertex specific to this candidate
2796 
2797  Int_t nTrks = trkArray->GetEntriesFast();
2798  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2799 
2801  // recalculating the vertex
2802 
2803  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2804  Float_t diamondcovxy[3];
2805  event->GetDiamondCovXY(diamondcovxy);
2806  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2807  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2808  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2809  vertexer->SetVtxStart(diamond);
2810  delete diamond; diamond=NULL;
2811  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2812  vertexer->SetOnlyFitter();
2813  }
2814  Int_t skipped[1000];
2815  Int_t nTrksToSkip=0,id;
2816  AliExternalTrackParam *t = 0;
2817  for(Int_t i=0; i<nTrks; i++) {
2818  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2819  id = (Int_t)t->GetID();
2820  if(id<0) continue;
2821  skipped[nTrksToSkip++] = id;
2822  }
2823  // TEMPORARY FIX
2824  // For AOD, skip also tracks without covariance matrix
2825  if(fInputAOD) {
2826  Double_t covtest[21];
2827  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2828  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2829  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2830  id = (Int_t)vtrack->GetID();
2831  if(id<0) continue;
2832  skipped[nTrksToSkip++] = id;
2833  }
2834  }
2835  }
2836  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2837  //
2838  vertexer->SetSkipTracks(nTrksToSkip,skipped);
2839  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2840 
2841  } else if(fRmTrksFromPrimVtx && nTrks>0) {
2842  // removing the prongs tracks
2843 
2844  TObjArray rmArray(nTrks);
2845  UShort_t *rmId = new UShort_t[nTrks];
2846  AliESDtrack *esdTrack = 0;
2847  AliESDtrack *t = 0;
2848  for(Int_t i=0; i<nTrks; i++) {
2849  t = (AliESDtrack*)trkArray->UncheckedAt(i);
2850  esdTrack = new AliESDtrack(*t);
2851  rmArray.AddLast(esdTrack);
2852  if(esdTrack->GetID()>=0) {
2853  rmId[i]=(UShort_t)esdTrack->GetID();
2854  } else {
2855  rmId[i]=9999;
2856  }
2857  }
2858  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2859  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2860  delete [] rmId; rmId=NULL;
2861  rmArray.Delete();
2862 
2863  }
2864 
2865  if(!vertexESD) return vertexAOD;
2866  if(vertexESD->GetNContributors()<=0) {
2867  //AliDebug(2,"vertexing failed");
2868  delete vertexESD; vertexESD=NULL;
2869  return vertexAOD;
2870  }
2871 
2872  delete vertexer; vertexer=NULL;
2873 
2874  }
2875 
2876  // convert to AliAODVertex
2877  Double_t pos[3],cov[6],chi2perNDF;
2878  vertexESD->GetXYZ(pos); // position
2879  vertexESD->GetCovMatrix(cov); //covariance matrix
2880  chi2perNDF = vertexESD->GetChi2toNDF();
2881  delete vertexESD; vertexESD=NULL;
2882 
2883  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2884 
2885  return vertexAOD;
2886 }
2887 //-----------------------------------------------------------------------------
2890 
2891  //printf("Preselections:\n");
2892  // fTrackFilter->Dump();
2893  if(fSecVtxWithKF) {
2894  printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2895  } else {
2896  printf("Secondary vertex with AliVertexerTracks\n");
2897  }
2898  if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2899  if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2900  if(fD0toKpi) {
2901  printf("Reconstruct D0->Kpi candidates with cuts:\n");
2903  }
2904  if(fDstar) {
2905  printf("Reconstruct D*->D0pi candidates with cuts:\n");
2906  if(fFindVertexForDstar) {
2907  printf(" Reconstruct a secondary vertex for the D*\n");
2908  } else {
2909  printf(" Assume the D* comes from the primary vertex\n");
2910  }
2912  }
2913  if(fJPSItoEle) {
2914  printf("Reconstruct J/psi from B candidates with cuts:\n");
2916  }
2917  if(f3Prong) {
2918  printf("Reconstruct 3 prong candidates.\n");
2919  printf(" D+->Kpipi cuts:\n");
2921  printf(" Ds->KKpi cuts:\n");
2923  printf(" Lc->pKpi cuts:\n");
2925  }
2926  if(f4Prong) {
2927  printf("Reconstruct 4 prong candidates.\n");
2928  printf(" D0->Kpipipi cuts:\n");
2930  }
2931  if(fCascades) {
2932  printf("Reconstruct cascade candidates formed with v0s.\n");
2933  printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
2935  printf(" D+ -> K0s pi cuts:\n");
2937  printf(" Ds -> K0s K cuts:\n");
2939  }
2940 
2941  return;
2942 }
2943 //-----------------------------------------------------------------------------
2945  Double_t &dispersion,Bool_t useTRefArray) const
2946 {
2948  //AliCodeTimerAuto("",0);
2949 
2950  AliESDVertex *vertexESD = 0;
2951  AliAODVertex *vertexAOD = 0;
2952 
2953  if(!fSecVtxWithKF) { // AliVertexerTracks
2954 
2955  fVertexerTracks->SetVtxStart(fV1);
2956  vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2957 
2958  if(!vertexESD) return vertexAOD;
2959 
2960  if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2961  //AliDebug(2,"vertexing failed");
2962  delete vertexESD; vertexESD=NULL;
2963  return vertexAOD;
2964  }
2965 
2966  Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
2967  if(vertRadius2>8.){
2968  // vertex outside beam pipe, reject candidate to avoid propagation through material
2969  delete vertexESD; vertexESD=NULL;
2970  return vertexAOD;
2971  }
2972 
2973  } else { // Kalman Filter vertexer (AliKFParticle)
2974 
2975  AliKFParticle::SetField(fBzkG);
2976 
2977  AliKFVertex vertexKF;
2978 
2979  Int_t nTrks = trkArray->GetEntriesFast();
2980  for(Int_t i=0; i<nTrks; i++) {
2981  AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2982  AliKFParticle daughterKF(*esdTrack,211);
2983  vertexKF.AddDaughter(daughterKF);
2984  }
2985  vertexESD = new AliESDVertex(vertexKF.Parameters(),
2986  vertexKF.CovarianceMatrix(),
2987  vertexKF.GetChi2(),
2988  vertexKF.GetNContributors());
2989 
2990  }
2991 
2992  // convert to AliAODVertex
2993  Double_t pos[3],cov[6],chi2perNDF;
2994  vertexESD->GetXYZ(pos); // position
2995  vertexESD->GetCovMatrix(cov); //covariance matrix
2996  chi2perNDF = vertexESD->GetChi2toNDF();
2997  dispersion = vertexESD->GetDispersion();
2998  delete vertexESD; vertexESD=NULL;
2999 
3000  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
3001  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
3002 
3003  return vertexAOD;
3004 }
3005 //-----------------------------------------------------------------------------
3008  //AliCodeTimerAuto("",0);
3009 
3010  Int_t retval=kFALSE;
3011  Double_t momentum[3];
3012  Double_t px[3],py[3],pz[3];
3013  for(Int_t iTrack=0; iTrack<3; iTrack++){
3014  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3015  track->GetPxPyPz(momentum);
3016  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3017  }
3018  retval = SelectInvMassAndPt3prong(px,py,pz);
3019 
3020  return retval;
3021 }
3022 
3023 //-----------------------------------------------------------------------------
3026  //AliCodeTimerAuto("",0);
3027 
3028  Int_t retval=kFALSE;
3029  Double_t momentum[3];
3030  Double_t px[4],py[4],pz[4];
3031 
3032  for(Int_t iTrack=0; iTrack<4; iTrack++){
3033  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3034  track->GetPxPyPz(momentum);
3035  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3036  }
3037 
3038  retval = SelectInvMassAndPt4prong(px,py,pz);
3039 
3040  return retval;
3041 }
3042 //-----------------------------------------------------------------------------
3045  //AliCodeTimerAuto("",0);
3046 
3047  Int_t retval=kFALSE;
3048  Double_t momentum[3];
3049  Double_t px[2],py[2],pz[2];
3050 
3051  for(Int_t iTrack=0; iTrack<2; iTrack++){
3052  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3053  track->GetPxPyPz(momentum);
3054  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3055  }
3056  retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
3057 
3058  return retval;
3059 }
3060 //-----------------------------------------------------------------------------
3062  Double_t *py,
3063  Double_t *pz){
3065  //AliCodeTimerAuto("",0);
3066 
3067  UInt_t pdg2[2];
3068  Int_t nprongs=2;
3069  Double_t minv2,mrange;
3070  Double_t lolim,hilim;
3071  Double_t minPt=0;
3072  Bool_t retval=kFALSE;
3073 
3074  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3075  fOKInvMassD0=kFALSE;
3076  // pt cut
3078  if(minPt>0.1)
3079  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3080  // mass cut
3081  mrange=fCutsD0toKpi->GetMassCut();
3082  lolim=fMassDzero-mrange;
3083  hilim=fMassDzero+mrange;
3084  pdg2[0]=211; pdg2[1]=321;
3085  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3086  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3087  retval=kTRUE;
3088  fOKInvMassD0=kTRUE;
3089  }
3090  pdg2[0]=321; pdg2[1]=211;
3091  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3092  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3093  retval=kTRUE;
3094  fOKInvMassD0=kTRUE;
3095  }
3096  return retval;
3097 }
3098 
3099 //-----------------------------------------------------------------------------
3101  Double_t *py,
3102  Double_t *pz){
3104  //AliCodeTimerAuto("",0);
3105 
3106  UInt_t pdg2[2];
3107  Int_t nprongs=2;
3108  Double_t minv2,mrange;
3109  Double_t lolim,hilim;
3110  Double_t minPt=0;
3111  Bool_t retval=kFALSE;
3112 
3113  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3114  fOKInvMassJpsi=kFALSE;
3115  // pt cut
3117  if(minPt>0.1)
3118  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3119  // mass cut
3120  mrange=fCutsJpsitoee->GetMassCut();
3121  lolim=fMassJpsi-mrange;
3122  hilim=fMassJpsi+mrange;
3123 
3124  pdg2[0]=11; pdg2[1]=11;
3125  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3126  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3127  retval=kTRUE;
3128  fOKInvMassJpsi=kTRUE;
3129  }
3130 
3131  return retval;
3132 }
3133 //-----------------------------------------------------------------------------
3135  Double_t *py,
3136  Double_t *pz,
3137  Int_t pidLcStatus){
3139  //AliCodeTimerAuto("",0);
3140 
3141  UInt_t pdg3[3];
3142  Int_t nprongs=3;
3143  Double_t minv2,mrange;
3144  Double_t lolim,hilim;
3145  Double_t minPt=0;
3146  Bool_t retval=kFALSE;
3147 
3148 
3149  fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
3150  fOKInvMassDplus=kFALSE;
3151  fOKInvMassDs=kFALSE;
3152  fOKInvMassLc=kFALSE;
3153  // pt cut
3155  minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
3156  if(minPt>0.1)
3157  if(fMassCalc3->Pt2() < minPt*minPt) return retval;
3158  // D+->Kpipi
3159  mrange=fCutsDplustoKpipi->GetMassCut();
3160  lolim=fMassDplus-mrange;
3161  hilim=fMassDplus+mrange;
3162  pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
3163  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3164  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3165  retval=kTRUE;
3166  fOKInvMassDplus=kTRUE;
3167  }
3168  // Ds+->KKpi
3169  mrange=fCutsDstoKKpi->GetMassCut();
3170  lolim=fMassDs-mrange;
3171  hilim=fMassDs+mrange;
3172  pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
3173  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3174  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3175  retval=kTRUE;
3176  fOKInvMassDs=kTRUE;
3177  }
3178  pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
3179  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3180  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3181  retval=kTRUE;
3182  fOKInvMassDs=kTRUE;
3183  }
3184  // Lc->pKpi
3185  mrange=fCutsLctopKpi->GetMassCut();
3186  lolim=fMassLambdaC-mrange;
3187  hilim=fMassLambdaC+mrange;
3188  if(pidLcStatus&1){
3189  pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
3190  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3191  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3192  retval=kTRUE;
3193  fOKInvMassLc=kTRUE;
3194  }
3195  }
3196  if(pidLcStatus&2){
3197  pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
3198  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3199  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3200  retval=kTRUE;
3201  fOKInvMassLc=kTRUE;
3202  }
3203  }
3204 
3205  return retval;
3206 }
3207 
3208 //-----------------------------------------------------------------------------
3210  Double_t *py,
3211  Double_t *pz){
3213  //AliCodeTimerAuto("",0);
3214 
3215  UInt_t pdg2[2];
3216  Int_t nprongs=2;
3217  Double_t minv2,mrange;
3218  Double_t lolim,hilim;
3219  Double_t minPt=0;
3220  Bool_t retval=kFALSE;
3221 
3222  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3223  fOKInvMassDstar=kFALSE;
3224  // pt cut
3226  if(minPt>0.1)
3227  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3228  // mass cut
3229  pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
3230  mrange=fCutsDStartoKpipi->GetMassCut();
3231  lolim=fMassDstar-mrange;
3232  hilim=fMassDstar+mrange;
3233  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3234  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3235  retval=kTRUE;
3236  fOKInvMassDstar=kTRUE;
3237  }
3238 
3239  return retval;
3240 }
3241 
3242 //-----------------------------------------------------------------------------
3244  Double_t *py,
3245  Double_t *pz){
3247  //AliCodeTimerAuto("",0);
3248 
3249  UInt_t pdg4[4];
3250  Int_t nprongs=4;
3251  Double_t minv2,mrange;
3252  Double_t lolim,hilim;
3253  Double_t minPt=0;
3254  Bool_t retval=kFALSE;
3255 
3256  // D0->Kpipipi without PID
3257  fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
3258  fOKInvMassD0to4p=kFALSE;
3259  // pt cut
3261  if(minPt>0.1)
3262  if(fMassCalc4->Pt2() < minPt*minPt) return retval;
3263  // mass cut
3264  mrange=fCutsD0toKpipipi->GetMassCut();
3265  lolim=fMassDzero-mrange;
3266  hilim=fMassDzero+mrange;
3267 
3268  pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
3269  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3270  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3271  retval=kTRUE;
3272  fOKInvMassD0to4p=kTRUE;
3273  }
3274 
3275  pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
3276  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3277  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3278  retval=kTRUE;
3279  fOKInvMassD0to4p=kTRUE;
3280  }
3281 
3282  pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
3283  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3284  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3285  retval=kTRUE;
3286  fOKInvMassD0to4p=kTRUE;
3287  }
3288 
3289  pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
3290  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3291  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3292  retval=kTRUE;
3293  fOKInvMassD0to4p=kTRUE;
3294  }
3295 
3296  return retval;
3297 }
3298 //-----------------------------------------------------------------------------
3300  Double_t *py,
3301  Double_t *pz){
3303  //AliCodeTimerAuto("",0);
3304 
3305  UInt_t pdg2[2];
3306  Int_t nprongs=2;
3307  Double_t minv2,mrange;
3308  Double_t lolim,hilim;
3309  Double_t minPt=0;
3310  Bool_t retval=kFALSE;
3311 
3312  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3313  minPt=fCutsLctoV0->GetMinPtCandidate();
3314  if(minPt>0.1)
3315  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3316 
3317  fOKInvMassLctoV0 = kFALSE;
3318 
3319  // LambdaC candidate
3320  if (fCutsLctoV0) {
3321  mrange = fCutsLctoV0->GetMassCut();
3322  lolim = fMassLambdaC - mrange;
3323  hilim = fMassLambdaC + mrange;
3324  pdg2[0] = 2212; pdg2[1] = 310;
3325  minv2 = fMassCalc2->InvMass2(2, pdg2);
3326  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3327  retval = kTRUE;
3328  fOKInvMassLctoV0 = kTRUE;
3329  }
3330  pdg2[0] = 211; pdg2[1] = 3122;
3331  minv2 = fMassCalc2->InvMass2(2, pdg2);
3332  if ((minv2>lolim*lolim) && (minv2<hilim*hilim)) {
3333  retval = kTRUE;
3334  fOKInvMassLctoV0 = kTRUE;
3335  }
3336  }
3337 
3338  // D+ cascade candidate
3339  if (fCutsDplustoK0spi) {
3340  mrange = fCutsDplustoK0spi->GetMassCut();
3341  lolim = fMassDplus - mrange;
3342  hilim = fMassDplus + mrange;
3343  pdg2[0] = 211; pdg2[1] = 310;
3344  minv2 = fMassCalc2->InvMass2(2, pdg2);
3345  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3346  retval = kTRUE;
3347  }
3348  }
3349 
3350  // Ds cascade candidate
3351  if (fCutsDstoK0sK) {
3352  mrange = fCutsDstoK0sK->GetMassCut();
3353  lolim = fMassDs - mrange;
3354  hilim = fMassDs + mrange;
3355  pdg2[0] = 321; pdg2[1] = 310;
3356  minv2 = fMassCalc2->InvMass2(2, pdg2);
3357  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3358  retval = kTRUE;
3359  }
3360  }
3361 
3362  return retval;
3363 }
3364 //-----------------------------------------------------------------------------
3366  Int_t trkEntries,
3367  TObjArray &seleTrksArray,
3368  TObjArray &tracksAtVertex,
3369  Int_t &nSeleTrks,
3370  UChar_t *seleFlags,Int_t *evtNumber)
3371 {
3377  //AliCodeTimerAuto("",0);
3378 
3379  const AliVVertex *vprimary = event->GetPrimaryVertex();
3380 
3381  if(fV1) { delete fV1; fV1=NULL; }
3382  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
3383 
3384  Int_t nindices=0;
3385  UShort_t *indices = 0;
3386  Double_t pos[3],cov[6];
3387  const Int_t entries = event->GetNumberOfTracks();
3388  AliCentrality* cent;
3389 
3390  if(!fInputAOD) { // ESD
3391  fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
3392  cent=((AliESDEvent*)event)->GetCentrality();
3393  } else { // AOD
3394  vprimary->GetXYZ(pos);
3395  vprimary->GetCovarianceMatrix(cov);
3396  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
3397  if(entries<=0) return;
3398  indices = new UShort_t[entries];
3399  memset(indices,0,sizeof(UShort_t)*entries);
3400  fAODMapSize = 100000;
3401  fAODMap = new Int_t[fAODMapSize];
3402  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
3403  cent=((AliAODEvent*)event)->GetCentrality();
3404  }
3405  Float_t centperc=0.1;
3406  if(event->GetRunNumber()<244824){
3407  centperc=cent->GetCentralityPercentile("V0M");
3408  }else{
3409  AliMultSelection *multSelection = (AliMultSelection * ) event->FindListObject("MultSelection");
3410  if(multSelection){
3411  centperc=multSelection->GetMultiplicityPercentile("V0M");
3412  Int_t qual = multSelection->GetEvSelCode();
3413  if(qual == 199 ) centperc=0.1; // use central cuts as default
3414  }
3415  }
3416  Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE,okBachelor=kFALSE;
3417  nSeleTrks=0;
3418 
3419  // transfer ITS tracks from event to arrays
3420  for(Int_t i=0; i<entries; i++) {
3421  AliVTrack *track;
3422  track = (AliVTrack*)event->GetTrack(i);
3423 
3424  // skip pure ITS SA tracks
3425  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
3426 
3427  // skip tracks without ITS
3428  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
3429 
3430  // skip tracks with negative ID
3431  // (these are duplicated TPC-only AOD tracks, for jet analysis...)
3432  if(track->GetID()<0) continue;
3433 
3434  // TEMPORARY: check that the cov matrix is there
3435  Double_t covtest[21];
3436  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
3437  //
3438 
3439  if(fInputAOD) {
3440  AliAODTrack *aodt = (AliAODTrack*)track;
3441  if(aodt->GetUsedForPrimVtxFit()) {
3442  indices[nindices]=aodt->GetID(); nindices++;
3443  }
3444  Int_t ind = (Int_t)aodt->GetID();
3445  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
3446  }
3447 
3448  AliESDtrack *esdt = 0;
3449 
3450  if(!fInputAOD) {
3451  esdt = (AliESDtrack*)track;
3452  } else {
3453  esdt = new AliESDtrack(track);
3454  }
3455 
3456  // single track selection
3457  okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE; okBachelor=kFALSE;
3458  if(fMixEvent && i<trkEntries){
3459  evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
3460  const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
3461  Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
3462  eventVtx->GetXYZ(vtxPos);
3463  vprimary->GetXYZ(primPos);
3464  eventVtx->GetCovarianceMatrix(primCov);
3465  for(Int_t ind=0;ind<3;ind++){
3466  trasl[ind]=vtxPos[ind]-primPos[ind];
3467  }
3468 
3469  Bool_t isTransl=esdt->Translate(trasl,primCov);
3470  if(!isTransl) {
3471  delete esdt;
3472  esdt = NULL;
3473  continue;
3474  }
3475  }
3476 
3477  if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong,okBachelor) && nSeleTrks<trkEntries) {
3478  esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
3479  seleTrksArray.AddLast(esdt);
3480  tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
3481  seleFlags[nSeleTrks]=0;
3482  if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
3483  if(okFor3Prong) SETBIT(seleFlags[nSeleTrks],kBit3Prong);
3484  if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
3485  if(okBachelor) SETBIT(seleFlags[nSeleTrks],kBitBachelor);
3486  // Check the PID
3487  SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
3488  SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3489  SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3490  Bool_t useTPC=kTRUE;
3491  if(fUseTOFPID){
3492  Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
3493  if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
3494  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3495  }
3496  Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
3497  if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
3498  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3499  }
3500  Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
3501  if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
3502  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3503  }
3504  if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
3505  }
3506  if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
3507  Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
3508  if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
3509  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3510  }
3511  Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
3512  if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
3513  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3514  }
3515  Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
3516  if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
3517  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3518  }
3519  }
3520  nSeleTrks++;
3521  } else {
3522  if(fInputAOD) delete esdt;
3523  esdt = NULL;
3524  continue;
3525  }
3526 
3527  } // end loop on tracks
3528 
3529  // primary vertex from AOD
3530  if(fInputAOD) {
3531  delete fV1; fV1=NULL;
3532  vprimary->GetXYZ(pos);
3533  vprimary->GetCovarianceMatrix(cov);
3534  Double_t chi2toNDF = vprimary->GetChi2perNDF();
3535  Int_t ncontr=nindices;
3536  if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
3537  Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
3538  fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
3539  fV1->SetTitle(vprimary->GetTitle());
3540  fV1->SetIndices(nindices,indices);
3541  }
3542  if(indices) { delete [] indices; indices=NULL; }
3543 
3544 
3545  return;
3546 }
3547 //-----------------------------------------------------------------------------
3549  //
3551  //
3552  if(fUsePidTag && cuts->GetPidHF()) {
3553  Bool_t usepid=cuts->GetIsUsePID();
3554  cuts->SetUsePID(kTRUE);
3555  if(cuts->IsSelectedPID(rd))
3556  rd->SetSelectionBit(bit);
3557  cuts->SetUsePID(usepid);
3558  }
3559  return;
3560 }
3561 //-----------------------------------------------------------------------------
3563  Float_t centralityperc,
3564  Bool_t &okDisplaced,
3565  Bool_t &okSoftPi,
3566  Bool_t &okFor3Prong,
3567  Bool_t &okBachelor) const
3568 {
3570 
3571  // this is needed to store the impact parameters
3572  //AliCodeTimerAuto("",0);
3573 
3574  trk->RelateToVertex(fV1,fBzkG,kVeryBig);
3575 
3576  UInt_t selectInfo;
3577  //
3578  // Track selection, displaced tracks -- 2 prong
3579  selectInfo = 0;
3580  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3582  // central PbPb events, tighter cuts
3583  selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
3584  }else{
3585  // standard cuts
3586  if(fTrackFilter) {
3587  selectInfo = fTrackFilter->IsSelected(trk);
3588  }
3589  }
3590  if(selectInfo) okDisplaced=kTRUE;
3591 
3592  // Track selection, displaced tracks -- 3 prong
3593  selectInfo = 0;
3594  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3596  // central PbPb events, tighter cuts
3597  selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
3598  }else{
3599  // standard cuts
3600  if(fTrackFilter) {
3601  selectInfo = fTrackFilter->IsSelected(trk);
3602  }
3603  }
3604  if(selectInfo) okFor3Prong=kTRUE;
3605 
3606  // Track selection, soft pions
3607  selectInfo = 0;
3608  if(fDstar && fTrackFilterSoftPi) {
3609  selectInfo = fTrackFilterSoftPi->IsSelected(trk);
3610  }
3611  if(selectInfo) okSoftPi=kTRUE;
3612 
3613  // Track selection, bachelor
3614  selectInfo = 0;
3615  if(fCascades){
3617  selectInfo = fTrackFilterBachelor->IsSelected(trk);
3618  }else{
3619  if(okDisplaced) selectInfo=1;
3620  }
3621  }
3622  if(selectInfo) okBachelor=kTRUE;
3623 
3624  if(okDisplaced || okSoftPi || okFor3Prong || okBachelor) return kTRUE;
3625 
3626  return kFALSE;
3627 }
3628 
3629 
3630 //-----------------------------------------------------------------------------
3631 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
3637  //
3638  //AliCodeTimerAuto("",0);
3639  Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
3640  AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
3641 
3642  // create the v0 neutral track to compute the DCA to the primary vertex
3643  Double_t xyz[3], pxpypz[3];
3644  esdV0->XvYvZv(xyz);
3645  esdV0->PxPyPz(pxpypz);
3646  Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
3647  AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
3648  if(!trackesdV0) {
3649  delete vertexV0;
3650  return 0;
3651  }
3652  Double_t d0z0[2],covd0z0[3];
3653  AliAODVertex *primVertexAOD = PrimaryVertex();
3654  trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3655  Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
3656  // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
3657  Double_t dcaV0DaughterToPrimVertex[2];
3658  AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
3659  AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
3660  if( !posV0track || !negV0track) {
3661  if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
3662  delete vertexV0;
3663  delete primVertexAOD;
3664  return 0;
3665  }
3666  posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3667  // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
3668  // else
3669  dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
3670  negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3671  // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
3672  // else
3673  dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
3674  Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
3675  Double_t pmom[3],nmom[3];
3676  esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
3677  esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
3678 
3679  AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
3680  aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
3681 
3682  delete trackesdV0;
3683  delete primVertexAOD;
3684 
3685  return aodV0;
3686 }
3687 //-----------------------------------------------------------------------------
3688 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
3690  //AliCodeTimerAuto("",0);
3691 
3692  const Double_t *par=extpar->GetParameter();
3693  const Double_t *cov=extpar->GetCovariance();
3694  Double_t alpha=extpar->GetAlpha();
3695  Double_t x=extpar->GetX();
3696  esdt->Set(x,alpha,par,cov);
3697  return;
3698 }
3699 //-----------------------------------------------------------------------------
3702 
3703  fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
3704  fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
3705  fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
3706  fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
3707  fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
3708  fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
3709 }
3710 //-----------------------------------------------------------------------------
3712  //
3714  //
3715 
3716 
3717  //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
3718 
3719 
3721  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
3722  return kFALSE;
3723  }
3724 
3726  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHfCutsDplustoK0spi (%d)\n",fV0TypeForCascadeVertex,fCutsDplustoK0spi->GetV0Type());
3727  return kFALSE;
3728  }
3729 
3731  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsDstoK0sK (%d)\n",fV0TypeForCascadeVertex,fCutsDstoK0sK->GetV0Type());
3732  return kFALSE;
3733  }
3734 
3735  return kTRUE;
3736 }
3737 //-----------------------------------------------------------------------------
3738 
Int_t charge
AliAODRecoDecay * fMassCalc4
for 3 prong
Bool_t fUsePidTag
upper momentum limit to apply TPC PID
void SetSelectionBit(Int_t i)
selection map
double Double_t
Definition: External.C:58
void SetupPID(AliVEvent *event)
AliAnalysisFilter * fTrackFilterSoftPi
Track Filter for displaced vertices in PbPb central events (tighter cuts) for 3 prong (D+...
AliVertexerTracks * fVertexerTracks
Float_t GetMassCut(Int_t iPtBin=0) const
AliPIDResponse * fPidResponse
event mixing
AliRDHFCutsDstoK0sK * fCutsDstoK0sK
D+->Kpipi cuts.
AliAODVertex * ReconstructSecondaryVertex(TObjArray *trkArray, Double_t &dispersion, Bool_t useTRefArray=kTRUE) const
Double_t fnSigmaTOFPionHi
Low cut value on n. of sigmas for pi TOF PID.
Double_t fnSigmaTOFProtonHi
Low cut value on n. of sigmas for p TOF PID.
AliAODVertex * PrimaryVertex(const TObjArray *trkArray=0x0, AliVEvent *event=0x0) const
AliAnalysisFilter * fTrackFilter3prongCentral
Track Filter for displaced vertices in PbPb central events (tighter cuts) for 2 prong (D0->Kpi) ...
AliAnalysisVertexingHF & operator=(const AliAnalysisVertexingHF &source)
virtual void PrintAll() const
Float_t GetMassCut(Int_t iPtBin=0) const
AliAODv0 * Getv0() const
AliRDHFCutsDstoKKpi * fCutsDstoKKpi
Ds->K0s+K.
Int_t IsD0FromDStarSelected(Double_t pt, TObject *obj, Int_t selectionLevel, AliAODEvent *aod) const
Bool_t fSecVtxWithKF
z componenent of field in kG
TList * list
TDirectory file where lists per trigger are stored in train ouput.
AliRDHFCutsLctopKpi * fCutsLctopKpi
Ds->KKpi cuts.
Bool_t SelectInvMassAndPt3prong(Double_t *px, Double_t *py, Double_t *pz, Int_t pidLcStatus=3)
Int_t fnTrksTotal
triplet fullfilling Lc inv mass selection
Double_t fnSigmaTOFProtonLow
High cut value on n. of sigmas for p TPC PID.
Bool_t fOKInvMassDplus
pair fullfilling Jpsi inv mass selection
Double_t fnSigmaTPCKaonLow
High cut value on n. of sigmas for pi TOF PID.
Bool_t SingleTrkCuts(AliESDtrack *trk, Float_t centralityperc, Bool_t &okDisplaced, Bool_t &okSoftPi, Bool_t &ok3prong, Bool_t &okBachelor) const
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Bool_t fUseTPCPIDOnlyIfNoTOF
switch use/not use TOF PID
AliRDHFCutsD0toKpi * fCutsD0toKpi
Track Filter for bachelor.
Bool_t fMixEvent
Like-sign triplets.
virtual void DeleteRecoD()
Float_t GetMassCut(Int_t iPtBin=0) const
Bool_t fUseTOFPID
switch use/not use TPC PID
Bool_t fOKInvMassD0to4p
combination fullfilling D* inv mass selection
Bool_t fOKInvMassLc
triplet fullfilling Ds inv mass selection
Float_t GetMassCut(Int_t iPtBin=0) const
Int_t * fAODMap
map between index and ID for AOD tracks
Double_t fnSigmaTPCPionLow
flag to control usage of PID tagging
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *)
Definition: AliRDHFCuts.h:288
Bool_t fFindVertexForCascades
reconstruct a secondary vertex or assume it's from the primary vertex
Double_t fnSigmaTPCKaonHi
Low cut value on n. of sigmas for K TPC PID.
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:237
Bool_t fUseKaonPIDfor3Prong
PID response.
AliRDHFCutsDplustoK0spi * fCutsDplustoK0spi
J/psi->ee cuts.
AliRDHFCutsDplustoKpipi * fCutsDplustoKpipi
D+->K0s+pi.
AliRDHFCutsD0toKpipipi * fCutsD0toKpipipi
Lc –> v0 + bachelor cuts.
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t fLikeSign3prong
Like-sign pairs.
Bool_t SelectInvMassAndPtCascade(Double_t *px, Double_t *py, Double_t *pz)
AliAODRecoDecayHF4Prong * Make4Prong(TObjArray *fourTrackArray, AliVEvent *event, AliAODVertex *secVert, const AliAODVertex *vertexp1n1, const AliAODVertex *vertexp1n1p2, Double_t dcap1n1, Double_t dcap1n2, Double_t dcap2n1, Double_t dcap2n2, Bool_t &ok4Prong)
Class for cuts on AOD reconstructed D+->Kpipi.
Float_t GetDCACut(Int_t iPtBin=0) const
Double_t fnSigmaTPCProtonLow
High cut value on n. of sigmas for K TOF PID.
void SetParametersAtVertex(AliESDtrack *esdt, const AliExternalTrackParam *extpar) const
int Int_t
Definition: External.C:63
Float_t GetMassCut(Int_t iPtBin=0) const
Float_t GetDCACut(Int_t iPtBin=0) const
Double_t fMaxMomForTPCPid
use TPC PID only for tracks that without TOF
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
unsigned int UInt_t
Definition: External.C:33
Int_t GetIsFilled() const
float Float_t
Definition: External.C:68
Class for cuts on AOD reconstructed Ds->K0S+K.
Float_t GetMassCut(Int_t iPtBin=0) const
Float_t GetMassCut(Int_t iPtBin=0) const
AliAnalysisFilter * fTrackFilter2prongCentral
Track Filter for displaced vertices.
Double_t fnSigmaTOFKaonLow
High cut value on n. of sigmas for K TPC PID.
Float_t fMaxCentPercentileForTightCuts
High cut value on n. of sigmas for p TOF PID.
AliAODTrack * GetBachelor() const
Bool_t fOKInvMassDstar
triplet fullfilling Lc inv mass selection
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Double_t fBzkG
vertexer, to compute secondary vertices
Float_t GetDCACut(Int_t iPtBin=0) const
Bool_t fUsePIDforLc2V0
PID for Lambdac: 0=no, 1=proton, 2=p and pi.
void SetProngIDs(Int_t nIDs, UShort_t *id)
Class for cuts on AOD reconstructed D+->K0S+pi.
virtual void PrintAll() const
UShort_t GetProngID(Int_t ip) const
Float_t GetMassCut(Int_t iPtBin=0) const
Bool_t fUseKaonPIDforDs
PID for Lambdac 2 V0: 0=no, 1=proton,.
Bool_t fUseTPCPID
Kaon PID usage for Ds.
Bool_t fOKInvMassLctoV0
4tracks fullfilling D0 inv mass selection
void SetPrimaryVtxRef(TObject *vtx)
primary vertex
short Short_t
Definition: External.C:23
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Bool_t fFindVertexForDstar
pointer to list of cuts for output file
Double_t fnSigmaTOFKaonHi
Low cut value on n. of sigmas for K TOF PID.
AliAODRecoDecay * fMassCalc2
to go faster in PbPb
Float_t GetDCACut(Int_t iPtBin=0) const
void FindCandidates(AliVEvent *event, TClonesArray *aodVerticesHFTClArr, TClonesArray *aodD0toKpiTClArr, TClonesArray *aodJPSItoEleTClArr, TClonesArray *aodCharm3ProngTClArr, TClonesArray *aodCharm4ProngTClArr, TClonesArray *aodDstarTClArr, TClonesArray *aodCascadesTClArr, TClonesArray *aodLikeSign2ProngTClArr, TClonesArray *aodLikeSign3ProngTClArr)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
void AddDaughterRefs(AliAODVertex *v, const AliVEvent *event, const TObjArray *trkArray) const
Bool_t PreSelect(TObject *obj, AliAODv0 *v0, AliVTrack *bachelorTrack)
Int_t fUsePIDforLc
Kaon PID usage for 3 prongs.
Bool_t IsEventSelected(AliVEvent *event)
Bool_t fLikeSign
cascades, Lc –> v0+track, D+ –> K0s+Pi, Ds –> K0s+K
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:206
virtual void PrintAll() const
Float_t GetDCACut(Int_t iPtBin=0) const
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
AliRDHFCutsJpsitoee * fCutsJpsitoee
D0->Kpi cuts.
Bool_t SelectInvMassAndPtJpsiee(Double_t *px, Double_t *py, Double_t *pz)
Bool_t fOKInvMassD0
for 4 prong
void Setd0errProngs(Int_t nprongs, Double_t *d0)
Bool_t fOKInvMassJpsi
pair fullfilling D0 inv mass selection
Bool_t SelectInvMassAndPtDstarD0pi(Double_t *px, Double_t *py, Double_t *pz)
Int_t fV0TypeForCascadeVertex
reconstruct a secondary vertex or assume it's from the primary vertex
unsigned short UShort_t
Definition: External.C:28
AliAnalysisFilter * fTrackFilterBachelor
Track Filter for D* soft pion.
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:260
Bool_t fMassCutBeforeVertexing
Select which V0 type we want to use for the cascas.
Int_t fAODMapSize
input from AOD (kTRUE) or ESD (kFALSE)
void FixReferences(AliAODEvent *aod)
void SetSelectionBitForPID(AliRDHFCuts *cuts, AliAODRecoDecayHF *rd, Int_t bit)
AliAODRecoCascadeHF * MakeCascade(TObjArray *twoTrackArray, AliVEvent *event, AliAODVertex *secVert, AliAODRecoDecayHF2Prong *rd2Prong, Double_t dca, Bool_t &okDstar)
void SelectTracksAndCopyVertex(const AliVEvent *event, Int_t trkEntries, TObjArray &seleTrksArray, TObjArray &tracksAtVertex, Int_t &nSeleTrks, UChar_t *seleFlags, Int_t *evtNumber)
bool Bool_t
Definition: External.C:53
void SetSigmaVert(Double_t sigmaVert)
Bool_t fRecoPrimVtxSkippingTrks
if kTRUE use KF vertexer, else AliVertexerTracks
AliAODRecoDecay * fMassCalc3
for 2 prong
TList * fListOfCuts
Dstar->D0pi cuts.
Bool_t RecoSecondaryVertexForCascades(AliVEvent *event, AliAODRecoCascadeHF *rc)
Float_t GetMassCut(Int_t iPtBin=0) const
Double_t fnSigmaTPCProtonHi
Low cut value on n. of sigmas for p TPC PID.
Double_t GetMinPtCandidate() const
Definition: AliRDHFCuts.h:274
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Float_t GetMassCut(Int_t iPtBin=0) const
AliRDHFCutsLctoV0 * fCutsLctoV0
Lc->pKpi cuts.
Float_t GetDCACut(Int_t iPtBin=0) const
void MapAODtracks(AliVEvent *aod)
AliAODRecoDecayHF3Prong * Make3Prong(TObjArray *threeTrackArray, AliVEvent *event, AliAODVertex *secVert, Double_t dispersion, const AliAODVertex *vertexp1n1, const AliAODVertex *vertexp2n1, Double_t dcap1n1, Double_t dcap2n1, Double_t dcap1p2, Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
AliRDHFCutsDStartoKpipi * fCutsDStartoKpipi
D0->Kpipipi cuts.
AliAnalysisFilter * fTrackFilter
max. centrality percentile for using tight cuts
Float_t GetDCACut(Int_t iPtBin=0) const
Bool_t SelectInvMassAndPtD0Kpi(Double_t *px, Double_t *py, Double_t *pz)
Bool_t GetUsePID(Int_t iPtBin=0) const
AliAODv0 * TransformESDv0toAODv0(AliESDv0 *esdv0, TObjArray *twoTrackArrayV0)
Bool_t SelectInvMassAndPt4prong(Double_t *px, Double_t *py, Double_t *pz)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Double_t fnSigmaTPCPionHi
Low cut value on n. of sigmas for pi TPC PID.
Bool_t fD0toKpi
primary vertex
Bool_t fOKInvMassDs
triplet fullfilling D+ inv mass selection
void SetIsFilled(Int_t filled)
void AddRefs(AliAODVertex *v, AliAODRecoDecayHF *rd, const AliVEvent *event, const TObjArray *trkArray) const
AliAODRecoDecayHF2Prong * Make2Prong(TObjArray *twoTrackArray1, AliVEvent *event, AliAODVertex *secVert, Double_t dcap1n1, Bool_t &okD0, Bool_t &okJPSI, Bool_t &okD0fromDstar, Bool_t refill=kFALSE, AliAODRecoDecayHF2Prong *rd=0x0)
Double_t fnSigmaTOFPionLow
High cut value on n. of sigmas for pi TPC PID.