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