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