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