Possibility to apply a cut on kaon PID in the 3 prong filtering
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisVertexingHF.cxx
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 "AliLog.h"
35 #include "AliVEvent.h"
36 #include "AliVVertex.h"
37 #include "AliVTrack.h"
38 #include "AliVertexerTracks.h"
39 #include "AliKFVertex.h"
40 #include "AliESDEvent.h"
41 #include "AliESDVertex.h"
42 #include "AliExternalTrackParam.h"
43 #include "AliNeutralTrackParam.h"
44 #include "AliESDtrack.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliAODEvent.h"
47 #include "AliPIDResponse.h"
48 #include "AliAODRecoDecay.h"
49 #include "AliAODRecoDecayHF.h"
50 #include "AliAODRecoDecayHF2Prong.h"
51 #include "AliAODRecoDecayHF3Prong.h"
52 #include "AliAODRecoDecayHF4Prong.h"
53 #include "AliAODRecoCascadeHF.h"
54 #include "AliRDHFCutsD0toKpi.h"
55 #include "AliRDHFCutsJpsitoee.h"
56 #include "AliRDHFCutsDplustoKpipi.h"
57 #include "AliRDHFCutsDstoKKpi.h"
58 #include "AliRDHFCutsLctopKpi.h"
59 #include "AliRDHFCutsLctoV0.h"
60 #include "AliRDHFCutsD0toKpipipi.h"
61 #include "AliRDHFCutsDStartoKpipi.h"
62 #include "AliAnalysisFilter.h"
63 #include "AliAnalysisVertexingHF.h"
64 #include "AliMixedEvent.h"
65 #include "AliESDv0.h"
66 #include "AliAODv0.h"
67 #include "AliCodeTimer.h"
68 #include <cstring>
69
70 ClassImp(AliAnalysisVertexingHF)
71
72 //----------------------------------------------------------------------------
73 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
74 fInputAOD(kFALSE),
75 fAODMapSize(0),
76 fAODMap(0),
77 fVertexerTracks(0x0),
78 fBzkG(0.),
79 fSecVtxWithKF(kFALSE),
80 fRecoPrimVtxSkippingTrks(kFALSE),
81 fRmTrksFromPrimVtx(kFALSE),
82 fV1(0x0),
83 fD0toKpi(kTRUE),
84 fJPSItoEle(kTRUE),
85 f3Prong(kTRUE),
86 f4Prong(kTRUE),
87 fDstar(kTRUE),
88 fCascades(kTRUE),
89 fLikeSign(kFALSE),
90 fLikeSign3prong(kFALSE),
91 fMixEvent(kFALSE),
92 fPidResponse(0x0),
93 fUseKaonPIDfor3Prong(kFALSE),
94 fnSigmaTOFKaonLow(5.),
95 fnSigmaTOFKaonHi(5.),
96 fMaxCentPercentileForTightCuts(-9999),
97 fTrackFilter(0x0),
98 fTrackFilter2prongCentral(0x0),
99 fTrackFilter3prongCentral(0x0),
100 fTrackFilterSoftPi(0x0),
101 fCutsD0toKpi(0x0),
102 fCutsJpsitoee(0x0),
103 fCutsDplustoKpipi(0x0),
104 fCutsDstoKKpi(0x0),
105 fCutsLctopKpi(0x0),
106 fCutsLctoV0(0x0),
107 fCutsD0toKpipipi(0x0),
108 fCutsDStartoKpipi(0x0),
109 fListOfCuts(0x0),
110 fFindVertexForDstar(kTRUE),
111 fFindVertexForCascades(kTRUE),
112 fMassCutBeforeVertexing(kFALSE),
113 fMassCalc2(0),
114 fMassCalc3(0),
115 fMassCalc4(0),
116 fOKInvMassD0(kFALSE),
117 fOKInvMassJpsi(kFALSE),
118 fOKInvMassDplus(kFALSE),
119 fOKInvMassDs(kFALSE),
120 fOKInvMassLc(kFALSE),
121 fOKInvMassDstar(kFALSE),
122 fOKInvMassD0to4p(kFALSE),
123 fOKInvMassLctoV0(kFALSE),
124 fnTrksTotal(0),
125 fnSeleTrksTotal(0),
126 fMassDzero(0.),
127 fMassDplus(0.),
128 fMassDs(0.),
129 fMassLambdaC(0.),
130 fMassDstar(0.),
131 fMassJpsi(0.)
132 {
133   // Default constructor
134
135   Double_t d02[2]={0.,0.};
136   Double_t d03[3]={0.,0.,0.};
137   Double_t d04[4]={0.,0.,0.,0.};
138   fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
139   fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
140   fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
141   SetMasses();
142 }
143 //--------------------------------------------------------------------------
144 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) : 
145 TNamed(source),
146 fInputAOD(source.fInputAOD),
147 fAODMapSize(source.fAODMapSize),
148 fAODMap(source.fAODMap),
149 fVertexerTracks(source.fVertexerTracks),
150 fBzkG(source.fBzkG),
151 fSecVtxWithKF(source.fSecVtxWithKF),
152 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
153 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
154 fV1(source.fV1),
155 fD0toKpi(source.fD0toKpi),
156 fJPSItoEle(source.fJPSItoEle),
157 f3Prong(source.f3Prong),
158 f4Prong(source.f4Prong),
159 fDstar(source.fDstar),
160 fCascades(source.fCascades),
161 fLikeSign(source.fLikeSign),
162 fLikeSign3prong(source.fLikeSign3prong),
163 fMixEvent(source.fMixEvent),
164 fPidResponse(source.fPidResponse),
165 fUseKaonPIDfor3Prong(source.fUseKaonPIDfor3Prong),
166 fnSigmaTOFKaonLow(source.fnSigmaTOFKaonLow),
167 fnSigmaTOFKaonHi(source.fnSigmaTOFKaonHi),
168 fMaxCentPercentileForTightCuts(source.fMaxCentPercentileForTightCuts),
169 fTrackFilter(source.fTrackFilter),
170 fTrackFilter2prongCentral(source.fTrackFilter2prongCentral),
171 fTrackFilter3prongCentral(source.fTrackFilter3prongCentral),
172 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
173 fCutsD0toKpi(source.fCutsD0toKpi),
174 fCutsJpsitoee(source.fCutsJpsitoee),
175 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
176 fCutsDstoKKpi(source.fCutsDstoKKpi),
177 fCutsLctopKpi(source.fCutsLctopKpi),
178 fCutsLctoV0(source.fCutsLctoV0),
179 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
180 fCutsDStartoKpipi(source.fCutsDStartoKpipi),
181 fListOfCuts(source.fListOfCuts),
182 fFindVertexForDstar(source.fFindVertexForDstar),
183 fFindVertexForCascades(source.fFindVertexForCascades),
184 fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
185 fMassCalc2(source.fMassCalc2),
186 fMassCalc3(source.fMassCalc3),
187 fMassCalc4(source.fMassCalc4),
188 fOKInvMassD0(source.fOKInvMassD0),
189 fOKInvMassJpsi(source.fOKInvMassJpsi),
190 fOKInvMassDplus(source.fOKInvMassDplus),
191 fOKInvMassDs(source.fOKInvMassDs),
192 fOKInvMassLc(source.fOKInvMassLc),
193 fOKInvMassDstar(source.fOKInvMassDstar),
194 fOKInvMassD0to4p(source.fOKInvMassD0to4p),
195 fOKInvMassLctoV0(source.fOKInvMassLctoV0),
196 fnTrksTotal(0),
197 fnSeleTrksTotal(0),
198 fMassDzero(source.fMassDzero),
199 fMassDplus(source.fMassDplus),
200 fMassDs(source.fMassDs),
201 fMassLambdaC(source.fMassLambdaC),
202 fMassDstar(source.fMassDstar),
203 fMassJpsi(source.fMassJpsi)
204 {
205   //
206   // Copy constructor
207   //
208 }
209 //--------------------------------------------------------------------------
210 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
211 {
212   //
213   // assignment operator
214   //
215   if(&source == this) return *this;
216   fInputAOD = source.fInputAOD;
217   fAODMapSize = source.fAODMapSize;
218   fVertexerTracks = source.fVertexerTracks;
219   fBzkG = source.fBzkG;
220   fSecVtxWithKF = source.fSecVtxWithKF;
221   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
222   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
223   fV1 = source.fV1;
224   fD0toKpi = source.fD0toKpi;
225   fJPSItoEle = source.fJPSItoEle;
226   f3Prong = source.f3Prong;
227   f4Prong = source.f4Prong;
228   fDstar = source.fDstar;
229   fCascades = source.fCascades;
230   fLikeSign = source.fLikeSign;
231   fLikeSign3prong = source.fLikeSign3prong;
232   fMixEvent = source.fMixEvent;
233   fPidResponse = source.fPidResponse;
234   fUseKaonPIDfor3Prong = source.fUseKaonPIDfor3Prong;
235   fnSigmaTOFKaonLow = source.fnSigmaTOFKaonLow;
236   fnSigmaTOFKaonHi = source.fnSigmaTOFKaonHi;
237   fMaxCentPercentileForTightCuts = source.fMaxCentPercentileForTightCuts;
238   fTrackFilter = source.fTrackFilter;
239   fTrackFilter2prongCentral = source.fTrackFilter2prongCentral;
240   fTrackFilter3prongCentral = source.fTrackFilter3prongCentral;
241   fTrackFilterSoftPi = source.fTrackFilterSoftPi;
242   fCutsD0toKpi = source.fCutsD0toKpi;
243   fCutsJpsitoee = source.fCutsJpsitoee;
244   fCutsDplustoKpipi = source.fCutsDplustoKpipi;
245   fCutsDstoKKpi = source.fCutsDstoKKpi;
246   fCutsLctopKpi = source.fCutsLctopKpi;
247   fCutsLctoV0 = source.fCutsLctoV0;
248   fCutsD0toKpipipi = source.fCutsD0toKpipipi;
249   fCutsDStartoKpipi = source.fCutsDStartoKpipi;
250   fListOfCuts = source.fListOfCuts;
251   fFindVertexForDstar = source.fFindVertexForDstar;
252   fFindVertexForCascades = source.fFindVertexForCascades;
253   fMassCutBeforeVertexing = source.fMassCutBeforeVertexing;
254   fMassCalc2 = source.fMassCalc2;
255   fMassCalc3 = source.fMassCalc3;
256   fMassCalc4 = source.fMassCalc4;
257   fOKInvMassD0 = source.fOKInvMassD0;
258   fOKInvMassJpsi = source.fOKInvMassJpsi;
259   fOKInvMassDplus = source.fOKInvMassDplus;
260   fOKInvMassDs = source.fOKInvMassDs;
261   fOKInvMassLc = source.fOKInvMassLc;
262   fOKInvMassDstar = source.fOKInvMassDstar;
263   fOKInvMassD0to4p = source.fOKInvMassD0to4p;
264   fOKInvMassLctoV0 = source.fOKInvMassLctoV0;
265   fMassDzero = source.fMassDzero;
266   fMassDplus = source.fMassDplus;
267   fMassDs = source.fMassDs;
268   fMassLambdaC = source.fMassLambdaC;
269   fMassDstar = source.fMassDstar;
270   fMassJpsi = source.fMassJpsi;
271
272   return *this;
273 }
274 //----------------------------------------------------------------------------
275 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
276   // Destructor
277   if(fV1) { delete fV1; fV1=0; }
278   delete fVertexerTracks;
279   if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
280   if(fTrackFilter2prongCentral) { delete fTrackFilter2prongCentral; fTrackFilter2prongCentral=0; }
281   if(fTrackFilter3prongCentral) { delete fTrackFilter3prongCentral; fTrackFilter3prongCentral=0; }
282   if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
283   if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
284   if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
285   if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
286   if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
287   if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
288   if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
289   if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
290   if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
291   if(fAODMap) { delete [] fAODMap; fAODMap=0; }
292   if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
293   if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
294   if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
295 }
296 //----------------------------------------------------------------------------
297 TList *AliAnalysisVertexingHF::FillListOfCuts() {
298   // Fill list of analysis cuts
299
300   TList *list = new TList();
301   list->SetOwner();
302   list->SetName("ListOfCuts");
303   
304   if(fCutsD0toKpi) {
305     AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
306     list->Add(cutsD0toKpi);
307   }
308   if(fCutsJpsitoee) {
309     AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
310     list->Add(cutsJpsitoee);
311   }
312   if(fCutsDplustoKpipi) {
313     AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
314     list->Add(cutsDplustoKpipi);
315   }
316   if(fCutsDstoKKpi) {
317     AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
318     list->Add(cutsDstoKKpi);
319   }
320   if(fCutsLctopKpi) {
321     AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
322     list->Add(cutsLctopKpi);
323   }
324   if(fCutsLctoV0){
325     AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
326     list->Add(cutsLctoV0);
327   }
328   if(fCutsD0toKpipipi) {
329     AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
330     list->Add(cutsD0toKpipipi);
331   }
332   if(fCutsDStartoKpipi) {
333     AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
334     list->Add(cutsDStartoKpipi);
335   }
336   
337   // keep a pointer to the list
338   fListOfCuts = list;
339
340   return list;
341 }
342 //----------------------------------------------------------------------------
343 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
344                                             TClonesArray *aodVerticesHFTClArr,
345                                             TClonesArray *aodD0toKpiTClArr,
346                                             TClonesArray *aodJPSItoEleTClArr,
347                                             TClonesArray *aodCharm3ProngTClArr,
348                                             TClonesArray *aodCharm4ProngTClArr,
349                                             TClonesArray *aodDstarTClArr,
350                                             TClonesArray *aodCascadesTClArr,
351                                             TClonesArray *aodLikeSign2ProngTClArr,
352                                             TClonesArray *aodLikeSign3ProngTClArr)
353 {
354   // Find heavy-flavour vertex candidates
355   // Input:  ESD or AOD
356   // Output: AOD (additional branches added)
357   //AliCodeTimerAuto("",0);
358
359   if(!fMixEvent){
360     TString evtype = event->IsA()->GetName();
361     fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
362   } // if we do mixing AliVEvent is a AliMixedEvent
363
364   if(fInputAOD) {
365     AliDebug(2,"Creating HF candidates from AOD");
366   } else {
367     AliDebug(2,"Creating HF candidates from ESD");
368   }
369
370   if(!aodVerticesHFTClArr) {
371     printf("ERROR: no aodVerticesHFTClArr");
372     return;
373   }
374   if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
375     printf("ERROR: no aodD0toKpiTClArr");
376     return;
377   }
378   if(fJPSItoEle && !aodJPSItoEleTClArr) {
379     printf("ERROR: no aodJPSItoEleTClArr");
380     return;
381   }
382   if(f3Prong && !aodCharm3ProngTClArr) {
383     printf("ERROR: no aodCharm3ProngTClArr");
384     return;
385   }
386   if(f4Prong && !aodCharm4ProngTClArr) {
387     printf("ERROR: no aodCharm4ProngTClArr");
388     return;
389   }
390   if(fDstar && !aodDstarTClArr) {
391     printf("ERROR: no aodDstarTClArr");
392     return;
393   }
394   if(fCascades && !aodCascadesTClArr){
395     printf("ERROR: no aodCascadesTClArr ");
396     return;
397   }
398   if(fLikeSign && !aodLikeSign2ProngTClArr) {
399     printf("ERROR: no aodLikeSign2ProngTClArr");
400     return;
401   }
402   if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
403     printf("ERROR: no aodLikeSign3ProngTClArr");
404     return;
405   }
406
407   // delete candidates from previous event and create references
408   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
409   aodVerticesHFTClArr->Delete();
410   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
411   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
412   if(fD0toKpi || fDstar)   {
413     aodD0toKpiTClArr->Delete();
414     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
415   }
416   if(fJPSItoEle) {
417     aodJPSItoEleTClArr->Delete();
418     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
419   }
420   if(f3Prong) {   
421     aodCharm3ProngTClArr->Delete();
422     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
423   }
424   if(f4Prong) {
425     aodCharm4ProngTClArr->Delete();
426     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
427   }
428   if(fDstar) {
429     aodDstarTClArr->Delete();
430     iDstar = aodDstarTClArr->GetEntriesFast();
431   }
432   if(fCascades) {
433     aodCascadesTClArr->Delete();
434     iCascades = aodCascadesTClArr->GetEntriesFast();
435   }
436   if(fLikeSign) {                                
437     aodLikeSign2ProngTClArr->Delete();                     
438     iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast(); 
439   }  
440   if(fLikeSign3prong && f3Prong) {                                
441     aodLikeSign3ProngTClArr->Delete();                     
442     iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast(); 
443   }  
444
445   TClonesArray &aodD0toKpiRef        = *aodD0toKpiTClArr;
446   TClonesArray &aodJPSItoEleRef      = *aodJPSItoEleTClArr;
447   TClonesArray &aodCharm3ProngRef    = *aodCharm3ProngTClArr;
448   TClonesArray &aodCharm4ProngRef    = *aodCharm4ProngTClArr;
449   TClonesArray &aodDstarRef          = *aodDstarTClArr;
450   TClonesArray &aodCascadesRef       = *aodCascadesTClArr;
451   TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
452   TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
453
454
455   AliAODRecoDecayHF2Prong *io2Prong  = 0;
456   AliAODRecoDecayHF3Prong *io3Prong  = 0;
457   AliAODRecoDecayHF4Prong *io4Prong  = 0;
458   AliAODRecoCascadeHF     *ioCascade = 0;
459
460   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
461   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
462   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
463   Bool_t   okDstar=kFALSE,okD0fromDstar=kFALSE;
464   Bool_t   okCascades=kFALSE;
465   AliESDtrack *postrack1 = 0;
466   AliESDtrack *postrack2 = 0;
467   AliESDtrack *negtrack1 = 0;
468   AliESDtrack *negtrack2 = 0;
469   AliESDtrack *trackPi   = 0;
470   Double_t mompos1[3],mompos2[3],momneg1[3],momneg2[3];
471   //   AliESDtrack *posV0track = 0;
472   //   AliESDtrack *negV0track = 0;
473   Float_t dcaMax = fCutsD0toKpi->GetDCACut();
474   if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
475   if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
476   if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
477   if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
478   if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
479   if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
480   
481   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
482
483
484   // get Bz
485   fBzkG = (Double_t)event->GetMagneticField(); 
486   if(!fVertexerTracks){
487     fVertexerTracks=new AliVertexerTracks(fBzkG);
488   }else{
489     Double_t oldField=fVertexerTracks->GetFieldkG();
490     if(oldField!=fBzkG) fVertexerTracks->SetFieldkG(fBzkG);
491   }
492
493   trkEntries = (Int_t)event->GetNumberOfTracks();
494   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
495   fnTrksTotal += trkEntries;
496
497   nv0 = (Int_t)event->GetNumberOfV0s();
498   AliDebug(1,Form(" Number of V0s: %d",nv0));
499
500   if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
501     AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
502     return;
503   }
504
505   // event selection + PID configuration
506   if(!fCutsD0toKpi->IsEventSelected(event)) return;
507   if(fCutsJpsitoee) fCutsJpsitoee->SetupPID(event);
508   if(fCutsDplustoKpipi) fCutsDplustoKpipi->SetupPID(event);
509   if(fCutsDstoKKpi) fCutsDstoKKpi->SetupPID(event);
510   if(fCutsLctopKpi) fCutsLctopKpi->SetupPID(event);
511   if(fCutsLctoV0) fCutsLctoV0->SetupPID(event);
512   if(fCutsD0toKpipipi) fCutsD0toKpipipi->SetupPID(event);
513   if(fCutsDStartoKpipi) fCutsDStartoKpipi->SetupPID(event);
514
515   // call function that applies sigle-track selection,
516   // for displaced tracks and soft pions (both charges) for D*,
517   // and retrieves primary vertex
518   TObjArray seleTrksArray(trkEntries);
519   TObjArray tracksAtVertex(trkEntries);
520   UChar_t  *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi, bit 2: 3 prong
521   Int_t     nSeleTrks=0;
522   Int_t *evtNumber    = new Int_t[trkEntries];
523   SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
524
525   AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
526   fnSeleTrksTotal += nSeleTrks;
527     
528
529   TObjArray *twoTrackArray1    = new TObjArray(2);
530   TObjArray *twoTrackArray2    = new TObjArray(2);
531   TObjArray *twoTrackArrayV0   = new TObjArray(2);
532   TObjArray *twoTrackArrayCasc = new TObjArray(2);
533   TObjArray *threeTrackArray   = new TObjArray(3);
534   TObjArray *fourTrackArray    = new TObjArray(4);
535
536   Double_t dispersion;
537   Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
538
539   AliAODRecoDecayHF   *rd = 0;
540   AliAODRecoCascadeHF *rc = 0;
541   AliAODv0            *v0 = 0;
542   AliESDv0         *esdV0 = 0;
543
544   Bool_t massCutOK=kTRUE;
545
546   // LOOP ON  POSITIVE  TRACKS
547   for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
548
549     //if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));  
550     //if(iTrkP1%1==0) printf("  1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);  
551
552     // get track from tracks array
553     postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
554     if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
555     postrack1->GetPxPyPz(mompos1);
556
557     // Make cascades with V0+track
558     // 
559     if(fCascades) {
560       // loop on V0's
561       for(iv0=0; iv0<nv0; iv0++){
562
563         //AliDebug(1,Form("   loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));   
564
565         // Get the V0 
566         if(fInputAOD) {
567           v0 = ((AliAODEvent*)event)->GetV0(iv0);
568         } else {
569           esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
570         }
571         if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) && 
572              (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
573         
574
575         // Get the tracks that form the V0
576         //  ( parameters at primary vertex )
577         //   and define an AliExternalTrackParam out of them
578         AliExternalTrackParam * posV0track;
579         AliExternalTrackParam * negV0track;
580
581         if(fInputAOD){
582           AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
583           AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
584           if( !posVV0track || !negVV0track ) continue;
585           //
586           // Apply some basic V0 daughter criteria
587           //
588           // bachelor must not be a v0-track                                                                  
589           if (posVV0track->GetID() == postrack1->GetID() ||
590               negVV0track->GetID() == postrack1->GetID()) continue;
591           // reject like-sign v0                                                                              
592           if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
593           // avoid ghost TPC tracks                                                                           
594           if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
595              !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
596           // Get AliExternalTrackParam out of the AliAODTracks
597           Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
598           posVV0track->PxPyPz(pxpypz);                    posVV0track->XvYvZv(xyz);
599           posVV0track->GetCovarianceXYZPxPyPz(cv);        sign=posVV0track->Charge();
600           posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
601           negVV0track->PxPyPz(pxpypz);                    negVV0track->XvYvZv(xyz);
602           negVV0track->GetCovarianceXYZPxPyPz(cv);        sign=negVV0track->Charge();
603           negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
604         }  else {
605           AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
606           AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
607           if( !posVV0track || !negVV0track ) continue;
608           //
609           // Apply some basic V0 daughter criteria
610           //
611           // bachelor must not be a v0-track                                                                  
612           if (posVV0track->GetID() == postrack1->GetID() ||
613               negVV0track->GetID() == postrack1->GetID()) continue;
614           // reject like-sign v0                                                                              
615           if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
616           // avoid ghost TPC tracks                                                                           
617           if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
618              !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
619           //  reject kinks (only necessary on AliESDtracks)
620           if (posVV0track->GetKinkIndex(0)>0  || negVV0track->GetKinkIndex(0)>0) continue;
621           // Get AliExternalTrackParam out of the AliESDtracks  
622           posV0track = new AliExternalTrackParam(*posVV0track);
623           negV0track = new AliExternalTrackParam(*negVV0track);
624
625           // Define the AODv0 from ESDv0 if reading ESDs
626           v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
627         }
628         if( !posV0track || !negV0track ){
629           AliDebug(1,Form(" Couldn't get the V0 daughters"));
630           continue;
631         }
632
633         // fill in the v0 two-external-track-param array
634         twoTrackArrayV0->AddAt(posV0track,0);
635         twoTrackArrayV0->AddAt(negV0track,1);
636
637         // Get the V0 dca
638         dcaV0 = v0->DcaV0Daughters();
639
640         // Define the V0 (neutral) track
641         AliNeutralTrackParam *trackV0;
642         if(fInputAOD) {
643           const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
644           if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
645         } else {  
646           Double_t xyz[3], pxpypz[3];
647           esdV0->XvYvZv(xyz);
648           esdV0->PxPyPz(pxpypz);
649           Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
650           trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
651         }
652         // Fill in the object array to create the cascade
653         twoTrackArrayCasc->AddAt(postrack1,0);
654         twoTrackArrayCasc->AddAt(trackV0,1);
655         // Compute the cascade vertex
656         AliAODVertex *vertexCasc = 0;
657         if(fFindVertexForCascades) {  
658           // DCA between the two tracks
659           dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
660           // Vertexing+   
661           vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
662         } else {
663           // assume Cascade decays at the primary vertex
664           Double_t pos[3],cov[6],chi2perNDF;
665           fV1->GetXYZ(pos);
666           fV1->GetCovMatrix(cov);
667           chi2perNDF = fV1->GetChi2toNDF();
668           vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
669           dcaCasc = 0.;
670         }
671         if(!vertexCasc) {
672           delete posV0track; posV0track=NULL;
673           delete negV0track; negV0track=NULL;
674           delete trackV0; trackV0=NULL;
675           if(!fInputAOD) {delete v0; v0=NULL;}
676           twoTrackArrayV0->Clear();
677           twoTrackArrayCasc->Clear();
678           continue; 
679         }
680
681         // Create and store the Cascade if passed the cuts
682         ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
683         if(okCascades && ioCascade) {
684           //AliDebug(1,Form("Storing a cascade object... "));
685           // add the vertex and the cascade to the AOD
686           AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
687           rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
688           rc->SetSecondaryVtx(vCasc);
689           vCasc->SetParent(rc);
690           rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
691           if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
692           AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
693           vCasc->AddDaughter(v0); // fill the 2prong V0 
694         }
695
696         // Clean up 
697         delete posV0track; posV0track=NULL;
698         delete negV0track; negV0track=NULL;
699         delete trackV0; trackV0=NULL;
700         twoTrackArrayV0->Clear();
701         twoTrackArrayCasc->Clear();
702         if(ioCascade) { delete ioCascade; ioCascade=NULL; }
703         if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
704         if(!fInputAOD) {delete v0; v0=NULL;}
705
706       } // end loop on V0's
707     } 
708   
709     // If there is less than 2 particles continue
710     if(trkEntries<2) {
711       AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
712       continue;
713     }
714
715     if(postrack1->Charge()<0 && !fLikeSign) continue;
716
717     // LOOP ON  NEGATIVE  TRACKS
718     for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
719
720       //if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));  
721       //if(iTrkN1%1==0) printf("    1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);  
722
723       if(iTrkN1==iTrkP1) continue;
724
725       // get track from tracks array
726       negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
727
728       if(negtrack1->Charge()>0 && !fLikeSign) continue;
729
730       if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
731
732       if(fMixEvent) {
733         if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
734       }
735
736       if(postrack1->Charge()==negtrack1->Charge()) { // like-sign 
737         isLikeSign2Prong=kTRUE;
738         if(!fLikeSign)    continue;
739         if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
740       } else { // unlike-sign
741         isLikeSign2Prong=kFALSE;
742         if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue;  // this is needed to avoid double-counting of unlike-sign
743         if(fMixEvent) {
744           if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
745         }
746        
747       }
748
749       // back to primary vertex
750       //      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
751       //      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
752       SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
753       SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
754       negtrack1->GetPxPyPz(momneg1);
755
756       // DCA between the two tracks
757       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
758       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
759
760       // Vertexing
761       twoTrackArray1->AddAt(postrack1,0);
762       twoTrackArray1->AddAt(negtrack1,1);
763       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
764       if(!vertexp1n1) {
765         twoTrackArray1->Clear();
766         negtrack1=0; 
767         continue; 
768       }
769
770       // 2 prong candidate
771       if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) { 
772       
773         io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
774       
775         if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
776           // add the vertex and the decay to the AOD
777           AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
778           if(!isLikeSign2Prong) {
779             if(okD0) {  
780               rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
781               rd->SetSecondaryVtx(v2Prong);
782               v2Prong->SetParent(rd);
783               AddRefs(v2Prong,rd,event,twoTrackArray1);
784             }
785             if(okJPSI) {
786               rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
787               rd->SetSecondaryVtx(v2Prong);
788               if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
789               AddRefs(v2Prong,rd,event,twoTrackArray1);
790             }
791           } else { // isLikeSign2Prong
792             rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
793             rd->SetSecondaryVtx(v2Prong);
794             v2Prong->SetParent(rd);
795             AddRefs(v2Prong,rd,event,twoTrackArray1);
796           }
797           // Set selection bit for PID
798           if(okD0) SetSelectionBitForPID(fCutsD0toKpi,rd,AliRDHFCuts::kD0toKpiPID);
799         }
800         // D* candidates
801         if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
802           // write references in io2Prong
803           if(fInputAOD) {
804             AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
805           } else {
806             vertexp1n1->AddDaughter(postrack1);
807             vertexp1n1->AddDaughter(negtrack1);
808           }
809           io2Prong->SetSecondaryVtx(vertexp1n1);
810           //printf("--->  %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
811           // create a track from the D0
812           AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
813
814           // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
815           for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
816
817             if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
818
819             if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
820
821             if(fMixEvent) {
822               if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] || 
823                  evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] || 
824                  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
825             }
826
827             //if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));  
828
829             trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
830             if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
831
832             // get track from tracks array
833             trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
834             //      trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
835             SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
836             twoTrackArrayCasc->AddAt(trackPi,0);
837             twoTrackArrayCasc->AddAt(trackD0,1);
838             if(!SelectInvMassAndPtDstarD0pi(twoTrackArrayCasc)){
839               twoTrackArrayCasc->Clear();
840               trackPi=0; 
841               continue; 
842             }
843
844             AliAODVertex *vertexCasc = 0;
845
846             if(fFindVertexForDstar) {
847               // DCA between the two tracks
848               dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
849               // Vertexing
850               vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
851             } else {
852               // assume Dstar decays at the primary vertex
853               Double_t pos[3],cov[6],chi2perNDF;
854               fV1->GetXYZ(pos);
855               fV1->GetCovMatrix(cov);
856               chi2perNDF = fV1->GetChi2toNDF();
857               vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
858               dcaCasc = 0.;
859             }
860             if(!vertexCasc) { 
861               twoTrackArrayCasc->Clear();
862               trackPi=0; 
863               continue; 
864             }
865
866             ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
867             if(okDstar) {
868               // add the D0 to the AOD (if not already done)
869               if(!okD0) {
870                 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
871                 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
872                 rd->SetSecondaryVtx(v2Prong);
873                 v2Prong->SetParent(rd);
874                 AddRefs(v2Prong,rd,event,twoTrackArray1);
875                 okD0=kTRUE; // this is done to add it only once
876               }
877               // add the vertex and the cascade to the AOD
878               AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc); 
879               rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
880               rc->SetSecondaryVtx(vCasc);
881               vCasc->SetParent(rc);
882               if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0 
883               AddRefs(vCasc,rc,event,twoTrackArrayCasc);
884               vCasc->AddDaughter(rd); // add the D0 (in ref #1)
885               // Set selection bit for PID
886               SetSelectionBitForPID(fCutsDStartoKpipi,rc,AliRDHFCuts::kDstarPID);
887             }
888             twoTrackArrayCasc->Clear();
889             trackPi=0; 
890             if(ioCascade) {delete ioCascade; ioCascade=NULL;}
891             delete vertexCasc; vertexCasc=NULL;
892           } // end loop on soft pi tracks
893
894           if(trackD0) {delete trackD0; trackD0=NULL;}
895
896         }
897         if(io2Prong) {delete io2Prong; io2Prong=NULL;}
898       }      
899
900       twoTrackArray1->Clear(); 
901       if( (!f3Prong && !f4Prong) || 
902           (isLikeSign2Prong && !f3Prong) ) { 
903         negtrack1=0; 
904         delete vertexp1n1; 
905         continue; 
906       }
907
908         
909       // 2nd LOOP  ON  POSITIVE  TRACKS 
910       for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
911
912         if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
913
914         //if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));  
915
916         // get track from tracks array
917         postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
918
919         if(postrack2->Charge()<0) continue; 
920
921         if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
922
923         // Check single tracks cuts specific for 3 prongs
924         if(!TESTBIT(seleFlags[iTrkP2],kBit3Prong)) continue;
925         if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
926         if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
927
928         if(fMixEvent) {
929           if(evtNumber[iTrkP1]==evtNumber[iTrkP2] || 
930              evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
931              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
932         }
933
934         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet 
935           if(!fLikeSign3prong) continue;
936           if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
937             isLikeSign3Prong=kTRUE;
938           } else { // not ok
939             continue;
940           }
941         } else { // normal triplet (+-+)
942           isLikeSign3Prong=kFALSE; 
943           if(fMixEvent) {
944             if(evtNumber[iTrkP1]==evtNumber[iTrkP2] || 
945                evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
946                evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
947           }
948         }
949
950         if(fUseKaonPIDfor3Prong){
951           Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(negtrack1,AliPID::kKaon);
952           if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)) continue;
953         }
954         // back to primary vertex
955         //      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
956         //      postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
957         //      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
958         SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
959         SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
960         SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
961       
962         //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
963
964         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
965         if(dcap2n1>dcaMax) { postrack2=0; continue; }
966         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
967         if(dcap1p2>dcaMax) { postrack2=0; continue; }
968
969         // check invariant mass cuts for D+,Ds,Lc
970         massCutOK=kTRUE;
971         if(f3Prong) {
972           if(postrack2->Charge()>0) {
973             threeTrackArray->AddAt(postrack1,0);
974             threeTrackArray->AddAt(negtrack1,1);
975             threeTrackArray->AddAt(postrack2,2);
976           } else {
977             threeTrackArray->AddAt(negtrack1,0);
978             threeTrackArray->AddAt(postrack1,1);
979             threeTrackArray->AddAt(postrack2,2);
980           }
981           if(fMassCutBeforeVertexing){
982             postrack2->GetPxPyPz(mompos2);
983             Double_t pxDau[3]={mompos1[0],momneg1[0],mompos2[0]};
984             Double_t pyDau[3]={mompos1[1],momneg1[1],mompos2[1]};
985             Double_t pzDau[3]={mompos1[2],momneg1[2],mompos2[2]};
986             //      massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
987             massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau);
988           }
989         }
990
991         if(f3Prong && !massCutOK) {
992           threeTrackArray->Clear();
993           if(!f4Prong) { 
994             postrack2=0; 
995             continue;
996           } 
997         } 
998         
999         // Vertexing
1000         twoTrackArray2->AddAt(postrack2,0);
1001         twoTrackArray2->AddAt(negtrack1,1);
1002         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1003         if(!vertexp2n1) { 
1004           twoTrackArray2->Clear();
1005           postrack2=0; 
1006           continue;
1007         }
1008
1009         // 3 prong candidates
1010         if(f3Prong && massCutOK) { 
1011
1012           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1013           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
1014           if(ok3Prong) {
1015             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1016             if(!isLikeSign3Prong) {
1017               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1018               rd->SetSecondaryVtx(v3Prong);
1019               v3Prong->SetParent(rd);
1020               AddRefs(v3Prong,rd,event,threeTrackArray);
1021             } else { // isLikeSign3Prong
1022               if(fLikeSign3prong){
1023                 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1024                 rd->SetSecondaryVtx(v3Prong);
1025                 v3Prong->SetParent(rd);
1026                 AddRefs(v3Prong,rd,event,threeTrackArray);
1027               }
1028             }
1029             // Set selection bit for PID
1030             SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
1031             SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
1032             SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
1033           }
1034           if(io3Prong) {delete io3Prong; io3Prong=NULL;} 
1035           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;} 
1036         }
1037
1038         // 4 prong candidates
1039         if(f4Prong 
1040            // don't make 4 prong with like-sign pairs and triplets
1041            && !isLikeSign2Prong && !isLikeSign3Prong
1042            // track-to-track dca cuts already now
1043            && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
1044            && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
1045
1046           // back to primary vertex
1047           //      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1048           //      postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1049           //      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1050           SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1051           SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1052           SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1053
1054           // Vertexing for these 3 (can be taken from above?)
1055           threeTrackArray->AddAt(postrack1,0);
1056           threeTrackArray->AddAt(negtrack1,1);
1057           threeTrackArray->AddAt(postrack2,2);
1058           AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1059
1060           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong) 
1061           for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1062
1063             if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1064
1065             //if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
1066
1067             // get track from tracks array
1068             negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1069
1070             if(negtrack2->Charge()>0) continue;
1071
1072             if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1073             if(fMixEvent){ 
1074               if(evtNumber[iTrkP1]==evtNumber[iTrkN2] || 
1075                  evtNumber[iTrkN1]==evtNumber[iTrkN2] || 
1076                  evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
1077                  evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
1078                  evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1079                  evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
1080             }
1081
1082             // back to primary vertex
1083             // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1084             // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1085             // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1086             // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1087             SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1088             SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1089             SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1090             SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1091
1092             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1093             if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1094             dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1095             if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1096
1097
1098             fourTrackArray->AddAt(postrack1,0);
1099             fourTrackArray->AddAt(negtrack1,1);
1100             fourTrackArray->AddAt(postrack2,2);
1101             fourTrackArray->AddAt(negtrack2,3);
1102
1103             // check invariant mass cuts for D0
1104             massCutOK=kTRUE;
1105             if(fMassCutBeforeVertexing) 
1106               massCutOK = SelectInvMassAndPt4prong(fourTrackArray);
1107             
1108             if(!massCutOK) {
1109               fourTrackArray->Clear(); 
1110               negtrack2=0; 
1111               continue; 
1112             }
1113
1114             // Vertexing
1115             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
1116             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
1117             if(ok4Prong) {
1118               AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
1119               rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
1120               rd->SetSecondaryVtx(v4Prong);
1121               v4Prong->SetParent(rd);
1122               AddRefs(v4Prong,rd,event,fourTrackArray);
1123             }
1124
1125             if(io4Prong) {delete io4Prong; io4Prong=NULL;} 
1126             if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;} 
1127             fourTrackArray->Clear();
1128             negtrack2 = 0;
1129
1130           } // end loop on negative tracks
1131
1132           threeTrackArray->Clear();
1133           delete vertexp1n1p2;
1134
1135         }
1136
1137         postrack2 = 0;
1138         delete vertexp2n1;
1139
1140       } // end 2nd loop on positive tracks
1141
1142       twoTrackArray2->Clear();
1143       
1144       // 2nd LOOP  ON  NEGATIVE  TRACKS (for 3 prong -+-)
1145       for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1146
1147         if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1148
1149         //if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
1150
1151         // get track from tracks array
1152         negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1153
1154         if(negtrack2->Charge()>0) continue;
1155
1156         if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1157
1158         // Check single tracks cuts specific for 3 prongs
1159         if(!TESTBIT(seleFlags[iTrkN2],kBit3Prong)) continue;
1160         if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
1161         if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
1162
1163         if(fMixEvent) {
1164           if(evtNumber[iTrkP1]==evtNumber[iTrkN2] || 
1165              evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1166              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1167         }
1168
1169         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet 
1170           if(!fLikeSign3prong) continue;
1171           if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1172             isLikeSign3Prong=kTRUE;
1173           } else { // not ok
1174             continue;
1175           }
1176         } else { // normal triplet (-+-)
1177           isLikeSign3Prong=kFALSE;
1178         }
1179
1180         if(fUseKaonPIDfor3Prong){
1181           Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(postrack1,AliPID::kKaon);
1182           if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)) continue;
1183         }
1184
1185         // back to primary vertex
1186         // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1187         // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1188         // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1189         SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1190         SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1191         SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1192         //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1193
1194         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1195         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1196         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1197         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1198
1199         threeTrackArray->AddAt(negtrack1,0);
1200         threeTrackArray->AddAt(postrack1,1);
1201         threeTrackArray->AddAt(negtrack2,2);
1202
1203         // check invariant mass cuts for D+,Ds,Lc
1204         massCutOK=kTRUE;
1205         if(fMassCutBeforeVertexing && f3Prong){ 
1206           negtrack2->GetPxPyPz(momneg2);
1207           Double_t pxDau[3]={momneg1[0],mompos1[0],momneg2[0]};
1208           Double_t pyDau[3]={momneg1[1],mompos1[1],momneg2[1]};
1209           Double_t pzDau[3]={momneg1[2],mompos1[2],momneg2[2]};
1210           //      massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1211           massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau);
1212         }
1213         if(!massCutOK) { 
1214           threeTrackArray->Clear();
1215           negtrack2=0; 
1216           continue; 
1217         }
1218         
1219         // Vertexing
1220         twoTrackArray2->AddAt(postrack1,0);
1221         twoTrackArray2->AddAt(negtrack2,1);
1222
1223         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1224         if(!vertexp1n2) { 
1225           twoTrackArray2->Clear();
1226           negtrack2=0; 
1227           continue; 
1228         }
1229
1230         if(f3Prong) { 
1231           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1232           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1233           if(ok3Prong) {
1234             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1235             if(!isLikeSign3Prong) {
1236               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1237               rd->SetSecondaryVtx(v3Prong);
1238               v3Prong->SetParent(rd);
1239               AddRefs(v3Prong,rd,event,threeTrackArray);
1240             } else { // isLikeSign3Prong
1241               if(fLikeSign3prong){
1242                 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1243                 rd->SetSecondaryVtx(v3Prong);
1244                 v3Prong->SetParent(rd);
1245                 AddRefs(v3Prong,rd,event,threeTrackArray);
1246               }
1247             }
1248             // Set selection bit for PID
1249             SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
1250             SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
1251             SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
1252           }
1253           if(io3Prong) {delete io3Prong; io3Prong=NULL;} 
1254           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1255         }
1256         threeTrackArray->Clear();
1257         negtrack2 = 0;
1258         delete vertexp1n2;
1259
1260       } // end 2nd loop on negative tracks
1261       
1262       twoTrackArray2->Clear();
1263       
1264       negtrack1 = 0;
1265       delete vertexp1n1; 
1266     } // end 1st loop on negative tracks
1267     
1268     postrack1 = 0;
1269   }  // end 1st loop on positive tracks
1270
1271
1272   AliDebug(1,Form(" Total HF vertices in event = %d;",
1273                   (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1274   if(fD0toKpi) {
1275     AliDebug(1,Form(" D0->Kpi in event = %d;",
1276                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1277   }
1278   if(fJPSItoEle) {
1279     AliDebug(1,Form(" JPSI->ee in event = %d;",
1280                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1281   }
1282   if(f3Prong) {
1283     AliDebug(1,Form(" Charm->3Prong in event = %d;",
1284                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1285   }
1286   if(f4Prong) {
1287     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1288                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1289   }
1290   if(fDstar) {
1291     AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1292                     (Int_t)aodDstarTClArr->GetEntriesFast()));
1293   }
1294   if(fCascades){
1295     AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1296                     (Int_t)aodCascadesTClArr->GetEntriesFast()));
1297   }
1298   if(fLikeSign) {
1299     AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1300                     (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1301   }
1302   if(fLikeSign3prong && f3Prong) {
1303     AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1304                     (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1305   }
1306     
1307
1308   twoTrackArray1->Delete();  delete twoTrackArray1;
1309   twoTrackArray2->Delete();  delete twoTrackArray2;
1310   twoTrackArrayCasc->Delete();  delete twoTrackArrayCasc;
1311   twoTrackArrayV0->Delete();  delete twoTrackArrayV0;
1312   threeTrackArray->Clear(); 
1313   threeTrackArray->Delete(); delete threeTrackArray;
1314   fourTrackArray->Delete();  delete fourTrackArray;
1315   delete [] seleFlags; seleFlags=NULL;
1316   if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1317   tracksAtVertex.Delete();
1318
1319   if(fInputAOD) {
1320     seleTrksArray.Delete(); 
1321     if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
1322   }
1323   
1324
1325   //printf("Trks: total %d  sele %d\n",fnTrksTotal,fnSeleTrksTotal);
1326
1327   return;
1328 }
1329 //----------------------------------------------------------------------------
1330 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1331                                      const AliVEvent *event,
1332                                      const TObjArray *trkArray) const
1333 {
1334   // Add the AOD tracks as daughters of the vertex (TRef)
1335   // Also add the references to the primary vertex and to the cuts
1336   //AliCodeTimerAuto("",0);
1337
1338   if(fInputAOD) {
1339     AddDaughterRefs(v,event,trkArray);
1340     rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1341   }
1342
1343   /*
1344   rd->SetListOfCutsRef((TList*)fListOfCuts);
1345   //fListOfCuts->Print();
1346   cout<<fListOfCuts<<endl;
1347   TList *l=(TList*)rd->GetListOfCuts();
1348   cout<<l<<endl;
1349   if(l) {l->Print(); }else{printf("error\n");}
1350   */
1351
1352   return;
1353 }       
1354 //----------------------------------------------------------------------------
1355 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1356                                              const AliVEvent *event,
1357                                              const TObjArray *trkArray) const
1358 {
1359   // Add the AOD tracks as daughters of the vertex (TRef)
1360   //AliCodeTimerAuto("",0);
1361
1362   Int_t nDg = v->GetNDaughters();
1363   TObject *dg = 0;
1364   if(nDg) dg = v->GetDaughter(0);
1365   
1366   if(dg) return; // daughters already added
1367
1368   Int_t nTrks = trkArray->GetEntriesFast();
1369
1370   AliExternalTrackParam *track = 0;
1371   AliAODTrack *aodTrack = 0;
1372   Int_t id;
1373
1374   for(Int_t i=0; i<nTrks; i++) {
1375     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1376     id = (Int_t)track->GetID();
1377     //printf("---> %d\n",id);
1378     if(id<0) continue; // this track is a AliAODRecoDecay
1379     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1380     v->AddDaughter(aodTrack);
1381   }
1382
1383   return;
1384 }       
1385 //---------------------------------------------------------------------------
1386 void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)  
1387 {
1388   // Checks that the references to the daughter tracks are properly
1389   // assigned and reassigns them if needed
1390   //
1391   //AliCodeTimerAuto("",0);
1392
1393
1394   TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1395   if(!inputArray) return;
1396
1397   AliAODTrack *track = 0;
1398   AliAODVertex *vertex = 0;
1399
1400   Bool_t needtofix=kFALSE;
1401   for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1402     vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1403     for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1404       track = (AliAODTrack*)vertex->GetDaughter(id);
1405       if(!track->GetStatus()) needtofix=kTRUE;
1406     }
1407     if(needtofix) break;
1408   }
1409
1410   if(!needtofix) return;
1411
1412
1413   printf("Fixing references\n");
1414
1415   fAODMapSize = 100000;
1416   fAODMap = new Int_t[fAODMapSize];
1417   memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
1418
1419   for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1420     track = aod->GetTrack(i);
1421
1422     // skip pure ITS SA tracks
1423     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1424
1425     // skip tracks without ITS
1426     if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1427
1428     // TEMPORARY: check that the cov matrix is there
1429     Double_t covtest[21];
1430     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1431     //
1432
1433     Int_t ind = (Int_t)track->GetID();
1434     if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
1435   }
1436
1437
1438   Int_t ids[4]={-1,-1,-1,-1};
1439   for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1440     Bool_t cascade=kFALSE;
1441     vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1442     Int_t id=0;
1443     Int_t nDgs = vertex->GetNDaughters();
1444     for(id=0; id<nDgs; id++) {
1445       track = (AliAODTrack*)vertex->GetDaughter(id);
1446       if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1447       ids[id]=(Int_t)track->GetID();
1448       vertex->RemoveDaughter(track);
1449     }
1450     if(cascade) continue;
1451     for(id=0; id<nDgs; id++) {
1452       if (ids[id]>-1 && ids[id] < fAODMapSize) {
1453         track = aod->GetTrack(fAODMap[ids[id]]);
1454         vertex->AddDaughter(track);
1455       }
1456     }
1457     
1458   }
1459
1460   return;
1461 }
1462 //----------------------------------------------------------------------------
1463 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1464                                    TObjArray *twoTrackArray,AliVEvent *event,
1465                                    AliAODVertex *secVert,
1466                                    AliAODRecoDecayHF2Prong *rd2Prong,
1467                                    Double_t dca,
1468                                    Bool_t &okDstar) 
1469 {
1470   // Make the cascade as a 2Prong decay and check if it passes Dstar
1471   // reconstruction cuts
1472   //AliCodeTimerAuto("",0);
1473
1474   okDstar = kFALSE;
1475
1476   Bool_t dummy1,dummy2,dummy3;
1477
1478   // We use Make2Prong to construct the AliAODRecoCascadeHF
1479   // (which inherits from AliAODRecoDecayHF2Prong) 
1480   AliAODRecoCascadeHF *theCascade = 
1481     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1482                                      dummy1,dummy2,dummy3);
1483   if(!theCascade) return 0x0;
1484
1485   // charge
1486   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1487   theCascade->SetCharge(trackPi->Charge());
1488
1489   //--- selection cuts
1490   //
1491   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1492   if(fInputAOD){
1493     Int_t idSoftPi=(Int_t)trackPi->GetID();
1494     if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
1495       AliAODTrack* trackPiAOD=(AliAODTrack*)event->GetTrack(fAODMap[idSoftPi]);
1496       tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
1497     }
1498   }else{
1499     tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1500   }
1501   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1502
1503   AliAODVertex *primVertexAOD=0;
1504   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1505     // take event primary vertex
1506     primVertexAOD = PrimaryVertex(); 
1507     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1508     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1509   }
1510   // select D*->D0pi
1511   if(fDstar) {
1512     okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1513     if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1514   }
1515   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1516   tmpCascade->UnsetOwnPrimaryVtx(); 
1517   delete tmpCascade; tmpCascade=NULL;
1518   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1519     rd2Prong->UnsetOwnPrimaryVtx();
1520   }
1521   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1522   //---
1523
1524   
1525   return theCascade;
1526 }
1527
1528
1529 //----------------------------------------------------------------------------
1530 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1531                                    TObjArray *twoTrackArray,AliVEvent *event,
1532                                    AliAODVertex *secVert,
1533                                    AliAODv0 *v0,
1534                                    Double_t dca,
1535                                    Bool_t &okCascades) 
1536 {
1537   //
1538   // Make the cascade as a 2Prong decay and check if it passes 
1539   // cascades reconstruction cuts
1540   //AliCodeTimerAuto("",0);
1541   
1542   //  AliDebug(2,Form("         building the cascade"));
1543   okCascades= kFALSE; 
1544   Bool_t dummy1,dummy2,dummy3;
1545
1546   // We use Make2Prong to construct the AliAODRecoCascadeHF
1547   // (which inherits from AliAODRecoDecayHF2Prong) 
1548   AliAODRecoCascadeHF *theCascade = 
1549     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1550                                      dummy1,dummy2,dummy3);
1551   if(!theCascade) return 0x0;
1552
1553   // bachelor track and charge
1554   AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1555   theCascade->SetCharge(trackBachelor->Charge());
1556
1557   //--- selection cuts
1558   //
1559
1560   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);  
1561   if(fInputAOD){
1562     Int_t idBachelor=(Int_t)trackBachelor->GetID();
1563     if (idBachelor > -1 && idBachelor < fAODMapSize) {
1564       AliAODTrack* trackBachelorAOD=(AliAODTrack*)event->GetTrack(fAODMap[idBachelor]);
1565       tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
1566     }
1567   }else{
1568     tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1569   }
1570   tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1571   
1572   AliAODVertex *primVertexAOD=0;
1573   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1574     // take event primary vertex
1575     primVertexAOD = PrimaryVertex(); 
1576     if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex(); 
1577     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1578   }
1579   
1580   // select Cascades
1581   if(fCascades && fInputAOD){
1582     okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1583   }
1584   else { 
1585     //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); 
1586     okCascades=kTRUE; 
1587   }// no cuts implemented from ESDs
1588   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1589   tmpCascade->UnsetOwnPrimaryVtx(); 
1590   delete tmpCascade; tmpCascade=NULL;
1591   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1592   //---
1593   
1594   return theCascade;
1595 }
1596
1597 //-----------------------------------------------------------------------------
1598 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1599                                    TObjArray *twoTrackArray,AliVEvent *event,
1600                                    AliAODVertex *secVert,Double_t dca,
1601                                    Bool_t &okD0,Bool_t &okJPSI,
1602                                    Bool_t &okD0fromDstar) 
1603 {
1604   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1605   // reconstruction cuts
1606   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1607   //AliCodeTimerAuto("",0);
1608
1609   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1610
1611   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1612   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1613   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1614
1615   // propagate tracks to secondary vertex, to compute inv. mass
1616   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1617   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1618
1619   Double_t momentum[3];
1620   postrack->GetPxPyPz(momentum);
1621   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
1622   negtrack->GetPxPyPz(momentum);
1623   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
1624
1625
1626   // invariant mass cut (try to improve coding here..)
1627   Bool_t okMassCut=kFALSE;
1628   if(!okMassCut && fD0toKpi)   if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
1629   if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
1630   if(!okMassCut && fDstar)     if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
1631   if(!okMassCut && fCascades)  if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
1632   if(!okMassCut) {
1633     //AliDebug(2," candidate didn't pass mass cut");
1634     return 0x0;    
1635   }
1636   // primary vertex to be used by this candidate
1637   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
1638   if(!primVertexAOD) return 0x0;
1639
1640
1641   Double_t d0z0[2],covd0z0[3];
1642   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1643   d0[0] = d0z0[0];
1644   d0err[0] = TMath::Sqrt(covd0z0[0]);
1645   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1646   d0[1] = d0z0[0];
1647   d0err[1] = TMath::Sqrt(covd0z0[0]);
1648   
1649   // create the object AliAODRecoDecayHF2Prong
1650   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1651   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1652   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1653   the2Prong->SetProngIDs(2,id);
1654   delete primVertexAOD; primVertexAOD=NULL;
1655
1656  
1657   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar 
1658
1659     // Add daughter references already here
1660     if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
1661
1662     // select D0->Kpi
1663     if(fD0toKpi)   {
1664       okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
1665       if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1666     }
1667     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
1668     // select J/psi from B
1669     if(fJPSItoEle)   {
1670       okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1671     }
1672     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
1673     // select D0->Kpi from Dstar
1674     if(fDstar)   {
1675       okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1676       if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1677     }
1678     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
1679   }
1680
1681
1682   // remove the primary vertex (was used only for selection)
1683   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1684     the2Prong->UnsetOwnPrimaryVtx();
1685   }
1686   
1687   // get PID info from ESD
1688   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1689   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1690   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1691   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1692   Double_t esdpid[10];
1693   for(Int_t i=0;i<5;i++) {
1694     esdpid[i]   = esdpid0[i];
1695     esdpid[5+i] = esdpid1[i];
1696   }
1697   the2Prong->SetPID(2,esdpid);
1698
1699   return the2Prong;  
1700 }
1701 //----------------------------------------------------------------------------
1702 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1703                              TObjArray *threeTrackArray,AliVEvent *event,
1704                              AliAODVertex *secVert,Double_t dispersion,
1705                              const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1706                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1707                              Bool_t &ok3Prong) 
1708 {
1709   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1710   // reconstruction cuts 
1711   // E.Bruna, F.Prino
1712
1713   //AliCodeTimerAuto("",0);
1714
1715   ok3Prong=kFALSE;
1716   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0; 
1717
1718   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1719   Double_t momentum[3];
1720
1721
1722   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1723   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1724   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1725
1726   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1727   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1728   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1729   postrack1->GetPxPyPz(momentum);
1730   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
1731   negtrack->GetPxPyPz(momentum);
1732   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
1733   postrack2->GetPxPyPz(momentum);
1734   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2]; 
1735
1736   // invariant mass cut for D+, Ds, Lc
1737   Bool_t okMassCut=kFALSE;
1738   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed 
1739   if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
1740   if(!okMassCut) {
1741     //AliDebug(2," candidate didn't pass mass cut");
1742     return 0x0;    
1743   }
1744
1745   // primary vertex to be used by this candidate
1746   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
1747   if(!primVertexAOD) return 0x0;
1748
1749   Double_t d0z0[2],covd0z0[3];
1750   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1751   d0[0]=d0z0[0];
1752   d0err[0] = TMath::Sqrt(covd0z0[0]);
1753   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1754   d0[1]=d0z0[0];
1755   d0err[1] = TMath::Sqrt(covd0z0[0]);
1756   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1757   d0[2]=d0z0[0];
1758   d0err[2] = TMath::Sqrt(covd0z0[0]);
1759
1760
1761   // create the object AliAODRecoDecayHF3Prong
1762   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1763   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1764   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]));
1765   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]));
1766   Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1767   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1768   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1769   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1770   the3Prong->SetProngIDs(3,id);
1771
1772   delete primVertexAOD; primVertexAOD=NULL;
1773
1774   // Add daughter references already here
1775   if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
1776
1777   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1778   if(f3Prong) {
1779     ok3Prong = kFALSE;
1780     
1781     if(fOKInvMassDplus && fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1782       ok3Prong = kTRUE;
1783       the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1784     }
1785     if(fOKInvMassDs && fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1786       ok3Prong = kTRUE;
1787       the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1788     }
1789     if(fOKInvMassLc && fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1790       ok3Prong = kTRUE;
1791       the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1792     } 
1793   }
1794   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1795
1796   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1797     the3Prong->UnsetOwnPrimaryVtx();
1798   }
1799
1800   // get PID info from ESD
1801   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1802   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1803   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1804   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1805   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1806   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1807   
1808   Double_t esdpid[15];
1809   for(Int_t i=0;i<5;i++) {
1810     esdpid[i]    = esdpid0[i];
1811     esdpid[5+i]  = esdpid1[i];
1812     esdpid[10+i] = esdpid2[i];
1813   }
1814   the3Prong->SetPID(3,esdpid);
1815
1816   return the3Prong;
1817 }
1818 //----------------------------------------------------------------------------
1819 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1820                              TObjArray *fourTrackArray,AliVEvent *event,
1821                              AliAODVertex *secVert,
1822                              const AliAODVertex *vertexp1n1,
1823                              const AliAODVertex *vertexp1n1p2,
1824                              Double_t dcap1n1,Double_t dcap1n2,
1825                              Double_t dcap2n1,Double_t dcap2n2,
1826                              Bool_t &ok4Prong) 
1827 {
1828   // Make 4Prong candidates and check if they pass D0toKpipipi
1829   // reconstruction cuts
1830   // G.E.Bruno, R.Romita
1831   //AliCodeTimerAuto("",0);
1832
1833   ok4Prong=kFALSE;
1834   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0; 
1835
1836   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1837
1838   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1839   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1840   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1841   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1842
1843   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1844   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1845   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1846   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1847
1848   Double_t momentum[3];
1849   postrack1->GetPxPyPz(momentum);
1850   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1851   negtrack1->GetPxPyPz(momentum);
1852   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1853   postrack2->GetPxPyPz(momentum);
1854   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1855   negtrack2->GetPxPyPz(momentum);
1856   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1857
1858   // invariant mass cut for rho or D0 (try to improve coding here..)
1859   Bool_t okMassCut=kFALSE;
1860   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed 
1861   if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) {      //no PID, to be implemented with PID
1862     if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
1863   }
1864   if(!okMassCut) {
1865     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1866     //printf(" candidate didn't pass mass cut\n");
1867     return 0x0;
1868   }
1869
1870   // primary vertex to be used by this candidate
1871   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1872   if(!primVertexAOD) return 0x0;
1873
1874   Double_t d0z0[2],covd0z0[3];
1875   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1876   d0[0]=d0z0[0];
1877   d0err[0] = TMath::Sqrt(covd0z0[0]);
1878   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1879   d0[1]=d0z0[0];
1880   d0err[1] = TMath::Sqrt(covd0z0[0]);
1881   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1882   d0[2]=d0z0[0];
1883   d0err[2] = TMath::Sqrt(covd0z0[0]);
1884   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1885   d0[3]=d0z0[0];
1886   d0err[3] = TMath::Sqrt(covd0z0[0]);
1887
1888
1889   // create the object AliAODRecoDecayHF4Prong
1890   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1891   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1892   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]));
1893   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]));
1894   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]));
1895   Short_t charge=0;
1896   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1897   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1898   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1899   the4Prong->SetProngIDs(4,id);
1900
1901   delete primVertexAOD; primVertexAOD=NULL;
1902
1903   ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1904
1905
1906   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1907     the4Prong->UnsetOwnPrimaryVtx();
1908   }
1909
1910  
1911   // get PID info from ESD
1912   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1913   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1914   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1915   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1916   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1917   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1918   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1919   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1920
1921   Double_t esdpid[20];
1922   for(Int_t i=0;i<5;i++) {
1923     esdpid[i]    = esdpid0[i];
1924     esdpid[5+i]  = esdpid1[i];
1925     esdpid[10+i] = esdpid2[i];
1926     esdpid[15+i] = esdpid3[i];
1927   }
1928   the4Prong->SetPID(4,esdpid);
1929   
1930   return the4Prong;
1931 }
1932 //-----------------------------------------------------------------------------
1933 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1934                                                     AliVEvent *event) const
1935 {
1936   // Returns primary vertex to be used for this candidate
1937   //AliCodeTimerAuto("",0);
1938
1939   AliESDVertex *vertexESD = 0;
1940   AliAODVertex *vertexAOD = 0;
1941
1942
1943   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) { 
1944     // primary vertex from the input event
1945     
1946     vertexESD = new AliESDVertex(*fV1);
1947
1948   } else {
1949     // primary vertex specific to this candidate
1950
1951     Int_t nTrks = trkArray->GetEntriesFast();
1952     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1953
1954     if(fRecoPrimVtxSkippingTrks) { 
1955       // recalculating the vertex
1956       
1957       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1958         Float_t diamondcovxy[3];
1959         event->GetDiamondCovXY(diamondcovxy);
1960         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1961         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1962         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1963         vertexer->SetVtxStart(diamond);
1964         delete diamond; diamond=NULL;
1965         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) 
1966           vertexer->SetOnlyFitter();
1967       }
1968       Int_t skipped[1000];
1969       Int_t nTrksToSkip=0,id;
1970       AliExternalTrackParam *t = 0;
1971       for(Int_t i=0; i<nTrks; i++) {
1972         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1973         id = (Int_t)t->GetID();
1974         if(id<0) continue;
1975         skipped[nTrksToSkip++] = id;
1976       }
1977       // TEMPORARY FIX
1978       // For AOD, skip also tracks without covariance matrix
1979       if(fInputAOD) {
1980         Double_t covtest[21];
1981         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1982           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1983           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1984             id = (Int_t)vtrack->GetID();
1985             if(id<0) continue;
1986             skipped[nTrksToSkip++] = id;
1987           }
1988         }
1989       }
1990       for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
1991       //
1992       vertexer->SetSkipTracks(nTrksToSkip,skipped);
1993       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
1994       
1995     } else if(fRmTrksFromPrimVtx && nTrks>0) { 
1996       // removing the prongs tracks
1997       
1998       TObjArray rmArray(nTrks);
1999       UShort_t *rmId = new UShort_t[nTrks];
2000       AliESDtrack *esdTrack = 0;
2001       AliESDtrack *t = 0;
2002       for(Int_t i=0; i<nTrks; i++) {
2003         t = (AliESDtrack*)trkArray->UncheckedAt(i);
2004         esdTrack = new AliESDtrack(*t);
2005         rmArray.AddLast(esdTrack);
2006         if(esdTrack->GetID()>=0) {
2007           rmId[i]=(UShort_t)esdTrack->GetID();
2008         } else {
2009           rmId[i]=9999;
2010         }
2011       }
2012       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
2013       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2014       delete [] rmId; rmId=NULL;
2015       rmArray.Delete();
2016       
2017     }
2018
2019     if(!vertexESD) return vertexAOD;
2020     if(vertexESD->GetNContributors()<=0) { 
2021       //AliDebug(2,"vertexing failed"); 
2022       delete vertexESD; vertexESD=NULL;
2023       return vertexAOD;
2024     }
2025
2026     delete vertexer; vertexer=NULL;
2027
2028   }
2029
2030   // convert to AliAODVertex
2031   Double_t pos[3],cov[6],chi2perNDF;
2032   vertexESD->GetXYZ(pos); // position
2033   vertexESD->GetCovMatrix(cov); //covariance matrix
2034   chi2perNDF = vertexESD->GetChi2toNDF();
2035   delete vertexESD; vertexESD=NULL;
2036
2037   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2038
2039   return vertexAOD;
2040 }
2041 //-----------------------------------------------------------------------------
2042 void AliAnalysisVertexingHF::PrintStatus() const {
2043   // Print parameters being used
2044
2045   //printf("Preselections:\n");
2046   //   fTrackFilter->Dump();
2047   if(fSecVtxWithKF) {
2048     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2049   } else {
2050     printf("Secondary vertex with AliVertexerTracks\n");
2051   }
2052   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2053   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2054   if(fD0toKpi) {
2055     printf("Reconstruct D0->Kpi candidates with cuts:\n");
2056     if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
2057   }
2058   if(fDstar) {
2059     printf("Reconstruct D*->D0pi candidates with cuts:\n");
2060     if(fFindVertexForDstar) {
2061       printf("    Reconstruct a secondary vertex for the D*\n");
2062     } else {
2063       printf("    Assume the D* comes from the primary vertex\n");
2064     }
2065     if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
2066   }
2067   if(fJPSItoEle) {
2068     printf("Reconstruct J/psi from B candidates with cuts:\n");
2069     if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
2070   }
2071   if(f3Prong) {
2072     printf("Reconstruct 3 prong candidates.\n");
2073     printf("  D+->Kpipi cuts:\n");
2074     if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
2075     printf("  Ds->KKpi cuts:\n");
2076     if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
2077     printf("  Lc->pKpi cuts:\n");
2078     if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
2079   }
2080   if(f4Prong) {
2081     printf("Reconstruct 4 prong candidates.\n");
2082     printf("  D0->Kpipipi cuts:\n");
2083     if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
2084   }
2085   if(fCascades) {
2086     printf("Reconstruct cascades candidates formed with v0s.\n");
2087     printf("  Lc -> k0s P & Lc -> L Pi cuts:\n");
2088     if(fCutsLctoV0) fCutsLctoV0->PrintAll();
2089   }
2090
2091   return;
2092 }
2093 //-----------------------------------------------------------------------------
2094 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
2095                                                                  Double_t &dispersion,Bool_t useTRefArray) const
2096 {
2097   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
2098   //AliCodeTimerAuto("",0);
2099
2100   AliESDVertex *vertexESD = 0;
2101   AliAODVertex *vertexAOD = 0;
2102
2103   if(!fSecVtxWithKF) { // AliVertexerTracks
2104
2105     fVertexerTracks->SetVtxStart(fV1);
2106     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2107
2108     if(!vertexESD) return vertexAOD;
2109
2110     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
2111       //AliDebug(2,"vertexing failed"); 
2112       delete vertexESD; vertexESD=NULL;
2113       return vertexAOD;
2114     }
2115     
2116     Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
2117     if(vertRadius2>8.){
2118       // vertex outside beam pipe, reject candidate to avoid propagation through material
2119       delete vertexESD; vertexESD=NULL;
2120       return vertexAOD;
2121     }
2122
2123   } else { // Kalman Filter vertexer (AliKFParticle)
2124
2125     AliKFParticle::SetField(fBzkG);
2126
2127     AliKFVertex vertexKF;
2128
2129     Int_t nTrks = trkArray->GetEntriesFast();
2130     for(Int_t i=0; i<nTrks; i++) {
2131       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2132       AliKFParticle daughterKF(*esdTrack,211);
2133       vertexKF.AddDaughter(daughterKF);
2134     }
2135     vertexESD = new AliESDVertex(vertexKF.Parameters(),
2136                                  vertexKF.CovarianceMatrix(),
2137                                  vertexKF.GetChi2(),
2138                                  vertexKF.GetNContributors());
2139
2140   }
2141
2142   // convert to AliAODVertex
2143   Double_t pos[3],cov[6],chi2perNDF;
2144   vertexESD->GetXYZ(pos); // position
2145   vertexESD->GetCovMatrix(cov); //covariance matrix
2146   chi2perNDF = vertexESD->GetChi2toNDF();
2147   dispersion = vertexESD->GetDispersion();
2148   delete vertexESD; vertexESD=NULL;
2149
2150   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2151   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2152
2153   return vertexAOD;
2154 }
2155 //-----------------------------------------------------------------------------
2156 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(TObjArray *trkArray){
2157   // Invariant mass cut on tracks
2158   //AliCodeTimerAuto("",0);
2159
2160   Int_t retval=kFALSE;
2161   Double_t momentum[3];
2162   Double_t px[3],py[3],pz[3];
2163   for(Int_t iTrack=0; iTrack<3; iTrack++){
2164     AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2165     track->GetPxPyPz(momentum);
2166     px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2]; 
2167   }
2168   retval = SelectInvMassAndPt3prong(px,py,pz);
2169
2170   return retval;
2171 }
2172
2173 //-----------------------------------------------------------------------------
2174 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(TObjArray *trkArray){
2175   // Invariant mass cut on tracks
2176   //AliCodeTimerAuto("",0);
2177
2178   Int_t retval=kFALSE;
2179   Double_t momentum[3];
2180   Double_t px[4],py[4],pz[4];
2181
2182   for(Int_t iTrack=0; iTrack<4; iTrack++){
2183     AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2184     track->GetPxPyPz(momentum);
2185     px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2]; 
2186   }
2187
2188   retval = SelectInvMassAndPt4prong(px,py,pz);
2189
2190   return retval;
2191 }
2192 //-----------------------------------------------------------------------------
2193 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(TObjArray *trkArray){
2194   // Invariant mass cut on tracks
2195   //AliCodeTimerAuto("",0);
2196
2197   Int_t retval=kFALSE;
2198   Double_t momentum[3];
2199   Double_t px[2],py[2],pz[2];
2200
2201   for(Int_t iTrack=0; iTrack<2; iTrack++){
2202     AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2203     track->GetPxPyPz(momentum);
2204     px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2]; 
2205   }
2206   retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
2207
2208   return retval;
2209 }
2210 //-----------------------------------------------------------------------------
2211 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtD0Kpi(Double_t *px,
2212                                                        Double_t *py,
2213                                                        Double_t *pz){
2214   // Check invariant mass cut and pt candidate cut
2215   //AliCodeTimerAuto("",0);
2216
2217   UInt_t pdg2[2];
2218   Int_t nprongs=2;
2219   Double_t minv2,mrange;
2220   Double_t lolim,hilim;
2221   Double_t minPt=0;
2222   Bool_t retval=kFALSE;
2223
2224   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2225   fOKInvMassD0=kFALSE;
2226   // pt cut
2227   minPt=fCutsD0toKpi->GetMinPtCandidate();
2228   if(minPt>0.1) 
2229     if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2230   // mass cut
2231   mrange=fCutsD0toKpi->GetMassCut();
2232   lolim=fMassDzero-mrange;
2233   hilim=fMassDzero+mrange;
2234   pdg2[0]=211; pdg2[1]=321;
2235   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2236   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2237     retval=kTRUE;
2238     fOKInvMassD0=kTRUE;
2239   }
2240   pdg2[0]=321; pdg2[1]=211;
2241   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2242   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2243     retval=kTRUE;
2244     fOKInvMassD0=kTRUE;
2245   }
2246   return retval;
2247 }
2248
2249 //-----------------------------------------------------------------------------
2250 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtJpsiee(Double_t *px,
2251                                                         Double_t *py,
2252                                                         Double_t *pz){
2253   // Check invariant mass cut and pt candidate cut
2254   //AliCodeTimerAuto("",0);
2255
2256   UInt_t pdg2[2];
2257   Int_t nprongs=2;
2258   Double_t minv2,mrange;
2259   Double_t lolim,hilim;
2260   Double_t minPt=0;
2261   Bool_t retval=kFALSE;
2262
2263   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2264   fOKInvMassJpsi=kFALSE;
2265   // pt cut
2266   minPt=fCutsJpsitoee->GetMinPtCandidate();
2267   if(minPt>0.1) 
2268     if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2269   // mass cut
2270   mrange=fCutsJpsitoee->GetMassCut();
2271   lolim=fMassJpsi-mrange;
2272   hilim=fMassJpsi+mrange;
2273
2274   pdg2[0]=11; pdg2[1]=11;
2275   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2276   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2277     retval=kTRUE;
2278     fOKInvMassJpsi=kTRUE;
2279   }
2280
2281   return retval;
2282 }
2283 //-----------------------------------------------------------------------------
2284 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(Double_t *px,
2285                                                         Double_t *py,
2286                                                         Double_t *pz){
2287   // Check invariant mass cut and pt candidate cut
2288   //AliCodeTimerAuto("",0);
2289
2290   UInt_t pdg3[3];
2291   Int_t nprongs=3;
2292   Double_t minv2,mrange;
2293   Double_t lolim,hilim;
2294   Double_t minPt=0;
2295   Bool_t retval=kFALSE;
2296
2297
2298   fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2299   fOKInvMassDplus=kFALSE;
2300   fOKInvMassDs=kFALSE;
2301   fOKInvMassLc=kFALSE;
2302   // pt cut
2303   minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
2304   minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
2305   if(minPt>0.1) 
2306     if(fMassCalc3->Pt2() < minPt*minPt) return retval;
2307   // D+->Kpipi
2308   mrange=fCutsDplustoKpipi->GetMassCut();
2309   lolim=fMassDplus-mrange;
2310   hilim=fMassDplus+mrange;
2311   pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2312   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2313   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2314     retval=kTRUE;
2315     fOKInvMassDplus=kTRUE;
2316   }
2317   // Ds+->KKpi
2318   mrange=fCutsDstoKKpi->GetMassCut();
2319   lolim=fMassDs-mrange;
2320   hilim=fMassDs+mrange;
2321   pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2322   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2323   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2324     retval=kTRUE;
2325     fOKInvMassDs=kTRUE;
2326   }
2327   pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2328   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2329   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2330     retval=kTRUE;
2331     fOKInvMassDs=kTRUE;
2332   }
2333   // Lc->pKpi
2334   mrange=fCutsLctopKpi->GetMassCut();
2335   lolim=fMassLambdaC-mrange;
2336   hilim=fMassLambdaC+mrange;
2337   pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2338   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2339   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2340     retval=kTRUE;
2341     fOKInvMassLc=kTRUE;
2342   }
2343   pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2344   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2345   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2346     retval=kTRUE;
2347     fOKInvMassLc=kTRUE;
2348   }
2349
2350   return retval;
2351 }
2352
2353 //-----------------------------------------------------------------------------
2354 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(Double_t *px,
2355                                                            Double_t *py,
2356                                                            Double_t *pz){
2357   // Check invariant mass cut and pt candidate cut
2358   //AliCodeTimerAuto("",0);
2359
2360   UInt_t pdg2[2];
2361   Int_t nprongs=2;
2362   Double_t minv2,mrange;
2363   Double_t lolim,hilim;
2364   Double_t minPt=0;
2365   Bool_t retval=kFALSE;
2366
2367   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2368   fOKInvMassDstar=kFALSE;
2369   // pt cut
2370   minPt=fCutsDStartoKpipi->GetMinPtCandidate();
2371   if(minPt>0.1) 
2372     if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2373   // mass cut
2374   pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2375   mrange=fCutsDStartoKpipi->GetMassCut();
2376   lolim=fMassDstar-mrange;
2377   hilim=fMassDstar+mrange;
2378   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2379   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2380     retval=kTRUE;
2381     fOKInvMassDstar=kTRUE;
2382   }
2383
2384   return retval;
2385 }
2386
2387 //-----------------------------------------------------------------------------
2388 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(Double_t *px,
2389                                                         Double_t *py,
2390                                                         Double_t *pz){
2391   // Check invariant mass cut and pt candidate cut
2392   //AliCodeTimerAuto("",0);
2393
2394   UInt_t pdg4[4];
2395   Int_t nprongs=4;
2396   Double_t minv2,mrange;
2397   Double_t lolim,hilim;
2398   Double_t minPt=0;
2399   Bool_t retval=kFALSE;
2400
2401   // D0->Kpipipi without PID
2402   fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2403   fOKInvMassD0to4p=kFALSE;
2404   // pt cut
2405   minPt=fCutsD0toKpipipi->GetMinPtCandidate();
2406   if(minPt>0.1) 
2407     if(fMassCalc4->Pt2() < minPt*minPt) return retval;
2408   // mass cut
2409   mrange=fCutsD0toKpipipi->GetMassCut();
2410   lolim=fMassDzero-mrange;
2411   hilim=fMassDzero+mrange;
2412
2413   pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2414   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2415   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2416     retval=kTRUE;
2417     fOKInvMassD0to4p=kTRUE;
2418   }
2419
2420   pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2421   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2422   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2423     retval=kTRUE;
2424     fOKInvMassD0to4p=kTRUE;
2425   }
2426
2427   pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2428   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2429   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2430     retval=kTRUE;
2431     fOKInvMassD0to4p=kTRUE;
2432   }
2433
2434   pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2435   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2436   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2437     retval=kTRUE;
2438     fOKInvMassD0to4p=kTRUE;
2439   }
2440
2441   return retval;
2442 }
2443 //-----------------------------------------------------------------------------
2444 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtCascade(Double_t *px,
2445                                                          Double_t *py,
2446                                                          Double_t *pz){
2447   // Check invariant mass cut and pt candidate cut
2448   //AliCodeTimerAuto("",0);
2449
2450   UInt_t pdg2[2];
2451   Int_t nprongs=2;
2452   Double_t minv2,mrange;
2453   Double_t lolim,hilim;
2454   Double_t minPt=0;
2455   Bool_t retval=kFALSE;
2456   
2457   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2458   minPt=fCutsLctoV0->GetMinPtCandidate();
2459   fOKInvMassLctoV0=kFALSE; 
2460   mrange=fCutsLctoV0->GetMassCut();
2461   lolim=fMassLambdaC-mrange;
2462   hilim=fMassLambdaC+mrange;
2463   pdg2[0]=2212;pdg2[1]=310;
2464   minv2=fMassCalc2->InvMass2(2,pdg2);
2465   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2466     retval=kTRUE;
2467     fOKInvMassLctoV0=kTRUE;
2468   }
2469   pdg2[0]=211;pdg2[1]=3122;
2470   minv2=fMassCalc2->InvMass2(2,pdg2);
2471   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2472     retval=kTRUE;
2473     fOKInvMassLctoV0=kTRUE;
2474   }
2475
2476   return retval;
2477 }
2478 //-----------------------------------------------------------------------------
2479 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2480                                                        Int_t trkEntries,
2481                                                        TObjArray &seleTrksArray,
2482                                                        TObjArray &tracksAtVertex,
2483                                                        Int_t &nSeleTrks,
2484                                                        UChar_t *seleFlags,Int_t *evtNumber)
2485 {
2486   // Apply single-track preselection.
2487   // Fill a TObjArray with selected tracks (for displaced vertices or
2488   // soft pion from D*). Selection flag stored in seleFlags.
2489   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2490   // In case of AOD input, also fill fAODMap for track index<->ID
2491   //AliCodeTimerAuto("",0);
2492
2493   const AliVVertex *vprimary = event->GetPrimaryVertex();
2494
2495   if(fV1) { delete fV1; fV1=NULL; }
2496   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
2497
2498   Int_t nindices=0;
2499   UShort_t *indices = 0;
2500   Double_t pos[3],cov[6];
2501   const Int_t entries = event->GetNumberOfTracks();
2502   AliCentrality* cent;
2503  
2504   if(!fInputAOD) { // ESD
2505     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2506     cent=((AliESDEvent*)event)->GetCentrality();
2507  } else {         // AOD
2508     vprimary->GetXYZ(pos);
2509     vprimary->GetCovarianceMatrix(cov);
2510     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2511     if(entries<=0) return;
2512     indices = new UShort_t[entries];
2513     memset(indices,0,sizeof(UShort_t)*entries);
2514     fAODMapSize = 100000;
2515     fAODMap = new Int_t[fAODMapSize];
2516     memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2517     cent=((AliAODEvent*)event)->GetCentrality();
2518   }
2519   Float_t centperc=cent->GetCentralityPercentile("V0M");
2520
2521   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE;
2522   nSeleTrks=0;
2523  
2524   // transfer ITS tracks from event to arrays
2525   for(Int_t i=0; i<entries; i++) {
2526     AliVTrack *track;
2527     track = (AliVTrack*)event->GetTrack(i);
2528
2529     // skip pure ITS SA tracks
2530     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2531
2532     // skip tracks without ITS
2533     if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2534
2535     // skip tracks with negative ID 
2536     // (these are duplicated TPC-only AOD tracks, for jet analysis...)
2537     if(track->GetID()<0) continue;
2538
2539     // TEMPORARY: check that the cov matrix is there
2540     Double_t covtest[21];
2541     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2542     //
2543
2544     if(fInputAOD) {
2545       AliAODTrack *aodt = (AliAODTrack*)track;
2546       if(aodt->GetUsedForPrimVtxFit()) { 
2547         indices[nindices]=aodt->GetID(); nindices++; 
2548       }
2549       Int_t ind = (Int_t)aodt->GetID();
2550       if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2551     }
2552
2553     AliESDtrack *esdt = 0;
2554
2555     if(!fInputAOD) {
2556       esdt = (AliESDtrack*)track;
2557     } else {
2558       esdt = new AliESDtrack(track);
2559     }
2560
2561     // single track selection
2562     okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE;
2563     if(fMixEvent && i<trkEntries){
2564       evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2565       const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2566       Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2567       eventVtx->GetXYZ(vtxPos);
2568       vprimary->GetXYZ(primPos);
2569       eventVtx->GetCovarianceMatrix(primCov);
2570       for(Int_t ind=0;ind<3;ind++){
2571         trasl[ind]=vtxPos[ind]-primPos[ind];
2572       }
2573       
2574       Bool_t isTransl=esdt->Translate(trasl,primCov);
2575       if(!isTransl) {
2576         delete esdt;
2577         esdt = NULL;
2578         continue;
2579       }
2580     }
2581
2582     if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong) && nSeleTrks<trkEntries) {
2583       esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
2584       seleTrksArray.AddLast(esdt);
2585       tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
2586       seleFlags[nSeleTrks]=0;
2587       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2588       if(okFor3Prong)    SETBIT(seleFlags[nSeleTrks],kBit3Prong);
2589       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2590       nSeleTrks++;
2591     } else {
2592       if(fInputAOD) delete esdt; 
2593       esdt = NULL;
2594       continue;
2595     } 
2596
2597   } // end loop on tracks
2598
2599   // primary vertex from AOD
2600   if(fInputAOD) {
2601     delete fV1; fV1=NULL;
2602     vprimary->GetXYZ(pos);
2603     vprimary->GetCovarianceMatrix(cov);
2604     Double_t chi2toNDF = vprimary->GetChi2perNDF();
2605     Int_t ncontr=nindices;
2606     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2607     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.); 
2608     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2609     fV1->SetTitle(vprimary->GetTitle());
2610     fV1->SetIndices(nindices,indices);
2611   }
2612   if(indices) { delete [] indices; indices=NULL; }
2613
2614
2615   return;
2616 }
2617 //-----------------------------------------------------------------------------
2618 void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
2619   //
2620   // Set the selection bit for PID
2621   //
2622   //AliCodeTimerAuto("",0);
2623   if(cuts->GetPidHF()) {
2624     Bool_t usepid=cuts->GetIsUsePID();
2625     cuts->SetUsePID(kTRUE);
2626     if(cuts->IsSelectedPID(rd))
2627       rd->SetSelectionBit(bit);
2628     cuts->SetUsePID(usepid);
2629   }
2630   return;
2631 }
2632 //-----------------------------------------------------------------------------
2633 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk, 
2634                                              Float_t centralityperc,
2635                                              Bool_t &okDisplaced,
2636                                              Bool_t &okSoftPi,
2637                                              Bool_t &okFor3Prong) const 
2638 {
2639   // Check if track passes some kinematical cuts  
2640
2641   // this is needed to store the impact parameters
2642   //AliCodeTimerAuto("",0);
2643
2644   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2645
2646   UInt_t selectInfo;
2647   //
2648   // Track selection, displaced tracks -- 2 prong
2649   selectInfo = 0; 
2650   if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0 
2651      && fTrackFilter2prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2652     // central PbPb events, tighter cuts
2653     selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
2654   }else{
2655     // standard cuts
2656     if(fTrackFilter) {
2657       selectInfo = fTrackFilter->IsSelected(trk);
2658     }
2659   }
2660   if(selectInfo) okDisplaced=kTRUE;
2661
2662   // Track selection, displaced tracks -- 3 prong
2663   selectInfo = 0; 
2664   if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0 
2665      && fTrackFilter3prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2666     // central PbPb events, tighter cuts
2667     selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
2668   }else{
2669     // standard cuts
2670     if(fTrackFilter) {
2671       selectInfo = fTrackFilter->IsSelected(trk);
2672     }
2673   }
2674   if(selectInfo) okFor3Prong=kTRUE;
2675
2676   // Track selection, soft pions
2677   selectInfo = 0; 
2678   if(fDstar && fTrackFilterSoftPi) {
2679     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2680   }
2681   if(selectInfo) okSoftPi=kTRUE;
2682
2683   if(okDisplaced || okSoftPi || okFor3Prong) return kTRUE;
2684
2685   return kFALSE;
2686 }
2687
2688
2689 //-----------------------------------------------------------------------------
2690 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2691   //
2692   // Transform ESDv0 to AODv0
2693   //
2694   //  this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2695   //  and creates an AODv0 out of them
2696   //
2697   //AliCodeTimerAuto("",0);
2698   Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2699   AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2700
2701   // create the v0 neutral track to compute the DCA to the primary vertex
2702   Double_t xyz[3], pxpypz[3];
2703   esdV0->XvYvZv(xyz);
2704   esdV0->PxPyPz(pxpypz);
2705   Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2706   AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2707   if(!trackesdV0) {
2708     delete vertexV0;
2709     return 0;
2710   }
2711   Double_t d0z0[2],covd0z0[3];
2712   AliAODVertex *primVertexAOD = PrimaryVertex(); 
2713   trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2714   Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2715   // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2716   Double_t dcaV0DaughterToPrimVertex[2];  
2717   AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2718   AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2719   if( !posV0track || !negV0track) {
2720     if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2721     delete vertexV0;
2722     delete primVertexAOD;
2723     return 0;
2724   }
2725   posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2726   //  if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2727   //  else 
2728   dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2729   negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);  
2730   //  if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2731   //  else 
2732   dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2733   Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2734   Double_t pmom[3],nmom[3];
2735   esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2736   esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2737
2738   AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2739   aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2740
2741   delete trackesdV0;
2742   delete primVertexAOD;
2743
2744   return aodV0;
2745 }
2746 //-----------------------------------------------------------------------------
2747 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
2748   // Set the stored track parameters at primary vertex into AliESDtrack
2749   //AliCodeTimerAuto("",0);
2750
2751   const Double_t *par=extpar->GetParameter();
2752   const Double_t *cov=extpar->GetCovariance();
2753   Double_t alpha=extpar->GetAlpha();
2754   Double_t x=extpar->GetX();
2755   esdt->Set(x,alpha,par,cov);
2756   return;
2757 }
2758 //-----------------------------------------------------------------------------
2759 void AliAnalysisVertexingHF::SetMasses(){
2760   // Set the hadron mass values from TDatabasePDG 
2761
2762   fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2763   fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2764   fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2765   fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2766   fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2767   fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2768 }