memory consumption measures only memory used by preprocessor
[u/mrichter/AliRoot.git] / PWG3 / vertexingOld / AliBtoJPSItoEle.cxx
CommitLineData
27e190fb 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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 BtoJPSItoEle class
18// for pp and PbPb interactions
19// Note: the two decay tracks are labelled: 0 (positive electron)
20// 1 (negative electron)
21// Origin: G.E. Bruno giuseppe.bruno@ba.infn.it
22// based on Class for charm golden channel (D0->Kpi)
23//----------------------------------------------------------------------------
24
25// #include <Riostream.h> // for debugging
26
27#include <TH1.h>
28#include <TH2.h>
29#include <TCanvas.h>
30#include <TPaveLabel.h>
31#include <TVector3.h>
0ac6b6e8 32#include <TString.h>
27e190fb 33
34#include "AliBtoJPSItoEle.h"
35
36ClassImp(AliBtoJPSItoEle)
37
38//----------------------------------------------------------------------------
cbd5890c 39AliBtoJPSItoEle::AliBtoJPSItoEle():
40fSignal(kFALSE),
41fJpsiPrimary(kFALSE),
42fEvent(0),
43fV1x(0.),
44fV1y(0.),
45fV1z(0.),
46fV2x(0.),
47fV2y(0.),
48fV2z(0.),
49fDCA(0.),
50fWgtJPsi(0.)
51{
27e190fb 52 // Default constructor
53
0ac6b6e8 54
27e190fb 55 fTrkNum[0] = 0;
56 fTrkNum[1] = 0;
57
27e190fb 58 fPx[0] = 0.;
59 fPy[0] = 0.;
60 fPz[0] = 0.;
61 fPx[1] = 0.;
62 fPy[1] = 0.;
63 fPz[1] = 0.;
64
65 fd0[0] = 0.;
66 fd0[1] = 0.;
67
68 fPdg[0] = 0;
69 fPdg[1] = 0;
70 fMum[0] = 0;
71 fMum[1] = 0;
72 fGMum[0] = 0;
73 fGMum[1] = 0;
74
75 fTagPi[0] = 0.;
76 fTagPi[1] = 0.;
77 fTagKa[0] = 0.;
78 fTagKa[1] = 0.;
79 fTagNid[0] = 0.;
80 fTagNid[1] = 0.;
81
27e190fb 82}
83//----------------------------------------------------------------------------
84AliBtoJPSItoEle::AliBtoJPSItoEle(Int_t ev,Int_t trkNum[2],
85 Double_t v1[3],Double_t v2[3],
86 Double_t dca,
cbd5890c 87 Double_t mom[6],Double_t d0[2]) :
88fSignal(kFALSE),
89fJpsiPrimary(kFALSE),
90fEvent(ev),
91fV1x(v1[0]),
92fV1y(v1[1]),
93fV1z(v1[2]),
94fV2x(v2[0]),
95fV2y(v2[1]),
96fV2z(v2[2]),
97fDCA(dca),
98fWgtJPsi(0.)
99{
27e190fb 100 // Constructor
101
27e190fb 102 fTrkNum[0] = trkNum[0];
103 fTrkNum[1] = trkNum[1];
104
0ac6b6e8 105 fV1x = v1[0];
106 fV1y = v1[1];
107 fV1z = v1[2];
108 fV2x = v2[0];
109 fV2y = v2[1];
110 fV2z = v2[2];
0ac6b6e8 111
27e190fb 112 fPx[0] = mom[0];
113 fPy[0] = mom[1];
114 fPz[0] = mom[2];
115 fPx[1] = mom[3];
116 fPy[1] = mom[4];
117 fPz[1] = mom[5];
118
119 fd0[0] = d0[0];
120 fd0[1] = d0[1];
121
122 fPdg[0] = 0;
123 fPdg[1] = 0;
124 fMum[0] = 0;
125 fMum[1] = 0;
126 fGMum[0] = 0;
127 fGMum[1] = 0;
128
129 fTagPi[0] = 0.;
130 fTagPi[1] = 0.;
131 fTagKa[0] = 0.;
132 fTagKa[1] = 0.;
133 fTagNid[0] = 0.;
134 fTagNid[1] = 0.;
0ac6b6e8 135
27e190fb 136}
137//----------------------------------------------------------------------------
138AliBtoJPSItoEle::~AliBtoJPSItoEle() {}
27e190fb 139//----------------------------------------------------------------------------
140void AliBtoJPSItoEle::ApplyPID(TString pidScheme) {
141 // Applies particle identification
142
143 if(!pidScheme.CompareTo("TRDTPCparam") && fPdg[0]==0) {
144 printf("AliBtoJPSItoEle::ApplyPID :\n Warning: TRD-TPC parameterized PID can be used only for simulation!\n");
145 return;
146 }
147
148 if(!pidScheme.CompareTo("TRDTPCparam")) {
149 // tagging of the positive track
150 if(TMath::Abs(fPdg[0])==11) { // electron
151 fTagEl[0] = 0.81;
152 fTagNid[0] = 1.-fTagEl[0];
153 }
154 else if(TMath::Abs(fPdg[0])==211) { // pion
155 fTagEl[0] = TRDTPCCombinedPIDParametrization(PChild(0));
156 fTagNid[0] = 1.-fTagEl[0];
157 }
158 else { // all the others
159 fTagEl[0] = 0.;
160 fTagNid[0] = 1.;
161 }
162 // tagging of the negative track
163 if(TMath::Abs(fPdg[1])==11) { // electron
164 fTagEl[1] = 0.81;
165 fTagNid[1] = 1.-fTagEl[1];
166 }
167 else if(TMath::Abs(fPdg[1])==211) { // pion
168 fTagEl[1] = TRDTPCCombinedPIDParametrization(PChild(1));
169 fTagNid[1] = 1.-fTagEl[1];
170 }
171 else { // all the others
172 fTagEl[1] = 0.;
173 fTagNid[1] = 1.;
174 }
175 }
176
177 if(!pidScheme.CompareTo("ESDCombinedPID")) {
178 fTagEl[0]=fPIDrespEl[0];
179 fTagEl[1]=fPIDrespEl[1];
180 fTagNid[0] = 1.-fTagEl[0];
181 fTagNid[1] = 1.-fTagEl[1];
182 }
183 return;
184}
185//----------------------------------------------------------------------------
186Double_t AliBtoJPSItoEle::ChildrenRelAngle() const {
187 // relative angle between K and pi
188
189 TVector3 mom0(fPx[0],fPy[0],fPz[0]);
190 TVector3 mom1(fPx[1],fPy[1],fPz[1]);
191
192 Double_t angle = mom0.Angle(mom1);
193
194 return angle;
195}
196//----------------------------------------------------------------------------
197void AliBtoJPSItoEle::ComputeWgts() {
198 // calculate the weights for PID
199
200
201 // assignement of the weights from PID
202 fWgtJPsi = fTagEl[0]*fTagEl[1]; // both assumed to be electrons
203
204
205 // if(fWgtJPsi<0.) cerr<<"AliBtoJPSItoEle::ComputeWgts() Negative weight!!!\n";
206
207
208 return;
209}
210//----------------------------------------------------------------------------
211void AliBtoJPSItoEle::CorrectWgt4BR(Double_t factor) {
212 // correct weights of background from charm
213
214 fWgtJPsi *= factor;
215
216 return;
217}
218//----------------------------------------------------------------------------
219Double_t AliBtoJPSItoEle::CosPointing() const {
220 // cosine of pointing angle in space
221
222 TVector3 mom(Px(),Py(),Pz());
223 TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
224
225 Double_t pta = mom.Angle(flight);
226
227 return TMath::Cos(pta);
228}
229//----------------------------------------------------------------------------
230Double_t AliBtoJPSItoEle::CosPointingXY() const {
231 // cosine of pointing angle in transverse plane
232
233 TVector3 momXY(Px(),Py(),0.);
234 TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
235
236 Double_t ptaXY = momXY.Angle(flightXY);
237
238 return TMath::Cos(ptaXY);
239}
240//----------------------------------------------------------------------------
241void AliBtoJPSItoEle::CosThetaStar(Double_t &ctsJPsi) const {
242 // cosine of decay angle in the J/Psi rest frame (of the negative electron)
243
244 Double_t pStar = TMath::Sqrt(TMath::Power(kMJPsi*kMJPsi-2.*kMe*kMe,2.)-4.*kMe*kMe*kMe*kMe)/(2.*kMJPsi);
245
246 Double_t beta = P()/Energy();
247 Double_t gamma = Energy()/kMJPsi;
248
249 ctsJPsi = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMe*kMe))/pStar;
250 // if(ctsJPsi > 1.) { cerr<<"AliBtoJPSItoEle::CosThetaStar: > 1 "<<ctsJPsi<<"!\n"; }
251 // if(ctsJPsi < -1.) { cerr<<"AliBtoJPSItoEle::CosThetaStar: < -1 "<<ctsJPsi<<"!\n"; }
252
253 return;
254}
255//----------------------------------------------------------------------------
256Double_t AliBtoJPSItoEle::Eta() const {
257 // pseudorapidity of the J/Psi
258
259 Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
260 Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
261 return eta;
262}
263//----------------------------------------------------------------------------
264Double_t AliBtoJPSItoEle::EtaChild(Int_t child) const {
265 // pseudorapidity of the decay tracks
266
267 Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
268 Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
269 return eta;
270}
271//----------------------------------------------------------------------------
272void AliBtoJPSItoEle::GetWgts(Double_t &WgtJPsi)
273 const {
274 // returns the weights for pid
275
276 WgtJPsi = fWgtJPsi;
277
278 return;
279}
280//----------------------------------------------------------------------------
281Double_t AliBtoJPSItoEle::ImpPar() const {
282 // J/Psi impact parameter in the bending plane
283
284 Double_t k = -(fV2x-fV1x)*Px()-(fV2y-fV1y)*Py();
285 k /= Pt()*Pt();
286 Double_t dx = fV2x-fV1x+k*Px();
287 Double_t dy = fV2y-fV1y+k*Py();
288 Double_t absDD = TMath::Sqrt(dx*dx+dy*dy);
289 TVector3 mom(Px(),Py(),Pz());
290 TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
291 TVector3 cross = mom.Cross(flight);
292 return (cross.Z()>0. ? absDD : -absDD);
293}
294//----------------------------------------------------------------------------
295void AliBtoJPSItoEle::InvMass(Double_t &mJPsi) const {
296 // invariant mass as J/Psi
297
298 Double_t energy[2];
299
300 // J/psi -> e- e+
301 energy[1] = TMath::Sqrt(kMe*kMe+PChild(1)*PChild(1));
302 energy[0] = TMath::Sqrt(kMe*kMe+PChild(0)*PChild(0));
303
304 mJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
305
306 return;
307
308}
309//----------------------------------------------------------------------------
310Double_t AliBtoJPSItoEle::Ql(Int_t child) const {
311 // longitudinal momentum of decay tracks w.r.t. to J/Psi momentum
312
313 Double_t qL;
314 TVector3 mom(fPx[child],fPy[child],fPz[child]);
315 TVector3 momJPsi(Px(),Py(),Pz());
316
317 qL = mom.Dot(momJPsi)/momJPsi.Mag();
318
319 return qL ;
320}
321//----------------------------------------------------------------------------
322Double_t AliBtoJPSItoEle::Qt() const {
323 // transverse momentum of decay tracks w.r.t. to JPsi momentum
324
325 TVector3 mom0(fPx[0],fPy[0],fPz[0]);
326 TVector3 momJPsi(Px(),Py(),Pz());
327
328 return mom0.Perp(momJPsi);
329}
330//----------------------------------------------------------------------------
331Bool_t AliBtoJPSItoEle::Select(const Double_t* cuts,Int_t& okB)
332 const {
333//
334// This function compares the B candidates with a set of cuts:
335//
336// cuts[0] = inv. mass half width [GeV]
337// cuts[1] = dca [micron]
338// cuts[2] = cosThetaStar
339// cuts[3] = pTP [GeV/c]
340// cuts[4] = pTN [GeV/c]
341// cuts[5] = d0P [micron] upper limit!
342// cuts[6] = d0N [micron] upper limit!
343// cuts[7] = d0d0 [micron^2]
344// cuts[8] = cosThetaPoint
345//
346// If the candidate doesn't pass the cuts it sets the weight to 0
347// and return kFALSE
348//
349 Double_t mJPsi,ctsJPsi;
350 okB=1;
351
352 if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okB = 0;
353 if(!okB) return kFALSE;
354
355 if(TMath::Abs(Getd0Child(1)) > cuts[5] ||
356 TMath::Abs(Getd0Child(0)) > cuts[6]) okB = 0;
357 if(!okB) return kFALSE;
358
359 if(GetDCA() > cuts[1]) { okB = 0; return kFALSE; }
360
361 InvMass(mJPsi);
362 if(TMath::Abs(mJPsi-kMJPsi) > cuts[0]) okB = 0;
363 if(!okB) return kFALSE;
364
365 CosThetaStar(ctsJPsi);
366 if(TMath::Abs(ctsJPsi) > cuts[2]) okB = 0;
367 if(!okB) return kFALSE;
368
369 if(ProdImpParams() > cuts[7]) { okB = 0; return kFALSE; }
370
371 if(CosPointing() < cuts[8]) { okB = 0; return kFALSE; }
372
373 return kTRUE;
374}
375//-----------------------------------------------------------------------------
376void AliBtoJPSItoEle::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
377 // Set combined PID detector response probabilities
378
379 fPIDrespEl[0] = resp0[0];
380 fPIDrespEl[1] = resp1[0];
381 fPIDrespMu[0] = resp0[1];
382 fPIDrespMu[1] = resp1[1];
383 fPIDrespPi[0] = resp0[2];
384 fPIDrespPi[1] = resp1[2];
385 fPIDrespKa[0] = resp0[3];
386 fPIDrespKa[1] = resp1[3];
387 fPIDrespPr[0] = resp0[4];
388 fPIDrespPr[1] = resp1[4];
389
390 return;
391}
392//-----------------------------------------------------------------------------
0ac6b6e8 393Double_t AliBtoJPSItoEle::GetTagEl(Int_t child) const {
394 // Get tag probabilities for electrons
395
396 if(child!=0 && child !=1) return -1;
397
398 return fTagEl[child];
399}
400//-----------------------------------------------------------------------------
27e190fb 401void AliBtoJPSItoEle::GetPIDresponse(Double_t resp0[5],Double_t resp1[5]) const {
402 // Get combined PID detector response probabilities
403
404 resp0[0] = fPIDrespEl[0];
405 resp1[0] = fPIDrespEl[1];
406 resp0[1] = fPIDrespMu[0];
407 resp1[1] = fPIDrespMu[1];
408 resp0[2] = fPIDrespPi[0];
409 resp1[2] = fPIDrespPi[1];
410 resp0[3] = fPIDrespKa[0];
411 resp1[3] = fPIDrespKa[1];
412 resp0[4] = fPIDrespPr[0];
413 resp1[4] = fPIDrespPr[1];
414
415 return;
416}
417//----------------------------------------------------------------------------
418Double_t AliBtoJPSItoEle::TRDTPCCombinedPIDParametrization(Double_t p) const {
419
420
421 // a first raw parametrization of the probability to misidentify a charged pion as electron as a
422 // function of the momentum, as given by the combined TPC and TRD response.
423 // PID cuts are set such that the probability for correct electron id is 90% in each of the two
424 // detectors
425
426// first estimate based on parameterization in the B-> single electron analysis
427 Double_t value=0;
428 Double_t p1 =11.;
429 Double_t p2=0.00007;
430 Double_t p3=0.007;
431 value=p2+p3*(1.-exp(-TMath::Power(p/p1,4.)));
432
27e190fb 433 value/=0.01; // here remove from TPC+TRD the TRD contribution estimated to be 0.01
0ac6b6e8 434// Better estimation based on TRD test beam (as presented by Andrea at Munster)
435// if (p<10.) value*=(1.32-0.18*p+0.076*p*p-0.0037*p*p*p)/100.;
436// if (p>10.) value*=(0.48+0.287*p)/100.;
437// From Silvia MASCIOCCHI (Darmstadt) at PWG3 AliceWeek October 2007
438 if (p<10.) value*=(0.44-0.06*p+0.031*p*p-0.0008*p*p*p)/100.;
439 if (p>10.) value*=(-0.67+0.28*p)/100.;
27e190fb 440
441 return value;
442}
443//----------------------------------------------------------------------------
444
445
446
447