1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
23 // Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita
24 //----------------------------------------------------------------------------
26 #include <TDatabasePDG.h>
29 #include "AliVertexerTracks.h"
30 #include "AliESDVertex.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"
39 ClassImp(AliAnalysisVertexingHF)
41 //----------------------------------------------------------------------------
42 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
43 fRecoPrimVtxSkippingTrks(kFALSE),
44 fRmTrksFromPrimVtx(kFALSE),
47 fD0toKpi(kTRUE),fJPSItoEle(kTRUE),f3Prong(kTRUE),f4Prong(kTRUE),
48 fITSrefit(kFALSE),fBothSPD(kTRUE),fMinITSCls(5),
49 fMinPtCut(0.),fMind0rphiCut(0.)
51 // Default constructor
57 //----------------------------------------------------------------------------
58 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
62 //----------------------------------------------------------------------------
63 void AliAnalysisVertexingHF::FindCandidates(AliESD *esd,TTree treeout[])
65 // Find heavy-flavour vertex candidates
68 AliAODRecoDecayHF2Prong *io2Prong = new AliAODRecoDecayHF2Prong;
69 AliAODRecoDecayHF3Prong *io3Prong = new AliAODRecoDecayHF3Prong;
70 AliAODRecoDecayHF4Prong *io4Prong = new AliAODRecoDecayHF4Prong;
72 Int_t itreeD0toKpi=-1,itreeJPSItoEle=-1,itree3Prong=-1,itree4Prong=-1;
73 Int_t initEntriesD0toKpi=0,initEntriesJPSItoEle=0,initEntries3Prong=0,initEntries4Prong=0;
76 //treeout[itree].Print();
77 treeout[itree].SetBranchAddress("D0toKpi",&io2Prong);
79 initEntriesD0toKpi = treeout[itreeD0toKpi].GetEntries();
83 //treeout[itree].Print();
84 treeout[itree].SetBranchAddress("JPSItoEle",&io2Prong);
86 initEntriesJPSItoEle = treeout[itreeJPSItoEle].GetEntries();
90 treeout[itree].SetBranchAddress("Charmto3Prong",&io3Prong);
92 initEntries3Prong = treeout[itree3Prong].GetEntries();
96 treeout[itree].SetBranchAddress("D0to4Prong",&io4Prong);
98 initEntries4Prong = treeout[itree4Prong].GetEntries();
100 delete io2Prong; io2Prong = NULL;
101 delete io3Prong; io3Prong = NULL;
102 delete io4Prong; io4Prong = NULL;
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);
117 // create the AliVertexerTracks object for the secondary vertex
118 AliVertexerTracks *vertexer2 = new AliVertexerTracks(esd->GetMagneticField());
119 vertexer2->SetDebug(0);
121 Int_t ev = (Int_t)esd->GetEventNumberInFile();
122 printf("--- Finding candidates in event %d\n",ev);
124 Double_t bfieldkG=(Double_t)esd->GetMagneticField();
126 trkEntries = (Int_t)esd->GetNumberOfTracks();
127 printf(" Number of tracks: %d\n",trkEntries);
128 if(trkEntries<2) return;
130 // retrieve primary vertex from the AliESD
131 if(!esd->GetPrimaryVertex()) {
132 printf(" No vertex in AliESD\n");
135 AliESDVertex copy(*(esd->GetPrimaryVertex()));
136 SetPrimaryVertex(©);
137 vertexer2->SetVtxStart(©);
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,
146 printf(" Pos. tracks: %d Neg. tracks: %d\n",nTrksP,nTrksN);
148 TObjArray *twoTrackArray1 = new TObjArray(2);
149 TObjArray *twoTrackArray2 = new TObjArray(2);
150 TObjArray *threeTrackArray = new TObjArray(3);
151 TObjArray *fourTrackArray = new TObjArray(4);
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);
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;
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;
182 twoTrackArray1->Clear();
183 if(!f3Prong && !f4Prong) {
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; }
198 // Vertexing with AliVertexerTracks
199 twoTrackArray2->AddAt(postrack2,0);
200 twoTrackArray2->AddAt(negtrack1,1);
201 AliESDVertex *vertexp2n1 = vertexer2->VertexForSelectedTracks(twoTrackArray2);
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;
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();
227 } // end loop on negative tracks
231 } // end 2nd loop on positive tracks
232 twoTrackArray2->Clear();
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; }
243 // Vertexing with AliVertexerTracks
244 twoTrackArray2->AddAt(postrack1,0);
245 twoTrackArray2->AddAt(negtrack2,1);
246 AliESDVertex *vertexp1n2 = vertexer2->VertexForSelectedTracks(twoTrackArray2);
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;
257 } // end 2nd loop on negative tracks
258 twoTrackArray2->Clear();
262 } // end 1st loop on negative tracks
265 } // end 1st loop on positive tracks
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;
281 // create a copy of this class to be written to output file
282 //AliAnalysisVertexingHF *copy = (AliAnalysisVertexingHF*)this->Clone("AnalysisVertexingHF");
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());
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());
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());
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());
313 //----------------------------------------------------------------------------
314 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
315 TObjArray *twoTrackArray1,AliESD *esd,
316 AliESDVertex *secVertexESD,Double_t dca,
317 Bool_t &okD0,Bool_t &okJPSI) const
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)
323 okD0=kFALSE; okJPSI=kFALSE;
325 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
326 Double_t bfieldkG=(Double_t)esd->GetMagneticField();
328 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0);
329 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1);
331 // propagate tracks to secondary vertex, to compute inv. mass
332 postrack->RelateToVertex(secVertexESD,bfieldkG,10.);
333 negtrack->RelateToVertex(secVertexESD,bfieldkG,10.);
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];
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;
347 if(fDebug) printf(" candidate didn't pass mass cut\n");
352 AliESDVertex *primVertex = fV1;
353 AliESDVertex *ownPrimVertex=0;
355 // primary vertex from *other* tracks in the event
356 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
357 ownPrimVertex = OwnPrimaryVertex(2,twoTrackArray1,esd);
361 if(ownPrimVertex->GetNContributors()<2) {
362 delete ownPrimVertex;
365 primVertex = ownPrimVertex;
370 Float_t d0z0[2],covd0z0[3];
371 postrack->RelateToVertex(primVertex,bfieldkG,10.);
372 postrack->GetImpactParameters(d0z0,covd0z0);
374 d0err[0] = TMath::Sqrt(covd0z0[0]);
375 negtrack->RelateToVertex(primVertex,bfieldkG,10.);
376 negtrack->GetImpactParameters(d0z0,covd0z0);
378 d0err[1] = TMath::Sqrt(covd0z0[0]);
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);
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
397 if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
398 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
402 // get PID info from ESD
404 postrack->GetESDpid(esdpid0);
406 negtrack->GetESDpid(esdpid1);
408 for(Int_t i=0;i<5;i++) {
409 esdpid[i] = esdpid0[i];
410 esdpid[5+i] = esdpid1[i];
412 the2Prong->SetPID(2,esdpid);
415 if(ownPrimVertex) delete ownPrimVertex;
419 //----------------------------------------------------------------------------
420 AliAODRecoDecayHF3Prong* 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
426 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
427 // reconstruction cuts
431 Double_t px[3],py[3],pz[3],d0[3],d0err[3];//d0z[3];
432 Float_t d0z0[2],covd0z0[3];
434 Double_t bfieldkG=(Double_t)esd->GetMagneticField();
436 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
437 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
438 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
440 AliESDVertex *primVertex = fV1;
442 postrack1->RelateToVertex(primVertex,bfieldkG,10.);
443 negtrack->RelateToVertex(primVertex,bfieldkG,10.);
444 postrack2->RelateToVertex(primVertex,bfieldkG,10.);
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];
454 postrack1->GetImpactParameters(d0z0,covd0z0);
456 d0err[0] = TMath::Sqrt(covd0z0[0]);
457 negtrack->GetImpactParameters(d0z0,covd0z0);
459 d0err[1] = TMath::Sqrt(covd0z0[0]);
460 postrack2->GetImpactParameters(d0z0,covd0z0);
462 d0err[2] = TMath::Sqrt(covd0z0[0]);
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;
469 if(fDebug) printf(" candidate didn't pass mass cut\n");
474 Short_t charge=(Short_t)(postrack1->GetSign()*postrack2->GetSign()*negtrack->GetSign());
476 AliESDVertex *ownPrimVertex = 0;
477 // primary vertex from *other* tracks in the event
478 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
479 ownPrimVertex = OwnPrimaryVertex(3,threeTrackArray,esd);
483 if(ownPrimVertex->GetNContributors()<2) {
484 delete ownPrimVertex;
487 primVertex = ownPrimVertex;
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();
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};
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]));
510 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert3PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
511 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
515 if(f3Prong) ok3Prong = the3Prong->SelectDplus(fDplusCuts);
516 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
518 // get PID info from ESD
520 postrack1->GetESDpid(esdpid0);
522 negtrack->GetESDpid(esdpid1);
524 postrack2->GetESDpid(esdpid2);
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];
533 the3Prong->SetPID(3,esdpid);
536 if(ownPrimVertex) delete ownPrimVertex;
540 //----------------------------------------------------------------------------
541 AliAODRecoDecayHF4Prong* 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
547 // Make 4Prong candidates and check if they pass D0toKpipipi
548 // reconstruction cuts
549 // G.E.Bruno, R.Romita
553 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
554 //Float_t d0z0[2],covd0z0[3];
556 Double_t bfieldkG=(Double_t)esd->GetMagneticField();
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);
566 AliESDVertex *primVertex = fV1;
568 postrack1->RelateToVertex(primVertex,bfieldkG,10.);
569 negtrack1->RelateToVertex(primVertex,bfieldkG,10.);
570 postrack2->RelateToVertex(primVertex,bfieldkG,10.);
571 negtrack2->RelateToVertex(primVertex,bfieldkG,10.);
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];
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;
587 // if(fDebug) printf(" candidate didn't pass mass cut\n");
591 AliESDVertex *ownPrimVertex = 0;
592 // primary vertex from *other* tracks in the event
593 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
594 ownPrimVertex = OwnPrimaryVertex(4,fourTrackArray,esd);
598 if(ownPrimVertex->GetNContributors()<2) {
599 delete ownPrimVertex;
602 primVertex = ownPrimVertex;
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();
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
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
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);
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 ...
640 // get PID info from ESD
642 postrack1->GetESDpid(esdpid0);
644 negtrack1->GetESDpid(esdpid1);
646 postrack2->GetESDpid(esdpid2);
648 negtrack2->GetESDpid(esdpid3);
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];
657 the4Prong->SetPID(4,esdpid);
659 if(ownPrimVertex) delete ownPrimVertex;
663 //-----------------------------------------------------------------------------
664 AliESDVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(Int_t ntrks,
668 // Returns primary vertex specific to this candidate
670 AliVertexerTracks *vertexer1 = new AliVertexerTracks(esd->GetMagneticField());
671 AliESDVertex *ownPrimVertex = 0;
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();
686 Int_t skipped[ntrks];
688 for(Int_t i=0; i<ntrks; i++) {
689 t = (AliESDtrack*)trkArray->UncheckedAt(i);
690 skipped[i] = t->GetID();
692 vertexer1->SetSkipTracks(ntrks,skipped);
693 ownPrimVertex = (AliESDVertex*)vertexer1->FindPrimaryVertex(esd);
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);
702 for(Int_t i=0; i<ntrks; i++) {
703 t = (AliESDtrack*)trkArray->UncheckedAt(i);
704 esdTrack = new AliESDtrack(*t);
708 Float_t diamondxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
709 ownPrimVertex = vertexer1->RemoveTracksFromVertex(fV1,rmTree,diamondxy);
710 delete rmTree; rmTree=NULL;
713 delete vertexer1; vertexer1=NULL;
715 return ownPrimVertex;
717 //-----------------------------------------------------------------------------
718 void AliAnalysisVertexingHF::PrintStatus() const {
719 // Print parameters being used
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");
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]);
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]);
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]);
772 //-----------------------------------------------------------------------------
773 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
777 Double_t *pz) const {
778 // Check invariant mass cut
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);
785 UInt_t pdg2[2],pdg3[3];
788 Bool_t retval=kFALSE;
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;
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;
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;
820 //-----------------------------------------------------------------------------
821 void AliAnalysisVertexingHF::SelectTracks(AliESD *esd,
822 TObjArray &trksP,Int_t &nTrksP,
823 TObjArray &trksN,Int_t &nTrksN) const
825 // Fill two TObjArrays with positive and negative tracks and
826 // apply single-track preselection
830 Int_t entries = (Int_t)esd->GetNumberOfTracks();
832 // transfer ITS tracks from ESD to arrays and to a tree
833 for(Int_t i=0; i<entries; i++) {
835 AliESDtrack *esdtrack = esd->GetTrack(i);
836 UInt_t status = esdtrack->GetStatus();
838 // require refit in ITS
839 if(fITSrefit && !(status&AliESDtrack::kITSrefit)) {
840 if(fDebug) printf("track %d is not kITSrefit\n",i);
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));
849 // require points on the 2 pixel layers
851 if(!TESTBIT(esdtrack->GetITSClusterMap(),0) ||
852 !TESTBIT(esdtrack->GetITSClusterMap(),1)) continue;
854 // single track selection
855 if(!SingleTrkCuts(*esdtrack,esd->GetMagneticField())) continue;
857 if(esdtrack->GetSign()<0) { // negative track
858 trksN.AddLast(esdtrack);
860 } else { // positive track
861 trksP.AddLast(esdtrack);
865 } // loop on ESD tracks
869 //-----------------------------------------------------------------------------
870 void 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)
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;
888 //-----------------------------------------------------------------------------
889 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
891 // Set the cuts for D0 selection
893 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
897 //-----------------------------------------------------------------------------
898 void 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)
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;
916 //-----------------------------------------------------------------------------
917 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
919 // Set the cuts for J/psi from B selection
921 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
925 //-----------------------------------------------------------------------------
926 void 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)
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;
948 //-----------------------------------------------------------------------------
949 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
951 // Set the cuts for Dplus selection
953 for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
957 //-----------------------------------------------------------------------------
958 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk,
959 Double_t bfieldkG) const
961 // Check if track passes some kinematical cuts
962 // Magnetic field "bfieldkG" (kG)
964 if(TMath::Abs(1./trk.GetParameter()[4]) < fMinPtCut) {
965 printf("pt %f\n",1./trk.GetParameter()[4]);
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]));
978 //-----------------------------------------------------------------------------