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