memory consumption measures only memory used by preprocessor
[u/mrichter/AliRoot.git] / PWG3 / vertexingOld / 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 #include <TString.h>
33
34 #include "AliBtoJPSItoEle.h"
35
36 ClassImp(AliBtoJPSItoEle)
37
38 //----------------------------------------------------------------------------
39 AliBtoJPSItoEle::AliBtoJPSItoEle():
40 fSignal(kFALSE),
41 fJpsiPrimary(kFALSE),
42 fEvent(0),
43 fV1x(0.),
44 fV1y(0.),
45 fV1z(0.),
46 fV2x(0.),
47 fV2y(0.),
48 fV2z(0.),
49 fDCA(0.),
50 fWgtJPsi(0.)
51 {
52   // Default constructor
53   
54
55   fTrkNum[0] = 0;
56   fTrkNum[1] = 0;
57
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
82 }
83 //----------------------------------------------------------------------------
84 AliBtoJPSItoEle::AliBtoJPSItoEle(Int_t ev,Int_t trkNum[2],
85                        Double_t v1[3],Double_t v2[3], 
86                        Double_t dca,
87                        Double_t mom[6],Double_t d0[2]) :
88 fSignal(kFALSE),
89 fJpsiPrimary(kFALSE),
90 fEvent(ev),
91 fV1x(v1[0]),
92 fV1y(v1[1]),
93 fV1z(v1[2]),
94 fV2x(v2[0]),
95 fV2y(v2[1]),
96 fV2z(v2[2]),
97 fDCA(dca),
98 fWgtJPsi(0.)
99 {
100   // Constructor
101
102   fTrkNum[0] = trkNum[0];
103   fTrkNum[1] = trkNum[1];
104
105   fV1x = v1[0];
106   fV1y = v1[1];
107   fV1z = v1[2];
108   fV2x = v2[0];
109   fV2y = v2[1];
110   fV2z = v2[2];
111
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.;
135
136 }
137 //----------------------------------------------------------------------------
138 AliBtoJPSItoEle::~AliBtoJPSItoEle() {}
139 //----------------------------------------------------------------------------
140 void 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 //----------------------------------------------------------------------------
186 Double_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 //----------------------------------------------------------------------------
197 void 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 //----------------------------------------------------------------------------
211 void AliBtoJPSItoEle::CorrectWgt4BR(Double_t factor) {
212   // correct weights of background from charm 
213
214   fWgtJPsi    *= factor;
215
216   return;
217 }
218 //----------------------------------------------------------------------------
219 Double_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 //----------------------------------------------------------------------------
230 Double_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 //----------------------------------------------------------------------------
241 void 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 //----------------------------------------------------------------------------
256 Double_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 //----------------------------------------------------------------------------
264 Double_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 //----------------------------------------------------------------------------
272 void AliBtoJPSItoEle::GetWgts(Double_t &WgtJPsi) 
273   const {
274   // returns the weights for pid
275
276     WgtJPsi    = fWgtJPsi; 
277
278   return;
279 }
280 //----------------------------------------------------------------------------
281 Double_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 //----------------------------------------------------------------------------
295 void 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 //----------------------------------------------------------------------------
310 Double_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 //----------------------------------------------------------------------------
322 Double_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 //----------------------------------------------------------------------------
331 Bool_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 //-----------------------------------------------------------------------------
376 void 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 //-----------------------------------------------------------------------------
393 Double_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 //-----------------------------------------------------------------------------
401 void 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 //----------------------------------------------------------------------------
418 Double_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
433   value/=0.01; // here remove from TPC+TRD the TRD contribution estimated to be 0.01
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.;
440
441   return value;
442 }
443 //----------------------------------------------------------------------------
444
445
446
447