Added prong IDs (Andrea)
[u/mrichter/AliRoot.git] / PWG3 / AliAnalysisVertexingHF.cxx
CommitLineData
6a213b59 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 store as objects deriving from AliAODRecoDecay.
19// An example of usage can be found in the macro AliVertexingHFTest.C.
20// Can be used as a task of AliAnalysisManager by means of the interface
21// class AliAnalysisTaskVertexingHF.
22//
23// Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita
24//----------------------------------------------------------------------------
25#include <TTree.h>
b82f6d67 26#include <TFile.h>
6a213b59 27#include <TDatabasePDG.h>
28#include <TString.h>
5e6f195c 29#include "AliESDEvent.h"
6a213b59 30#include "AliVertexerTracks.h"
b82f6d67 31#include "AliKFParticle.h"
6a213b59 32#include "AliESDVertex.h"
33#include "AliESDv0.h"
34#include "AliAODRecoDecay.h"
35#include "AliAODRecoDecayHF.h"
36#include "AliAODRecoDecayHF2Prong.h"
37#include "AliAODRecoDecayHF3Prong.h"
38#include "AliAODRecoDecayHF4Prong.h"
39#include "AliAnalysisVertexingHF.h"
40
41ClassImp(AliAnalysisVertexingHF)
42
43//----------------------------------------------------------------------------
44AliAnalysisVertexingHF::AliAnalysisVertexingHF():
b82f6d67 45fBzkG(0.),
46fSecVtxWithKF(kFALSE),
6a213b59 47fRecoPrimVtxSkippingTrks(kFALSE),
48fRmTrksFromPrimVtx(kFALSE),
49fV1(0x0),
50fDebug(0),
b82f6d67 51fD0toKpi(kTRUE),
52fJPSItoEle(kTRUE),
53f3Prong(kTRUE),
54f4Prong(kTRUE),
55fITSrefit(kFALSE),
56fBothSPD(kTRUE),
57fMinITSCls(5),
58fMinPtCut(0.),
59fMind0rphiCut(0.)
6a213b59 60{
61 // Default constructor
62
63 SetD0toKpiCuts();
64 SetBtoJPSICuts();
65 SetDplusCuts();
b82f6d67 66 SetDsCuts();
67 SetLcCuts();
6a213b59 68}
13977a79 69//--------------------------------------------------------------------------
70AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
71TNamed(source),
b82f6d67 72fBzkG(source.fBzkG),
73fSecVtxWithKF(source.fSecVtxWithKF),
13977a79 74fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
75fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
76fV1(source.fV1),
77fDebug(source.fDebug),
78fD0toKpi(source.fD0toKpi),
79fJPSItoEle(source.fJPSItoEle),
80f3Prong(source.f3Prong),
81f4Prong(source.f4Prong),
82fITSrefit(source.fITSrefit),
83fBothSPD(source.fBothSPD),
84fMinITSCls(source.fMinITSCls),
85fMinPtCut(source.fMinPtCut),
86fMind0rphiCut(source.fMind0rphiCut)
87{
88 //
89 // Copy constructor
90 //
91 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
92 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
93 for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
b82f6d67 94 for(Int_t i=0; i<1; i++) fDsCuts[i]=source.fDsCuts[i];
95 for(Int_t i=0; i<1; i++) fLcCuts[i]=source.fLcCuts[i];
13977a79 96}
97//--------------------------------------------------------------------------
98AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
99{
100 //
101 // assignment operator
102 //
103 if(&source == this) return *this;
b82f6d67 104 fBzkG = source.fBzkG;
105 fSecVtxWithKF = source.fSecVtxWithKF;
13977a79 106 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
107 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
108 fV1 = source.fV1;
109 fDebug = source.fDebug;
110 fD0toKpi = source.fD0toKpi;
111 fJPSItoEle = source.fJPSItoEle;
112 f3Prong = source.f3Prong;
113 f4Prong = source.f4Prong;
114 fITSrefit = source.fITSrefit;
115 fBothSPD = source.fBothSPD;
116 fMinITSCls = source.fMinITSCls;
117 fMinPtCut = source.fMinPtCut;
118 fMind0rphiCut = source.fMind0rphiCut;
119
120 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
121 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
122 for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
b82f6d67 123 for(Int_t i=0; i<1; i++) fDsCuts[i]=source.fDsCuts[i];
124 for(Int_t i=0; i<1; i++) fLcCuts[i]=source.fLcCuts[i];
13977a79 125
126 return *this;
127}
6a213b59 128//----------------------------------------------------------------------------
129AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
130 // Destructor
13977a79 131 if(fV1) { delete fV1; fV1=0; }
6a213b59 132}
133//----------------------------------------------------------------------------
260836a9 134void AliAnalysisVertexingHF::FindCandidates(AliESDEvent *esd,TTree *treeout[])
6a213b59 135{
136 // Find heavy-flavour vertex candidates
6a213b59 137
b82f6d67 138 AliAODRecoDecayHF2Prong *io2Prong = new AliAODRecoDecayHF2Prong();
139 AliAODRecoDecayHF3Prong *io3Prong = new AliAODRecoDecayHF3Prong();
140 AliAODRecoDecayHF4Prong *io4Prong = new AliAODRecoDecayHF4Prong();
6a213b59 141 Int_t itree=0;
142 Int_t itreeD0toKpi=-1,itreeJPSItoEle=-1,itree3Prong=-1,itree4Prong=-1;
143 Int_t initEntriesD0toKpi=0,initEntriesJPSItoEle=0,initEntries3Prong=0,initEntries4Prong=0;
144 if(fD0toKpi) {
145 itreeD0toKpi=itree;
260836a9 146 treeout[itree]->SetBranchAddress("D0toKpi",&io2Prong);
6a213b59 147 itree++;
260836a9 148 initEntriesD0toKpi = treeout[itreeD0toKpi]->GetEntries();
6a213b59 149 }
150 if(fJPSItoEle) {
151 itreeJPSItoEle=itree;
260836a9 152 treeout[itree]->SetBranchAddress("JPSItoEle",&io2Prong);
6a213b59 153 itree++;
260836a9 154 initEntriesJPSItoEle = treeout[itreeJPSItoEle]->GetEntries();
6a213b59 155 }
156 if(f3Prong) {
157 itree3Prong=itree;
260836a9 158 treeout[itree]->SetBranchAddress("Charmto3Prong",&io3Prong);
6a213b59 159 itree++;
260836a9 160 initEntries3Prong = treeout[itree3Prong]->GetEntries();
6a213b59 161 }
162 if(f4Prong) {
163 itree4Prong=itree;
260836a9 164 treeout[itree]->SetBranchAddress("D0to4Prong",&io4Prong);
6a213b59 165 itree++;
260836a9 166 initEntries4Prong = treeout[itree4Prong]->GetEntries();
6a213b59 167 }
168 delete io2Prong; io2Prong = NULL;
169 delete io3Prong; io3Prong = NULL;
170 delete io4Prong; io4Prong = NULL;
171
172 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries;
173 Int_t nTrksP=0,nTrksN=0;
174 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2;
175 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
176 AliESDtrack *postrack1 = 0;
177 AliESDtrack *postrack2 = 0;
178 AliESDtrack *negtrack1 = 0;
179 AliESDtrack *negtrack2 = 0;
180 Double_t dcaMax = fD0toKpiCuts[1];
181 if(dcaMax<fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
182 if(dcaMax<fDplusCuts[11]) dcaMax=fDplusCuts[11];
183 if(fDebug) printf(" dca cut set to %f cm\n",dcaMax);
184
6a213b59 185 Int_t ev = (Int_t)esd->GetEventNumberInFile();
186 printf("--- Finding candidates in event %d\n",ev);
187
b82f6d67 188 fBzkG = (Double_t)esd->GetMagneticField();
6a213b59 189
190 trkEntries = (Int_t)esd->GetNumberOfTracks();
191 printf(" Number of tracks: %d\n",trkEntries);
192 if(trkEntries<2) return;
193
b82f6d67 194 // retrieve primary vertex from the AliESDEvent
6a213b59 195 if(!esd->GetPrimaryVertex()) {
196 printf(" No vertex in AliESD\n");
197 return;
198 }
199 AliESDVertex copy(*(esd->GetPrimaryVertex()));
200 SetPrimaryVertex(&copy);
6a213b59 201
202 // call function which applies sigle-track selection and
203 // separetes positives and negatives
204 TObjArray trksP(trkEntries/2); // will become TClonesArray
205 TObjArray trksN(trkEntries/2); // will become TClonesArray
206 SelectTracks(esd,trksP,nTrksP,
207 trksN,nTrksN);
208
209 printf(" Pos. tracks: %d Neg. tracks: %d\n",nTrksP,nTrksN);
210
211 TObjArray *twoTrackArray1 = new TObjArray(2);
212 TObjArray *twoTrackArray2 = new TObjArray(2);
213 TObjArray *threeTrackArray = new TObjArray(3);
214 TObjArray *fourTrackArray = new TObjArray(4);
215
216 // LOOP ON POSITIVE TRACKS
217 for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) {
218 if(fDebug) if(iTrkP1%1==0) printf(" Processing positive track number %d of %d\n",iTrkP1,nTrksP);
219 // get track from track array
220 postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1);
221
222 // LOOP ON NEGATIVE TRACKS
223 for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) {
224 if(fDebug) if(iTrkN1%1==0) printf(" Processing negative track number %d of %d\n",iTrkN1,nTrksN);
225 // get track from tracks array
226 negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1);
227 // DCA between the two tracks
b82f6d67 228 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
6a213b59 229 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
b82f6d67 230 // Vertexing
6a213b59 231 twoTrackArray1->AddAt(postrack1,0);
232 twoTrackArray1->AddAt(negtrack1,1);
b82f6d67 233 AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1);
6a213b59 234 if(vertexp1n1->GetNContributors()!=2) {
235 if(fDebug) printf("two-track vertexing failed\n");
806873d4 236 twoTrackArray1->Clear();
237 delete vertexp1n1;
238 negtrack1=0;
239 continue;
6a213b59 240 }
241 if(fD0toKpi || fJPSItoEle) {
242 io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI);
260836a9 243 if(okD0) treeout[itreeD0toKpi]->Fill();
244 if(okJPSI) treeout[itreeJPSItoEle]->Fill();
6a213b59 245 delete io2Prong; io2Prong=NULL;
246 }
247
248 twoTrackArray1->Clear();
249 if(!f3Prong && !f4Prong) {
250 negtrack1=0;
251 delete vertexp1n1;
252 continue;
253 }
254
255 // 2nd LOOP ON POSITIVE TRACKS
256 for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) {
257 // get track from tracks array
258 postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2);
b82f6d67 259 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
6a213b59 260 if(dcap2n1>dcaMax) { postrack2=0; continue; }
b82f6d67 261 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
6a213b59 262 if(dcap1p2>dcaMax) { postrack2=0; continue; }
263
b82f6d67 264 // Vertexing
6a213b59 265 twoTrackArray2->AddAt(postrack2,0);
266 twoTrackArray2->AddAt(negtrack1,1);
b82f6d67 267 AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2);
6a213b59 268 if(f3Prong) {
269 threeTrackArray->AddAt(postrack1,0);
270 threeTrackArray->AddAt(negtrack1,1);
271 threeTrackArray->AddAt(postrack2,2);
272 io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
260836a9 273 if(ok3Prong) treeout[itree3Prong]->Fill();
6a213b59 274 if(io3Prong) delete io3Prong; io3Prong=NULL;
275 }
276 if(f4Prong) {
277 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
278 for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
279 // get track from tracks array
280 negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
b82f6d67 281 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
6a213b59 282 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
b82f6d67 283 // Vertexing
6a213b59 284 fourTrackArray->AddAt(postrack1,0);
285 fourTrackArray->AddAt(negtrack1,1);
286 fourTrackArray->AddAt(postrack2,2);
287 fourTrackArray->AddAt(negtrack2,3);
288 io4Prong = Make4Prong(fourTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
260836a9 289 if(ok4Prong) treeout[itree4Prong]->Fill();
6a213b59 290 delete io4Prong; io4Prong=NULL;
291 fourTrackArray->Clear();
292 negtrack2 = 0;
293 } // end loop on negative tracks
294 }
295 postrack2 = 0;
296 delete vertexp2n1;
297 } // end 2nd loop on positive tracks
298 twoTrackArray2->Clear();
299
300 // 2nd LOOP ON NEGATIVE TRACKS
301 for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
302 // get track from tracks array
303 negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
b82f6d67 304 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
6a213b59 305 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
b82f6d67 306 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
6a213b59 307 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
308
b82f6d67 309 // Vertexing
6a213b59 310 twoTrackArray2->AddAt(postrack1,0);
311 twoTrackArray2->AddAt(negtrack2,1);
b82f6d67 312 AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2);
6a213b59 313 if(f3Prong) {
314 threeTrackArray->AddAt(negtrack1,0);
315 threeTrackArray->AddAt(postrack1,1);
316 threeTrackArray->AddAt(negtrack2,2);
317 io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
260836a9 318 if(ok3Prong) treeout[itree3Prong]->Fill();
6a213b59 319 if(io3Prong) delete io3Prong; io3Prong=NULL;
320 }
321 negtrack2 = 0;
322 delete vertexp1n2;
323 } // end 2nd loop on negative tracks
324 twoTrackArray2->Clear();
325
326 negtrack1 = 0;
327 delete vertexp1n1;
328 } // end 1st loop on negative tracks
329
330 postrack1 = 0;
331 } // end 1st loop on positive tracks
332
333
334
6a213b59 335 //printf("delete twoTr 1\n");
336 twoTrackArray1->Delete(); delete twoTrackArray1;
337 //printf("delete twoTr 2\n");
338 twoTrackArray2->Delete(); delete twoTrackArray2;
339 //printf("delete threeTr 1\n");
340 threeTrackArray->Clear();
341 threeTrackArray->Delete(); delete threeTrackArray;
342 //printf("delete fourTr 1\n");
343 fourTrackArray->Delete(); delete fourTrackArray;
344
345
346 // create a copy of this class to be written to output file
347 //AliAnalysisVertexingHF *copy = (AliAnalysisVertexingHF*)this->Clone("AnalysisVertexingHF");
348
349 // print statistics
350 if(fD0toKpi) {
351 printf(" D0->Kpi: event %d = %d; total = %d;\n",
352 (Int_t)esd->GetEventNumberInFile(),
260836a9 353 (Int_t)treeout[itreeD0toKpi]->GetEntries()-initEntriesD0toKpi,
354 (Int_t)treeout[itreeD0toKpi]->GetEntries());
6a213b59 355 }
356 if(fJPSItoEle) {
357 printf(" JPSI->ee: event %d = %d; total = %d;\n",
358 (Int_t)esd->GetEventNumberInFile(),
260836a9 359 (Int_t)treeout[itreeJPSItoEle]->GetEntries()-initEntriesJPSItoEle,
360 (Int_t)treeout[itreeJPSItoEle]->GetEntries());
6a213b59 361 }
362 if(f3Prong) {
363 printf(" Charm->3Prong: event %d = %d; total = %d;\n",
364 (Int_t)esd->GetEventNumberInFile(),
260836a9 365 (Int_t)treeout[itree3Prong]->GetEntries()-initEntries3Prong,
366 (Int_t)treeout[itree3Prong]->GetEntries());
6a213b59 367 }
368 if(f4Prong) {
369 printf(" Charm->4Prong: event %d = %d; total = %d;\n",
370 (Int_t)esd->GetEventNumberInFile(),
260836a9 371 (Int_t)treeout[itree4Prong]->GetEntries()-initEntries4Prong,
372 (Int_t)treeout[itree4Prong]->GetEntries());
6a213b59 373 }
374
375
376 return;
377}
378//----------------------------------------------------------------------------
379AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
5e6f195c 380 TObjArray *twoTrackArray1,AliESDEvent *esd,
6a213b59 381 AliESDVertex *secVertexESD,Double_t dca,
382 Bool_t &okD0,Bool_t &okJPSI) const
383{
384 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
385 // reconstruction cuts
386 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
387
388 okD0=kFALSE; okJPSI=kFALSE;
389
390 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
6a213b59 391
392 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0);
393 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1);
394
395 // propagate tracks to secondary vertex, to compute inv. mass
b82f6d67 396 postrack->RelateToVertex(secVertexESD,fBzkG,10.);
397 negtrack->RelateToVertex(secVertexESD,fBzkG,10.);
6a213b59 398
399 Double_t momentum[3];
400 postrack->GetPxPyPz(momentum);
401 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
402 negtrack->GetPxPyPz(momentum);
403 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
404
405
406 // invariant mass cut (try to improve coding here..)
407 Bool_t okMassCut=kFALSE;
408 if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
409 if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
410 if(!okMassCut) {
411 if(fDebug) printf(" candidate didn't pass mass cut\n");
412 return 0x0;
413 }
414
415
416 AliESDVertex *primVertex = fV1;
417 AliESDVertex *ownPrimVertex=0;
418
419 // primary vertex from *other* tracks in the event
420 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
421 ownPrimVertex = OwnPrimaryVertex(2,twoTrackArray1,esd);
422 if(!ownPrimVertex) {
423 return 0x0;
424 } else {
425 if(ownPrimVertex->GetNContributors()<2) {
426 delete ownPrimVertex;
427 return 0x0;
428 } else {
429 primVertex = ownPrimVertex;
430 }
431 }
432 }
433
434 Float_t d0z0[2],covd0z0[3];
b82f6d67 435 postrack->RelateToVertex(primVertex,fBzkG,10.);
6a213b59 436 postrack->GetImpactParameters(d0z0,covd0z0);
437 d0[0] = d0z0[0];
438 d0err[0] = TMath::Sqrt(covd0z0[0]);
b82f6d67 439 negtrack->RelateToVertex(primVertex,fBzkG,10.);
6a213b59 440 negtrack->GetImpactParameters(d0z0,covd0z0);
441 d0[1] = d0z0[0];
442 d0err[1] = TMath::Sqrt(covd0z0[0]);
443
444 // create the object AliAODRecoDecayHF2Prong
445 Double_t pos[3],cov[6];
446 secVertexESD->GetXYZ(pos); // position
447 secVertexESD->GetCovMatrix(cov); //covariance matrix
448 AliAODVertex *secVertexAOD = new AliAODVertex(pos,cov,secVertexESD->GetChi2toNDF());
449 primVertex->GetXYZ(pos); // position
450 primVertex->GetCovMatrix(cov); //covariance matrix
451 AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
452 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVertexAOD,px,py,pz,d0,d0err,dca);
b82f6d67 453 the2Prong->SetOwnSecondaryVtx(secVertexAOD);//temporary
6a213b59 454 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
455
456 // select D0->Kpi
457 Int_t checkD0,checkD0bar;
458 if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
459 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
460 // select J/psi from B
461 Int_t checkJPSI;
462 if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
463 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
464
465
466 if(okD0 || okJPSI) {
467 // get PID info from ESD
468 Double_t esdpid0[5];
469 postrack->GetESDpid(esdpid0);
470 Double_t esdpid1[5];
471 negtrack->GetESDpid(esdpid1);
472 Double_t esdpid[10];
473 for(Int_t i=0;i<5;i++) {
474 esdpid[i] = esdpid0[i];
475 esdpid[5+i] = esdpid1[i];
476 }
477 the2Prong->SetPID(2,esdpid);
478 }
479
480 if(ownPrimVertex) delete ownPrimVertex;
b82f6d67 481
6a213b59 482 return the2Prong;
483}
484//----------------------------------------------------------------------------
485AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
5e6f195c 486 TObjArray *threeTrackArray,AliESDEvent *esd,
6a213b59 487 AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
488 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
489 Bool_t &ok3Prong) const
490{
491 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
492 // reconstruction cuts
493 // E.Bruna, F.Prino
494
495 ok3Prong=kFALSE;
496 Double_t px[3],py[3],pz[3],d0[3],d0err[3];//d0z[3];
497 Float_t d0z0[2],covd0z0[3];
498
6a213b59 499
500 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
501 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
502 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
503
504 AliESDVertex *primVertex = fV1;
505
b82f6d67 506 postrack1->RelateToVertex(primVertex,fBzkG,10.);
507 negtrack->RelateToVertex(primVertex,fBzkG,10.);
508 postrack2->RelateToVertex(primVertex,fBzkG,10.);
6a213b59 509
510 Double_t momentum[3];
511 postrack1->GetPxPyPz(momentum);
512 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
513 negtrack->GetPxPyPz(momentum);
514 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
515 postrack2->GetPxPyPz(momentum);
516 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
517
518 postrack1->GetImpactParameters(d0z0,covd0z0);
519 d0[0]=d0z0[0];
520 d0err[0] = TMath::Sqrt(covd0z0[0]);
521 negtrack->GetImpactParameters(d0z0,covd0z0);
522 d0[1]=d0z0[0];
523 d0err[1] = TMath::Sqrt(covd0z0[0]);
524 postrack2->GetImpactParameters(d0z0,covd0z0);
525 d0[2]=d0z0[0];
526 d0err[2] = TMath::Sqrt(covd0z0[0]);
527
528
b82f6d67 529 // invariant mass cut for D+, Ds, Lc
6a213b59 530 Bool_t okMassCut=kFALSE;
531 if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
532 if(!okMassCut) {
533 if(fDebug) printf(" candidate didn't pass mass cut\n");
534 return 0x0;
535 }
536
537 //charge
538 Short_t charge=(Short_t)(postrack1->GetSign()*postrack2->GetSign()*negtrack->GetSign());
539
540 AliESDVertex *ownPrimVertex = 0;
541 // primary vertex from *other* tracks in the event
542 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
543 ownPrimVertex = OwnPrimaryVertex(3,threeTrackArray,esd);
544 if(!ownPrimVertex) {
545 return 0x0;
546 } else {
547 if(ownPrimVertex->GetNContributors()<2) {
548 delete ownPrimVertex;
549 return 0x0;
550 } else {
551 primVertex = ownPrimVertex;
552 }
553 }
554 }
555
556 // create the object AliAODRecoDecayHF3Prong
b82f6d67 557 AliESDVertex* secVert3Prong = ReconstructSecondaryVertex(threeTrackArray);
6a213b59 558 Double_t pos[3],cov[6],sigmavert;
559 secVert3Prong->GetXYZ(pos); // position
560 secVert3Prong->GetCovMatrix(cov); //covariance matrix
561 sigmavert=secVert3Prong->GetDispersion();
562
563 AliAODVertex *secVert3PrAOD = new AliAODVertex(pos,cov,secVert3Prong->GetChi2toNDF());
564 primVertex->GetXYZ(pos); // position
565 primVertex->GetCovMatrix(cov); //covariance matrix
566 AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
567 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
568
569 Double_t dist12=TMath::Sqrt((vertexp1n1->GetXv()-pos[0])*(vertexp1n1->GetXv()-pos[0])+(vertexp1n1->GetYv()-pos[1])*(vertexp1n1->GetYv()-pos[1])+(vertexp1n1->GetZv()-pos[2])*(vertexp1n1->GetZv()-pos[2]));
570 Double_t dist23=TMath::Sqrt((vertexp2n1->GetXv()-pos[0])*(vertexp2n1->GetXv()-pos[0])+(vertexp2n1->GetYv()-pos[1])*(vertexp2n1->GetYv()-pos[1])+(vertexp2n1->GetZv()-pos[2])*(vertexp2n1->GetZv()-pos[2]));
571
572 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert3PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
573 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
b82f6d67 574 the3Prong->SetOwnSecondaryVtx(secVert3PrAOD);//temporary
6a213b59 575
576
b82f6d67 577 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
578 if(f3Prong) {
579 ok3Prong = kFALSE;
580 Int_t ok1,ok2;
581 if(the3Prong->SelectDplus(fDplusCuts)) ok3Prong = kTRUE;
582 if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE;
583 if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
584 }
6a213b59 585 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
586 if(ok3Prong) {
587 // get PID info from ESD
588 Double_t esdpid0[5];
589 postrack1->GetESDpid(esdpid0);
590 Double_t esdpid1[5];
591 negtrack->GetESDpid(esdpid1);
592 Double_t esdpid2[5];
593 postrack2->GetESDpid(esdpid2);
594
595
596 Double_t esdpid[15];
597 for(Int_t i=0;i<5;i++) {
598 esdpid[i] = esdpid0[i];
599 esdpid[5+i] = esdpid1[i];
600 esdpid[10+i] = esdpid2[i];
601 }
602 the3Prong->SetPID(3,esdpid);
603 }
604
605 if(ownPrimVertex) delete ownPrimVertex;
606
607 return the3Prong;
608}
609//----------------------------------------------------------------------------
610AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
5e6f195c 611 TObjArray *fourTrackArray,AliESDEvent *esd,
6a213b59 612 AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
613 Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
614 Bool_t &ok4Prong) const
615{
616 // Make 4Prong candidates and check if they pass D0toKpipipi
617 // reconstruction cuts
618 // G.E.Bruno, R.Romita
619
620 ok4Prong=kFALSE;
621
622 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
623 //Float_t d0z0[2],covd0z0[3];
624
b82f6d67 625 px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
6a213b59 626
627 //charge
628 Short_t charge=0;
629
630 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
631 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
632 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
633 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
634
635 AliESDVertex *primVertex = fV1;
636
b82f6d67 637 postrack1->RelateToVertex(primVertex,fBzkG,10.);
638 negtrack1->RelateToVertex(primVertex,fBzkG,10.);
639 postrack2->RelateToVertex(primVertex,fBzkG,10.);
640 negtrack2->RelateToVertex(primVertex,fBzkG,10.);
6a213b59 641
642 Double_t momentum[3];
643 postrack1->GetPxPyPz(momentum);
644 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
645 negtrack1->GetPxPyPz(momentum);
646 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
647 postrack2->GetPxPyPz(momentum);
648 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
649 negtrack2->GetPxPyPz(momentum);
650 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
651
652 // invariant mass cut for rho or D0 (try to improve coding here..)
653 //Bool_t okMassCut=kFALSE;
654 //if(!okMassCut) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
655 //if(!okMassCut) {
656 // if(fDebug) printf(" candidate didn't pass mass cut\n");
657 // return 0x0;
658 //}
659
660 AliESDVertex *ownPrimVertex = 0;
661 // primary vertex from *other* tracks in the event
662 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
663 ownPrimVertex = OwnPrimaryVertex(4,fourTrackArray,esd);
664 if(!ownPrimVertex) {
665 return 0x0;
666 } else {
667 if(ownPrimVertex->GetNContributors()<2) {
668 delete ownPrimVertex;
669 return 0x0;
670 } else {
671 primVertex = ownPrimVertex;
672 }
673 }
674 }
675
676 // create the object AliAODRecoDecayHF4Prong
b82f6d67 677 AliESDVertex* secVert4Prong = ReconstructSecondaryVertex(fourTrackArray);
6a213b59 678 Double_t pos[3],cov[6],sigmavert;
679 secVert4Prong->GetXYZ(pos); // position
680 secVert4Prong->GetCovMatrix(cov); //covariance matrix
681 sigmavert=secVert4Prong->GetDispersion();
682
683 AliAODVertex *secVert4PrAOD = new AliAODVertex(pos,cov,secVert4Prong->GetChi2toNDF());
684 primVertex->GetXYZ(pos); // position
685 primVertex->GetCovMatrix(cov); //covariance matrix
686 AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
687 //Double_t dca[6]={dcap1n1,dcap2n1,dcap1p2,0.,0.,0.}; //
688 Double_t dca[6]={0.,0.,0.,0.,0.,0.}; // modify it
689
690 Double_t dist12=TMath::Sqrt((vertexp1n1->GetXv()-pos[0])*(vertexp1n1->GetXv()-pos[0])+(vertexp1n1->GetYv()-pos[1])*(vertexp1n1->GetYv()-pos[1])+(vertexp1n1->GetZv()-pos[2])*(vertexp1n1->GetZv()-pos[2]));
691 Double_t dist23=TMath::Sqrt((vertexp2n1->GetXv()-pos[0])*(vertexp2n1->GetXv()-pos[0])+(vertexp2n1->GetYv()-pos[1])*(vertexp2n1->GetYv()-pos[1])+(vertexp2n1->GetZv()-pos[2])*(vertexp2n1->GetZv()-pos[2]));
692 Double_t dist14=0.; // to be implemented
693 Double_t dist34=0.; // to be implemented
694
695 //AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
696 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge);
697 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
b82f6d67 698 the4Prong->SetOwnSecondaryVtx(secVert4PrAOD);//temporary
6a213b59 699
700
701 // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available
702 // select D0->Kpipipi
703 //Int_t checkD0,checkD0bar;
704 // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar);
705 ok4Prong=kFALSE; //for the time being ...
706
707
708 // get PID info from ESD
709 Double_t esdpid0[5];
710 postrack1->GetESDpid(esdpid0);
711 Double_t esdpid1[5];
712 negtrack1->GetESDpid(esdpid1);
713 Double_t esdpid2[5];
714 postrack2->GetESDpid(esdpid2);
715 Double_t esdpid3[5];
716 negtrack2->GetESDpid(esdpid3);
717
718 Double_t esdpid[20];
719 for(Int_t i=0;i<5;i++) {
720 esdpid[i] = esdpid0[i];
721 esdpid[5+i] = esdpid1[i];
722 esdpid[10+i] = esdpid2[i];
723 esdpid[15+i] = esdpid3[i];
724 }
725 the4Prong->SetPID(4,esdpid);
726
727 if(ownPrimVertex) delete ownPrimVertex;
728
729 return the4Prong;
730}
731//-----------------------------------------------------------------------------
732AliESDVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(Int_t ntrks,
733 TObjArray *trkArray,
5e6f195c 734 AliESDEvent *esd) const
6a213b59 735{
736 // Returns primary vertex specific to this candidate
737
738 AliVertexerTracks *vertexer1 = new AliVertexerTracks(esd->GetMagneticField());
739 AliESDVertex *ownPrimVertex = 0;
740
741 // recalculating the vertex
742 if(fRecoPrimVtxSkippingTrks) {
743 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
744 Float_t diamondcovxy[3];
745 esd->GetDiamondCovXY(diamondcovxy);
746 Double_t pos[3]={esd->GetDiamondX(),esd->GetDiamondY(),0.};
747 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.};
748 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
749 vertexer1->SetVtxStart(diamond);
750 delete diamond; diamond=NULL;
751 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
752 vertexer1->SetOnlyFitter();
753 }
d70e06c7 754 Int_t skipped[10];
6a213b59 755 AliESDtrack *t = 0;
756 for(Int_t i=0; i<ntrks; i++) {
757 t = (AliESDtrack*)trkArray->UncheckedAt(i);
d70e06c7 758 skipped[i] = (Int_t)t->GetID();
6a213b59 759 }
760 vertexer1->SetSkipTracks(ntrks,skipped);
761 ownPrimVertex = (AliESDVertex*)vertexer1->FindPrimaryVertex(esd);
762 }
763
764 // removing the prongs tracks
765 if(fRmTrksFromPrimVtx) {
f65f93f3 766 TObjArray rmArray(ntrks);
767 UShort_t *rmId = new UShort_t[ntrks];
6a213b59 768 AliESDtrack *esdTrack = 0;
6a213b59 769 AliESDtrack *t = 0;
770 for(Int_t i=0; i<ntrks; i++) {
771 t = (AliESDtrack*)trkArray->UncheckedAt(i);
772 esdTrack = new AliESDtrack(*t);
f65f93f3 773 rmArray.AddLast(esdTrack);
774 rmId[i]=(UShort_t)esdTrack->GetID();
6a213b59 775 delete esdTrack;
776 }
777 Float_t diamondxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
f65f93f3 778 ownPrimVertex = vertexer1->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
779 delete [] rmId; rmId=NULL;
780 rmArray.Delete();
6a213b59 781 }
782
783 delete vertexer1; vertexer1=NULL;
784
785 return ownPrimVertex;
786}
787//-----------------------------------------------------------------------------
788void AliAnalysisVertexingHF::PrintStatus() const {
789 // Print parameters being used
790
791 printf("Preselections:\n");
792 printf(" fITSrefit = %d\n",(Int_t)fITSrefit);
793 printf(" fBothSPD = %d\n",(Int_t)fBothSPD);
794 printf(" fMinITSCls = %d\n",fMinITSCls);
795 printf(" fMinPtCut = %f GeV/c\n",fMinPtCut);
796 printf(" fMind0rphiCut = %f cm\n",fMind0rphiCut);
4d5c9633 797 if(fSecVtxWithKF) {
798 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
799 } else {
800 printf("Secondary vertex with AliVertexerTracks\n");
801 }
6a213b59 802 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
803 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
804 if(fD0toKpi) {
805 printf("Reconstruct D0->Kpi candidates with cuts:\n");
806 printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
807 printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
808 printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
809 printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
810 printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
811 printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
812 printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
813 printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
814 printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
815 }
816 if(fJPSItoEle) {
817 printf("Reconstruct J/psi from B candidates with cuts:\n");
818 printf(" |M-MJPSI| [GeV] < %f\n",fBtoJPSICuts[0]);
819 printf(" dca [cm] < %f\n",fBtoJPSICuts[1]);
820 printf(" cosThetaStar < %f\n",fBtoJPSICuts[2]);
821 printf(" pTP [GeV/c] > %f\n",fBtoJPSICuts[3]);
822 printf(" pTN [GeV/c] > %f\n",fBtoJPSICuts[4]);
823 printf(" |d0P| [cm] < %f\n",fBtoJPSICuts[5]);
824 printf(" |d0N| [cm] < %f\n",fBtoJPSICuts[6]);
825 printf(" d0d0 [cm^2] < %f\n",fBtoJPSICuts[7]);
826 printf(" cosThetaPoint > %f\n",fBtoJPSICuts[8]);
827 }
828 if(f3Prong) {
829 printf("Reconstruct 3 prong candidates.\n");
4d5c9633 830 printf(" D+->Kpipi cuts:\n");
6a213b59 831 printf(" |M-MD+| [GeV] < %f\n",fDplusCuts[0]);
832 printf(" pTK [GeV/c] > %f\n",fDplusCuts[1]);
833 printf(" pTPi [GeV/c] > %f\n",fDplusCuts[2]);
834 printf(" |d0K| [cm] > %f\n",fDplusCuts[3]);
835 printf(" |d0Pi| [cm] > %f\n",fDplusCuts[4]);
836 printf(" dist12 [cm] < %f\n",fDplusCuts[5]);
837 printf(" sigmavert [cm] < %f\n",fDplusCuts[6]);
838 printf(" dist prim-sec [cm] > %f\n",fDplusCuts[7]);
839 printf(" pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
840 printf(" cosThetaPoint > %f\n",fDplusCuts[9]);
841 printf(" Sum d0^2 [cm^2] > %f\n",fDplusCuts[10]);
842 printf(" dca cut [cm] < %f\n",fDplusCuts[11]);
4d5c9633 843 printf(" Ds->KKpi cuts:\n");
844 printf(" |M-MDs| [GeV] < %f\n",fDsCuts[0]);
845 printf(" Lc->pKpi cuts:\n");
846 printf(" |M-MD+| [GeV] < %f\n",fLcCuts[0]);
6a213b59 847 }
848
849 return;
850}
851//-----------------------------------------------------------------------------
852Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
853 Int_t nprongs,
854 Double_t *px,
855 Double_t *py,
856 Double_t *pz) const {
857 // Check invariant mass cut
858
859 Short_t dummycharge=0;
860 Double_t *dummyd0 = new Double_t[nprongs];
861 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
862 delete [] dummyd0;
863
864 UInt_t pdg2[2],pdg3[3];
865 Double_t mPDG,minv;
866
867 Bool_t retval=kFALSE;
868 switch (decay)
869 {
870 case 0: // D0->Kpi
871 pdg2[0]=211; pdg2[1]=321;
872 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
873 minv = rd->InvMass(nprongs,pdg2);
874 if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
875 pdg2[0]=321; pdg2[1]=211;
876 minv = rd->InvMass(nprongs,pdg2);
877 if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
878 break;
879 case 1: // JPSI->ee
880 pdg2[0]=11; pdg2[1]=11;
881 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
882 minv = rd->InvMass(nprongs,pdg2);
883 if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
884 break;
885 case 2: // D+->Kpipi
886 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
887 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
888 minv = rd->InvMass(nprongs,pdg3);
889 if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
b82f6d67 890 // Ds+->KKpi
891 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
892 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
893 minv = rd->InvMass(nprongs,pdg3);
894 if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
895 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
896 minv = rd->InvMass(nprongs,pdg3);
897 if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
898 // Lc->pKpi
899 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
900 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
901 minv = rd->InvMass(nprongs,pdg3);
902 if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
903 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
904 minv = rd->InvMass(nprongs,pdg3);
905 if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
6a213b59 906 break;
907 default:
908 break;
909 }
910
911 delete rd;
912
913 return retval;
914}
915//-----------------------------------------------------------------------------
5e6f195c 916void AliAnalysisVertexingHF::SelectTracks(AliESDEvent *esd,
6a213b59 917 TObjArray &trksP,Int_t &nTrksP,
918 TObjArray &trksN,Int_t &nTrksN) const
919{
920 // Fill two TObjArrays with positive and negative tracks and
921 // apply single-track preselection
922
923 nTrksP=0,nTrksN=0;
924
925 Int_t entries = (Int_t)esd->GetNumberOfTracks();
926
927 // transfer ITS tracks from ESD to arrays and to a tree
928 for(Int_t i=0; i<entries; i++) {
929
930 AliESDtrack *esdtrack = esd->GetTrack(i);
931 UInt_t status = esdtrack->GetStatus();
932
933 // require refit in ITS
934 if(fITSrefit && !(status&AliESDtrack::kITSrefit)) {
935 if(fDebug) printf("track %d is not kITSrefit\n",i);
936 continue;
937 }
938
939 // require minimum # of ITS points
940 if(esdtrack->GetNcls(0)<fMinITSCls) {
941 if(fDebug) printf("track %d has %d ITS cls\n",i,esdtrack->GetNcls(0));
942 continue;
943 }
944 // require points on the 2 pixel layers
945 if(fBothSPD)
946 if(!TESTBIT(esdtrack->GetITSClusterMap(),0) ||
947 !TESTBIT(esdtrack->GetITSClusterMap(),1)) continue;
948
949 // single track selection
b82f6d67 950 if(!SingleTrkCuts(*esdtrack)) continue;
6a213b59 951
952 if(esdtrack->GetSign()<0) { // negative track
953 trksN.AddLast(esdtrack);
954 nTrksN++;
955 } else { // positive track
956 trksP.AddLast(esdtrack);
957 nTrksP++;
958 }
959
960 } // loop on ESD tracks
961
962 return;
963}
964//-----------------------------------------------------------------------------
965void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
966 Double_t cut2,Double_t cut3,Double_t cut4,
967 Double_t cut5,Double_t cut6,
968 Double_t cut7,Double_t cut8)
969{
970 // Set the cuts for D0 selection
971 fD0toKpiCuts[0] = cut0;
972 fD0toKpiCuts[1] = cut1;
973 fD0toKpiCuts[2] = cut2;
974 fD0toKpiCuts[3] = cut3;
975 fD0toKpiCuts[4] = cut4;
976 fD0toKpiCuts[5] = cut5;
977 fD0toKpiCuts[6] = cut6;
978 fD0toKpiCuts[7] = cut7;
979 fD0toKpiCuts[8] = cut8;
980
981 return;
982}
983//-----------------------------------------------------------------------------
984void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
985{
986 // Set the cuts for D0 selection
987
988 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
989
990 return;
991}
992//-----------------------------------------------------------------------------
993void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
994 Double_t cut2,Double_t cut3,Double_t cut4,
995 Double_t cut5,Double_t cut6,
996 Double_t cut7,Double_t cut8)
997{
998 // Set the cuts for J/psi from B selection
999 fBtoJPSICuts[0] = cut0;
1000 fBtoJPSICuts[1] = cut1;
1001 fBtoJPSICuts[2] = cut2;
1002 fBtoJPSICuts[3] = cut3;
1003 fBtoJPSICuts[4] = cut4;
1004 fBtoJPSICuts[5] = cut5;
1005 fBtoJPSICuts[6] = cut6;
1006 fBtoJPSICuts[7] = cut7;
1007 fBtoJPSICuts[8] = cut8;
1008
1009 return;
1010}
1011//-----------------------------------------------------------------------------
1012void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1013{
1014 // Set the cuts for J/psi from B selection
1015
1016 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1017
1018 return;
1019}
1020//-----------------------------------------------------------------------------
1021void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1022 Double_t cut2,Double_t cut3,Double_t cut4,
1023 Double_t cut5,Double_t cut6,
1024 Double_t cut7,Double_t cut8,
1025 Double_t cut9,Double_t cut10,Double_t cut11)
1026{
1027 // Set the cuts for Dplus->Kpipi selection
1028 fDplusCuts[0] = cut0;
1029 fDplusCuts[1] = cut1;
1030 fDplusCuts[2] = cut2;
1031 fDplusCuts[3] = cut3;
1032 fDplusCuts[4] = cut4;
1033 fDplusCuts[5] = cut5;
1034 fDplusCuts[6] = cut6;
1035 fDplusCuts[7] = cut7;
1036 fDplusCuts[8] = cut8;
1037 fDplusCuts[9] = cut9;
1038 fDplusCuts[10] = cut10;
1039 fDplusCuts[11] = cut11;
1040
1041 return;
1042}
1043//-----------------------------------------------------------------------------
1044void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1045{
b82f6d67 1046 // Set the cuts for Dplus->Kpipi selection
6a213b59 1047
1048 for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1049
1050 return;
1051}
1052//-----------------------------------------------------------------------------
b82f6d67 1053void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0)
1054{
1055 // Set the cuts for Ds->KKpi selection
1056 fDsCuts[0] = cut0;
1057
1058 return;
1059}
1060//-----------------------------------------------------------------------------
1061void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[1])
1062{
1063 // Set the cuts for Ds->KKpi selection
1064
1065 for(Int_t i=0; i<1; i++) fDsCuts[i] = cuts[i];
1066
1067 return;
1068}
1069//-----------------------------------------------------------------------------
1070void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0)
1071{
1072 // Set the cuts for Lc->pKpi selection
1073 fLcCuts[0] = cut0;
1074
1075 return;
1076}
1077//-----------------------------------------------------------------------------
1078void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[1])
1079{
1080 // Set the cuts for Lc->pKpi selection
1081
1082 for(Int_t i=0; i<1; i++) fLcCuts[i] = cuts[i];
1083
1084 return;
1085}
1086//-----------------------------------------------------------------------------
1087Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk) const
6a213b59 1088{
1089 // Check if track passes some kinematical cuts
6a213b59 1090
1091 if(TMath::Abs(1./trk.GetParameter()[4]) < fMinPtCut) {
1092 printf("pt %f\n",1./trk.GetParameter()[4]);
1093 return kFALSE;
1094 }
b82f6d67 1095 trk.RelateToVertex(fV1,fBzkG,10.);
6a213b59 1096 Float_t d0z0[2],covd0z0[3];
1097 trk.GetImpactParameters(d0z0,covd0z0);
1098 if(TMath::Abs(d0z0[0]) < fMind0rphiCut) {
1099 printf("d0rphi %f\n",TMath::Abs(d0z0[0]));
1100 return kFALSE;
1101 }
1102
1103 return kTRUE;
1104}
1105//-----------------------------------------------------------------------------
b82f6d67 1106AliESDVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray) const
1107{
1108 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1109
1110 AliESDVertex *vertex = 0;
1111
1112 if(!fSecVtxWithKF) { // AliVertexerTracks
1113
1114 AliVertexerTracks *vertexer2 = new AliVertexerTracks(fBzkG);
1115 vertexer2->SetDebug(0);
1116 vertexer2->SetVtxStart(fV1);
f65f93f3 1117 vertex = (AliESDVertex*)vertexer2->VertexForSelectedESDTracks(trkArray);
b82f6d67 1118 delete vertexer2;
1119
1120 } else { // Kalman Filter vertexer (AliKFParticle)
1121
1122 AliKFParticle::SetField(fBzkG);
1123
1124 AliKFParticle vertexKF;
1125
1126 Int_t nTrks = trkArray->GetEntriesFast();
1127 for(Int_t i=0; i<nTrks; i++) {
1128 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1129 AliKFParticle daughterKF(*esdTrack,211);
1130 vertexKF.AddDaughter(daughterKF);
1131 }
1132 vertex = new AliESDVertex();
1133 vertexKF.CopyToESDVertex(*vertex);
1134
1135 }
1136
1137 return vertex;
1138}
1139//-----------------------------------------------------------------------------
6a213b59 1140
1141
1142