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