Removing PWG3TrackExtrapInterface from PWG3base library (Gines)
[u/mrichter/AliRoot.git] / PWG3 / AliBtoJPSItoEle.cxx
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
35 ClassImp(AliBtoJPSItoEle)
36
37 //----------------------------------------------------------------------------
38 AliBtoJPSItoEle::AliBtoJPSItoEle():
39 fSignal(kFALSE),
40 fJpsiPrimary(kFALSE),
41 fEvent(0),
42 fV1x(0.),
43 fV1y(0.),
44 fV1z(0.),
45 fV2x(0.),
46 fV2y(0.),
47 fV2z(0.),
48 fDCA(0.),
49 fWgtJPsi(0.)
50 {
51   // Default constructor
52   
53   fTrkNum[0] = 0;
54   fTrkNum[1] = 0;
55
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
80 }
81 //----------------------------------------------------------------------------
82 AliBtoJPSItoEle::AliBtoJPSItoEle(Int_t ev,Int_t trkNum[2],
83                        Double_t v1[3],Double_t v2[3], 
84                        Double_t dca,
85                        Double_t mom[6],Double_t d0[2]):
86 fSignal(kFALSE),
87 fJpsiPrimary(kFALSE),
88 fEvent(ev),
89 fV1x(v1[0]),
90 fV1y(v1[1]),
91 fV1z(v1[2]),
92 fV2x(v2[0]),
93 fV2y(v2[1]),
94 fV2z(v2[2]),
95 fDCA(dca),
96 fWgtJPsi(0.)
97 {
98   // Constructor
99
100   fTrkNum[0] = trkNum[0];
101   fTrkNum[1] = trkNum[1];
102
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.;
126 }
127 //----------------------------------------------------------------------------
128 AliBtoJPSItoEle::~AliBtoJPSItoEle() {}
129 //----------------------------------------------------------------------------
130 void 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 //----------------------------------------------------------------------------
176 Double_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 //----------------------------------------------------------------------------
187 void 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 //----------------------------------------------------------------------------
201 void AliBtoJPSItoEle::CorrectWgt4BR(Double_t factor) {
202   // correct weights of background from charm 
203
204   fWgtJPsi    *= factor;
205
206   return;
207 }
208 //----------------------------------------------------------------------------
209 Double_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 //----------------------------------------------------------------------------
220 Double_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 //----------------------------------------------------------------------------
231 void 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 //----------------------------------------------------------------------------
246 Double_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 //----------------------------------------------------------------------------
254 Double_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 //----------------------------------------------------------------------------
262 void AliBtoJPSItoEle::GetWgts(Double_t &WgtJPsi) 
263   const {
264   // returns the weights for pid
265
266     WgtJPsi    = fWgtJPsi; 
267
268   return;
269 }
270 //----------------------------------------------------------------------------
271 Double_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 //----------------------------------------------------------------------------
285 void 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 //----------------------------------------------------------------------------
300 Double_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 //----------------------------------------------------------------------------
312 Double_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 //----------------------------------------------------------------------------
321 Bool_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 //-----------------------------------------------------------------------------
366 void 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 //-----------------------------------------------------------------------------
383 void 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 //----------------------------------------------------------------------------
400 Double_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