Ref to the 2Prong moved from AliAODRecoCascadeHF to AliAODVertex
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisVertexingHF.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //----------------------------------------------------------------------------
17 //    Implementation of the heavy-flavour vertexing analysis class
18 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
19 // To be used as a task of AliAnalysisManager by means of the interface
20 // class AliAnalysisTaskSEVertexingHF. 
21 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
22 //
23 //  Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita, X.M.Zhang
24 //  Contact: andrea.dainese@lnl.infn.it
25 //----------------------------------------------------------------------------
26 #include <TFile.h>
27 #include <TDatabasePDG.h>
28 #include <TString.h>
29 #include "AliLog.h"
30 #include "AliVEvent.h"
31 #include "AliVVertex.h"
32 #include "AliVTrack.h"
33 #include "AliVertexerTracks.h"
34 #include "AliKFVertex.h"
35 #include "AliESDEvent.h"
36 #include "AliESDVertex.h"
37 #include "AliExternalTrackParam.h"
38 #include "AliNeutralTrackParam.h"
39 #include "AliESDtrack.h"
40 #include "AliESDtrackCuts.h"
41 #include "AliAODEvent.h"
42 #include "AliAODRecoDecay.h"
43 #include "AliAODRecoDecayHF.h"
44 #include "AliAODRecoDecayHF2Prong.h"
45 #include "AliAODRecoDecayHF3Prong.h"
46 #include "AliAODRecoDecayHF4Prong.h"
47 #include "AliAODRecoCascadeHF.h"
48 #include "AliAnalysisFilter.h"
49 #include "AliAnalysisVertexingHF.h"
50
51 ClassImp(AliAnalysisVertexingHF)
52
53 //----------------------------------------------------------------------------
54 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
55 fInputAOD(kFALSE),
56 fAODMap(0),
57 fBzkG(0.),
58 fSecVtxWithKF(kFALSE),
59 fRecoPrimVtxSkippingTrks(kFALSE),
60 fRmTrksFromPrimVtx(kFALSE),
61 fV1(0x0),
62 fD0toKpi(kTRUE),
63 fJPSItoEle(kTRUE),
64 f3Prong(kTRUE),
65 f4Prong(kTRUE),
66 fDstar(kTRUE),
67 fTrackFilter(0x0),
68 fTrackFilterSoftPi(0x0)
69 {
70   // Default constructor
71
72   SetD0toKpiCuts();
73   SetBtoJPSICuts();
74   SetDplusCuts();
75   SetDsCuts();
76   SetLcCuts();
77   SetDstarCuts();
78
79 }
80 //--------------------------------------------------------------------------
81 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) : 
82 TNamed(source),
83 fInputAOD(source.fInputAOD),
84 fAODMap(source.fAODMap),
85 fBzkG(source.fBzkG),
86 fSecVtxWithKF(source.fSecVtxWithKF),
87 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
88 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
89 fV1(source.fV1),
90 fD0toKpi(source.fD0toKpi),
91 fJPSItoEle(source.fJPSItoEle),
92 f3Prong(source.f3Prong),
93 f4Prong(source.f4Prong),
94 fDstar(source.fDstar),
95 fTrackFilter(source.fTrackFilter),
96 fTrackFilterSoftPi(source.fTrackFilterSoftPi)
97 {
98   //
99   // Copy constructor
100   //
101   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
102   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
103   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
104   for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
105   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
106   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
107 }
108 //--------------------------------------------------------------------------
109 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
110 {
111   //
112   // assignment operator
113   //
114   if(&source == this) return *this;
115   fInputAOD = source.fInputAOD;
116   fBzkG = source.fBzkG;
117   fSecVtxWithKF = source.fSecVtxWithKF;
118   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
119   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
120   fV1 = source.fV1;
121   fD0toKpi = source.fD0toKpi;
122   fJPSItoEle = source.fJPSItoEle;
123   f3Prong = source.f3Prong;
124   f4Prong = source.f4Prong;
125   fDstar = source.fDstar;
126   fTrackFilter = source.fTrackFilter;
127   fTrackFilterSoftPi = source.fTrackFilterSoftPi;
128
129   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
130   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
131   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
132   for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
133   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
134   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
135
136   return *this;
137 }
138 //----------------------------------------------------------------------------
139 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
140   // Destructor
141   if(fV1) { delete fV1; fV1=0; }
142   if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
143   if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
144 }
145 //----------------------------------------------------------------------------
146 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
147                                             TClonesArray *aodVerticesHFTClArr,
148                                             TClonesArray *aodD0toKpiTClArr,
149                                             TClonesArray *aodJPSItoEleTClArr,
150                                             TClonesArray *aodCharm3ProngTClArr,
151                                             TClonesArray *aodCharm4ProngTClArr,
152                                             TClonesArray *aodDstarTClArr)
153 {
154   // Find heavy-flavour vertex candidates
155   // Input:  ESD or AOD
156   // Output: AOD (additional branches added)
157
158
159   TString evtype = event->IsA()->GetName();
160   fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
161
162   if(!aodVerticesHFTClArr) {
163     printf("ERROR: no aodVerticesHFTClArr");
164     return;
165   }
166   if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
167     printf("ERROR: no aodD0toKpiTClArr");
168     return;
169   }
170   if(fJPSItoEle && !aodJPSItoEleTClArr) {
171     printf("ERROR: no aodJPSItoEleTClArr");
172     return;
173   }
174   if(f3Prong && !aodCharm3ProngTClArr) {
175     printf("ERROR: no aodCharm3ProngTClArr");
176     return;
177   }
178   if(f4Prong && !aodCharm4ProngTClArr) {
179     printf("ERROR: no aodCharm4ProngTClArr");
180     return;
181   }
182   if(fDstar && !aodDstarTClArr) {
183     printf("ERROR: no aodDstarTClArr");
184     return;
185   }
186
187   // delete candidates from previous event and create references
188   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0;
189   aodVerticesHFTClArr->Delete();
190   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
191   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
192   if(fD0toKpi || fDstar)   {
193     aodD0toKpiTClArr->Delete();
194     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
195   }
196   if(fJPSItoEle) {
197     aodJPSItoEleTClArr->Delete();
198     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
199   }
200   if(f3Prong) {   
201     aodCharm3ProngTClArr->Delete();
202     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
203   }
204   if(f4Prong) {
205     aodCharm4ProngTClArr->Delete();
206     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
207   }
208   if(fDstar) {
209     aodDstarTClArr->Delete();
210     iDstar = aodDstarTClArr->GetEntriesFast();
211   }
212   TClonesArray &aodD0toKpiRef      = *aodD0toKpiTClArr;
213   TClonesArray &aodJPSItoEleRef    = *aodJPSItoEleTClArr;
214   TClonesArray &aodCharm3ProngRef  = *aodCharm3ProngTClArr;
215   TClonesArray &aodCharm4ProngRef  = *aodCharm4ProngTClArr;
216   TClonesArray &aodDstarRef        = *aodDstarTClArr;
217   
218
219   AliAODRecoDecayHF2Prong *io2Prong  = 0;
220   AliAODRecoDecayHF3Prong *io3Prong  = 0;
221   AliAODRecoDecayHF4Prong *io4Prong  = 0;
222   AliAODRecoCascadeHF     *ioCascade = 0;
223
224   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries;
225   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcaCasc;
226   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
227   Bool_t   okDstar=kFALSE,okD0fromDstar=kFALSE;
228   AliESDtrack *postrack1 = 0;
229   AliESDtrack *postrack2 = 0;
230   AliESDtrack *negtrack1 = 0;
231   AliESDtrack *negtrack2 = 0;
232   AliESDtrack *trackPi   = 0;
233   Double_t dcaMax = fD0toKpiCuts[1];
234   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
235   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
236   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
237
238
239   // get Bz
240   fBzkG = (Double_t)event->GetMagneticField(); 
241
242   trkEntries = (Int_t)event->GetNumberOfTracks();
243   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
244
245   if(trkEntries<2) {
246     AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
247     return;
248   }
249     
250   const AliVVertex *primary = event->GetPrimaryVertex();
251   if(!primary) {
252     AliDebug(1," No primary vertex from tracks");
253     return;
254   }
255   TString primTitle = primary->GetTitle();
256   if(!primTitle.Contains("VertexerTracks")) {
257     AliDebug(1," No primary vertex from tracks");
258     return;
259   }
260
261   // call function that applies sigle-track selection,
262   // for displaced tracks and soft pions (both charges) for D*,
263   // and retrieves primary vertex
264   TObjArray seleTrksArray(trkEntries);
265   UChar_t  *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
266   Int_t     nSeleTrks=0;
267   SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags);
268     
269   AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
270     
271   TObjArray *twoTrackArray1    = new TObjArray(2);
272   TObjArray *twoTrackArray2    = new TObjArray(2);
273   TObjArray *twoTrackArrayCasc = new TObjArray(2);
274   TObjArray *threeTrackArray   = new TObjArray(3);
275   TObjArray *fourTrackArray    = new TObjArray(4);
276   
277   Double_t dispersion;
278
279   AliAODRecoDecayHF   *rd = 0;
280   AliAODRecoCascadeHF *rc = 0;
281
282   // LOOP ON  POSITIVE  TRACKS
283   for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
284     if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));  
285     // get track from tracks array
286     postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
287     if(postrack1->Charge()<0 || !TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
288     // LOOP ON  NEGATIVE  TRACKS
289     for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
290       if(iTrkN1==iTrkP1) continue;
291       if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));  
292       // get track from tracks array
293       negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
294       if(negtrack1->Charge()>0 || !TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
295       // back to primary vertex
296       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
297       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
298       // DCA between the two tracks
299       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
300       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
301       // Vertexing
302       twoTrackArray1->AddAt(postrack1,0);
303       twoTrackArray1->AddAt(negtrack1,1);
304       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
305       if(!vertexp1n1) { 
306         twoTrackArray1->Clear();
307         negtrack1=0; 
308         continue; 
309       }
310       if(fD0toKpi || fJPSItoEle || fDstar) { 
311         io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
312         if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI)) {
313           // add the vertex and the decay to the AOD
314           AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
315           if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
316           if(okD0) {  
317             rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
318             rd->SetSecondaryVtx(v2Prong);
319             v2Prong->SetParent(rd);
320           }
321           if(okJPSI) {
322             rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
323             rd->SetSecondaryVtx(v2Prong);
324             if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
325           }
326         }
327         if(fDstar && okD0fromDstar) {
328           if(fInputAOD) { // write references in io2Prong
329             AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
330             io2Prong->SetSecondaryVtx(vertexp1n1);
331           }
332           // create a track from the D0
333           AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
334           // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
335           for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
336             if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
337             if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
338             if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));  
339             // get track from tracks array
340             trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
341
342             trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
343             trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
344             // DCA between the two tracks
345             dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
346             // Vertexing
347             twoTrackArrayCasc->AddAt(trackPi,0);
348             twoTrackArrayCasc->AddAt(trackD0,1);
349             AliAODVertex *vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
350             if(!vertexCasc) { 
351               twoTrackArrayCasc->Clear();
352               trackPi=0; 
353               continue; 
354             }
355             ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
356             if(okDstar) {
357               if(!okD0) {
358                 // add the D0 to the AOD (if not already done)
359                 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
360                 if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
361                 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
362                 rd->SetSecondaryVtx(v2Prong);
363                 v2Prong->SetParent(rd);
364                 okD0=kTRUE; // this is done to add it only once
365               }
366               // add the vertex and the cascade to the AOD
367               AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc); 
368               if(fInputAOD) {
369                 AddDaughterRefs(vCasc,event,twoTrackArrayCasc); // add the pion
370               } else {
371                 vCasc->AddDaughter(rd); // just to fill ref #0 
372               }
373               vCasc->AddDaughter(rd); // add the D0 (in ref #1)
374               rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
375               rc->SetSecondaryVtx(vCasc);
376               vCasc->SetParent(rc);
377             }
378             twoTrackArrayCasc->Clear();
379             trackPi=0; 
380             ioCascade=NULL;
381             delete vertexCasc;
382           } // end loop on soft pi tracks
383           if(trackD0) {delete trackD0; trackD0=NULL;}
384         }
385         io2Prong=NULL;
386       }      
387       twoTrackArray1->Clear(); 
388       if(!f3Prong && !f4Prong)  { 
389         negtrack1=0; 
390         delete vertexp1n1; 
391         continue; 
392       }
393
394         
395       // 2nd LOOP  ON  POSITIVE  TRACKS 
396       for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
397         if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
398         if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));  
399         // get track from tracks array
400         postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
401         if(postrack2->Charge()<0 || !TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
402         // back to primary vertex
403         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
404         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
405         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
406         //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
407         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
408         if(dcap2n1>dcaMax) { postrack2=0; continue; }
409         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
410         if(dcap1p2>dcaMax) { postrack2=0; continue; }
411         
412         // Vertexing
413         twoTrackArray2->AddAt(postrack2,0);
414         twoTrackArray2->AddAt(negtrack1,1);
415         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
416         if(!vertexp2n1) { 
417           twoTrackArray2->Clear();
418           postrack2=0; 
419           continue; 
420         }
421         if(f3Prong) { 
422           threeTrackArray->AddAt(postrack1,0);
423           threeTrackArray->AddAt(negtrack1,1);
424           threeTrackArray->AddAt(postrack2,2);
425           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
426           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
427           if(ok3Prong) {
428             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
429             if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
430             rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
431             rd->SetSecondaryVtx(v3Prong);
432             v3Prong->SetParent(rd);
433           }
434           if(io3Prong) io3Prong=NULL; 
435         }
436         if(f4Prong) {
437           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong) 
438           for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
439             if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
440             if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
441             // get track from tracks array
442             negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
443             if(negtrack2->Charge()>0 || !TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
444             // back to primary vertex
445             postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
446             postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
447             negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
448             negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
449             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
450             if(dcap1n2>dcaMax) { negtrack2=0; continue; }
451             // Vertexing
452             fourTrackArray->AddAt(postrack1,0);
453             fourTrackArray->AddAt(negtrack1,1);
454             fourTrackArray->AddAt(postrack2,2);
455             fourTrackArray->AddAt(negtrack2,3);
456             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
457             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
458             if(ok4Prong) {
459               AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
460               if(fInputAOD) AddDaughterRefs(v4Prong,event,fourTrackArray);
461               rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
462               rd->SetSecondaryVtx(v4Prong);
463               v4Prong->SetParent(rd);
464             }
465             if(io4Prong) io4Prong=NULL; 
466             fourTrackArray->Clear();
467             negtrack2 = 0;
468           } // end loop on negative tracks
469         }
470         postrack2 = 0;
471         delete vertexp2n1;
472       } // end 2nd loop on positive tracks
473       twoTrackArray2->Clear();
474       
475       // 2nd LOOP  ON  NEGATIVE  TRACKS 
476       for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
477         if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
478         if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
479         // get track from tracks array
480         negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
481         if(negtrack2->Charge()>0 || !TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
482         // back to primary vertex
483         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
484         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
485         negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
486         //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
487         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
488         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
489         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
490         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
491         
492         // Vertexing
493         twoTrackArray2->AddAt(postrack1,0);
494         twoTrackArray2->AddAt(negtrack2,1);
495
496         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
497         if(!vertexp1n2) { 
498           twoTrackArray2->Clear();
499           negtrack2=0; 
500           continue; 
501         }
502         if(f3Prong) { 
503           threeTrackArray->AddAt(negtrack1,0);
504           threeTrackArray->AddAt(postrack1,1);
505           threeTrackArray->AddAt(negtrack2,2);
506           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
507           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
508           if(ok3Prong) {
509             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
510             if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
511             rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
512             rd->SetSecondaryVtx(v3Prong);
513             v3Prong->SetParent(rd);
514           }
515           if(io3Prong) io3Prong=NULL; 
516         }
517         negtrack2 = 0;
518         delete vertexp1n2;
519       } // end 2nd loop on negative tracks
520       twoTrackArray2->Clear();
521       
522       negtrack1 = 0;
523       delete vertexp1n1; 
524     } // end 1st loop on negative tracks
525     
526     postrack1 = 0;
527   }  // end 1st loop on positive tracks
528
529
530   AliDebug(1,Form(" Total HF vertices in event = %d;",
531                   (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
532   if(fD0toKpi) {
533     AliDebug(1,Form(" D0->Kpi in event = %d;",
534                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
535   }
536   if(fJPSItoEle) {
537     AliDebug(1,Form(" JPSI->ee in event = %d;",
538                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
539   }
540   if(f3Prong) {
541     AliDebug(1,Form(" Charm->3Prong in event = %d;",
542                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
543   }
544   if(f4Prong) {
545     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
546                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
547   }
548   if(fDstar) {
549     AliDebug(1,Form(" D*->D0pi in event = %d;\n",
550                     (Int_t)aodDstarTClArr->GetEntriesFast()));
551   }
552     
553
554   twoTrackArray1->Delete();  delete twoTrackArray1;
555   twoTrackArray2->Delete();  delete twoTrackArray2;
556   twoTrackArrayCasc->Delete();  delete twoTrackArrayCasc;
557   threeTrackArray->Clear(); 
558   threeTrackArray->Delete(); delete threeTrackArray;
559   fourTrackArray->Delete();  delete fourTrackArray;
560   delete [] seleFlags; seleFlags=NULL;
561
562   if(fInputAOD) {
563     seleTrksArray.Delete();
564   }
565
566   return;
567 }
568 //----------------------------------------------------------------------------
569 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,AliVEvent *event,
570                                              TObjArray *trkArray) const
571 {
572   // Add the AOD tracks as daughters of the vertex (TRef)
573
574   Int_t nTrks = trkArray->GetEntriesFast();
575
576   AliExternalTrackParam *track = 0;
577   AliAODTrack *aodTrack = 0;
578   Int_t id;
579
580   for(Int_t i=0; i<nTrks; i++) {
581     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
582     id = (Int_t)track->GetID();
583     //printf("---> %d\n",id);
584     if(id<0) continue; // this track is a AliAODRecoDecay
585     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
586     v->AddDaughter(aodTrack);
587   }
588
589   return;
590 }       
591 //----------------------------------------------------------------------------
592 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
593                                    TObjArray *twoTrackArray,AliVEvent *event,
594                                    AliAODVertex *secVert,
595                                    AliAODRecoDecayHF2Prong *rd2Prong,
596                                    Double_t dca,
597                                    Bool_t &okDstar) const
598 {
599   // Make the cascade as a 2Prong decay and check if it passes Dstar
600   // reconstruction cuts
601
602   okDstar = kFALSE;
603
604   Bool_t dummy1,dummy2,dummy3;
605
606   // We use Make2Prong to construct the AliAODRecoCascadeHF
607   // (which inherits from AliAODRecoDecayHF2Prong) 
608   AliAODRecoCascadeHF *theCascade = 
609     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
610                                      dummy1,dummy2,dummy3);
611   if(!theCascade) return 0x0;
612
613   // charge
614   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
615   theCascade->SetCharge(trackPi->Charge());
616
617   //--- selection cuts
618   //
619   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
620   tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
621   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
622   AliAODVertex *primVertexAOD=0;
623   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
624     // take event primary vertex
625     primVertexAOD = PrimaryVertex(); 
626     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
627     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
628   }
629   // select D*->D0pi
630   if(fDstar) {
631     Bool_t testD0=kTRUE;
632     okDstar = tmpCascade->SelectDstar(fDstarCuts,fD0fromDstarCuts,testD0);
633   }
634   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
635   tmpCascade->UnsetOwnPrimaryVtx(); 
636   delete tmpCascade; tmpCascade=NULL;
637   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
638     rd2Prong->UnsetOwnPrimaryVtx();
639     delete primVertexAOD; primVertexAOD=NULL;
640   }
641   //---
642   
643   return theCascade;
644 }
645 //-----------------------------------------------------------------------------
646 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
647                                    TObjArray *twoTrackArray,AliVEvent *event,
648                                    AliAODVertex *secVert,Double_t dca,
649                                    Bool_t &okD0,Bool_t &okJPSI,
650                                    Bool_t &okD0fromDstar) const
651 {
652   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
653   // reconstruction cuts
654   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
655
656   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
657
658   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
659
660   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
661   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
662
663   // propagate tracks to secondary vertex, to compute inv. mass
664   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
665   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
666
667   Double_t momentum[3];
668   postrack->GetPxPyPz(momentum);
669   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
670   negtrack->GetPxPyPz(momentum);
671   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
672
673
674   // invariant mass cut (try to improve coding here..)
675   Bool_t okMassCut=kFALSE;
676   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
677   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
678   if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
679   if(!okMassCut) {
680     AliDebug(2," candidate didn't pass mass cut");
681     return 0x0;    
682   }
683
684   // primary vertex to be used by this candidate
685   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
686   if(!primVertexAOD) return 0x0;
687
688   Double_t d0z0[2],covd0z0[3];
689   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
690   d0[0] = d0z0[0];
691   d0err[0] = TMath::Sqrt(covd0z0[0]);
692   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
693   d0[1] = d0z0[0];
694   d0err[1] = TMath::Sqrt(covd0z0[0]);
695
696   // create the object AliAODRecoDecayHF2Prong
697   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
698   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
699   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
700   the2Prong->SetProngIDs(2,id);
701
702
703   // select D0->Kpi
704   Int_t checkD0,checkD0bar;
705   if(fD0toKpi)   okD0          = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
706   //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
707   // select J/psi from B
708   Int_t checkJPSI;
709   if(fJPSItoEle) okJPSI        = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
710   //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
711   // select D0->Kpi from Dstar
712   if(fDstar)     okD0fromDstar = the2Prong->SelectD0(fD0fromDstarCuts,checkD0,checkD0bar);
713   //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
714
715   // remove the primary vertex (was used only for selection)
716   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the2Prong->UnsetOwnPrimaryVtx();
717
718   // get PID info from ESD
719   Double_t esdpid0[5];
720   postrack->GetESDpid(esdpid0);
721   Double_t esdpid1[5];
722   negtrack->GetESDpid(esdpid1);
723   Double_t esdpid[10];
724   for(Int_t i=0;i<5;i++) {
725     esdpid[i]   = esdpid0[i];
726     esdpid[5+i] = esdpid1[i];
727   }
728   the2Prong->SetPID(2,esdpid);
729
730   return the2Prong;  
731 }
732 //----------------------------------------------------------------------------
733 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
734                              TObjArray *threeTrackArray,AliVEvent *event,
735                              AliAODVertex *secVert,Double_t dispersion,
736                              AliAODVertex *vertexp1n1,AliAODVertex *vertexp2n1,
737                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
738                              Bool_t &ok3Prong) const
739 {
740   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
741   // reconstruction cuts 
742   // E.Bruna, F.Prino
743
744
745   ok3Prong=kFALSE;
746   if(!secVert) return 0x0; 
747
748   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
749   Double_t momentum[3];
750
751
752   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
753   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
754   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
755
756   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
757   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
758   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
759   postrack1->GetPxPyPz(momentum);
760   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
761   negtrack->GetPxPyPz(momentum);
762   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
763   postrack2->GetPxPyPz(momentum);
764   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2]; 
765
766   // invariant mass cut for D+, Ds, Lc
767   Bool_t okMassCut=kFALSE;
768   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
769   if(!okMassCut) {
770     AliDebug(2," candidate didn't pass mass cut");
771     return 0x0;    
772   }
773
774   // primary vertex to be used by this candidate
775   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
776   if(!primVertexAOD) return 0x0;
777
778   Double_t d0z0[2],covd0z0[3];
779   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
780   d0[0]=d0z0[0];
781   d0err[0] = TMath::Sqrt(covd0z0[0]);
782   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
783   d0[1]=d0z0[0];
784   d0err[1] = TMath::Sqrt(covd0z0[0]);
785   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
786   d0[2]=d0z0[0];
787   d0err[2] = TMath::Sqrt(covd0z0[0]);
788
789
790   // create the object AliAODRecoDecayHF3Prong
791   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
792   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
793   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]));
794   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]));
795   Short_t charge=(Short_t)(postrack1->Charge()*postrack2->Charge()*negtrack->Charge());
796   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
797   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
798   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
799   the3Prong->SetProngIDs(3,id);
800
801
802   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
803   if(f3Prong) {
804     ok3Prong = kFALSE;
805     Int_t ok1,ok2;
806     if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
807     if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE;
808     if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
809   }
810   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
811
812   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the3Prong->UnsetOwnPrimaryVtx();
813
814   // get PID info from ESD
815   Double_t esdpid0[5];
816   postrack1->GetESDpid(esdpid0);
817   Double_t esdpid1[5];
818   negtrack->GetESDpid(esdpid1);
819   Double_t esdpid2[5];
820   postrack2->GetESDpid(esdpid2);
821   
822   Double_t esdpid[15];
823   for(Int_t i=0;i<5;i++) {
824     esdpid[i]    = esdpid0[i];
825     esdpid[5+i]  = esdpid1[i];
826     esdpid[10+i] = esdpid2[i];
827   }
828   the3Prong->SetPID(3,esdpid);
829
830   return the3Prong;
831 }
832 //----------------------------------------------------------------------------
833 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
834                              TObjArray *fourTrackArray,AliVEvent *event,
835                              AliAODVertex *secVert,
836                              AliAODVertex *vertexp1n1,AliAODVertex *vertexp2n1,
837                              Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
838                              Bool_t &ok4Prong) const
839 {
840   // Make 4Prong candidates and check if they pass D0toKpipipi
841   // reconstruction cuts
842   // G.E.Bruno, R.Romita
843
844   ok4Prong=kFALSE;
845   if(!secVert) return 0x0; 
846
847   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
848
849   px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
850
851   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
852   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
853   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
854   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
855
856   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
857   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
858   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
859   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
860
861   Double_t momentum[3];
862   postrack1->GetPxPyPz(momentum);
863   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
864   negtrack1->GetPxPyPz(momentum);
865   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
866   postrack2->GetPxPyPz(momentum);
867   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
868   negtrack2->GetPxPyPz(momentum);
869   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
870
871   // invariant mass cut for rho or D0 (try to improve coding here..)
872   //Bool_t okMassCut=kFALSE;
873   //if(!okMassCut) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
874   //if(!okMassCut) {
875   //  if(fDebug) printf(" candidate didn't pass mass cut\n");
876   //  return 0x0;
877   //}
878
879   // primary vertex to be used by this candidate
880   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
881   if(!primVertexAOD) return 0x0;
882
883   /*
884     Double_t d0z0[2],covd0z0[3];
885     postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
886     negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
887     postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
888     negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
889   */
890
891   // create the object AliAODRecoDecayHF4Prong
892   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
893   Double_t dca[6]={0.,0.,0.,0.,0.,0.}; //  modify it
894   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]));
895   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]));
896   Double_t dist14=0.; // to be implemented
897   Double_t dist34=0.; // to be implemented
898   Short_t charge=0;
899   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge);
900   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
901   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
902   the4Prong->SetProngIDs(4,id);
903
904
905   // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available
906   // select D0->Kpipipi
907   //Int_t checkD0,checkD0bar;   
908   // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar); 
909   ok4Prong=kFALSE;  //for the time being ...
910
911   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the4Prong->UnsetOwnPrimaryVtx();
912
913   // get PID info from ESD
914   Double_t esdpid0[5];
915   postrack1->GetESDpid(esdpid0);
916   Double_t esdpid1[5];
917   negtrack1->GetESDpid(esdpid1);
918   Double_t esdpid2[5];
919   postrack2->GetESDpid(esdpid2);
920   Double_t esdpid3[5];
921   negtrack2->GetESDpid(esdpid3);
922
923   Double_t esdpid[20];
924   for(Int_t i=0;i<5;i++) {
925     esdpid[i]    = esdpid0[i];
926     esdpid[5+i]  = esdpid1[i];
927     esdpid[10+i] = esdpid2[i];
928     esdpid[15+i] = esdpid3[i];
929   }
930   the4Prong->SetPID(4,esdpid);
931
932   return the4Prong;
933 }
934 //-----------------------------------------------------------------------------
935 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(TObjArray *trkArray,
936                                                     AliVEvent *event) const
937 {
938   // Returns primary vertex to be used for this candidate
939  
940   AliESDVertex *vertexESD = 0;
941   AliAODVertex *vertexAOD = 0;
942
943
944   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) { 
945     // primary vertex from the input event
946     
947     vertexESD = new AliESDVertex(*fV1);
948
949   } else {
950     // primary vertex specific to this candidate
951
952     Int_t nTrks = trkArray->GetEntriesFast();
953     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
954
955     if(fRecoPrimVtxSkippingTrks) { 
956       // recalculating the vertex
957       
958       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
959         Float_t diamondcovxy[3];
960         event->GetDiamondCovXY(diamondcovxy);
961         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
962         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.};
963         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
964         vertexer->SetVtxStart(diamond);
965         delete diamond; diamond=NULL;
966         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) 
967           vertexer->SetOnlyFitter();
968       }
969       Int_t skipped[10];
970       Int_t nTrksToSkip=0,id;
971       AliExternalTrackParam *t = 0;
972       for(Int_t i=0; i<nTrks; i++) {
973         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
974         id = (Int_t)t->GetID();
975         if(id<0) continue;
976         skipped[nTrksToSkip++] = id;
977       }
978       vertexer->SetSkipTracks(nTrksToSkip,skipped);
979       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
980       
981     } else if(fRmTrksFromPrimVtx) { 
982       // removing the prongs tracks
983       
984       TObjArray rmArray(nTrks);
985       UShort_t *rmId = new UShort_t[nTrks];
986       AliESDtrack *esdTrack = 0;
987       AliESDtrack *t = 0;
988       for(Int_t i=0; i<nTrks; i++) {
989         t = (AliESDtrack*)trkArray->UncheckedAt(i);
990         esdTrack = new AliESDtrack(*t);
991         rmArray.AddLast(esdTrack);
992         if(esdTrack->GetID()>=0) {
993           rmId[i]=(UShort_t)esdTrack->GetID();
994         } else {
995           rmId[i]=9999;
996         }
997       }
998       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
999       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1000       delete [] rmId; rmId=NULL;
1001       rmArray.Delete();
1002       
1003     }
1004
1005     if(!vertexESD) return vertexAOD;
1006     if(vertexESD->GetNContributors()<=0) { 
1007       AliDebug(2,"vertexing failed"); 
1008       delete vertexESD; vertexESD=NULL;
1009       return vertexAOD;
1010     }
1011
1012     delete vertexer; vertexer=NULL;
1013
1014   }
1015
1016   // convert to AliAODVertex
1017   Double_t pos[3],cov[6],chi2perNDF;
1018   vertexESD->GetXYZ(pos); // position
1019   vertexESD->GetCovMatrix(cov); //covariance matrix
1020   chi2perNDF = vertexESD->GetChi2toNDF();
1021   delete vertexESD; vertexESD=NULL;
1022
1023   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1024
1025   return vertexAOD;
1026 }
1027 //-----------------------------------------------------------------------------
1028 void AliAnalysisVertexingHF::PrintStatus() const {
1029   // Print parameters being used
1030
1031   //printf("Preselections:\n");
1032   //   fTrackFilter->Dump();
1033   if(fSecVtxWithKF) {
1034     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1035   } else {
1036     printf("Secondary vertex with AliVertexerTracks\n");
1037   }
1038   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1039   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1040   if(fD0toKpi) {
1041     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1042     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
1043     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
1044     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
1045     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
1046     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
1047     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
1048     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
1049     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
1050     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
1051   }
1052   if(fDstar) {
1053     printf("    |M-MD*| [GeV]    < %f\n",fDstarCuts[0]);
1054     printf("    |M_Kpipi-M_Kpi-(MD*-MD0)| [GeV]  < %f\n",fDstarCuts[1]);
1055     printf("    pTpisoft [GeV/c]    > %f\n",fDstarCuts[2]);
1056     printf("    pTpisoft [GeV/c]    < %f\n",fDstarCuts[3]);
1057     printf("    Theta(pisoft,D0plane) < %f\n",fDstarCuts[4]);
1058     printf("Reconstruct D*->D0pi candidates with cuts:\n");
1059     printf("   D0 from D* cuts:\n");
1060     printf("    |M-MD0| [GeV]    < %f\n",fD0fromDstarCuts[0]);
1061     printf("    dca    [cm]  < %f\n",fD0fromDstarCuts[1]);
1062     printf("    cosThetaStar     < %f\n",fD0fromDstarCuts[2]);
1063     printf("    pTK     [GeV/c]    > %f\n",fD0fromDstarCuts[3]);
1064     printf("    pTpi    [GeV/c]    > %f\n",fD0fromDstarCuts[4]);
1065     printf("    |d0K|  [cm]  < %f\n",fD0fromDstarCuts[5]);
1066     printf("    |d0pi| [cm]  < %f\n",fD0fromDstarCuts[6]);
1067     printf("    d0d0  [cm^2] < %f\n",fD0fromDstarCuts[7]);
1068     printf("    cosThetaPoint    > %f\n",fD0fromDstarCuts[8]);
1069   }
1070   if(fJPSItoEle) {
1071     printf("Reconstruct J/psi from B candidates with cuts:\n");
1072     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
1073     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
1074     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
1075     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
1076     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
1077     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
1078     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
1079     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
1080     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
1081   }
1082   if(f3Prong) {
1083     printf("Reconstruct 3 prong candidates.\n");
1084     printf("  D+->Kpipi cuts:\n");
1085     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
1086     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
1087     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
1088     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
1089     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
1090     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
1091     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
1092     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1093     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1094     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
1095     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
1096     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
1097     printf("  Ds->KKpi cuts:\n");
1098     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
1099     printf("    pTK     [GeV/c]    > %f\n",fDsCuts[1]);
1100     printf("    pTPi    [GeV/c]    > %f\n",fDsCuts[2]);
1101     printf("    |d0K|  [cm]  > %f\n",fDsCuts[3]);
1102     printf("    |d0Pi| [cm]  > %f\n",fDsCuts[4]);
1103     printf("    dist12    [cm]  < %f\n",fDsCuts[5]);
1104     printf("    sigmavert [cm]   < %f\n",fDsCuts[6]);
1105     printf("    dist prim-sec [cm] > %f\n",fDsCuts[7]);
1106     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
1107     printf("    cosThetaPoint    > %f\n",fDsCuts[9]);
1108     printf("    Sum d0^2 [cm^2]  > %f\n",fDsCuts[10]);
1109     printf("    dca cut [cm]  < %f\n",fDsCuts[11]);
1110     printf("    Inv. Mass  phi/K0* [GeV]  < %f\n",fDsCuts[12]);
1111     printf("  Lc->pKpi cuts:\n");
1112     printf("    |M-MLc| [GeV]    < %f\n",fLcCuts[0]);
1113     printf("    pTP     [GeV/c]    > %f\n",fLcCuts[1]);
1114     printf("    pTPi and pTK [GeV/c]    > %f\n",fLcCuts[2]);
1115     printf("    |d0P|  [cm]  > %f\n",fLcCuts[3]);
1116     printf("    |d0Pi| and |d0K| [cm]  > %f\n",fLcCuts[4]);
1117     printf("    dist12    [cm]  < %f\n",fLcCuts[5]);
1118     printf("    sigmavert [cm]   < %f\n",fLcCuts[6]);
1119     printf("    dist prim-sec [cm] > %f\n",fLcCuts[7]);
1120     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
1121     printf("    cosThetaPoint    > %f\n",fLcCuts[9]);
1122     printf("    Sum d0^2 [cm^2]  > %f\n",fLcCuts[10]);
1123     printf("    dca cut [cm]  < %f\n",fLcCuts[11]);
1124     printf("  Ds->KKpi cuts:\n");
1125   }
1126
1127   return;
1128 }
1129 //-----------------------------------------------------------------------------
1130 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1131                                                                  Double_t &dispersion,Bool_t useTRefArray) const
1132 {
1133   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1134
1135   AliESDVertex *vertexESD = 0;
1136   AliAODVertex *vertexAOD = 0;
1137
1138   if(!fSecVtxWithKF) { // AliVertexerTracks
1139
1140     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1141     vertexer->SetVtxStart(fV1);
1142     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1143     delete vertexer; vertexer=NULL;
1144
1145     if(!vertexESD) return vertexAOD;
1146
1147     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
1148       AliDebug(2,"vertexing failed"); 
1149       delete vertexESD; vertexESD=NULL;
1150       return vertexAOD;
1151     }
1152
1153   } else { // Kalman Filter vertexer (AliKFParticle)
1154
1155     AliKFParticle::SetField(fBzkG);
1156
1157     AliKFVertex vertexKF;
1158
1159     Int_t nTrks = trkArray->GetEntriesFast();
1160     for(Int_t i=0; i<nTrks; i++) {
1161       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1162       AliKFParticle daughterKF(*esdTrack,211);
1163       vertexKF.AddDaughter(daughterKF);
1164     }
1165     vertexESD = new AliESDVertex(vertexKF.Parameters(),
1166                                  vertexKF.CovarianceMatrix(),
1167                                  vertexKF.GetChi2(),
1168                                  vertexKF.GetNContributors());
1169
1170   }
1171
1172   // convert to AliAODVertex
1173   Double_t pos[3],cov[6],chi2perNDF;
1174   vertexESD->GetXYZ(pos); // position
1175   vertexESD->GetCovMatrix(cov); //covariance matrix
1176   chi2perNDF = vertexESD->GetChi2toNDF();
1177   dispersion = vertexESD->GetDispersion();
1178   delete vertexESD; vertexESD=NULL;
1179
1180   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1181   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1182
1183   return vertexAOD;
1184 }
1185 //-----------------------------------------------------------------------------
1186 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1187                                              Int_t nprongs,
1188                                              Double_t *px,
1189                                              Double_t *py,
1190                                              Double_t *pz) const {
1191   // Check invariant mass cut
1192
1193   Short_t dummycharge=0;
1194   Double_t *dummyd0 = new Double_t[nprongs];
1195   for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1196   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1197   delete [] dummyd0;
1198
1199   UInt_t pdg2[2],pdg3[3];
1200   Double_t mPDG,minv;
1201
1202   Bool_t retval=kFALSE;
1203   switch (decay) 
1204     { 
1205     case 0:                  // D0->Kpi
1206       pdg2[0]=211; pdg2[1]=321;
1207       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1208       minv = rd->InvMass(nprongs,pdg2);
1209       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1210       pdg2[0]=321; pdg2[1]=211;
1211       minv = rd->InvMass(nprongs,pdg2);
1212       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1213       break;
1214     case 1:                  // JPSI->ee
1215       pdg2[0]=11; pdg2[1]=11;
1216       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1217       minv = rd->InvMass(nprongs,pdg2);
1218       if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1219       break;
1220     case 2:                  // D+->Kpipi
1221       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1222       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1223       minv = rd->InvMass(nprongs,pdg3);
1224       if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1225                             // Ds+->KKpi
1226       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1227       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1228       minv = rd->InvMass(nprongs,pdg3);
1229       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1230       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1231       minv = rd->InvMass(nprongs,pdg3);
1232       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1233                             // Lc->pKpi
1234       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1235       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1236       minv = rd->InvMass(nprongs,pdg3);
1237       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1238       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1239       minv = rd->InvMass(nprongs,pdg3);
1240       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE; 
1241       break;
1242     case 3:                  // D*->D0pi
1243       pdg2[0]=421; pdg2[1]=211;
1244       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1245       minv = rd->InvMass(nprongs,pdg2);
1246       if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
1247       break;
1248     default:
1249       printf("SelectInvMass(): wrong decay selection\n");
1250       break;
1251     }
1252
1253   delete rd;
1254
1255   return retval;
1256 }
1257 //-----------------------------------------------------------------------------
1258 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(AliVEvent *event,
1259                                    TObjArray &seleTrksArray,Int_t &nSeleTrks,
1260                                                        UChar_t *seleFlags)
1261 {
1262   // Apply single-track preselection.
1263   // Fill a TObjArray with selected tracks (for displaced vertices or
1264   // soft pion from D*). Selection flag stored in seleFlags.
1265   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1266   // In case of AOD input, also fill fAODMap for track index<->ID
1267
1268   const AliVVertex *vprimary = event->GetPrimaryVertex();
1269
1270   if(fV1) { delete fV1; fV1=NULL; }
1271   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
1272
1273   Int_t nindices=0;
1274   UShort_t *indices = 0;
1275   Double_t pos[3],cov[6];
1276
1277   if(!fInputAOD) { // ESD
1278     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
1279   } else {         // AOD
1280     vprimary->GetXYZ(pos);
1281     vprimary->GetCovarianceMatrix(cov);
1282     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1283     indices = new UShort_t[event->GetNumberOfTracks()];
1284     fAODMap = new Int_t[100000];
1285   }
1286
1287
1288   Int_t entries = (Int_t)event->GetNumberOfTracks();
1289   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
1290   nSeleTrks=0;
1291  
1292   // transfer ITS tracks from event to arrays
1293   for(Int_t i=0; i<entries; i++) {
1294     AliVTrack *track = (AliVTrack*)event->GetTrack(i);
1295
1296     if(fInputAOD) {
1297       AliAODTrack *aodt = (AliAODTrack*)track;
1298       if(aodt->GetUsedForPrimVtxFit()) { 
1299         indices[nindices]=aodt->GetID(); nindices++; 
1300       }
1301       fAODMap[(Int_t)aodt->GetID()] = i;
1302     }
1303
1304     AliESDtrack *esdt = 0;
1305     if(!fInputAOD) {
1306       esdt = (AliESDtrack*)track;
1307     } else {
1308       esdt = new AliESDtrack(track);
1309     }
1310
1311     // single track selection
1312     okDisplaced=kFALSE; okSoftPi=kFALSE;
1313     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
1314       seleTrksArray.AddLast(esdt);
1315       seleFlags[nSeleTrks]=0;
1316       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
1317       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
1318       nSeleTrks++;
1319     } else {
1320       if(fInputAOD) delete esdt; 
1321       esdt = NULL;
1322       continue;
1323     } 
1324
1325   } // end loop on tracks
1326
1327   // primary vertex from AOD
1328   if(fInputAOD) {
1329     delete fV1; fV1=NULL;
1330     vprimary->GetXYZ(pos);
1331     vprimary->GetCovarianceMatrix(cov);
1332     Double_t chi2toNDF = vprimary->GetChi2perNDF();
1333     Int_t ncontr=nindices;
1334     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
1335     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.); 
1336     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
1337     fV1->SetTitle(vprimary->GetTitle());
1338     fV1->SetIndices(nindices,indices);
1339   }
1340   if(indices) { delete [] indices; indices=NULL; }
1341
1342
1343   return;
1344 }
1345 //-----------------------------------------------------------------------------
1346 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1347                                    Double_t cut2,Double_t cut3,Double_t cut4,
1348                                    Double_t cut5,Double_t cut6,
1349                                    Double_t cut7,Double_t cut8) 
1350 {
1351   // Set the cuts for D0 selection
1352   fD0toKpiCuts[0] = cut0;
1353   fD0toKpiCuts[1] = cut1;
1354   fD0toKpiCuts[2] = cut2;
1355   fD0toKpiCuts[3] = cut3;
1356   fD0toKpiCuts[4] = cut4;
1357   fD0toKpiCuts[5] = cut5;
1358   fD0toKpiCuts[6] = cut6;
1359   fD0toKpiCuts[7] = cut7;
1360   fD0toKpiCuts[8] = cut8;
1361
1362   return;
1363 }
1364 //-----------------------------------------------------------------------------
1365 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9]) 
1366 {
1367   // Set the cuts for D0 selection
1368
1369   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1370
1371   return;
1372 }
1373 //-----------------------------------------------------------------------------
1374 void AliAnalysisVertexingHF::SetD0fromDstarCuts(Double_t cut0,Double_t cut1,
1375                                    Double_t cut2,Double_t cut3,Double_t cut4,
1376                                    Double_t cut5,Double_t cut6,
1377                                    Double_t cut7,Double_t cut8) 
1378 {
1379   // Set the cuts for D0 from D* selection
1380   fD0fromDstarCuts[0] = cut0;
1381   fD0fromDstarCuts[1] = cut1;
1382   fD0fromDstarCuts[2] = cut2;
1383   fD0fromDstarCuts[3] = cut3;
1384   fD0fromDstarCuts[4] = cut4;
1385   fD0fromDstarCuts[5] = cut5;
1386   fD0fromDstarCuts[6] = cut6;
1387   fD0fromDstarCuts[7] = cut7;
1388   fD0fromDstarCuts[8] = cut8;
1389
1390   return;
1391 }
1392 //-----------------------------------------------------------------------------
1393 void AliAnalysisVertexingHF::SetD0fromDstarCuts(const Double_t cuts[9]) 
1394 {
1395   // Set the cuts for D0 from D* selection
1396
1397   for(Int_t i=0; i<9; i++) fD0fromDstarCuts[i] = cuts[i];
1398
1399   return;
1400 }
1401 //-----------------------------------------------------------------------------
1402 void AliAnalysisVertexingHF::SetDstarCuts(Double_t cut0,Double_t cut1,
1403                                           Double_t cut2,Double_t cut3,
1404                                           Double_t cut4)
1405 {
1406   // Set the cuts for D* selection
1407   fDstarCuts[0] = cut0;
1408   fDstarCuts[1] = cut1;
1409   fDstarCuts[2] = cut2;
1410   fDstarCuts[3] = cut3;
1411   fDstarCuts[4] = cut4;
1412
1413   return;
1414 }
1415 //-----------------------------------------------------------------------------
1416 void AliAnalysisVertexingHF::SetDstarCuts(const Double_t cuts[5])
1417 {
1418   // Set the cuts for D* selection
1419
1420   for(Int_t i=0; i<5; i++) fDstarCuts[i] = cuts[i];
1421
1422   return;
1423 }
1424 //-----------------------------------------------------------------------------
1425 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1426                                    Double_t cut2,Double_t cut3,Double_t cut4,
1427                                    Double_t cut5,Double_t cut6,
1428                                    Double_t cut7,Double_t cut8) 
1429 {
1430   // Set the cuts for J/psi from B selection
1431   fBtoJPSICuts[0] = cut0;
1432   fBtoJPSICuts[1] = cut1;
1433   fBtoJPSICuts[2] = cut2;
1434   fBtoJPSICuts[3] = cut3;
1435   fBtoJPSICuts[4] = cut4;
1436   fBtoJPSICuts[5] = cut5;
1437   fBtoJPSICuts[6] = cut6;
1438   fBtoJPSICuts[7] = cut7;
1439   fBtoJPSICuts[8] = cut8;
1440
1441   return;
1442 }
1443 //-----------------------------------------------------------------------------
1444 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9]) 
1445 {
1446   // Set the cuts for J/psi from B selection
1447
1448   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1449
1450   return;
1451 }
1452 //-----------------------------------------------------------------------------
1453 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1454                                    Double_t cut2,Double_t cut3,Double_t cut4,
1455                                    Double_t cut5,Double_t cut6,
1456                                    Double_t cut7,Double_t cut8,
1457                                    Double_t cut9,Double_t cut10,Double_t cut11)
1458 {
1459   // Set the cuts for Dplus->Kpipi selection
1460   fDplusCuts[0] = cut0;
1461   fDplusCuts[1] = cut1;
1462   fDplusCuts[2] = cut2;
1463   fDplusCuts[3] = cut3;
1464   fDplusCuts[4] = cut4;
1465   fDplusCuts[5] = cut5;
1466   fDplusCuts[6] = cut6;
1467   fDplusCuts[7] = cut7;
1468   fDplusCuts[8] = cut8;
1469   fDplusCuts[9] = cut9;
1470   fDplusCuts[10] = cut10;
1471   fDplusCuts[11] = cut11;
1472
1473   return;
1474 }
1475 //-----------------------------------------------------------------------------
1476 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12]) 
1477 {
1478   // Set the cuts for Dplus->Kpipi selection
1479
1480   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1481
1482   return;
1483 }
1484 //-----------------------------------------------------------------------------
1485 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
1486                                    Double_t cut2,Double_t cut3,Double_t cut4,
1487                                    Double_t cut5,Double_t cut6,
1488                                    Double_t cut7,Double_t cut8,
1489                                    Double_t cut9,Double_t cut10,
1490                                    Double_t cut11,Double_t cut12)
1491 {
1492   // Set the cuts for Ds->KKpi selection
1493   fDsCuts[0] = cut0;
1494   fDsCuts[1] = cut1;
1495   fDsCuts[2] = cut2;
1496   fDsCuts[3] = cut3;
1497   fDsCuts[4] = cut4;
1498   fDsCuts[5] = cut5;
1499   fDsCuts[6] = cut6;
1500   fDsCuts[7] = cut7;
1501   fDsCuts[8] = cut8;
1502   fDsCuts[9] = cut9;
1503   fDsCuts[10] = cut10;
1504   fDsCuts[11] = cut11;
1505   fDsCuts[12] = cut12;
1506
1507   return;
1508 }
1509 //-----------------------------------------------------------------------------
1510 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13]) 
1511 {
1512   // Set the cuts for Ds->KKpi selection
1513
1514   for(Int_t i=0; i<13; i++) fDsCuts[i] = cuts[i];
1515
1516   return;
1517 }
1518 //-----------------------------------------------------------------------------
1519 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
1520                                    Double_t cut2,Double_t cut3,Double_t cut4,
1521                                    Double_t cut5,Double_t cut6,
1522                                    Double_t cut7,Double_t cut8,
1523                                    Double_t cut9,Double_t cut10,Double_t cut11)
1524 {
1525   // Set the cuts for Lc->pKpi selection
1526   fLcCuts[0] = cut0;
1527   fLcCuts[1] = cut1;
1528   fLcCuts[2] = cut2;
1529   fLcCuts[3] = cut3;
1530   fLcCuts[4] = cut4;
1531   fLcCuts[5] = cut5;
1532   fLcCuts[6] = cut6;
1533   fLcCuts[7] = cut7;
1534   fLcCuts[8] = cut8;
1535   fLcCuts[9] = cut9;
1536   fLcCuts[10] = cut10;
1537   fLcCuts[11] = cut11;
1538
1539   return;
1540 }
1541 //-----------------------------------------------------------------------------
1542 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12]) 
1543 {
1544   // Set the cuts for Lc->pKpi selection
1545
1546   for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
1547
1548   return;
1549 }
1550 //-----------------------------------------------------------------------------
1551 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
1552                                              Bool_t &okDisplaced,Bool_t &okSoftPi) const 
1553 {
1554   // Check if track passes some kinematical cuts  
1555
1556   // this is needed to store the impact parameters
1557   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
1558
1559   UInt_t selectInfo;
1560   //
1561   // Track selection, displaced tracks
1562   selectInfo = 0; 
1563   if(fTrackFilter) {
1564     selectInfo = fTrackFilter->IsSelected(trk);
1565   }
1566   if(selectInfo) okDisplaced=kTRUE;
1567   // Track selection, soft pions
1568   selectInfo = 0; 
1569   if(fDstar && fTrackFilterSoftPi) {
1570     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
1571   }
1572   if(selectInfo) okSoftPi=kTRUE;
1573
1574   if(okDisplaced || okSoftPi) return kTRUE;
1575
1576   return kFALSE;
1577 }
1578 //-----------------------------------------------------------------------------