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