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