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