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