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