AliPhysics  251aa1e (251aa1e)
 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 
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 
1705  TObjArray *threeTrackArray = new TObjArray(3);
1706 
1707  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1708  if(!track1)return kFALSE;
1709  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1710  if(!track2)return kFALSE;
1711  AliESDtrack *postrack1 = 0;
1712  AliESDtrack *negtrack1 = 0;
1713  postrack1 = new AliESDtrack(track1);
1714  negtrack1 = new AliESDtrack(track2);
1715 
1716  // DCA between the two tracks
1717  Double_t xdummy, ydummy;
1718  fBzkG = (Double_t)event->GetMagneticField();
1719  Double_t dca12 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1720 
1721  const AliVVertex *vprimary = event->GetPrimaryVertex();
1722  Double_t pos[3];
1723  Double_t cov[6];
1724  vprimary->GetXYZ(pos);
1725  vprimary->GetCovarianceMatrix(cov);
1726  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1727  fV1->GetCovMatrix(cov);
1728  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1729 
1730  AliAODTrack *track3 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(2)]);
1731  if(!track3)return kFALSE;
1732  AliESDtrack *esdt3 = new AliESDtrack(track3);
1733 
1734  Double_t dca2;
1735  Double_t dca3;
1736  threeTrackArray->AddAt(postrack1,0);
1737  threeTrackArray->AddAt(negtrack1,1);
1738  threeTrackArray->AddAt(esdt3,2);
1739  dca2 = esdt3->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1740  dca3 = esdt3->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1741  Double_t dispersion;
1742 
1743  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray, dispersion);
1744  if (!secVert3PrAOD) {
1745  threeTrackArray->Clear();
1746  threeTrackArray->Delete(); delete threeTrackArray;
1747  delete fV1; fV1=0;
1748  delete postrack1; postrack1=NULL;
1749  delete negtrack1; negtrack1=NULL;
1750  delete esdt3; esdt3=NULL;
1751  return kFALSE;
1752  }
1753 
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= Make2Prong(twoTrackArray1, event, vtxRec, dca12, okD0, okJPSI, okD0FromDstar,refill,rd);
1817  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1818  rd->SetIsFilled(2);
1819  delete fV1; fV1=0;
1820  twoTrackArray1->Clear();
1821  twoTrackArray1->Delete(); delete twoTrackArray1;
1822  delete esdt1; esdt1=NULL;
1823  delete esdt2; esdt2=NULL;
1824  return kTRUE;
1825 }
1826 //----------------------------------------------------------------------------
1828  // method to retrieve daughters from trackID
1829  // and fill on-the-fly the data member of rCasc and their AliAODRecoDecayHF2Prong daughters
1830  if(rCasc->GetIsFilled()!=0) return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2: already refilled)
1831  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1832  TObjArray *twoTrackArrayCasc = new TObjArray(2);
1833 
1834  AliAODTrack *trackB =(AliAODTrack*)event->GetTrack(fAODMap[rCasc->GetProngID(0)]);//retrieve bachelor from the trackID through the AOD index map
1835  if(!trackB)return kFALSE;
1836 
1837  AliNeutralTrackParam *trackV0=NULL;
1838  AliAODv0 *v0 =NULL;
1839  AliAODRecoDecayHF2Prong *trackD0=NULL;
1840 
1841  if(DStar){
1842  TClonesArray *inputArrayD0=(TClonesArray*)event->GetList()->FindObject("D0toKpi");
1843  if(!inputArrayD0) return kFALSE;
1844  trackD0=(AliAODRecoDecayHF2Prong*)inputArrayD0->At(rCasc->GetProngID(1));
1845  if(!trackD0)return kFALSE;
1846  if(!FillRecoCand(event,trackD0)) return kFALSE; //fill missing info of the corresponding D0 daughter
1847 
1848  trackV0 = new AliNeutralTrackParam(trackD0);
1849 
1850  }else{//is a V0 candidate
1851  v0 = ((AliAODEvent*)event)->GetV0(rCasc->GetProngID(1));
1852  if(!v0) return kFALSE;
1853  // Define the V0 (neutral) track
1854  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1855  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
1856  }
1857 
1858  AliESDtrack *esdB = new AliESDtrack(trackB);
1859 
1860  twoTrackArrayCasc->AddAt(esdB,0);
1861  twoTrackArrayCasc->AddAt(trackV0,1);
1862 
1863  fBzkG = (Double_t)event->GetMagneticField();
1864  const AliVVertex *vprimary = event->GetPrimaryVertex();
1865 
1866  Double_t pos[3];
1867  Double_t cov[6];
1868  vprimary->GetXYZ(pos);
1869  vprimary->GetCovarianceMatrix(cov);
1870  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1871  fV1->GetCovMatrix(cov);
1872 
1873  Double_t dca = 0.;
1874  AliAODVertex *vtxCasc = 0x0;
1875  Double_t chi2perNDF = fV1->GetChi2toNDF();
1876  if (recoSecVtx) {
1877  Double_t dispersion, xdummy, ydummy;
1878  dca = esdB->GetDCA(trackV0,fBzkG,xdummy,ydummy);
1879  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
1880  vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
1881  } else {
1882  vtxCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
1883  }
1884  if(!vtxCasc) {
1885  twoTrackArrayCasc->Clear();
1886  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1887  delete fV1; fV1=0;
1888  delete esdB; esdB=NULL;
1889  delete vtxCasc;vtxCasc=NULL;
1890  trackB=NULL;
1891  delete trackV0; trackV0=NULL;
1892  if(!DStar){
1893  v0=NULL;
1894  }
1895  return kFALSE;
1896  }
1897  vtxCasc->SetParent(rCasc);
1898  rCasc->SetSecondaryVtx(vtxCasc);
1899  AddDaughterRefs(vtxCasc,(AliAODEvent*)event,twoTrackArrayCasc);
1900  if(DStar)vtxCasc->AddDaughter(trackD0);
1901  else vtxCasc->AddDaughter(v0);
1902  rCasc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1903 
1904  Bool_t refill =kTRUE;
1905 
1906  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1907  // propagate tracks to secondary vertex, to compute inv. mass
1908  esdB->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1909  trackV0->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1910  Double_t momentum[3];
1911  esdB->GetPxPyPz(momentum);
1912  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1913  trackV0->GetPxPyPz(momentum);
1914  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1915 
1916  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc,event);
1917  if(!primVertexAOD){
1918  delete fV1; fV1=0;
1919  delete vtxCasc; vtxCasc=NULL;
1920  twoTrackArrayCasc->Clear();
1921  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1922  delete esdB; esdB=NULL;
1923  delete trackV0; trackV0=NULL;
1924  if(!DStar)v0=NULL;
1925  return kFALSE;
1926  }
1927  Double_t d0z0[2],covd0z0[3];
1928  esdB->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1929  d0[0] = d0z0[0];
1930  d0err[0] = TMath::Sqrt(covd0z0[0]);
1931  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1932  d0[1] = d0z0[0];
1933  d0err[1] = TMath::Sqrt(covd0z0[0]);
1934  rCasc->SetPxPyPzProngs(2,px,py,pz);
1935  rCasc->SetDCA(dca);
1936  rCasc->Setd0Prongs(2,d0);
1937  rCasc->Setd0errProngs(2,d0err);
1938  rCasc->SetCharge(esdB->Charge());
1939  // get PID info from ESD
1940  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1941  if(esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
1942  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1943  Double_t esdpid[10];
1944  for(Int_t i=0;i<5;i++) {
1945  esdpid[i] = esdpid0[i];
1946  esdpid[5+i] = esdpid1[i];
1947  }
1948  rCasc->SetPID(2,esdpid);
1949  rCasc->SetIsFilled(2);
1950 
1951 
1952  delete fV1; fV1=0;
1953  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1954  twoTrackArrayCasc->Clear();
1955  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1956  delete esdB; esdB=NULL;
1957  delete trackV0; trackV0=NULL;
1958  if(!DStar)v0=NULL;
1959  return kTRUE;
1960 }
1961 //---------------------------------------------------------------------------
1963 {
1968 
1969 
1970  if (rc->GetIsFilled()==0) return kFALSE; // Reduced dAOD with candidates not yet refilled
1971  if (!fAODMap) MapAODtracks(event); // Fill the AOD index map if it is not yet done
1972 
1973 
1974  // - Get bachelor and V0 from the cascade
1975  AliAODTrack *trackB = dynamic_cast<AliAODTrack*> (rc->GetBachelor());
1976  if (!trackB) return kFALSE;
1977 
1978  AliAODv0 *v0 = dynamic_cast<AliAODv0*> (rc->Getv0());
1979  if (!v0) return kFALSE;
1980 
1981 
1982  // - Cast bachelor (AOD->ESD) and V0 (AliAODv0->AliNeutralTrackParam)
1983  AliESDtrack *esdB = new AliESDtrack(trackB);
1984  if (!esdB) return kFALSE;
1985 
1986  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1987  if (!trackVV0) return kFALSE;
1988 
1989  AliNeutralTrackParam *trackV0 = new AliNeutralTrackParam(trackVV0);
1990  if (!trackV0) return kFALSE;
1991 
1992  TObjArray *twoTrackArrayCasc = new TObjArray(2);
1993  twoTrackArrayCasc->AddAt(esdB, 0);
1994  twoTrackArrayCasc->AddAt(trackV0, 1);
1995 
1996 
1997  Double_t dispersion, xdummy, ydummy, pos[3], cov[6];
1998 
1999 
2000  // - Some ingredients will be needed:
2001  // magnetic field
2002  fBzkG = (Double_t) event->GetMagneticField();
2003 
2004  // primary vertex
2005  const AliVVertex *vprimary = event->GetPrimaryVertex();
2006  vprimary->GetXYZ(pos);
2007  vprimary->GetCovarianceMatrix(cov);
2008  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2009 
2010  // vertexer
2011  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
2012 
2013 
2014  // - Compute the DCA and the decay vertex
2015  Double_t dca = esdB->GetDCA(trackV0, fBzkG, xdummy, ydummy);
2016  AliAODVertex *vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc, dispersion, kFALSE);
2017 
2018  if (!vtxCasc) {
2019  twoTrackArrayCasc->Clear();
2020  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2021  delete fV1; fV1=0;
2022  delete esdB; esdB=0;
2023  delete vtxCasc; vtxCasc=0;
2024  delete trackV0; trackV0=0;
2025  trackB=0; v0=0;
2026  return kFALSE;
2027  }
2028 
2029 
2030  // - Linking the cascade with the new secondary vertex
2031  vtxCasc->SetParent(rc);
2032  rc->SetSecondaryVtx(vtxCasc);
2033  AddDaughterRefs(vtxCasc, (AliAODEvent*)event, twoTrackArrayCasc);
2034  vtxCasc->AddDaughter(v0);
2035  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
2036 
2037 
2038  // - Propagate tracks to secondary vertex, to get momenta
2039  Double_t momentum[3], px[2], py[2], pz[2], d0[2], d0err[2], d0z0[2], covd0z0[3];
2040  esdB->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2041  trackV0->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2042  esdB->GetPxPyPz(momentum);
2043  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2044  trackV0->GetPxPyPz(momentum);
2045  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2046 
2047 
2048  // - Propagate tracks to primary vertex, to get impact parameters
2049  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc, event);
2050  if (!primVertexAOD) {
2051  twoTrackArrayCasc->Clear();
2052  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2053  delete fV1; fV1=0;
2054  delete esdB; esdB=0;
2055  delete vtxCasc; vtxCasc=0;
2056  delete trackV0; trackV0=0;
2057  trackB=0; v0=0;
2058  return kFALSE;
2059  }
2060 
2061  esdB->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2062  d0[0] = d0z0[0];
2063  d0err[0] = TMath::Sqrt(covd0z0[0]);
2064  trackV0->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2065  d0[1] = d0z0[0];
2066  d0err[1] = TMath::Sqrt(covd0z0[0]);
2067 
2068 
2069  // - Get PID info from ESD
2070  Double_t esdpid[10];
2071  Double_t esdpid0[5] = {0., 0., 0., 0., 0.};
2072  Double_t esdpid1[5] = {0., 0., 0., 0., 0.};
2073  if (esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
2074  for (Int_t ipid=0; ipid<5; ipid++) {
2075  esdpid[ipid] = esdpid0[ipid];
2076  esdpid[5+ipid] = esdpid1[ipid];
2077  }
2078 
2079 
2080  // - Set data members
2081  rc->SetOwnPrimaryVtx(primVertexAOD);
2082  rc->SetDCA(dca);
2083  rc->SetPxPyPzProngs(2, px, py, pz);
2084  rc->Setd0Prongs(2, d0);
2085  rc->Setd0errProngs(2, d0err);
2086  rc->SetCharge(esdB->Charge());
2087  rc->SetPID(2, esdpid);
2088 
2089  rc->SetIsFilled(2); // To clean the newly reconstructed secondary vertex with the CleanUp task
2090 
2091 
2092  twoTrackArrayCasc->Clear();
2093  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2094  delete fV1; fV1=0;
2095  delete primVertexAOD; primVertexAOD=0;
2096  delete esdB; esdB=0;
2097  delete trackV0; trackV0=0;
2098  trackB=0; v0=0;
2099 
2100  return kTRUE;
2101 }
2102 
2103 //---------------------------------------------------------------------------
2105  TObjArray *twoTrackArray,AliVEvent *event,
2106  AliAODVertex *secVert,
2107  AliAODRecoDecayHF2Prong *rd2Prong,
2108  Double_t dca,
2109  Bool_t &okDstar)
2110 {
2113  //AliCodeTimerAuto("",0);
2114  UInt_t ntref=TProcessID::GetObjectCount();
2115  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2116  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2117  }
2118 
2119  okDstar = kFALSE;
2120 
2121  Bool_t dummy1,dummy2,dummy3;
2122 
2123  // We use Make2Prong to construct the AliAODRecoCascadeHF
2124  // (which inherits from AliAODRecoDecayHF2Prong)
2125  AliAODRecoCascadeHF *theCascade =
2126  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2127  dummy1,dummy2,dummy3);
2128  if(!theCascade) return 0x0;
2129 
2130  // charge
2131  AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2132  theCascade->SetCharge(trackPi->Charge());
2133 
2134  //--- selection cuts
2135  //
2136  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2137  if(fInputAOD){
2138  Int_t idSoftPi=(Int_t)trackPi->GetID();
2139  if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
2140  AliAODTrack* trackPiAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idSoftPi]));
2141  if(!trackPiAOD) AliFatal("Not a standard AOD");
2142  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
2143  }
2144  }else{
2145  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
2146  }
2147  tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
2148 
2149  AliAODVertex *primVertexAOD=0;
2151  // take event primary vertex
2152  primVertexAOD = PrimaryVertex();
2153  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2154  rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
2155  }
2156  // select D*->D0pi
2157  if(fDstar) {
2159  if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
2160  }
2161  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2162  tmpCascade->UnsetOwnPrimaryVtx();
2163  delete tmpCascade; tmpCascade=NULL;
2165  rd2Prong->UnsetOwnPrimaryVtx();
2166  }
2167  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2168  //---
2169 
2170 
2171  return theCascade;
2172 }
2173 
2174 
2175 //----------------------------------------------------------------------------
2177  TObjArray *twoTrackArray,AliVEvent *event,
2178  AliAODVertex *secVert,
2179  AliAODv0 *v0,
2180  Double_t dca,
2181  Bool_t &okCascades)
2182 {
2185  //AliCodeTimerAuto("",0);
2186  UInt_t ntref=TProcessID::GetObjectCount();
2187  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2188  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2189  }
2190 
2191  // preselection to reduce the number of TRefs
2192  Double_t px[2],py[2],pz[2];
2193  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2194  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2195  // propagate tracks to secondary vertex, to compute inv. mass
2196  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2197  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2198  Double_t momentum[3];
2199  postrack->GetPxPyPz(momentum);
2200  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2201  negtrack->GetPxPyPz(momentum);
2202  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2203  if(!SelectInvMassAndPtCascade(px,py,pz)) return 0x0;
2204  Double_t dummyd0[2]={0,0};
2205  Double_t dummyd0err[2]={0,0};
2206  AliAODRecoCascadeHF tmpCasc(0x0,postrack->Charge(),px, py, pz, dummyd0, dummyd0err,0.);
2207  Bool_t presel=fCutsLctoV0->PreSelect(&tmpCasc,v0,postrack);
2208  if(!presel) return 0x0;
2209 
2210  // AliDebug(2,Form(" building the cascade"));
2211  okCascades = kFALSE;
2212  Bool_t dummy1,dummy2,dummy3;
2213 
2214  // We use Make2Prong to construct the AliAODRecoCascadeHF
2215  // (which inherits from AliAODRecoDecayHF2Prong)
2216  AliAODRecoCascadeHF *theCascade =
2217  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2218  dummy1,dummy2,dummy3);
2219  if(!theCascade) return 0x0;
2220 
2221  // bachelor track and charge
2222  AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2223  theCascade->SetCharge(trackBachelor->Charge());
2224 
2225  //--- selection cuts
2226  //
2227 
2228  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2229  if(fInputAOD){
2230  Int_t idBachelor=(Int_t)trackBachelor->GetID();
2231  if (idBachelor > -1 && idBachelor < fAODMapSize) {
2232  AliAODTrack* trackBachelorAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idBachelor]));
2233  if(!trackBachelorAOD) AliFatal("Not a standard AOD");
2234  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
2235  }
2236  }else{
2237  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
2238  }
2239  tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
2240 
2241  AliAODVertex *primVertexAOD=0;
2243  // take event primary vertex
2244  primVertexAOD = PrimaryVertex();
2245  if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
2246  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2247  }
2248 
2249  // select Cascades
2250  if (fCascades && fInputAOD) {
2251  if (fCutsLctoV0->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2252  okCascades = kTRUE;
2254  }
2255  if (fCutsDplustoK0spi) {
2257  okCascades = kTRUE;
2259  }
2260  }
2261  if (fCutsDstoK0sK){
2262  if (fCutsDstoK0sK->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2263  okCascades = kTRUE;
2265  }
2266  }
2267  }
2268  else {
2269  //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
2270  okCascades=kTRUE;
2271  } // no cuts implemented from ESDs
2272  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2273  tmpCascade->UnsetOwnPrimaryVtx();
2274  delete tmpCascade; tmpCascade=NULL;
2275  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2276  //---
2277 
2278  return theCascade;
2279 }
2280 
2281 //-----------------------------------------------------------------------------
2283  TObjArray *twoTrackArray,AliVEvent *event,
2284  AliAODVertex *secVert,Double_t dca,
2285  Bool_t &okD0,Bool_t &okJPSI,
2286  Bool_t &okD0fromDstar, Bool_t refill, AliAODRecoDecayHF2Prong *rd)
2287 {
2290  // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
2291  // AliCodeTimerAuto("",0);
2292  okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
2293 
2294  UInt_t ntref=TProcessID::GetObjectCount();
2295  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2296  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2297  }
2298 
2299  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
2300  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2301  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2302 
2303  // propagate tracks to secondary vertex, to compute inv. mass
2304  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2305  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2306 
2307  Double_t momentum[3];
2308  postrack->GetPxPyPz(momentum);
2309  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2310  negtrack->GetPxPyPz(momentum);
2311  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2312 
2313  if(!refill){//skip if it is called in refill step because already checked
2314  // invariant mass cut (try to improve coding here..)
2315  Bool_t okMassCut=kFALSE;
2316  if(!okMassCut && fD0toKpi) if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
2317  if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
2318  if(!okMassCut && fDstar) if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
2319  if(!okMassCut && fCascades) if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
2320  if(!okMassCut) {
2321  //AliDebug(2," candidate didn't pass mass cut");
2322  return 0x0;
2323  }
2324  }
2325  // primary vertex to be used by this candidate
2326  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
2327  if(!primVertexAOD) return 0x0;
2328 
2329  Double_t d0z0[2],covd0z0[3];
2330  postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2331  d0[0] = d0z0[0];
2332  d0err[0] = TMath::Sqrt(covd0z0[0]);
2333  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2334  d0[1] = d0z0[0];
2335  d0err[1] = TMath::Sqrt(covd0z0[0]);
2336  AliAODRecoDecayHF2Prong *the2Prong;
2337  // create the object AliAODRecoDecayHF2Prong
2338  if(!refill){
2339  the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
2340  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2341  UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
2342  the2Prong->SetProngIDs(2,id);
2343  if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
2344  // Add daughter references already here
2345  if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2346 
2347  // select D0->Kpi
2348  if(fD0toKpi) {
2350  if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
2351  }
2352  //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
2353  // select J/psi from B
2354  if(fJPSItoEle) {
2356  }
2357  //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
2358  // select D0->Kpi from Dstar
2359  if(fDstar) {
2360  okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
2361  if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
2362  }
2363  //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
2364  }
2365  }else{
2366  the2Prong =rd;
2367  the2Prong->SetSecondaryVtx(secVert);
2368  secVert->SetParent(the2Prong);
2369  AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2370  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2371  the2Prong->SetPxPyPzProngs(2,px,py,pz);
2372  the2Prong->SetDCA(dca);
2373  the2Prong->Setd0Prongs(2,d0);
2374  the2Prong->Setd0errProngs(2,d0err);
2375  the2Prong->SetCharge(0);
2376  }
2377  delete primVertexAOD; primVertexAOD=NULL;
2378 
2379  // remove the primary vertex (was used only for selection)
2381  the2Prong->UnsetOwnPrimaryVtx();
2382  }
2383 
2384  // get PID info from ESD
2385  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2386  if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
2387  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2388  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2389  Double_t esdpid[10];
2390  for(Int_t i=0;i<5;i++) {
2391  esdpid[i] = esdpid0[i];
2392  esdpid[5+i] = esdpid1[i];
2393  }
2394  the2Prong->SetPID(2,esdpid);
2395  return the2Prong;
2396 }
2397 //----------------------------------------------------------------------------
2399  TObjArray *threeTrackArray,AliVEvent *event,
2400  AliAODVertex *secVert,Double_t dispersion,
2401  const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
2402  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2403  Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
2404 {
2407  // E.Bruna, F.Prino
2408  // AliCodeTimerAuto("",0);
2409 
2410 
2411  UInt_t ntref=TProcessID::GetObjectCount();
2412  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2413  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2414  }
2415 
2416  ok3Prong=kFALSE;
2417  if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
2418 
2419  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2420  Double_t momentum[3];
2421 
2422 
2423  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2424  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2425  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2426 
2427  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2428  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2429  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2430  postrack1->GetPxPyPz(momentum);
2431  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2432  negtrack->GetPxPyPz(momentum);
2433  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2434  postrack2->GetPxPyPz(momentum);
2435  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2436 
2437  // invariant mass cut for D+, Ds, Lc
2438  Bool_t okMassCut=kFALSE;
2439  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2440  if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
2441  if(!okMassCut) {
2442  //AliDebug(2," candidate didn't pass mass cut");
2443  return 0x0;
2444  }
2445 
2446  // primary vertex to be used by this candidate
2447  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2448  if(!primVertexAOD) return 0x0;
2449 
2450  Double_t d0z0[2],covd0z0[3];
2451  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2452  d0[0]=d0z0[0];
2453  d0err[0] = TMath::Sqrt(covd0z0[0]);
2454  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2455  d0[1]=d0z0[0];
2456  d0err[1] = TMath::Sqrt(covd0z0[0]);
2457  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2458  d0[2]=d0z0[0];
2459  d0err[2] = TMath::Sqrt(covd0z0[0]);
2460 
2461 
2462  // create the object AliAODRecoDecayHF3Prong
2463  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2464  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2465  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]));
2466  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]));
2467  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2468 
2469 
2470  // construct the candidate passing a NULL pointer for the secondary vertex to avoid creation of TRef
2471  AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(0x0,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
2472  the3Prong->SetOwnPrimaryVtx(primVertexAOD);
2473  // add a pointer to the secondary vertex via SetOwnSecondaryVtx (no TRef created)
2474  AliAODVertex* ownsecv=secVert->CloneWithoutRefs();
2475  the3Prong->SetOwnSecondaryVtx(ownsecv);
2476  UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
2477  the3Prong->SetProngIDs(3,id);
2478 
2479  delete primVertexAOD; primVertexAOD=NULL;
2480 
2481  // disable PID, which requires the TRefs to the daughter tracks
2482  fCutsDplustoKpipi->SetUsePID(kFALSE);
2483  fCutsDstoKKpi->SetUsePID(kFALSE);
2484  fCutsLctopKpi->SetUsePID(kFALSE);
2485 
2486  // select D+->Kpipi, Ds->KKpi, Lc->pKpi
2487  if(f3Prong) {
2488  ok3Prong = kFALSE;
2489 
2491  ok3Prong = kTRUE;
2493  }
2494  if(useForDs && fOKInvMassDs){
2496  ok3Prong = kTRUE;
2498  }
2499  }
2500  if(useForLc && fOKInvMassLc){
2502  ok3Prong = kTRUE;
2504  }
2505  }
2506  }
2507  //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
2508  the3Prong->UnsetOwnSecondaryVtx();
2510  the3Prong->UnsetOwnPrimaryVtx();
2511  }
2512 
2513  // Add TRefs to secondary vertex and daughter tracks only for candidates passing the filtering cuts
2514  if(ok3Prong && fInputAOD){
2515  the3Prong->SetSecondaryVtx(secVert);
2516  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2517  }
2518 
2519  // get PID info from ESD
2520  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2521  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2522  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2523  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2524  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2525  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2526 
2527  Double_t esdpid[15];
2528  for(Int_t i=0;i<5;i++) {
2529  esdpid[i] = esdpid0[i];
2530  esdpid[5+i] = esdpid1[i];
2531  esdpid[10+i] = esdpid2[i];
2532  }
2533  the3Prong->SetPID(3,esdpid);
2534 
2535  return the3Prong;
2536 }
2537 //----------------------------------------------------------------------------
2539  TObjArray *threeTrackArray,AliVEvent *event,
2540  AliAODVertex *secVert,Double_t dispersion,
2541  Double32_t dist12, Double32_t dist23,
2542  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2544 {
2545  // Fill on-the-fly the 3prong data member missing info
2546  // set TRef of vertex and daughters
2547  // do not recalculate the two-track secondary vertex dist12 and dist23
2548  // because they are stored in the AliAODRecoDecayHF3Prong candidate
2549  // reconstructed in the FindCandidate step
2550  // do not check if it is a Lambdac, Ds or Dplus because it is already check in the FindCandidate step
2551  // and the info is stored
2552  // AliCodeTimerAuto("",0);
2553 
2554  UInt_t ntref=TProcessID::GetObjectCount();
2555  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2556  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2557  }
2558 
2559  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2560  Double_t momentum[3];
2561 
2562 
2563  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2564  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2565  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2566 
2567  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2568  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2569  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2570  postrack1->GetPxPyPz(momentum);
2571  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2572  negtrack->GetPxPyPz(momentum);
2573  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2574  postrack2->GetPxPyPz(momentum);
2575  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2576  // primary vertex to be used by this candidate
2577  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2578  if(!primVertexAOD) return 0x0;
2579  Double_t d0z0[2],covd0z0[3];
2580  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2581  d0[0]=d0z0[0];
2582  d0err[0] = TMath::Sqrt(covd0z0[0]);
2583  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2584  d0[1]=d0z0[0];
2585  d0err[1] = TMath::Sqrt(covd0z0[0]);
2586  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2587  d0[2]=d0z0[0];
2588  d0err[2] = TMath::Sqrt(covd0z0[0]);
2589 
2590  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2591  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2592  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2593 
2594  rd->SetSecondaryVtx(secVert);
2595  secVert->SetParent(rd);
2596  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2597 
2598  rd->SetPxPyPzProngs(3,px,py,pz);
2599  rd->SetDCAs(3,dca);
2600  rd->Setd0Prongs(3,d0);
2601  rd->Setd0errProngs(3,d0err);
2602  rd->SetCharge(charge);
2603  rd->SetOwnPrimaryVtx(primVertexAOD);
2604  rd->SetSigmaVert(dispersion);
2605  delete primVertexAOD; primVertexAOD=NULL;
2606 
2608  rd->UnsetOwnPrimaryVtx();
2609  }
2610 
2611  // get PID info from ESD
2612  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2613  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2614  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2615  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2616  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2617  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2618 
2619  Double_t esdpid[15];
2620  for(Int_t i=0;i<5;i++) {
2621  esdpid[i] = esdpid0[i];
2622  esdpid[5+i] = esdpid1[i];
2623  esdpid[10+i] = esdpid2[i];
2624  }
2625  rd->SetPID(3,esdpid);
2626  return rd;
2627 }
2628 //----------------------------------------------------------------------------
2630  TObjArray *fourTrackArray,AliVEvent *event,
2631  AliAODVertex *secVert,
2632  const AliAODVertex *vertexp1n1,
2633  const AliAODVertex *vertexp1n1p2,
2634  Double_t dcap1n1,Double_t dcap1n2,
2635  Double_t dcap2n1,Double_t dcap2n2,
2636  Bool_t &ok4Prong)
2637 {
2640  // G.E.Bruno, R.Romita
2641  // AliCodeTimerAuto("",0);
2642 
2643  UInt_t ntref=TProcessID::GetObjectCount();
2644  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2645  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2646  }
2647 
2648  ok4Prong=kFALSE;
2649  if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
2650 
2651  Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
2652 
2653  AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
2654  AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
2655  AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
2656  AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
2657 
2658  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2659  negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2660  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2661  negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2662 
2663  Double_t momentum[3];
2664  postrack1->GetPxPyPz(momentum);
2665  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2666  negtrack1->GetPxPyPz(momentum);
2667  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2668  postrack2->GetPxPyPz(momentum);
2669  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2670  negtrack2->GetPxPyPz(momentum);
2671  px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
2672 
2673  // invariant mass cut for rho or D0 (try to improve coding here..)
2674  Bool_t okMassCut=kFALSE;
2675  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2676  if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
2677  if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
2678  }
2679  if(!okMassCut) {
2680  //if(fDebug) printf(" candidate didn't pass mass cut\n");
2681  //printf(" candidate didn't pass mass cut\n");
2682  return 0x0;
2683  }
2684 
2685  // primary vertex to be used by this candidate
2686  AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
2687  if(!primVertexAOD) return 0x0;
2688 
2689  Double_t d0z0[2],covd0z0[3];
2690  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2691  d0[0]=d0z0[0];
2692  d0err[0] = TMath::Sqrt(covd0z0[0]);
2693  negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2694  d0[1]=d0z0[0];
2695  d0err[1] = TMath::Sqrt(covd0z0[0]);
2696  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2697  d0[2]=d0z0[0];
2698  d0err[2] = TMath::Sqrt(covd0z0[0]);
2699  negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2700  d0[3]=d0z0[0];
2701  d0err[3] = TMath::Sqrt(covd0z0[0]);
2702 
2703 
2704  // create the object AliAODRecoDecayHF4Prong
2705  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2706  Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2707  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]));
2708  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]));
2709  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]));
2710  Short_t charge=0;
2711  AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2712  the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2713  UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2714  the4Prong->SetProngIDs(4,id);
2715 
2716  delete primVertexAOD; primVertexAOD=NULL;
2717 
2719 
2720 
2722  the4Prong->UnsetOwnPrimaryVtx();
2723  }
2724 
2725 
2726  // get PID info from ESD
2727  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2728  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2729  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2730  if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2731  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2732  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2733  Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2734  if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2735 
2736  Double_t esdpid[20];
2737  for(Int_t i=0;i<5;i++) {
2738  esdpid[i] = esdpid0[i];
2739  esdpid[5+i] = esdpid1[i];
2740  esdpid[10+i] = esdpid2[i];
2741  esdpid[15+i] = esdpid3[i];
2742  }
2743  the4Prong->SetPID(4,esdpid);
2744 
2745  return the4Prong;
2746 }
2747 //----------------------------------------------------------------------------------
2749  //assign and save in fAODMap the index of the AliAODTrack track
2750  //ordering them on the basis of selected criteria
2751 
2752  fAODMapSize = 100000;
2753  fAODMap = new Int_t[fAODMapSize];
2754  AliAODTrack *track=0;
2755  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2756  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
2757  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
2758  if(!track) AliFatal("Not a standard AOD");
2759  // skip pure ITS SA tracks
2760  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2761 
2762  // skip tracks without ITS
2763  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2764 
2765  // TEMPORARY: check that the cov matrix is there
2766  Double_t covtest[21];
2767  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2768  //
2769 
2770  Int_t ind = (Int_t)track->GetID();
2771  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2772  }
2773  return;
2774 }
2775 //-----------------------------------------------------------------------------
2776 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2777  AliVEvent *event) const
2778 {
2780  //AliCodeTimerAuto("",0);
2781 
2782  AliESDVertex *vertexESD = 0;
2783  AliAODVertex *vertexAOD = 0;
2784 
2785 
2787  // primary vertex from the input event
2788 
2789  vertexESD = new AliESDVertex(*fV1);
2790 
2791  } else {
2792  // primary vertex specific to this candidate
2793 
2794  Int_t nTrks = trkArray->GetEntriesFast();
2795  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2796 
2798  // recalculating the vertex
2799 
2800  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2801  Float_t diamondcovxy[3];
2802  event->GetDiamondCovXY(diamondcovxy);
2803  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2804  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2805  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2806  vertexer->SetVtxStart(diamond);
2807  delete diamond; diamond=NULL;
2808  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2809  vertexer->SetOnlyFitter();
2810  }
2811  Int_t skipped[1000];
2812  Int_t nTrksToSkip=0,id;
2813  AliExternalTrackParam *t = 0;
2814  for(Int_t i=0; i<nTrks; i++) {
2815  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2816  id = (Int_t)t->GetID();
2817  if(id<0) continue;
2818  skipped[nTrksToSkip++] = id;
2819  }
2820  // TEMPORARY FIX
2821  // For AOD, skip also tracks without covariance matrix
2822  if(fInputAOD) {
2823  Double_t covtest[21];
2824  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2825  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2826  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2827  id = (Int_t)vtrack->GetID();
2828  if(id<0) continue;
2829  skipped[nTrksToSkip++] = id;
2830  }
2831  }
2832  }
2833  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2834  //
2835  vertexer->SetSkipTracks(nTrksToSkip,skipped);
2836  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2837 
2838  } else if(fRmTrksFromPrimVtx && nTrks>0) {
2839  // removing the prongs tracks
2840 
2841  TObjArray rmArray(nTrks);
2842  UShort_t *rmId = new UShort_t[nTrks];
2843  AliESDtrack *esdTrack = 0;
2844  AliESDtrack *t = 0;
2845  for(Int_t i=0; i<nTrks; i++) {
2846  t = (AliESDtrack*)trkArray->UncheckedAt(i);
2847  esdTrack = new AliESDtrack(*t);
2848  rmArray.AddLast(esdTrack);
2849  if(esdTrack->GetID()>=0) {
2850  rmId[i]=(UShort_t)esdTrack->GetID();
2851  } else {
2852  rmId[i]=9999;
2853  }
2854  }
2855  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2856  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2857  delete [] rmId; rmId=NULL;
2858  rmArray.Delete();
2859 
2860  }
2861 
2862  if(!vertexESD) return vertexAOD;
2863  if(vertexESD->GetNContributors()<=0) {
2864  //AliDebug(2,"vertexing failed");
2865  delete vertexESD; vertexESD=NULL;
2866  return vertexAOD;
2867  }
2868 
2869  delete vertexer; vertexer=NULL;
2870 
2871  }
2872 
2873  // convert to AliAODVertex
2874  Double_t pos[3],cov[6],chi2perNDF;
2875  vertexESD->GetXYZ(pos); // position
2876  vertexESD->GetCovMatrix(cov); //covariance matrix
2877  chi2perNDF = vertexESD->GetChi2toNDF();
2878  delete vertexESD; vertexESD=NULL;
2879 
2880  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2881 
2882  return vertexAOD;
2883 }
2884 //-----------------------------------------------------------------------------
2887 
2888  //printf("Preselections:\n");
2889  // fTrackFilter->Dump();
2890  if(fSecVtxWithKF) {
2891  printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2892  } else {
2893  printf("Secondary vertex with AliVertexerTracks\n");
2894  }
2895  if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2896  if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2897  if(fD0toKpi) {
2898  printf("Reconstruct D0->Kpi candidates with cuts:\n");
2900  }
2901  if(fDstar) {
2902  printf("Reconstruct D*->D0pi candidates with cuts:\n");
2903  if(fFindVertexForDstar) {
2904  printf(" Reconstruct a secondary vertex for the D*\n");
2905  } else {
2906  printf(" Assume the D* comes from the primary vertex\n");
2907  }
2909  }
2910  if(fJPSItoEle) {
2911  printf("Reconstruct J/psi from B candidates with cuts:\n");
2913  }
2914  if(f3Prong) {
2915  printf("Reconstruct 3 prong candidates.\n");
2916  printf(" D+->Kpipi cuts:\n");
2918  printf(" Ds->KKpi cuts:\n");
2920  printf(" Lc->pKpi cuts:\n");
2922  }
2923  if(f4Prong) {
2924  printf("Reconstruct 4 prong candidates.\n");
2925  printf(" D0->Kpipipi cuts:\n");
2927  }
2928  if(fCascades) {
2929  printf("Reconstruct cascade candidates formed with v0s.\n");
2930  printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
2932  printf(" D+ -> K0s pi cuts:\n");
2934  printf(" Ds -> K0s K cuts:\n");
2936  }
2937 
2938  return;
2939 }
2940 //-----------------------------------------------------------------------------
2942  Double_t &dispersion,Bool_t useTRefArray) const
2943 {
2945  //AliCodeTimerAuto("",0);
2946 
2947  AliESDVertex *vertexESD = 0;
2948  AliAODVertex *vertexAOD = 0;
2949 
2950  if(!fSecVtxWithKF) { // AliVertexerTracks
2951 
2952  fVertexerTracks->SetVtxStart(fV1);
2953  vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2954 
2955  if(!vertexESD) return vertexAOD;
2956 
2957  if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2958  //AliDebug(2,"vertexing failed");
2959  delete vertexESD; vertexESD=NULL;
2960  return vertexAOD;
2961  }
2962 
2963  Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
2964  if(vertRadius2>8.){
2965  // vertex outside beam pipe, reject candidate to avoid propagation through material
2966  delete vertexESD; vertexESD=NULL;
2967  return vertexAOD;
2968  }
2969 
2970  } else { // Kalman Filter vertexer (AliKFParticle)
2971 
2972  AliKFParticle::SetField(fBzkG);
2973 
2974  AliKFVertex vertexKF;
2975 
2976  Int_t nTrks = trkArray->GetEntriesFast();
2977  for(Int_t i=0; i<nTrks; i++) {
2978  AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2979  AliKFParticle daughterKF(*esdTrack,211);
2980  vertexKF.AddDaughter(daughterKF);
2981  }
2982  vertexESD = new AliESDVertex(vertexKF.Parameters(),
2983  vertexKF.CovarianceMatrix(),
2984  vertexKF.GetChi2(),
2985  vertexKF.GetNContributors());
2986 
2987  }
2988 
2989  // convert to AliAODVertex
2990  Double_t pos[3],cov[6],chi2perNDF;
2991  vertexESD->GetXYZ(pos); // position
2992  vertexESD->GetCovMatrix(cov); //covariance matrix
2993  chi2perNDF = vertexESD->GetChi2toNDF();
2994  dispersion = vertexESD->GetDispersion();
2995  delete vertexESD; vertexESD=NULL;
2996 
2997  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2998  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2999 
3000  return vertexAOD;
3001 }
3002 //-----------------------------------------------------------------------------
3005  //AliCodeTimerAuto("",0);
3006 
3007  Int_t retval=kFALSE;
3008  Double_t momentum[3];
3009  Double_t px[3],py[3],pz[3];
3010  for(Int_t iTrack=0; iTrack<3; iTrack++){
3011  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3012  track->GetPxPyPz(momentum);
3013  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3014  }
3015  retval = SelectInvMassAndPt3prong(px,py,pz);
3016 
3017  return retval;
3018 }
3019 
3020 //-----------------------------------------------------------------------------
3023  //AliCodeTimerAuto("",0);
3024 
3025  Int_t retval=kFALSE;
3026  Double_t momentum[3];
3027  Double_t px[4],py[4],pz[4];
3028 
3029  for(Int_t iTrack=0; iTrack<4; iTrack++){
3030  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3031  track->GetPxPyPz(momentum);
3032  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3033  }
3034 
3035  retval = SelectInvMassAndPt4prong(px,py,pz);
3036 
3037  return retval;
3038 }
3039 //-----------------------------------------------------------------------------
3042  //AliCodeTimerAuto("",0);
3043 
3044  Int_t retval=kFALSE;
3045  Double_t momentum[3];
3046  Double_t px[2],py[2],pz[2];
3047 
3048  for(Int_t iTrack=0; iTrack<2; iTrack++){
3049  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3050  track->GetPxPyPz(momentum);
3051  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3052  }
3053  retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
3054 
3055  return retval;
3056 }
3057 //-----------------------------------------------------------------------------
3059  Double_t *py,
3060  Double_t *pz){
3062  //AliCodeTimerAuto("",0);
3063 
3064  UInt_t pdg2[2];
3065  Int_t nprongs=2;
3066  Double_t minv2,mrange;
3067  Double_t lolim,hilim;
3068  Double_t minPt=0;
3069  Bool_t retval=kFALSE;
3070 
3071  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3072  fOKInvMassD0=kFALSE;
3073  // pt cut
3075  if(minPt>0.1)
3076  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3077  // mass cut
3078  mrange=fCutsD0toKpi->GetMassCut();
3079  lolim=fMassDzero-mrange;
3080  hilim=fMassDzero+mrange;
3081  pdg2[0]=211; pdg2[1]=321;
3082  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3083  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3084  retval=kTRUE;
3085  fOKInvMassD0=kTRUE;
3086  }
3087  pdg2[0]=321; pdg2[1]=211;
3088  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3089  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3090  retval=kTRUE;
3091  fOKInvMassD0=kTRUE;
3092  }
3093  return retval;
3094 }
3095 
3096 //-----------------------------------------------------------------------------
3098  Double_t *py,
3099  Double_t *pz){
3101  //AliCodeTimerAuto("",0);
3102 
3103  UInt_t pdg2[2];
3104  Int_t nprongs=2;
3105  Double_t minv2,mrange;
3106  Double_t lolim,hilim;
3107  Double_t minPt=0;
3108  Bool_t retval=kFALSE;
3109 
3110  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3111  fOKInvMassJpsi=kFALSE;
3112  // pt cut
3114  if(minPt>0.1)
3115  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3116  // mass cut
3117  mrange=fCutsJpsitoee->GetMassCut();
3118  lolim=fMassJpsi-mrange;
3119  hilim=fMassJpsi+mrange;
3120 
3121  pdg2[0]=11; pdg2[1]=11;
3122  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3123  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3124  retval=kTRUE;
3125  fOKInvMassJpsi=kTRUE;
3126  }
3127 
3128  return retval;
3129 }
3130 //-----------------------------------------------------------------------------
3132  Double_t *py,
3133  Double_t *pz,
3134  Int_t pidLcStatus){
3136  //AliCodeTimerAuto("",0);
3137 
3138  UInt_t pdg3[3];
3139  Int_t nprongs=3;
3140  Double_t minv2,mrange;
3141  Double_t lolim,hilim;
3142  Double_t minPt=0;
3143  Bool_t retval=kFALSE;
3144 
3145 
3146  fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
3147  fOKInvMassDplus=kFALSE;
3148  fOKInvMassDs=kFALSE;
3149  fOKInvMassLc=kFALSE;
3150  // pt cut
3152  minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
3153  if(minPt>0.1)
3154  if(fMassCalc3->Pt2() < minPt*minPt) return retval;
3155  // D+->Kpipi
3156  mrange=fCutsDplustoKpipi->GetMassCut();
3157  lolim=fMassDplus-mrange;
3158  hilim=fMassDplus+mrange;
3159  pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
3160  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3161  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3162  retval=kTRUE;
3163  fOKInvMassDplus=kTRUE;
3164  }
3165  // Ds+->KKpi
3166  mrange=fCutsDstoKKpi->GetMassCut();
3167  lolim=fMassDs-mrange;
3168  hilim=fMassDs+mrange;
3169  pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
3170  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3171  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3172  retval=kTRUE;
3173  fOKInvMassDs=kTRUE;
3174  }
3175  pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
3176  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3177  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3178  retval=kTRUE;
3179  fOKInvMassDs=kTRUE;
3180  }
3181  // Lc->pKpi
3182  mrange=fCutsLctopKpi->GetMassCut();
3183  lolim=fMassLambdaC-mrange;
3184  hilim=fMassLambdaC+mrange;
3185  if(pidLcStatus&1){
3186  pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
3187  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3188  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3189  retval=kTRUE;
3190  fOKInvMassLc=kTRUE;
3191  }
3192  }
3193  if(pidLcStatus&2){
3194  pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
3195  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3196  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3197  retval=kTRUE;
3198  fOKInvMassLc=kTRUE;
3199  }
3200  }
3201 
3202  return retval;
3203 }
3204 
3205 //-----------------------------------------------------------------------------
3207  Double_t *py,
3208  Double_t *pz){
3210  //AliCodeTimerAuto("",0);
3211 
3212  UInt_t pdg2[2];
3213  Int_t nprongs=2;
3214  Double_t minv2,mrange;
3215  Double_t lolim,hilim;
3216  Double_t minPt=0;
3217  Bool_t retval=kFALSE;
3218 
3219  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3220  fOKInvMassDstar=kFALSE;
3221  // pt cut
3223  if(minPt>0.1)
3224  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3225  // mass cut
3226  pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
3227  mrange=fCutsDStartoKpipi->GetMassCut();
3228  lolim=fMassDstar-mrange;
3229  hilim=fMassDstar+mrange;
3230  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3231  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3232  retval=kTRUE;
3233  fOKInvMassDstar=kTRUE;
3234  }
3235 
3236  return retval;
3237 }
3238 
3239 //-----------------------------------------------------------------------------
3241  Double_t *py,
3242  Double_t *pz){
3244  //AliCodeTimerAuto("",0);
3245 
3246  UInt_t pdg4[4];
3247  Int_t nprongs=4;
3248  Double_t minv2,mrange;
3249  Double_t lolim,hilim;
3250  Double_t minPt=0;
3251  Bool_t retval=kFALSE;
3252 
3253  // D0->Kpipipi without PID
3254  fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
3255  fOKInvMassD0to4p=kFALSE;
3256  // pt cut
3258  if(minPt>0.1)
3259  if(fMassCalc4->Pt2() < minPt*minPt) return retval;
3260  // mass cut
3261  mrange=fCutsD0toKpipipi->GetMassCut();
3262  lolim=fMassDzero-mrange;
3263  hilim=fMassDzero+mrange;
3264 
3265  pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
3266  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3267  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3268  retval=kTRUE;
3269  fOKInvMassD0to4p=kTRUE;
3270  }
3271 
3272  pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
3273  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3274  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3275  retval=kTRUE;
3276  fOKInvMassD0to4p=kTRUE;
3277  }
3278 
3279  pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
3280  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3281  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3282  retval=kTRUE;
3283  fOKInvMassD0to4p=kTRUE;
3284  }
3285 
3286  pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
3287  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3288  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3289  retval=kTRUE;
3290  fOKInvMassD0to4p=kTRUE;
3291  }
3292 
3293  return retval;
3294 }
3295 //-----------------------------------------------------------------------------
3297  Double_t *py,
3298  Double_t *pz){
3300  //AliCodeTimerAuto("",0);
3301 
3302  UInt_t pdg2[2];
3303  Int_t nprongs=2;
3304  Double_t minv2,mrange;
3305  Double_t lolim,hilim;
3306  Double_t minPt=0;
3307  Bool_t retval=kFALSE;
3308 
3309  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3310  minPt=fCutsLctoV0->GetMinPtCandidate();
3311  if(minPt>0.1)
3312  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3313 
3314  fOKInvMassLctoV0 = kFALSE;
3315 
3316  // LambdaC candidate
3317  if (fCutsLctoV0) {
3318  mrange = fCutsLctoV0->GetMassCut();
3319  lolim = fMassLambdaC - mrange;
3320  hilim = fMassLambdaC + mrange;
3321  pdg2[0] = 2212; pdg2[1] = 310;
3322  minv2 = fMassCalc2->InvMass2(2, pdg2);
3323  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3324  retval = kTRUE;
3325  fOKInvMassLctoV0 = kTRUE;
3326  }
3327  pdg2[0] = 211; pdg2[1] = 3122;
3328  minv2 = fMassCalc2->InvMass2(2, pdg2);
3329  if ((minv2>lolim*lolim) && (minv2<hilim*hilim)) {
3330  retval = kTRUE;
3331  fOKInvMassLctoV0 = kTRUE;
3332  }
3333  }
3334 
3335  // D+ cascade candidate
3336  if (fCutsDplustoK0spi) {
3337  mrange = fCutsDplustoK0spi->GetMassCut();
3338  lolim = fMassDplus - mrange;
3339  hilim = fMassDplus + mrange;
3340  pdg2[0] = 211; pdg2[1] = 310;
3341  minv2 = fMassCalc2->InvMass2(2, pdg2);
3342  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3343  retval = kTRUE;
3344  }
3345  }
3346 
3347  // Ds cascade candidate
3348  if (fCutsDstoK0sK) {
3349  mrange = fCutsDstoK0sK->GetMassCut();
3350  lolim = fMassDs - mrange;
3351  hilim = fMassDs + mrange;
3352  pdg2[0] = 321; pdg2[1] = 310;
3353  minv2 = fMassCalc2->InvMass2(2, pdg2);
3354  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3355  retval = kTRUE;
3356  }
3357  }
3358 
3359  return retval;
3360 }
3361 //-----------------------------------------------------------------------------
3363  Int_t trkEntries,
3364  TObjArray &seleTrksArray,
3365  TObjArray &tracksAtVertex,
3366  Int_t &nSeleTrks,
3367  UChar_t *seleFlags,Int_t *evtNumber)
3368 {
3374  //AliCodeTimerAuto("",0);
3375 
3376  const AliVVertex *vprimary = event->GetPrimaryVertex();
3377 
3378  if(fV1) { delete fV1; fV1=NULL; }
3379  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
3380 
3381  Int_t nindices=0;
3382  UShort_t *indices = 0;
3383  Double_t pos[3],cov[6];
3384  const Int_t entries = event->GetNumberOfTracks();
3385  AliCentrality* cent;
3386 
3387  if(!fInputAOD) { // ESD
3388  fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
3389  cent=((AliESDEvent*)event)->GetCentrality();
3390  } else { // AOD
3391  vprimary->GetXYZ(pos);
3392  vprimary->GetCovarianceMatrix(cov);
3393  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
3394  if(entries<=0) return;
3395  indices = new UShort_t[entries];
3396  memset(indices,0,sizeof(UShort_t)*entries);
3397  fAODMapSize = 100000;
3398  fAODMap = new Int_t[fAODMapSize];
3399  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
3400  cent=((AliAODEvent*)event)->GetCentrality();
3401  }
3402  Float_t centperc=0.1;
3403  if(event->GetRunNumber()<244824){
3404  centperc=cent->GetCentralityPercentile("V0M");
3405  }else{
3406  AliMultSelection *multSelection = (AliMultSelection * ) event->FindListObject("MultSelection");
3407  if(multSelection){
3408  centperc=multSelection->GetMultiplicityPercentile("V0M");
3409  Int_t qual = multSelection->GetEvSelCode();
3410  if(qual == 199 ) centperc=0.1; // use central cuts as default
3411  }
3412  }
3413  Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE,okBachelor=kFALSE;
3414  nSeleTrks=0;
3415 
3416  // transfer ITS tracks from event to arrays
3417  for(Int_t i=0; i<entries; i++) {
3418  AliVTrack *track;
3419  track = (AliVTrack*)event->GetTrack(i);
3420 
3421  // skip pure ITS SA tracks
3422  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
3423 
3424  // skip tracks without ITS
3425  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
3426 
3427  // skip tracks with negative ID
3428  // (these are duplicated TPC-only AOD tracks, for jet analysis...)
3429  if(track->GetID()<0) continue;
3430 
3431  // TEMPORARY: check that the cov matrix is there
3432  Double_t covtest[21];
3433  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
3434  //
3435 
3436  if(fInputAOD) {
3437  AliAODTrack *aodt = (AliAODTrack*)track;
3438  if(aodt->GetUsedForPrimVtxFit()) {
3439  indices[nindices]=aodt->GetID(); nindices++;
3440  }
3441  Int_t ind = (Int_t)aodt->GetID();
3442  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
3443  }
3444 
3445  AliESDtrack *esdt = 0;
3446 
3447  if(!fInputAOD) {
3448  esdt = (AliESDtrack*)track;
3449  } else {
3450  esdt = new AliESDtrack(track);
3451  }
3452 
3453  // single track selection
3454  okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE; okBachelor=kFALSE;
3455  if(fMixEvent && i<trkEntries){
3456  evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
3457  const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
3458  Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
3459  eventVtx->GetXYZ(vtxPos);
3460  vprimary->GetXYZ(primPos);
3461  eventVtx->GetCovarianceMatrix(primCov);
3462  for(Int_t ind=0;ind<3;ind++){
3463  trasl[ind]=vtxPos[ind]-primPos[ind];
3464  }
3465 
3466  Bool_t isTransl=esdt->Translate(trasl,primCov);
3467  if(!isTransl) {
3468  delete esdt;
3469  esdt = NULL;
3470  continue;
3471  }
3472  }
3473 
3474  if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong,okBachelor) && nSeleTrks<trkEntries) {
3475  esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
3476  seleTrksArray.AddLast(esdt);
3477  tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
3478  seleFlags[nSeleTrks]=0;
3479  if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
3480  if(okFor3Prong) SETBIT(seleFlags[nSeleTrks],kBit3Prong);
3481  if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
3482  if(okBachelor) SETBIT(seleFlags[nSeleTrks],kBitBachelor);
3483  // Check the PID
3484  SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
3485  SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3486  SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3487  Bool_t useTPC=kTRUE;
3488  if(fUseTOFPID){
3489  Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
3490  if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
3491  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3492  }
3493  Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
3494  if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
3495  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3496  }
3497  Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
3498  if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
3499  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3500  }
3501  if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
3502  }
3503  if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
3504  Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
3505  if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
3506  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3507  }
3508  Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
3509  if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
3510  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3511  }
3512  Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
3513  if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
3514  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3515  }
3516  }
3517  nSeleTrks++;
3518  } else {
3519  if(fInputAOD) delete esdt;
3520  esdt = NULL;
3521  continue;
3522  }
3523 
3524  } // end loop on tracks
3525 
3526  // primary vertex from AOD
3527  if(fInputAOD) {
3528  delete fV1; fV1=NULL;
3529  vprimary->GetXYZ(pos);
3530  vprimary->GetCovarianceMatrix(cov);
3531  Double_t chi2toNDF = vprimary->GetChi2perNDF();
3532  Int_t ncontr=nindices;
3533  if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
3534  Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
3535  fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
3536  fV1->SetTitle(vprimary->GetTitle());
3537  fV1->SetIndices(nindices,indices);
3538  }
3539  if(indices) { delete [] indices; indices=NULL; }
3540 
3541 
3542  return;
3543 }
3544 //-----------------------------------------------------------------------------
3546  //
3548  //
3549  if(fUsePidTag && cuts->GetPidHF()) {
3550  Bool_t usepid=cuts->GetIsUsePID();
3551  cuts->SetUsePID(kTRUE);
3552  if(cuts->IsSelectedPID(rd))
3553  rd->SetSelectionBit(bit);
3554  cuts->SetUsePID(usepid);
3555  }
3556  return;
3557 }
3558 //-----------------------------------------------------------------------------
3560  Float_t centralityperc,
3561  Bool_t &okDisplaced,
3562  Bool_t &okSoftPi,
3563  Bool_t &okFor3Prong,
3564  Bool_t &okBachelor) const
3565 {
3567 
3568  // this is needed to store the impact parameters
3569  //AliCodeTimerAuto("",0);
3570 
3571  trk->RelateToVertex(fV1,fBzkG,kVeryBig);
3572 
3573  UInt_t selectInfo;
3574  //
3575  // Track selection, displaced tracks -- 2 prong
3576  selectInfo = 0;
3577  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3579  // central PbPb events, tighter cuts
3580  selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
3581  }else{
3582  // standard cuts
3583  if(fTrackFilter) {
3584  selectInfo = fTrackFilter->IsSelected(trk);
3585  }
3586  }
3587  if(selectInfo) okDisplaced=kTRUE;
3588 
3589  // Track selection, displaced tracks -- 3 prong
3590  selectInfo = 0;
3591  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3593  // central PbPb events, tighter cuts
3594  selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
3595  }else{
3596  // standard cuts
3597  if(fTrackFilter) {
3598  selectInfo = fTrackFilter->IsSelected(trk);
3599  }
3600  }
3601  if(selectInfo) okFor3Prong=kTRUE;
3602 
3603  // Track selection, soft pions
3604  selectInfo = 0;
3605  if(fDstar && fTrackFilterSoftPi) {
3606  selectInfo = fTrackFilterSoftPi->IsSelected(trk);
3607  }
3608  if(selectInfo) okSoftPi=kTRUE;
3609 
3610  // Track selection, bachelor
3611  selectInfo = 0;
3612  if(fCascades){
3614  selectInfo = fTrackFilterBachelor->IsSelected(trk);
3615  }else{
3616  if(okDisplaced) selectInfo=1;
3617  }
3618  }
3619  if(selectInfo) okBachelor=kTRUE;
3620 
3621  if(okDisplaced || okSoftPi || okFor3Prong || okBachelor) return kTRUE;
3622 
3623  return kFALSE;
3624 }
3625 
3626 
3627 //-----------------------------------------------------------------------------
3628 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
3634  //
3635  //AliCodeTimerAuto("",0);
3636  Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
3637  AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
3638 
3639  // create the v0 neutral track to compute the DCA to the primary vertex
3640  Double_t xyz[3], pxpypz[3];
3641  esdV0->XvYvZv(xyz);
3642  esdV0->PxPyPz(pxpypz);
3643  Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
3644  AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
3645  if(!trackesdV0) {
3646  delete vertexV0;
3647  return 0;
3648  }
3649  Double_t d0z0[2],covd0z0[3];
3650  AliAODVertex *primVertexAOD = PrimaryVertex();
3651  trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3652  Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
3653  // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
3654  Double_t dcaV0DaughterToPrimVertex[2];
3655  AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
3656  AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
3657  if( !posV0track || !negV0track) {
3658  if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
3659  delete vertexV0;
3660  delete primVertexAOD;
3661  return 0;
3662  }
3663  posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3664  // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
3665  // else
3666  dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
3667  negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3668  // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
3669  // else
3670  dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
3671  Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
3672  Double_t pmom[3],nmom[3];
3673  esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
3674  esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
3675 
3676  AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
3677  aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
3678 
3679  delete trackesdV0;
3680  delete primVertexAOD;
3681 
3682  return aodV0;
3683 }
3684 //-----------------------------------------------------------------------------
3685 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
3687  //AliCodeTimerAuto("",0);
3688 
3689  const Double_t *par=extpar->GetParameter();
3690  const Double_t *cov=extpar->GetCovariance();
3691  Double_t alpha=extpar->GetAlpha();
3692  Double_t x=extpar->GetX();
3693  esdt->Set(x,alpha,par,cov);
3694  return;
3695 }
3696 //-----------------------------------------------------------------------------
3699 
3700  fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
3701  fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
3702  fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
3703  fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
3704  fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
3705  fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
3706 }
3707 //-----------------------------------------------------------------------------
3709  //
3711  //
3712 
3713 
3714  //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
3715 
3716 
3718  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
3719  return kFALSE;
3720  }
3721 
3723  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHfCutsDplustoK0spi (%d)\n",fV0TypeForCascadeVertex,fCutsDplustoK0spi->GetV0Type());
3724  return kFALSE;
3725  }
3726 
3728  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsDstoK0sK (%d)\n",fV0TypeForCascadeVertex,fCutsDstoK0sK->GetV0Type());
3729  return kFALSE;
3730  }
3731 
3732  return kTRUE;
3733 }
3734 //-----------------------------------------------------------------------------
3735 
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
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:283
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:232
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:205
virtual void PrintAll() const
Float_t GetDCACut(Int_t iPtBin=0) const
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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:255
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:269
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.