790033dfc17ad77aa4d5246d75a261e1d6b75e8f
[u/mrichter/AliRoot.git] / PWG3 / 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 //----------------------------------------------------------------------------
17 //    Implementation of the heavy-flavour vertexing analysis class
18 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
19 // To be used as a task of AliAnalysisManager by means of the interface
20 // class AliAnalysisTaskSEVertexingHF. 
21 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
22 //
23 //  Contact: andrea.dainese@pd.infn.it
24 //  Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
25 //                F.Prino, R.Romita, X.M.Zhang
26 //----------------------------------------------------------------------------
27 #include <TFile.h>
28 #include <TDatabasePDG.h>
29 #include <TString.h>
30 #include "AliLog.h"
31 #include "AliVEvent.h"
32 #include "AliVVertex.h"
33 #include "AliVTrack.h"
34 #include "AliVertexerTracks.h"
35 #include "AliKFVertex.h"
36 #include "AliESDEvent.h"
37 #include "AliESDVertex.h"
38 #include "AliExternalTrackParam.h"
39 #include "AliNeutralTrackParam.h"
40 #include "AliESDtrack.h"
41 #include "AliESDtrackCuts.h"
42 #include "AliAODEvent.h"
43 #include "AliAODRecoDecay.h"
44 #include "AliAODRecoDecayHF.h"
45 #include "AliAODRecoDecayHF2Prong.h"
46 #include "AliAODRecoDecayHF3Prong.h"
47 #include "AliAODRecoDecayHF4Prong.h"
48 #include "AliAODRecoCascadeHF.h"
49 #include "AliRDHFCutsD0toKpi.h"
50 #include "AliRDHFCutsJpsitoee.h"
51 #include "AliRDHFCutsDplustoKpipi.h"
52 #include "AliRDHFCutsDstoKKpi.h"
53 #include "AliRDHFCutsLctopKpi.h"
54 #include "AliRDHFCutsD0toKpipipi.h"
55 #include "AliAnalysisFilter.h"
56 #include "AliAnalysisVertexingHF.h"
57 #include "AliMixedEvent.h"
58
59 ClassImp(AliAnalysisVertexingHF)
60
61 //----------------------------------------------------------------------------
62 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
63 fInputAOD(kFALSE),
64 fAODMapSize(0),
65 fAODMap(0),
66 fBzkG(0.),
67 fSecVtxWithKF(kFALSE),
68 fRecoPrimVtxSkippingTrks(kFALSE),
69 fRmTrksFromPrimVtx(kFALSE),
70 fV1(0x0),
71 fD0toKpi(kTRUE),
72 fJPSItoEle(kTRUE),
73 f3Prong(kTRUE),
74 f4Prong(kTRUE),
75 fDstar(kTRUE),
76 fLikeSign(kFALSE),
77 fMixEvent(kFALSE),
78 fTrackFilter(0x0),
79 fTrackFilterSoftPi(0x0),
80 fCutsD0toKpi(0x0),
81 fCutsJpsitoee(0x0),
82 fCutsDplustoKpipi(0x0),
83 fCutsDstoKKpi(0x0),
84 fCutsLctopKpi(0x0),
85 fCutsD0toKpipipi(0x0),
86 fCutsD0fromDstar(0x0),
87 fFindVertexForDstar(kTRUE)
88 {
89   // Default constructor
90
91   SetD0toKpiCuts();
92   SetBtoJPSICuts();
93   SetDplusCuts();
94   SetDsCuts();
95   SetLcCuts();
96   SetDstarCuts();
97   SetD0to4ProngsCuts();
98 }
99 //--------------------------------------------------------------------------
100 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) : 
101 TNamed(source),
102 fInputAOD(source.fInputAOD),
103 fAODMapSize(source.fAODMapSize),
104 fAODMap(source.fAODMap),
105 fBzkG(source.fBzkG),
106 fSecVtxWithKF(source.fSecVtxWithKF),
107 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
108 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
109 fV1(source.fV1),
110 fD0toKpi(source.fD0toKpi),
111 fJPSItoEle(source.fJPSItoEle),
112 f3Prong(source.f3Prong),
113 f4Prong(source.f4Prong),
114 fDstar(source.fDstar),
115 fLikeSign(source.fLikeSign),
116 fMixEvent(source.fMixEvent),
117 fTrackFilter(source.fTrackFilter),
118 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
119 fCutsD0toKpi(source.fCutsD0toKpi),
120 fCutsJpsitoee(source.fCutsJpsitoee),
121 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
122 fCutsDstoKKpi(source.fCutsDstoKKpi),
123 fCutsLctopKpi(source.fCutsLctopKpi),
124 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
125 fCutsD0fromDstar(source.fCutsD0fromDstar),
126 fFindVertexForDstar(source.fFindVertexForDstar)
127 {
128   //
129   // Copy constructor
130   //
131   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
132   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
133   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
134   for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
135   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
136   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
137   for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
138 }
139 //--------------------------------------------------------------------------
140 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
141 {
142   //
143   // assignment operator
144   //
145   if(&source == this) return *this;
146   fInputAOD = source.fInputAOD;
147   fAODMapSize = source.fAODMapSize;
148   fBzkG = source.fBzkG;
149   fSecVtxWithKF = source.fSecVtxWithKF;
150   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
151   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
152   fV1 = source.fV1;
153   fD0toKpi = source.fD0toKpi;
154   fJPSItoEle = source.fJPSItoEle;
155   f3Prong = source.f3Prong;
156   f4Prong = source.f4Prong;
157   fDstar = source.fDstar;
158   fLikeSign = source.fLikeSign;
159   fMixEvent = source.fMixEvent;
160   fTrackFilter = source.fTrackFilter;
161   fTrackFilterSoftPi = source.fTrackFilterSoftPi;
162   fCutsD0toKpi = source.fCutsD0toKpi;
163   fCutsJpsitoee = source.fCutsJpsitoee;
164   fCutsDplustoKpipi = source.fCutsDplustoKpipi;
165   fCutsDstoKKpi = source.fCutsDstoKKpi;
166   fCutsLctopKpi = source.fCutsLctopKpi;
167   fCutsD0toKpipipi = source.fCutsD0toKpipipi;
168   fCutsD0fromDstar = source.fCutsD0fromDstar;
169   fFindVertexForDstar = source.fFindVertexForDstar;
170
171   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
172   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
173   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
174   for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
175   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
176   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
177   for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
178
179   return *this;
180 }
181 //----------------------------------------------------------------------------
182 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
183   // Destructor
184   if(fV1) { delete fV1; fV1=0; }
185   if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
186   if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
187   if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
188   if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
189   if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
190   if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
191   if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
192   if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
193   if(fCutsD0fromDstar) { delete fCutsD0fromDstar; fCutsD0fromDstar=0; }
194   if(fAODMap) { delete fAODMap; fAODMap=0; }
195 }
196 //----------------------------------------------------------------------------
197 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
198                                             TClonesArray *aodVerticesHFTClArr,
199                                             TClonesArray *aodD0toKpiTClArr,
200                                             TClonesArray *aodJPSItoEleTClArr,
201                                             TClonesArray *aodCharm3ProngTClArr,
202                                             TClonesArray *aodCharm4ProngTClArr,
203                                             TClonesArray *aodDstarTClArr,
204                                             TClonesArray *aodLikeSign2ProngTClArr,
205                                             TClonesArray *aodLikeSign3ProngTClArr)
206 {
207   // Find heavy-flavour vertex candidates
208   // Input:  ESD or AOD
209   // Output: AOD (additional branches added)
210
211   if(!fMixEvent){
212     TString evtype = event->IsA()->GetName();
213     fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
214   } // if we do mixing AliVEvent is a AliMixedEvent
215
216   if(fInputAOD) {
217     AliDebug(2,"Creating HF candidates from AOD");
218   } else {
219     AliDebug(2,"Creating HF candidates from ESD");
220   }
221
222   if(!aodVerticesHFTClArr) {
223     printf("ERROR: no aodVerticesHFTClArr");
224     return;
225   }
226   if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
227     printf("ERROR: no aodD0toKpiTClArr");
228     return;
229   }
230   if(fJPSItoEle && !aodJPSItoEleTClArr) {
231     printf("ERROR: no aodJPSItoEleTClArr");
232     return;
233   }
234   if(f3Prong && !aodCharm3ProngTClArr) {
235     printf("ERROR: no aodCharm3ProngTClArr");
236     return;
237   }
238   if(f4Prong && !aodCharm4ProngTClArr) {
239     printf("ERROR: no aodCharm4ProngTClArr");
240     return;
241   }
242   if(fDstar && !aodDstarTClArr) {
243     printf("ERROR: no aodDstarTClArr");
244     return;
245   }
246   if(fLikeSign && !aodLikeSign2ProngTClArr) {
247     printf("ERROR: no aodLikeSign2ProngTClArr");
248     return;
249   }
250   if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
251     printf("ERROR: no aodLikeSign2ProngTClArr");
252     return;
253   }
254
255   // delete candidates from previous event and create references
256   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
257   aodVerticesHFTClArr->Delete();
258   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
259   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
260   if(fD0toKpi || fDstar)   {
261     aodD0toKpiTClArr->Delete();
262     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
263   }
264   if(fJPSItoEle) {
265     aodJPSItoEleTClArr->Delete();
266     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
267   }
268   if(f3Prong) {   
269     aodCharm3ProngTClArr->Delete();
270     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
271   }
272   if(f4Prong) {
273     aodCharm4ProngTClArr->Delete();
274     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
275   }
276   if(fDstar) {
277     aodDstarTClArr->Delete();
278     iDstar = aodDstarTClArr->GetEntriesFast();
279   }
280   if(fLikeSign) {                                
281     aodLikeSign2ProngTClArr->Delete();                     
282     iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast(); 
283   }  
284   if(fLikeSign && f3Prong) {                                
285     aodLikeSign3ProngTClArr->Delete();                     
286     iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast(); 
287   }  
288
289   TClonesArray &aodD0toKpiRef        = *aodD0toKpiTClArr;
290   TClonesArray &aodJPSItoEleRef      = *aodJPSItoEleTClArr;
291   TClonesArray &aodCharm3ProngRef    = *aodCharm3ProngTClArr;
292   TClonesArray &aodCharm4ProngRef    = *aodCharm4ProngTClArr;
293   TClonesArray &aodDstarRef          = *aodDstarTClArr;
294   TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
295   TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
296
297
298   AliAODRecoDecayHF2Prong *io2Prong  = 0;
299   AliAODRecoDecayHF3Prong *io3Prong  = 0;
300   AliAODRecoDecayHF4Prong *io4Prong  = 0;
301   AliAODRecoCascadeHF     *ioCascade = 0;
302
303   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries;
304   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaCasc;
305   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
306   Bool_t   okDstar=kFALSE,okD0fromDstar=kFALSE;
307   AliESDtrack *postrack1 = 0;
308   AliESDtrack *postrack2 = 0;
309   AliESDtrack *negtrack1 = 0;
310   AliESDtrack *negtrack2 = 0;
311   AliESDtrack *trackPi   = 0;
312   Double_t dcaMax = fD0toKpiCuts[1];
313   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
314   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
315   if(dcaMax < fD0to4ProngsCuts[1])  dcaMax=fD0to4ProngsCuts[1];
316   /*
317   Double_t dcaMax = fCutsD0toKpi->GetDCACut();
318   if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
319   if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
320   if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
321   if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
322   if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
323   */
324   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
325
326
327   // get Bz
328   fBzkG = (Double_t)event->GetMagneticField(); 
329
330   trkEntries = (Int_t)event->GetNumberOfTracks();
331   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
332
333   if(trkEntries<2) {
334     AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
335     return;
336   }
337     
338   // event selection
339   //if(!fCutsD0toKpi->IsEventSelected(event)) return;
340   const AliVVertex *primary = event->GetPrimaryVertex();
341   if(!primary) {
342     AliDebug(1," No primary vertex from tracks");
343     return;
344   }
345
346   TString primTitle=primary->GetTitle();
347   if(!primTitle.Contains("VertexerTracks") ||
348      primary->GetNContributors()<=0) {
349     AliDebug(1," No primary vertex from tracks");
350     return;
351   }
352   
353
354   // call function that applies sigle-track selection,
355   // for displaced tracks and soft pions (both charges) for D*,
356   // and retrieves primary vertex
357   TObjArray seleTrksArray(trkEntries);
358   UChar_t  *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
359   Int_t     nSeleTrks=0;
360   Int_t *evtNumber    = new Int_t[trkEntries];
361   SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
362     
363   AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
364     
365   TObjArray *twoTrackArray1    = new TObjArray(2);
366   TObjArray *twoTrackArray2    = new TObjArray(2);
367   TObjArray *twoTrackArrayCasc = new TObjArray(2);
368   TObjArray *threeTrackArray   = new TObjArray(3);
369   TObjArray *fourTrackArray    = new TObjArray(4);
370   
371   Double_t dispersion;
372   Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
373
374   AliAODRecoDecayHF   *rd = 0;
375   AliAODRecoCascadeHF *rc = 0;
376
377   // LOOP ON  POSITIVE  TRACKS
378   for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
379
380     if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));  
381
382     // get track from tracks array
383     postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
384
385     if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
386
387     if(postrack1->Charge()<0 && !fLikeSign) continue;
388
389     // LOOP ON  NEGATIVE  TRACKS
390     for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
391
392       if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));  
393
394       if(iTrkN1==iTrkP1) continue;
395
396       // get track from tracks array
397       negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
398
399       if(negtrack1->Charge()>0 && !fLikeSign) continue;
400
401       if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
402
403       if(fMixEvent) {
404         if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
405       }
406
407       if(postrack1->Charge()==negtrack1->Charge()) { // like-sign 
408         isLikeSign2Prong=kTRUE;
409         if(!fLikeSign)    continue;
410         if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
411       } else { // unlike-sign
412         isLikeSign2Prong=kFALSE;
413         if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue;  // this is needed to avoid double-counting of unlike-sign
414         if(fMixEvent) {
415           if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
416         }
417        
418       }
419
420       // back to primary vertex
421       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
422       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
423
424       // DCA between the two tracks
425       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
426       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
427
428       // Vertexing
429       twoTrackArray1->AddAt(postrack1,0);
430       twoTrackArray1->AddAt(negtrack1,1);
431       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
432       if(!vertexp1n1) { 
433         twoTrackArray1->Clear();
434         negtrack1=0; 
435         continue; 
436       }
437
438       // 2 prong candidate
439       if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) { 
440       
441         io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
442       
443         if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
444           // add the vertex and the decay to the AOD
445           AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
446           if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
447           if(!isLikeSign2Prong) {
448             if(okD0) {  
449               rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
450               rd->SetSecondaryVtx(v2Prong);
451               v2Prong->SetParent(rd);
452               if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
453             }
454             if(okJPSI) {
455               rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
456               rd->SetSecondaryVtx(v2Prong);
457               if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
458               if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
459             }
460           } else { // isLikeSign2Prong
461             rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
462             rd->SetSecondaryVtx(v2Prong);
463             v2Prong->SetParent(rd);
464             if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
465           }
466         }
467         // D* candidates
468         if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
469           // write references in io2Prong
470           if(fInputAOD) {
471             AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
472           } else {
473             vertexp1n1->AddDaughter(postrack1);
474             vertexp1n1->AddDaughter(negtrack1);
475           }
476           io2Prong->SetSecondaryVtx(vertexp1n1);
477           //printf("--->  %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
478           // create a track from the D0
479           AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
480
481           // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
482           for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
483
484             if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
485
486             if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
487
488             if(fMixEvent) {
489               if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] || 
490                  evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] || 
491                  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
492             }
493
494             if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));  
495
496             trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
497             if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
498
499             // get track from tracks array
500             trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
501             trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
502
503             twoTrackArrayCasc->AddAt(trackPi,0);
504             twoTrackArrayCasc->AddAt(trackD0,1);
505
506             AliAODVertex *vertexCasc = 0;
507
508             if(fFindVertexForDstar) {
509               // DCA between the two tracks
510               dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
511               // Vertexing
512               vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
513             } else {
514               // assume Dstar decays at the primary vertex
515               Double_t pos[3],cov[6],chi2perNDF;
516               fV1->GetXYZ(pos);
517               fV1->GetCovMatrix(cov);
518               chi2perNDF = fV1->GetChi2toNDF();
519               vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
520               dcaCasc = 0.;
521             }
522             if(!vertexCasc) { 
523               twoTrackArrayCasc->Clear();
524               trackPi=0; 
525               continue; 
526             }
527
528             ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
529             if(okDstar) {
530               // add the D0 to the AOD (if not already done)
531               if(!okD0) {
532                 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
533                 if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
534                 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
535                 rd->SetSecondaryVtx(v2Prong);
536                 v2Prong->SetParent(rd);
537                 if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
538                 okD0=kTRUE; // this is done to add it only once
539               }
540               // add the vertex and the cascade to the AOD
541               AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc); 
542               if(fInputAOD) {
543                 AddDaughterRefs(vCasc,event,twoTrackArrayCasc); // add the pion
544               } else {
545                 vCasc->AddDaughter(rd); // just to fill ref #0 
546               }
547               vCasc->AddDaughter(rd); // add the D0 (in ref #1)
548               rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
549               rc->SetSecondaryVtx(vCasc);
550               vCasc->SetParent(rc);
551               if(fInputAOD) rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
552             }
553             twoTrackArrayCasc->Clear();
554             trackPi=0; 
555             if(ioCascade) {delete ioCascade; ioCascade=NULL;}
556             delete vertexCasc; vertexCasc=NULL;
557           } // end loop on soft pi tracks
558
559           if(trackD0) {delete trackD0; trackD0=NULL;}
560
561         }
562         if(io2Prong) {delete io2Prong; io2Prong=NULL;}
563       }      
564
565       twoTrackArray1->Clear(); 
566       if( (!f3Prong && !f4Prong) || 
567           (isLikeSign2Prong && !f3Prong) ) { 
568         negtrack1=0; 
569         delete vertexp1n1; 
570         continue; 
571       }
572
573         
574       // 2nd LOOP  ON  POSITIVE  TRACKS 
575       for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
576
577         if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
578
579         if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));  
580
581         // get track from tracks array
582         postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
583
584         if(postrack2->Charge()<0) continue; 
585
586         if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
587
588         if(fMixEvent) {
589           if(evtNumber[iTrkP1]==evtNumber[iTrkP2] || 
590              evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
591              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
592         }
593
594         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet 
595           if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
596             isLikeSign3Prong=kTRUE;
597           } else { // not ok
598             continue;
599           }
600         } else { // normal triplet (+-+)
601           isLikeSign3Prong=kFALSE; 
602           if(fMixEvent) {
603             if(evtNumber[iTrkP1]==evtNumber[iTrkP2] || 
604                evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
605                evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
606           }
607         }
608
609         // back to primary vertex
610         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
611         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
612         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
613         //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
614
615         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
616         if(dcap2n1>dcaMax) { postrack2=0; continue; }
617         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
618         if(dcap1p2>dcaMax) { postrack2=0; continue; }
619         
620         // Vertexing
621         twoTrackArray2->AddAt(postrack2,0);
622         twoTrackArray2->AddAt(negtrack1,1);
623         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
624         if(!vertexp2n1) { 
625           twoTrackArray2->Clear();
626           postrack2=0; 
627           continue;
628         }
629
630         // 3 prong candidates
631         if(f3Prong) { 
632           if(postrack2->Charge()>0) {
633             threeTrackArray->AddAt(postrack1,0);
634             threeTrackArray->AddAt(negtrack1,1);
635             threeTrackArray->AddAt(postrack2,2);
636           } else {
637             threeTrackArray->AddAt(negtrack1,0);
638             threeTrackArray->AddAt(postrack1,1);
639             threeTrackArray->AddAt(postrack2,2);
640           }
641
642           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
643           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
644           if(ok3Prong) {
645             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
646             if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
647             if(!isLikeSign3Prong) {
648               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
649               rd->SetSecondaryVtx(v3Prong);
650               v3Prong->SetParent(rd);
651             } else { // isLikeSign3Prong
652               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
653               rd->SetSecondaryVtx(v3Prong);
654               v3Prong->SetParent(rd);
655               if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
656             }
657           }
658           if(io3Prong) {delete io3Prong; io3Prong=NULL;} 
659           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;} 
660         }
661
662         // 4 prong candidates
663         if(f4Prong 
664            // don't make 4 prong with like-sign pairs and triplets
665            && !isLikeSign2Prong && !isLikeSign3Prong
666            // track-to-track dca cuts already now
667            && dcap1n1 < fD0to4ProngsCuts[1]
668            && dcap2n1 < fD0to4ProngsCuts[1]) {
669
670           // back to primary vertex
671           postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
672           postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
673           negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
674           // Vertexing for these 3 (can be taken from above?)
675           threeTrackArray->AddAt(postrack1,0);
676           threeTrackArray->AddAt(negtrack1,1);
677           threeTrackArray->AddAt(postrack2,2);
678           AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
679
680           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong) 
681           for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
682
683             if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
684
685             if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
686
687             // get track from tracks array
688             negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
689
690             if(negtrack2->Charge()>0) continue;
691
692             if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
693             if(fMixEvent){ 
694               if(evtNumber[iTrkP1]==evtNumber[iTrkN2] || 
695                  evtNumber[iTrkN1]==evtNumber[iTrkN2] || 
696                  evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
697                  evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
698                  evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
699                  evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
700             }
701
702             // back to primary vertex
703             postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
704             postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
705             negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
706             negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
707             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
708             if(dcap1n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
709             dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
710             if(dcap2n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
711
712             // Vertexing
713             fourTrackArray->AddAt(postrack1,0);
714             fourTrackArray->AddAt(negtrack1,1);
715             fourTrackArray->AddAt(postrack2,2);
716             fourTrackArray->AddAt(negtrack2,3);
717
718             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
719             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
720             if(ok4Prong) {
721               AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
722               if(fInputAOD) AddDaughterRefs(v4Prong,event,fourTrackArray);
723               rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
724               rd->SetSecondaryVtx(v4Prong);
725               v4Prong->SetParent(rd);
726               if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
727             }
728
729             if(io4Prong) {delete io4Prong; io4Prong=NULL;} 
730             if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;} 
731             fourTrackArray->Clear();
732             negtrack2 = 0;
733
734           } // end loop on negative tracks
735
736           threeTrackArray->Clear();
737           delete vertexp1n1p2;
738
739         }
740
741         postrack2 = 0;
742         delete vertexp2n1;
743
744       } // end 2nd loop on positive tracks
745
746       twoTrackArray2->Clear();
747       
748       // 2nd LOOP  ON  NEGATIVE  TRACKS (for 3 prong -+-)
749       for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
750
751         if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
752
753         if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
754
755         // get track from tracks array
756         negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
757
758         if(negtrack2->Charge()>0) continue;
759
760         if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
761
762         if(fMixEvent) {
763           if(evtNumber[iTrkP1]==evtNumber[iTrkN2] || 
764              evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
765              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
766         }
767
768         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet 
769           if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
770             isLikeSign3Prong=kTRUE;
771           } else { // not ok
772             continue;
773           }
774         } else { // normal triplet (-+-)
775           isLikeSign3Prong=kFALSE;
776         }
777
778         // back to primary vertex
779         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
780         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
781         negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
782         //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
783
784         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
785         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
786         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
787         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
788         
789         // Vertexing
790         twoTrackArray2->AddAt(postrack1,0);
791         twoTrackArray2->AddAt(negtrack2,1);
792
793         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
794         if(!vertexp1n2) { 
795           twoTrackArray2->Clear();
796           negtrack2=0; 
797           continue; 
798         }
799
800         if(f3Prong) { 
801           threeTrackArray->AddAt(negtrack1,0);
802           threeTrackArray->AddAt(postrack1,1);
803           threeTrackArray->AddAt(negtrack2,2);
804           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
805           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
806           if(ok3Prong) {
807             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
808             if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
809             if(!isLikeSign3Prong) {
810               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
811               rd->SetSecondaryVtx(v3Prong);
812               v3Prong->SetParent(rd);
813             } else { // isLikeSign3Prong
814               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
815               rd->SetSecondaryVtx(v3Prong);
816               v3Prong->SetParent(rd);
817               if(fInputAOD) rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
818             }
819           }
820           if(io3Prong) {delete io3Prong; io3Prong=NULL;} 
821           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
822         }
823
824         negtrack2 = 0;
825         delete vertexp1n2;
826
827       } // end 2nd loop on negative tracks
828       
829       twoTrackArray2->Clear();
830       
831       negtrack1 = 0;
832       delete vertexp1n1; 
833     } // end 1st loop on negative tracks
834     
835     postrack1 = 0;
836   }  // end 1st loop on positive tracks
837
838
839   AliDebug(1,Form(" Total HF vertices in event = %d;",
840                   (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
841   if(fD0toKpi) {
842     AliDebug(1,Form(" D0->Kpi in event = %d;",
843                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
844   }
845   if(fJPSItoEle) {
846     AliDebug(1,Form(" JPSI->ee in event = %d;",
847                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
848   }
849   if(f3Prong) {
850     AliDebug(1,Form(" Charm->3Prong in event = %d;",
851                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
852   }
853   if(f4Prong) {
854     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
855                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
856   }
857   if(fDstar) {
858     AliDebug(1,Form(" D*->D0pi in event = %d;\n",
859                     (Int_t)aodDstarTClArr->GetEntriesFast()));
860   }
861   if(fLikeSign) {
862     AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
863                     (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
864   }
865   if(fLikeSign && f3Prong) {
866     AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
867                     (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
868   }
869     
870
871   twoTrackArray1->Delete();  delete twoTrackArray1;
872   twoTrackArray2->Delete();  delete twoTrackArray2;
873   twoTrackArrayCasc->Delete();  delete twoTrackArrayCasc;
874   threeTrackArray->Clear(); 
875   threeTrackArray->Delete(); delete threeTrackArray;
876   fourTrackArray->Delete();  delete fourTrackArray;
877   delete [] seleFlags; seleFlags=NULL;
878   if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
879
880   if(fInputAOD) {
881     seleTrksArray.Delete(); 
882     if(fAODMap) { delete fAODMap; fAODMap=NULL; }
883   }
884
885   return;
886 }
887 //----------------------------------------------------------------------------
888 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
889                                              const AliVEvent *event,
890                                              const TObjArray *trkArray) const
891 {
892   // Add the AOD tracks as daughters of the vertex (TRef)
893
894   Int_t nTrks = trkArray->GetEntriesFast();
895
896   AliExternalTrackParam *track = 0;
897   AliAODTrack *aodTrack = 0;
898   Int_t id;
899
900   for(Int_t i=0; i<nTrks; i++) {
901     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
902     id = (Int_t)track->GetID();
903     //printf("---> %d\n",id);
904     if(id<0) continue; // this track is a AliAODRecoDecay
905     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
906     v->AddDaughter(aodTrack);
907   }
908
909   return;
910 }       
911 //----------------------------------------------------------------------------
912 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
913                                    TObjArray *twoTrackArray,AliVEvent *event,
914                                    AliAODVertex *secVert,
915                                    AliAODRecoDecayHF2Prong *rd2Prong,
916                                    Double_t dca,
917                                    Bool_t &okDstar) const
918 {
919   // Make the cascade as a 2Prong decay and check if it passes Dstar
920   // reconstruction cuts
921
922   okDstar = kFALSE;
923
924   Bool_t dummy1,dummy2,dummy3;
925
926   // We use Make2Prong to construct the AliAODRecoCascadeHF
927   // (which inherits from AliAODRecoDecayHF2Prong) 
928   AliAODRecoCascadeHF *theCascade = 
929     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
930                                      dummy1,dummy2,dummy3);
931   if(!theCascade) return 0x0;
932
933   // charge
934   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
935   theCascade->SetCharge(trackPi->Charge());
936
937   //--- selection cuts
938   //
939   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
940   tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
941   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
942   AliAODVertex *primVertexAOD=0;
943   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
944     // take event primary vertex
945     primVertexAOD = PrimaryVertex(); 
946     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
947     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
948   }
949   // select D*->D0pi
950   if(fDstar) {
951     Bool_t testD0=kTRUE;
952     okDstar = tmpCascade->SelectDstar(fDstarCuts,fD0fromDstarCuts,testD0);
953   }
954   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
955   tmpCascade->UnsetOwnPrimaryVtx(); 
956   delete tmpCascade; tmpCascade=NULL;
957   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
958     rd2Prong->UnsetOwnPrimaryVtx();
959   }
960   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
961   //---
962   
963   return theCascade;
964 }
965 //-----------------------------------------------------------------------------
966 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
967                                    TObjArray *twoTrackArray,AliVEvent *event,
968                                    AliAODVertex *secVert,Double_t dca,
969                                    Bool_t &okD0,Bool_t &okJPSI,
970                                    Bool_t &okD0fromDstar) const
971 {
972   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
973   // reconstruction cuts
974   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
975
976   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
977
978   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
979   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
980   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
981
982   // propagate tracks to secondary vertex, to compute inv. mass
983   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
984   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
985
986   Double_t momentum[3];
987   postrack->GetPxPyPz(momentum);
988   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
989   negtrack->GetPxPyPz(momentum);
990   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
991
992
993   // invariant mass cut (try to improve coding here..)
994   Bool_t okMassCut=kFALSE;
995   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
996   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
997   if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
998   if(!okMassCut) {
999     AliDebug(2," candidate didn't pass mass cut");
1000     return 0x0;    
1001   }
1002   // primary vertex to be used by this candidate
1003   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
1004   if(!primVertexAOD) return 0x0;
1005
1006
1007   Double_t d0z0[2],covd0z0[3];
1008   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1009   d0[0] = d0z0[0];
1010   d0err[0] = TMath::Sqrt(covd0z0[0]);
1011   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1012   d0[1] = d0z0[0];
1013   d0err[1] = TMath::Sqrt(covd0z0[0]);
1014   
1015   // create the object AliAODRecoDecayHF2Prong
1016   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1017   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1018   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1019   the2Prong->SetProngIDs(2,id);
1020   delete primVertexAOD; primVertexAOD=NULL;
1021
1022  
1023   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar 
1024     // select D0->Kpi
1025     Int_t checkD0,checkD0bar;
1026     if(fD0toKpi)   okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
1027     //if(fD0toKpi)   okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1028     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
1029     // select J/psi from B
1030     Int_t checkJPSI;
1031     if(fJPSItoEle) okJPSI        = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
1032     //if(fJPSItpEle)   okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1033     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
1034     // select D0->Kpi from Dstar
1035     if(fDstar)     okD0fromDstar = the2Prong->SelectD0(fD0fromDstarCuts,checkD0,checkD0bar);
1036     //if(fDstar)   okD0fromDstar = (Bool_t)fCutsD0fromDstar->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1037     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
1038   }
1039
1040
1041   // remove the primary vertex (was used only for selection)
1042   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1043     the2Prong->UnsetOwnPrimaryVtx();
1044   }
1045   
1046   // get PID info from ESD
1047   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1048   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1049   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1050   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1051   Double_t esdpid[10];
1052   for(Int_t i=0;i<5;i++) {
1053     esdpid[i]   = esdpid0[i];
1054     esdpid[5+i] = esdpid1[i];
1055   }
1056   the2Prong->SetPID(2,esdpid);
1057
1058   return the2Prong;  
1059 }
1060 //----------------------------------------------------------------------------
1061 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1062                              TObjArray *threeTrackArray,AliVEvent *event,
1063                              AliAODVertex *secVert,Double_t dispersion,
1064                              const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1065                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1066                              Bool_t &ok3Prong) const
1067 {
1068   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1069   // reconstruction cuts 
1070   // E.Bruna, F.Prino
1071
1072
1073   ok3Prong=kFALSE;
1074   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0; 
1075
1076   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1077   Double_t momentum[3];
1078
1079
1080   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1081   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1082   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1083
1084   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1085   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1086   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1087   postrack1->GetPxPyPz(momentum);
1088   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
1089   negtrack->GetPxPyPz(momentum);
1090   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
1091   postrack2->GetPxPyPz(momentum);
1092   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2]; 
1093
1094   // invariant mass cut for D+, Ds, Lc
1095   Bool_t okMassCut=kFALSE;
1096   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1097   if(!okMassCut) {
1098     AliDebug(2," candidate didn't pass mass cut");
1099     return 0x0;    
1100   }
1101
1102   // primary vertex to be used by this candidate
1103   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
1104   if(!primVertexAOD) return 0x0;
1105
1106   Double_t d0z0[2],covd0z0[3];
1107   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1108   d0[0]=d0z0[0];
1109   d0err[0] = TMath::Sqrt(covd0z0[0]);
1110   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1111   d0[1]=d0z0[0];
1112   d0err[1] = TMath::Sqrt(covd0z0[0]);
1113   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1114   d0[2]=d0z0[0];
1115   d0err[2] = TMath::Sqrt(covd0z0[0]);
1116
1117
1118   // create the object AliAODRecoDecayHF3Prong
1119   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1120   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1121   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]));
1122   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]));
1123   Short_t charge=(Short_t)(postrack1->Charge()*postrack2->Charge()*negtrack->Charge());
1124   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1125   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1126   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1127   the3Prong->SetProngIDs(3,id);
1128
1129   delete primVertexAOD; primVertexAOD=NULL;
1130
1131   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1132   if(f3Prong) {
1133     ok3Prong = kFALSE;
1134     Int_t ok1,ok2;
1135     Int_t dum1,dum2;
1136     if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
1137     if(the3Prong->SelectDs(fDsCuts,ok1,ok2,dum1,dum2)) ok3Prong = kTRUE;
1138     if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
1139     /*
1140     if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1141     if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1142     if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1143     */
1144   }
1145   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1146
1147   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1148     the3Prong->UnsetOwnPrimaryVtx();
1149   }
1150
1151   // get PID info from ESD
1152   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1153   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1154   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1155   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1156   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1157   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1158   
1159   Double_t esdpid[15];
1160   for(Int_t i=0;i<5;i++) {
1161     esdpid[i]    = esdpid0[i];
1162     esdpid[5+i]  = esdpid1[i];
1163     esdpid[10+i] = esdpid2[i];
1164   }
1165   the3Prong->SetPID(3,esdpid);
1166
1167   return the3Prong;
1168 }
1169 //----------------------------------------------------------------------------
1170 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1171                              TObjArray *fourTrackArray,AliVEvent *event,
1172                              AliAODVertex *secVert,
1173                              const AliAODVertex *vertexp1n1,
1174                              const AliAODVertex *vertexp1n1p2,
1175                              Double_t dcap1n1,Double_t dcap1n2,
1176                              Double_t dcap2n1,Double_t dcap2n2,
1177                              Bool_t &ok4Prong) const
1178 {
1179   // Make 4Prong candidates and check if they pass D0toKpipipi
1180   // reconstruction cuts
1181   // G.E.Bruno, R.Romita
1182
1183   ok4Prong=kFALSE;
1184   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0; 
1185
1186   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1187
1188   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1189   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1190   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1191   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1192
1193   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1194   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1195   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1196   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1197
1198   Double_t momentum[3];
1199   postrack1->GetPxPyPz(momentum);
1200   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1201   negtrack1->GetPxPyPz(momentum);
1202   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1203   postrack2->GetPxPyPz(momentum);
1204   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1205   negtrack2->GetPxPyPz(momentum);
1206   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1207
1208   // invariant mass cut for rho or D0 (try to improve coding here..)
1209   Bool_t okMassCut=kFALSE;
1210   if(!okMassCut && TMath::Abs(fD0to4ProngsCuts[8])<1.e-13){      //no PID, to be implemented with PID
1211     if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1212   }
1213   if(!okMassCut) {
1214     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1215     //printf(" candidate didn't pass mass cut\n");
1216     return 0x0;
1217   }
1218
1219   // primary vertex to be used by this candidate
1220   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1221   if(!primVertexAOD) return 0x0;
1222
1223   Double_t d0z0[2],covd0z0[3];
1224   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1225   d0[0]=d0z0[0];
1226   d0err[0] = TMath::Sqrt(covd0z0[0]);
1227   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1228   d0[1]=d0z0[0];
1229   d0err[1] = TMath::Sqrt(covd0z0[0]);
1230   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1231   d0[2]=d0z0[0];
1232   d0err[2] = TMath::Sqrt(covd0z0[0]);
1233   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1234   d0[3]=d0z0[0];
1235   d0err[3] = TMath::Sqrt(covd0z0[0]);
1236
1237
1238   // create the object AliAODRecoDecayHF4Prong
1239   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1240   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1241   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]));
1242   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]));
1243   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]));
1244   Short_t charge=0;
1245   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1246   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1247   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1248   the4Prong->SetProngIDs(4,id);
1249
1250   delete primVertexAOD; primVertexAOD=NULL;
1251
1252   Int_t checkD0,checkD0bar;
1253   ok4Prong=the4Prong->SelectD0(fD0to4ProngsCuts,checkD0,checkD0bar);
1254   //ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1255
1256
1257   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1258     the4Prong->UnsetOwnPrimaryVtx();
1259   }
1260
1261  
1262   // get PID info from ESD
1263   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1264   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1265   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1266   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1267   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1268   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1269   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1270   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1271
1272   Double_t esdpid[20];
1273   for(Int_t i=0;i<5;i++) {
1274     esdpid[i]    = esdpid0[i];
1275     esdpid[5+i]  = esdpid1[i];
1276     esdpid[10+i] = esdpid2[i];
1277     esdpid[15+i] = esdpid3[i];
1278   }
1279   the4Prong->SetPID(4,esdpid);
1280   
1281   return the4Prong;
1282 }
1283 //-----------------------------------------------------------------------------
1284 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1285                                                     AliVEvent *event) const
1286 {
1287   // Returns primary vertex to be used for this candidate
1288  
1289   AliESDVertex *vertexESD = 0;
1290   AliAODVertex *vertexAOD = 0;
1291
1292
1293   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) { 
1294     // primary vertex from the input event
1295     
1296     vertexESD = new AliESDVertex(*fV1);
1297
1298   } else {
1299     // primary vertex specific to this candidate
1300
1301     Int_t nTrks = trkArray->GetEntriesFast();
1302     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1303
1304     if(fRecoPrimVtxSkippingTrks) { 
1305       // recalculating the vertex
1306       
1307       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1308         Float_t diamondcovxy[3];
1309         event->GetDiamondCovXY(diamondcovxy);
1310         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1311         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1312         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1313         vertexer->SetVtxStart(diamond);
1314         delete diamond; diamond=NULL;
1315         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) 
1316           vertexer->SetOnlyFitter();
1317       }
1318       Int_t skipped[1000];
1319       Int_t nTrksToSkip=0,id;
1320       AliExternalTrackParam *t = 0;
1321       for(Int_t i=0; i<nTrks; i++) {
1322         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1323         id = (Int_t)t->GetID();
1324         if(id<0) continue;
1325         skipped[nTrksToSkip++] = id;
1326       }
1327       // TEMPORARY FIX
1328       // For AOD, skip also tracks without covariance matrix
1329       if(fInputAOD) {
1330         Double_t covtest[21];
1331         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1332           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1333           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1334             id = (Int_t)t->GetID();
1335             if(id<0) continue;
1336             skipped[nTrksToSkip++] = id;
1337           }
1338         }
1339       }
1340       //
1341       vertexer->SetSkipTracks(nTrksToSkip,skipped);
1342       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
1343       
1344     } else if(fRmTrksFromPrimVtx) { 
1345       // removing the prongs tracks
1346       
1347       TObjArray rmArray(nTrks);
1348       UShort_t *rmId = new UShort_t[nTrks];
1349       AliESDtrack *esdTrack = 0;
1350       AliESDtrack *t = 0;
1351       for(Int_t i=0; i<nTrks; i++) {
1352         t = (AliESDtrack*)trkArray->UncheckedAt(i);
1353         esdTrack = new AliESDtrack(*t);
1354         rmArray.AddLast(esdTrack);
1355         if(esdTrack->GetID()>=0) {
1356           rmId[i]=(UShort_t)esdTrack->GetID();
1357         } else {
1358           rmId[i]=9999;
1359         }
1360       }
1361       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1362       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1363       delete [] rmId; rmId=NULL;
1364       rmArray.Delete();
1365       
1366     }
1367
1368     if(!vertexESD) return vertexAOD;
1369     if(vertexESD->GetNContributors()<=0) { 
1370       AliDebug(2,"vertexing failed"); 
1371       delete vertexESD; vertexESD=NULL;
1372       return vertexAOD;
1373     }
1374
1375     delete vertexer; vertexer=NULL;
1376
1377   }
1378
1379   // convert to AliAODVertex
1380   Double_t pos[3],cov[6],chi2perNDF;
1381   vertexESD->GetXYZ(pos); // position
1382   vertexESD->GetCovMatrix(cov); //covariance matrix
1383   chi2perNDF = vertexESD->GetChi2toNDF();
1384   delete vertexESD; vertexESD=NULL;
1385
1386   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1387
1388   return vertexAOD;
1389 }
1390 //-----------------------------------------------------------------------------
1391 void AliAnalysisVertexingHF::PrintStatus() const {
1392   // Print parameters being used
1393
1394   //printf("Preselections:\n");
1395   //   fTrackFilter->Dump();
1396   if(fSecVtxWithKF) {
1397     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1398   } else {
1399     printf("Secondary vertex with AliVertexerTracks\n");
1400   }
1401   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1402   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1403   if(fD0toKpi) {
1404     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1405     if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1406     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
1407     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
1408     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
1409     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
1410     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
1411     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
1412     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
1413     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
1414     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
1415   }
1416   if(fDstar) {
1417     if(fFindVertexForDstar) {
1418       printf("    Reconstruct a secondary vertex for the D*\n");
1419     } else {
1420       printf("    Assume the D* comes from the primary vertex\n");
1421     }
1422     printf("    |M-MD*| [GeV]    < %f\n",fDstarCuts[0]);
1423     printf("    |M_Kpipi-M_Kpi-(MD*-MD0)| [GeV]  < %f\n",fDstarCuts[1]);
1424     printf("    pTpisoft [GeV/c]    > %f\n",fDstarCuts[2]);
1425     printf("    pTpisoft [GeV/c]    < %f\n",fDstarCuts[3]);
1426     printf("    Theta(pisoft,D0plane) < %f\n",fDstarCuts[4]);
1427     printf("Reconstruct D*->D0pi candidates with cuts:\n");
1428     printf("   D0 from D* cuts:\n");
1429     if(fCutsD0fromDstar) fCutsD0fromDstar->PrintAll();
1430     printf("    |M-MD0| [GeV]    < %f\n",fD0fromDstarCuts[0]);
1431     printf("    dca    [cm]  < %f\n",fD0fromDstarCuts[1]);
1432     printf("    cosThetaStar     < %f\n",fD0fromDstarCuts[2]);
1433     printf("    pTK     [GeV/c]    > %f\n",fD0fromDstarCuts[3]);
1434     printf("    pTpi    [GeV/c]    > %f\n",fD0fromDstarCuts[4]);
1435     printf("    |d0K|  [cm]  < %f\n",fD0fromDstarCuts[5]);
1436     printf("    |d0pi| [cm]  < %f\n",fD0fromDstarCuts[6]);
1437     printf("    d0d0  [cm^2] < %f\n",fD0fromDstarCuts[7]);
1438     printf("    cosThetaPoint    > %f\n",fD0fromDstarCuts[8]);
1439   }
1440   if(fJPSItoEle) {
1441     printf("Reconstruct J/psi from B candidates with cuts:\n");
1442     if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1443     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
1444     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
1445     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
1446     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
1447     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
1448     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
1449     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
1450     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
1451     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
1452   }
1453   if(f3Prong) {
1454     printf("Reconstruct 3 prong candidates.\n");
1455     printf("  D+->Kpipi cuts:\n");
1456     if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1457     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
1458     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
1459     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
1460     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
1461     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
1462     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
1463     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
1464     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1465     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1466     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
1467     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
1468     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
1469     printf("  Ds->KKpi cuts:\n");
1470     if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1471     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
1472     printf("    pTK     [GeV/c]    > %f\n",fDsCuts[1]);
1473     printf("    pTPi    [GeV/c]    > %f\n",fDsCuts[2]);
1474     printf("    |d0K|  [cm]  > %f\n",fDsCuts[3]);
1475     printf("    |d0Pi| [cm]  > %f\n",fDsCuts[4]);
1476     printf("    dist12    [cm]  < %f\n",fDsCuts[5]);
1477     printf("    sigmavert [cm]   < %f\n",fDsCuts[6]);
1478     printf("    dist prim-sec [cm] > %f\n",fDsCuts[7]);
1479     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
1480     printf("    cosThetaPoint    > %f\n",fDsCuts[9]);
1481     printf("    Sum d0^2 [cm^2]  > %f\n",fDsCuts[10]);
1482     printf("    dca cut [cm]  < %f\n",fDsCuts[11]);
1483     printf("    Inv. Mass  phi [GeV]  < %f\n",fDsCuts[12]);
1484     printf("    Inv. Mass  K0* [GeV]  < %f\n",fDsCuts[13]);
1485     printf("  Lc->pKpi cuts:\n");
1486     if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1487     printf("    |M-MLc| [GeV]    < %f\n",fLcCuts[0]);
1488     printf("    pTP     [GeV/c]    > %f\n",fLcCuts[1]);
1489     printf("    pTPi and pTK [GeV/c]    > %f\n",fLcCuts[2]);
1490     printf("    |d0P|  [cm]  > %f\n",fLcCuts[3]);
1491     printf("    |d0Pi| and |d0K| [cm]  > %f\n",fLcCuts[4]);
1492     printf("    dist12    [cm]  < %f\n",fLcCuts[5]);
1493     printf("    sigmavert [cm]   < %f\n",fLcCuts[6]);
1494     printf("    dist prim-sec [cm] > %f\n",fLcCuts[7]);
1495     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
1496     printf("    cosThetaPoint    > %f\n",fLcCuts[9]);
1497     printf("    Sum d0^2 [cm^2]  > %f\n",fLcCuts[10]);
1498     printf("    dca cut [cm]  < %f\n",fLcCuts[11]);
1499   }
1500   if(f4Prong) {
1501     printf("Reconstruct 4 prong candidates.\n");
1502     printf("  D0->Kpipipi cuts:\n");
1503     if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1504   }
1505
1506   return;
1507 }
1508 //-----------------------------------------------------------------------------
1509 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1510                                                                  Double_t &dispersion,Bool_t useTRefArray) const
1511 {
1512   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1513
1514   AliESDVertex *vertexESD = 0;
1515   AliAODVertex *vertexAOD = 0;
1516
1517   if(!fSecVtxWithKF) { // AliVertexerTracks
1518
1519     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1520     vertexer->SetVtxStart(fV1);
1521     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1522     delete vertexer; vertexer=NULL;
1523
1524     if(!vertexESD) return vertexAOD;
1525
1526     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
1527       AliDebug(2,"vertexing failed"); 
1528       delete vertexESD; vertexESD=NULL;
1529       return vertexAOD;
1530     }
1531
1532   } else { // Kalman Filter vertexer (AliKFParticle)
1533
1534     AliKFParticle::SetField(fBzkG);
1535
1536     AliKFVertex vertexKF;
1537
1538     Int_t nTrks = trkArray->GetEntriesFast();
1539     for(Int_t i=0; i<nTrks; i++) {
1540       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1541       AliKFParticle daughterKF(*esdTrack,211);
1542       vertexKF.AddDaughter(daughterKF);
1543     }
1544     vertexESD = new AliESDVertex(vertexKF.Parameters(),
1545                                  vertexKF.CovarianceMatrix(),
1546                                  vertexKF.GetChi2(),
1547                                  vertexKF.GetNContributors());
1548
1549   }
1550
1551   // convert to AliAODVertex
1552   Double_t pos[3],cov[6],chi2perNDF;
1553   vertexESD->GetXYZ(pos); // position
1554   vertexESD->GetCovMatrix(cov); //covariance matrix
1555   chi2perNDF = vertexESD->GetChi2toNDF();
1556   dispersion = vertexESD->GetDispersion();
1557   delete vertexESD; vertexESD=NULL;
1558
1559   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1560   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1561
1562   return vertexAOD;
1563 }
1564 //-----------------------------------------------------------------------------
1565 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1566                                              Int_t nprongs,
1567                                              Double_t *px,
1568                                              Double_t *py,
1569                                              Double_t *pz) const {
1570   // Check invariant mass cut
1571
1572   Short_t dummycharge=0;
1573   Double_t *dummyd0 = new Double_t[nprongs];
1574   for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1575   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1576   delete [] dummyd0; dummyd0=NULL;
1577
1578   UInt_t pdg2[2],pdg3[3],pdg4[4];
1579   Double_t mPDG,minv;
1580
1581   Bool_t retval=kFALSE;
1582   switch (decay) 
1583     { 
1584     case 0:                  // D0->Kpi
1585       pdg2[0]=211; pdg2[1]=321;
1586       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1587       minv = rd->InvMass(nprongs,pdg2);
1588       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1589       //if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1590       pdg2[0]=321; pdg2[1]=211;
1591       minv = rd->InvMass(nprongs,pdg2);
1592       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1593       //if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1594       break;
1595     case 1:                  // JPSI->ee
1596       pdg2[0]=11; pdg2[1]=11;
1597       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1598       minv = rd->InvMass(nprongs,pdg2);
1599       if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1600       //if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
1601       break;
1602     case 2:                  // D+->Kpipi
1603       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1604       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1605       minv = rd->InvMass(nprongs,pdg3);
1606       if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1607       //if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
1608                             // Ds+->KKpi
1609       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1610       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1611       minv = rd->InvMass(nprongs,pdg3);
1612       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1613       //if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1614       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1615       minv = rd->InvMass(nprongs,pdg3);
1616       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1617       //if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1618                             // Lc->pKpi
1619       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1620       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1621       minv = rd->InvMass(nprongs,pdg3);
1622       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1623       //if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1624       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1625       minv = rd->InvMass(nprongs,pdg3);
1626       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE; 
1627       //if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1628       break;
1629     case 3:                  // D*->D0pi
1630       pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
1631       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1632       minv = rd->InvMass(nprongs,pdg2);
1633       if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
1634       break;
1635     case 4:                 // D0->Kpipipi without PID
1636       pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1637       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1638       minv = rd->InvMass(nprongs,pdg4);
1639       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1640       //if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1641       pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1642       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1643       minv = rd->InvMass(nprongs,pdg4);
1644       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1645       //if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1646       pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1647       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1648       minv = rd->InvMass(nprongs,pdg4);
1649       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1650       //if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1651       pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1652       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1653       minv = rd->InvMass(nprongs,pdg4);
1654       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1655       //if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1656       break;
1657     default:
1658       printf("SelectInvMass(): wrong decay selection\n");
1659       break;
1660     }
1661
1662   delete rd; rd=NULL;
1663
1664   return retval;
1665 }
1666 //-----------------------------------------------------------------------------
1667 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
1668                                    TObjArray &seleTrksArray,Int_t &nSeleTrks,
1669                                    UChar_t *seleFlags,Int_t *evtNumber)
1670 {
1671   // Apply single-track preselection.
1672   // Fill a TObjArray with selected tracks (for displaced vertices or
1673   // soft pion from D*). Selection flag stored in seleFlags.
1674   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1675   // In case of AOD input, also fill fAODMap for track index<->ID
1676
1677   const AliVVertex *vprimary = event->GetPrimaryVertex();
1678
1679   if(fV1) { delete fV1; fV1=NULL; }
1680   if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1681
1682   Int_t nindices=0;
1683   UShort_t *indices = 0;
1684   Double_t pos[3],cov[6];
1685
1686   if(!fInputAOD) { // ESD
1687     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
1688   } else {         // AOD
1689     vprimary->GetXYZ(pos);
1690     vprimary->GetCovarianceMatrix(cov);
1691     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1692     indices = new UShort_t[event->GetNumberOfTracks()];
1693     fAODMapSize = 100000;
1694     fAODMap = new Int_t[fAODMapSize];
1695   }
1696
1697
1698   Int_t entries = (Int_t)event->GetNumberOfTracks();
1699   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
1700   nSeleTrks=0;
1701  
1702   // transfer ITS tracks from event to arrays
1703   for(Int_t i=0; i<entries; i++) {
1704     AliVTrack *track;
1705     track = (AliVTrack*)event->GetTrack(i);
1706
1707     // skip pure ITS SA tracks
1708     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1709
1710     // TEMPORARY: check that the cov matrix is there
1711     Double_t covtest[21];
1712     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1713     //
1714
1715     if(fInputAOD) {
1716       AliAODTrack *aodt = (AliAODTrack*)track;
1717       if(aodt->GetUsedForPrimVtxFit()) { 
1718         indices[nindices]=aodt->GetID(); nindices++; 
1719       }
1720       fAODMap[(Int_t)aodt->GetID()] = i;
1721     }
1722
1723     AliESDtrack *esdt = 0;
1724
1725     if(!fInputAOD) {
1726       esdt = (AliESDtrack*)track;
1727     } else {
1728       esdt = new AliESDtrack(track);
1729     }
1730
1731     // single track selection
1732     okDisplaced=kFALSE; okSoftPi=kFALSE;
1733     if(fMixEvent){
1734       evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
1735       const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
1736       Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
1737       eventVtx->GetXYZ(vtxPos);
1738       vprimary->GetXYZ(primPos);
1739       eventVtx->GetCovarianceMatrix(primCov);
1740       for(Int_t ind=0;ind<3;ind++){
1741         trasl[ind]=vtxPos[ind]-primPos[ind];
1742       }
1743       
1744       Bool_t isTransl=esdt->Translate(trasl,primCov);
1745       if(!isTransl) {
1746         delete esdt;
1747         esdt = NULL;
1748         continue;
1749       }
1750     }
1751
1752     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
1753       seleTrksArray.AddLast(esdt);
1754       seleFlags[nSeleTrks]=0;
1755       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
1756       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
1757       nSeleTrks++;
1758     } else {
1759       if(fInputAOD) delete esdt; 
1760       esdt = NULL;
1761       continue;
1762     } 
1763
1764   } // end loop on tracks
1765
1766   // primary vertex from AOD
1767   if(fInputAOD) {
1768     delete fV1; fV1=NULL;
1769     vprimary->GetXYZ(pos);
1770     vprimary->GetCovarianceMatrix(cov);
1771     Double_t chi2toNDF = vprimary->GetChi2perNDF();
1772     Int_t ncontr=nindices;
1773     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
1774     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.); 
1775     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
1776     fV1->SetTitle(vprimary->GetTitle());
1777     fV1->SetIndices(nindices,indices);
1778   }
1779   if(indices) { delete [] indices; indices=NULL; }
1780
1781
1782   return;
1783 }
1784 //-----------------------------------------------------------------------------
1785 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1786                                    Double_t cut2,Double_t cut3,Double_t cut4,
1787                                    Double_t cut5,Double_t cut6,
1788                                    Double_t cut7,Double_t cut8) 
1789 {
1790   // Set the cuts for D0 selection
1791   fD0toKpiCuts[0] = cut0;
1792   fD0toKpiCuts[1] = cut1;
1793   fD0toKpiCuts[2] = cut2;
1794   fD0toKpiCuts[3] = cut3;
1795   fD0toKpiCuts[4] = cut4;
1796   fD0toKpiCuts[5] = cut5;
1797   fD0toKpiCuts[6] = cut6;
1798   fD0toKpiCuts[7] = cut7;
1799   fD0toKpiCuts[8] = cut8;
1800
1801   return;
1802 }
1803 //-----------------------------------------------------------------------------
1804 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9]) 
1805 {
1806   // Set the cuts for D0 selection
1807
1808   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1809
1810   return;
1811 }
1812 //-----------------------------------------------------------------------------
1813 void AliAnalysisVertexingHF::SetD0fromDstarCuts(Double_t cut0,Double_t cut1,
1814                                    Double_t cut2,Double_t cut3,Double_t cut4,
1815                                    Double_t cut5,Double_t cut6,
1816                                    Double_t cut7,Double_t cut8) 
1817 {
1818   // Set the cuts for D0 from D* selection
1819   fD0fromDstarCuts[0] = cut0;
1820   fD0fromDstarCuts[1] = cut1;
1821   fD0fromDstarCuts[2] = cut2;
1822   fD0fromDstarCuts[3] = cut3;
1823   fD0fromDstarCuts[4] = cut4;
1824   fD0fromDstarCuts[5] = cut5;
1825   fD0fromDstarCuts[6] = cut6;
1826   fD0fromDstarCuts[7] = cut7;
1827   fD0fromDstarCuts[8] = cut8;
1828
1829   return;
1830 }
1831 //-----------------------------------------------------------------------------
1832 void AliAnalysisVertexingHF::SetD0fromDstarCuts(const Double_t cuts[9]) 
1833 {
1834   // Set the cuts for D0 from D* selection
1835
1836   for(Int_t i=0; i<9; i++) fD0fromDstarCuts[i] = cuts[i];
1837
1838   return;
1839 }
1840 //-----------------------------------------------------------------------------
1841 void AliAnalysisVertexingHF::SetDstarCuts(Double_t cut0,Double_t cut1,
1842                                           Double_t cut2,Double_t cut3,
1843                                           Double_t cut4)
1844 {
1845   // Set the cuts for D* selection
1846   fDstarCuts[0] = cut0;
1847   fDstarCuts[1] = cut1;
1848   fDstarCuts[2] = cut2;
1849   fDstarCuts[3] = cut3;
1850   fDstarCuts[4] = cut4;
1851
1852   return;
1853 }
1854 //-----------------------------------------------------------------------------
1855 void AliAnalysisVertexingHF::SetDstarCuts(const Double_t cuts[5])
1856 {
1857   // Set the cuts for D* selection
1858
1859   for(Int_t i=0; i<5; i++) fDstarCuts[i] = cuts[i];
1860
1861   return;
1862 }
1863 //-----------------------------------------------------------------------------
1864 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1865                                    Double_t cut2,Double_t cut3,Double_t cut4,
1866                                    Double_t cut5,Double_t cut6,
1867                                    Double_t cut7,Double_t cut8) 
1868 {
1869   // Set the cuts for J/psi from B selection
1870   fBtoJPSICuts[0] = cut0;
1871   fBtoJPSICuts[1] = cut1;
1872   fBtoJPSICuts[2] = cut2;
1873   fBtoJPSICuts[3] = cut3;
1874   fBtoJPSICuts[4] = cut4;
1875   fBtoJPSICuts[5] = cut5;
1876   fBtoJPSICuts[6] = cut6;
1877   fBtoJPSICuts[7] = cut7;
1878   fBtoJPSICuts[8] = cut8;
1879
1880   return;
1881 }
1882 //-----------------------------------------------------------------------------
1883 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9]) 
1884 {
1885   // Set the cuts for J/psi from B selection
1886
1887   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1888
1889   return;
1890 }
1891 //-----------------------------------------------------------------------------
1892 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1893                                    Double_t cut2,Double_t cut3,Double_t cut4,
1894                                    Double_t cut5,Double_t cut6,
1895                                    Double_t cut7,Double_t cut8,
1896                                    Double_t cut9,Double_t cut10,Double_t cut11)
1897 {
1898   // Set the cuts for Dplus->Kpipi selection
1899   fDplusCuts[0] = cut0;
1900   fDplusCuts[1] = cut1;
1901   fDplusCuts[2] = cut2;
1902   fDplusCuts[3] = cut3;
1903   fDplusCuts[4] = cut4;
1904   fDplusCuts[5] = cut5;
1905   fDplusCuts[6] = cut6;
1906   fDplusCuts[7] = cut7;
1907   fDplusCuts[8] = cut8;
1908   fDplusCuts[9] = cut9;
1909   fDplusCuts[10] = cut10;
1910   fDplusCuts[11] = cut11;
1911
1912   return;
1913 }
1914 //-----------------------------------------------------------------------------
1915 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12]) 
1916 {
1917   // Set the cuts for Dplus->Kpipi selection
1918
1919   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1920
1921   return;
1922 }
1923 //-----------------------------------------------------------------------------
1924 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
1925                                        Double_t cut2,Double_t cut3,
1926                                        Double_t cut4,Double_t cut5,
1927                                        Double_t cut6,Double_t cut7,
1928                                        Double_t cut8,Double_t cut9,
1929                                        Double_t cut10,Double_t cut11,
1930                                        Double_t cut12,Double_t cut13)
1931 {
1932   // Set the cuts for Ds->KKpi selection
1933   fDsCuts[0] = cut0;
1934   fDsCuts[1] = cut1;
1935   fDsCuts[2] = cut2;
1936   fDsCuts[3] = cut3;
1937   fDsCuts[4] = cut4;
1938   fDsCuts[5] = cut5;
1939   fDsCuts[6] = cut6;
1940   fDsCuts[7] = cut7;
1941   fDsCuts[8] = cut8;
1942   fDsCuts[9] = cut9;
1943   fDsCuts[10] = cut10;
1944   fDsCuts[11] = cut11;
1945   fDsCuts[12] = cut12;
1946   fDsCuts[13] = cut13;
1947
1948   return;
1949 }
1950 //-----------------------------------------------------------------------------
1951 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13]) 
1952 {
1953   // Set the cuts for Ds->KKpi selection
1954
1955   for(Int_t i=0; i<14; i++) fDsCuts[i] = cuts[i];
1956
1957   return;
1958 }
1959 //-----------------------------------------------------------------------------
1960 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
1961                                    Double_t cut2,Double_t cut3,Double_t cut4,
1962                                    Double_t cut5,Double_t cut6,
1963                                    Double_t cut7,Double_t cut8,
1964                                    Double_t cut9,Double_t cut10,Double_t cut11)
1965 {
1966   // Set the cuts for Lc->pKpi selection
1967   fLcCuts[0] = cut0;
1968   fLcCuts[1] = cut1;
1969   fLcCuts[2] = cut2;
1970   fLcCuts[3] = cut3;
1971   fLcCuts[4] = cut4;
1972   fLcCuts[5] = cut5;
1973   fLcCuts[6] = cut6;
1974   fLcCuts[7] = cut7;
1975   fLcCuts[8] = cut8;
1976   fLcCuts[9] = cut9;
1977   fLcCuts[10] = cut10;
1978   fLcCuts[11] = cut11;
1979
1980   return;
1981 }
1982 //-----------------------------------------------------------------------------
1983 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12]) 
1984 {
1985   // Set the cuts for Lc->pKpi selection
1986
1987   for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
1988
1989   return;
1990 }
1991 //-----------------------------------------------------------------------------
1992 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(Double_t cut0,Double_t cut1,
1993                                    Double_t cut2,Double_t cut3,Double_t cut4,
1994                                    Double_t cut5,Double_t cut6,
1995                                    Double_t cut7,Double_t cut8)
1996 {
1997   // Set the cuts for D0->Kpipipi selection
1998
1999   fD0to4ProngsCuts[0] = cut0;
2000   fD0to4ProngsCuts[1] = cut1;
2001   fD0to4ProngsCuts[2] = cut2;
2002   fD0to4ProngsCuts[3] = cut3;
2003   fD0to4ProngsCuts[4] = cut4;
2004   fD0to4ProngsCuts[5] = cut5;
2005   fD0to4ProngsCuts[6] = cut6;
2006   fD0to4ProngsCuts[7] = cut7;
2007   fD0to4ProngsCuts[8] = cut8;
2008
2009   return;
2010 }
2011 //-----------------------------------------------------------------------------
2012 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(const Double_t cuts[9])
2013 {
2014   // Set the cuts for D0->Kpipipi selection
2015
2016   for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i] = cuts[i];
2017
2018   return;
2019 }
2020 //-----------------------------------------------------------------------------
2021 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2022                                              Bool_t &okDisplaced,Bool_t &okSoftPi) const 
2023 {
2024   // Check if track passes some kinematical cuts  
2025
2026   // this is needed to store the impact parameters
2027   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2028
2029   UInt_t selectInfo;
2030   //
2031   // Track selection, displaced tracks
2032   selectInfo = 0; 
2033   if(fTrackFilter) {
2034     selectInfo = fTrackFilter->IsSelected(trk);
2035   }
2036   if(selectInfo) okDisplaced=kTRUE;
2037   // Track selection, soft pions
2038   selectInfo = 0; 
2039   if(fDstar && fTrackFilterSoftPi) {
2040     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2041   }
2042   if(selectInfo) okSoftPi=kTRUE;
2043
2044   if(okDisplaced || okSoftPi) return kTRUE;
2045
2046   return kFALSE;
2047 }
2048 //-----------------------------------------------------------------------------