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