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 #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