[u/mrichter/AliRoot.git] / PWG3 / vertexingOld / AliBtoJPSItoEle.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//----------------------------------------------------------------------------
//               Implementation of the BtoJPSItoEle class
//                  for pp and PbPb interactions
// Note: the two decay tracks are labelled: 0 (positive electron)
//                                          1 (negative electron)
//            Origin: G.E. Bruno    giuseppe.bruno@ba.infn.it            
//  based on Class for charm golden channel (D0->Kpi)
//----------------------------------------------------------------------------

// #include <Riostream.h> // for debugging

#include <TH1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TPaveLabel.h>
#include <TVector3.h>
#include <TString.h>

#include "AliBtoJPSItoEle.h"

ClassImp(AliBtoJPSItoEle)

//----------------------------------------------------------------------------
AliBtoJPSItoEle::AliBtoJPSItoEle():
fSignal(kFALSE),
fJpsiPrimary(kFALSE),
fEvent(0),
fV1x(0.),
fV1y(0.),
fV1z(0.),
fV2x(0.),
fV2y(0.),
fV2z(0.),
fDCA(0.),
fWgtJPsi(0.)
{
  // Default constructor
  

  fTrkNum[0] = 0;
  fTrkNum[1] = 0;

  fPx[0] = 0.;
  fPy[0] = 0.;
  fPz[0] = 0.;
  fPx[1] = 0.;
  fPy[1] = 0.;
  fPz[1] = 0.;

  fd0[0] = 0.;
  fd0[1] = 0.;

  fPdg[0] = 0;
  fPdg[1] = 0;
  fMum[0] = 0;
  fMum[1] = 0;
  fGMum[0] = 0;
  fGMum[1] = 0;

  fTagPi[0] = 0.;
  fTagPi[1] = 0.;
  fTagKa[0] = 0.;
  fTagKa[1] = 0.;
  fTagNid[0] = 0.;
  fTagNid[1] = 0.;

}
//----------------------------------------------------------------------------
AliBtoJPSItoEle::AliBtoJPSItoEle(Int_t ev,Int_t trkNum[2],
		       Double_t v1[3],Double_t v2[3], 
		       Double_t dca,
		       Double_t mom[6],Double_t d0[2]) :
fSignal(kFALSE),
fJpsiPrimary(kFALSE),
fEvent(ev),
fV1x(v1[0]),
fV1y(v1[1]),
fV1z(v1[2]),
fV2x(v2[0]),
fV2y(v2[1]),
fV2z(v2[2]),
fDCA(dca),
fWgtJPsi(0.)
{
  // Constructor

  fTrkNum[0] = trkNum[0];
  fTrkNum[1] = trkNum[1];

  fV1x = v1[0];
  fV1y = v1[1];
  fV1z = v1[2];
  fV2x = v2[0];
  fV2y = v2[1];
  fV2z = v2[2];

  fPx[0] = mom[0];
  fPy[0] = mom[1];
  fPz[0] = mom[2];
  fPx[1] = mom[3];
  fPy[1] = mom[4];
  fPz[1] = mom[5];

  fd0[0] = d0[0];
  fd0[1] = d0[1];

  fPdg[0] = 0;
  fPdg[1] = 0;
  fMum[0] = 0;
  fMum[1] = 0;
  fGMum[0] = 0;
  fGMum[1] = 0;

  fTagPi[0]  = 0.;
  fTagPi[1]  = 0.;
  fTagKa[0]  = 0.;
  fTagKa[1]  = 0.;
  fTagNid[0] = 0.;
  fTagNid[1] = 0.;

}
//----------------------------------------------------------------------------
AliBtoJPSItoEle::~AliBtoJPSItoEle() {}
//----------------------------------------------------------------------------
void AliBtoJPSItoEle::ApplyPID(TString pidScheme) {
  // Applies particle identification

  if(!pidScheme.CompareTo("TRDTPCparam")  && fPdg[0]==0) {
    printf("AliBtoJPSItoEle::ApplyPID :\n Warning: TRD-TPC parameterized PID can be used only for simulation!\n"); 
    return;
  }

  if(!pidScheme.CompareTo("TRDTPCparam")) {
    // tagging of the positive track
    if(TMath::Abs(fPdg[0])==11) { // electron
      fTagEl[0] = 0.81;
      fTagNid[0] = 1.-fTagEl[0];
    }
    else if(TMath::Abs(fPdg[0])==211) { // pion
      fTagEl[0]   = TRDTPCCombinedPIDParametrization(PChild(0));
      fTagNid[0] = 1.-fTagEl[0];
    } 
    else { // all the others 
      fTagEl[0]  = 0.;
      fTagNid[0] = 1.;
    } 
    // tagging of the negative track
    if(TMath::Abs(fPdg[1])==11) { // electron
      fTagEl[1] = 0.81;
      fTagNid[1] = 1.-fTagEl[1];
    }
    else if(TMath::Abs(fPdg[1])==211) { // pion
      fTagEl[1]   = TRDTPCCombinedPIDParametrization(PChild(1));
      fTagNid[1] = 1.-fTagEl[1];
    }
    else { // all the others
      fTagEl[1]  = 0.;
      fTagNid[1] = 1.;
    }
  }

  if(!pidScheme.CompareTo("ESDCombinedPID")) {
    fTagEl[0]=fPIDrespEl[0];
    fTagEl[1]=fPIDrespEl[1];
    fTagNid[0] = 1.-fTagEl[0];
    fTagNid[1] = 1.-fTagEl[1];
  }
  return;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::ChildrenRelAngle() const {
  // relative angle between K and pi

  TVector3 mom0(fPx[0],fPy[0],fPz[0]);
  TVector3 mom1(fPx[1],fPy[1],fPz[1]);

  Double_t angle = mom0.Angle(mom1);

  return angle; 
}
//----------------------------------------------------------------------------
void AliBtoJPSItoEle::ComputeWgts() {
  // calculate the weights for PID


  // assignement of the weights from PID
  fWgtJPsi    = fTagEl[0]*fTagEl[1]; // both assumed to be electrons 

  
  // if(fWgtJPsi<0.) cerr<<"AliBtoJPSItoEle::ComputeWgts()  Negative weight!!!\n";
  

  return;
}
//----------------------------------------------------------------------------
void AliBtoJPSItoEle::CorrectWgt4BR(Double_t factor) {
  // correct weights of background from charm 

  fWgtJPsi    *= factor;

  return;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::CosPointing() const {
  // cosine of pointing angle in space

  TVector3 mom(Px(),Py(),Pz());
  TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);

  Double_t pta = mom.Angle(flight);

  return TMath::Cos(pta); 
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::CosPointingXY() const {
  // cosine of pointing angle in transverse plane

  TVector3 momXY(Px(),Py(),0.);
  TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);

  Double_t ptaXY = momXY.Angle(flightXY);

  return TMath::Cos(ptaXY); 
}
//----------------------------------------------------------------------------
void AliBtoJPSItoEle::CosThetaStar(Double_t &ctsJPsi) const {
  // cosine of decay angle in the J/Psi rest frame (of the negative electron)

  Double_t pStar = TMath::Sqrt(TMath::Power(kMJPsi*kMJPsi-2.*kMe*kMe,2.)-4.*kMe*kMe*kMe*kMe)/(2.*kMJPsi);

  Double_t beta = P()/Energy();
  Double_t gamma = Energy()/kMJPsi;

  ctsJPsi = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMe*kMe))/pStar;
  // if(ctsJPsi > 1.)  { cerr<<"AliBtoJPSItoEle::CosThetaStar: > 1 "<<ctsJPsi<<"!\n"; }
  // if(ctsJPsi < -1.) { cerr<<"AliBtoJPSItoEle::CosThetaStar: < -1 "<<ctsJPsi<<"!\n"; }

  return;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::Eta() const {
  // pseudorapidity of the J/Psi

  Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
  Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
  return eta;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::EtaChild(Int_t child) const {
  // pseudorapidity of the decay tracks

  Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
  Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
  return eta;
}
//----------------------------------------------------------------------------
void AliBtoJPSItoEle::GetWgts(Double_t &WgtJPsi) 
  const {
  // returns the weights for pid

    WgtJPsi    = fWgtJPsi; 

  return;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::ImpPar() const {
  // J/Psi impact parameter in the bending plane
  
    Double_t k = -(fV2x-fV1x)*Px()-(fV2y-fV1y)*Py();
    k /= Pt()*Pt();
    Double_t dx = fV2x-fV1x+k*Px();
    Double_t dy = fV2y-fV1y+k*Py();
    Double_t absDD = TMath::Sqrt(dx*dx+dy*dy);
    TVector3 mom(Px(),Py(),Pz());
    TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
    TVector3 cross = mom.Cross(flight);
    return (cross.Z()>0. ? absDD : -absDD);
}
//----------------------------------------------------------------------------
void AliBtoJPSItoEle::InvMass(Double_t &mJPsi) const {
  // invariant mass as J/Psi

  Double_t energy[2];

  // J/psi -> e- e+
  energy[1] = TMath::Sqrt(kMe*kMe+PChild(1)*PChild(1));
  energy[0] = TMath::Sqrt(kMe*kMe+PChild(0)*PChild(0));

  mJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
    
  return;

}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::Ql(Int_t child) const {
  // longitudinal momentum of decay tracks w.r.t. to J/Psi momentum

  Double_t qL;
  TVector3 mom(fPx[child],fPy[child],fPz[child]);
  TVector3 momJPsi(Px(),Py(),Pz());

  qL = mom.Dot(momJPsi)/momJPsi.Mag();

  return qL ;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::Qt() const {
  // transverse momentum of decay tracks w.r.t. to JPsi momentum  

  TVector3 mom0(fPx[0],fPy[0],fPz[0]);
  TVector3 momJPsi(Px(),Py(),Pz());

  return mom0.Perp(momJPsi);
}
//----------------------------------------------------------------------------
Bool_t AliBtoJPSItoEle::Select(const Double_t* cuts,Int_t& okB) 
  const {
//
// This function compares the B candidates with a set of cuts:
//
// cuts[0] = inv. mass half width [GeV]   
// cuts[1] = dca [micron]
// cuts[2] = cosThetaStar 
// cuts[3] = pTP [GeV/c]
// cuts[4] = pTN [GeV/c]
// cuts[5] = d0P [micron]   upper limit!
// cuts[6] = d0N [micron]  upper limit!
// cuts[7] = d0d0 [micron^2]
// cuts[8] = cosThetaPoint
//
// If the candidate doesn't pass the cuts it sets the weight to 0
// and return kFALSE
//
  Double_t mJPsi,ctsJPsi;
  okB=1; 

  if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okB = 0;
  if(!okB) return kFALSE;

  if(TMath::Abs(Getd0Child(1)) > cuts[5] || 
     TMath::Abs(Getd0Child(0)) > cuts[6]) okB = 0;
  if(!okB) return kFALSE;

  if(GetDCA() > cuts[1]) { okB = 0; return kFALSE; }

  InvMass(mJPsi);
  if(TMath::Abs(mJPsi-kMJPsi)    > cuts[0]) okB = 0;
  if(!okB) return kFALSE;

  CosThetaStar(ctsJPsi);
  if(TMath::Abs(ctsJPsi)    > cuts[2]) okB = 0;
  if(!okB) return kFALSE;

  if(ProdImpParams() > cuts[7]) { okB = 0; return kFALSE; }

  if(CosPointing()   < cuts[8]) { okB = 0; return kFALSE; }

  return kTRUE;
}
//-----------------------------------------------------------------------------
void AliBtoJPSItoEle::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
  // Set combined PID detector response probabilities

  fPIDrespEl[0] = resp0[0];
  fPIDrespEl[1] = resp1[0];
  fPIDrespMu[0] = resp0[1];
  fPIDrespMu[1] = resp1[1];
  fPIDrespPi[0] = resp0[2];
  fPIDrespPi[1] = resp1[2];
  fPIDrespKa[0] = resp0[3];
  fPIDrespKa[1] = resp1[3];
  fPIDrespPr[0] = resp0[4];
  fPIDrespPr[1] = resp1[4];

  return;
} 
//-----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::GetTagEl(Int_t child) const {
  // Get tag probabilities for electrons
                                                                                                                          
  if(child!=0 && child !=1) return -1;

  return fTagEl[child];
}
//-----------------------------------------------------------------------------
void AliBtoJPSItoEle::GetPIDresponse(Double_t resp0[5],Double_t resp1[5]) const {
  // Get combined PID detector response probabilities
                                                                                                                          
  resp0[0] = fPIDrespEl[0];
  resp1[0] = fPIDrespEl[1]; 
  resp0[1] = fPIDrespMu[0];
  resp1[1] = fPIDrespMu[1];
  resp0[2] = fPIDrespPi[0];
  resp1[2] = fPIDrespPi[1];
  resp0[3] = fPIDrespKa[0];
  resp1[3] = fPIDrespKa[1];
  resp0[4] = fPIDrespPr[0];
  resp1[4] = fPIDrespPr[1];
                                                                                                                          
  return;
}
//----------------------------------------------------------------------------
Double_t AliBtoJPSItoEle::TRDTPCCombinedPIDParametrization(Double_t p) const {


  // a first raw parametrization of the probability to misidentify a charged pion as electron as a 
  // function of the momentum, as given by the combined TPC and TRD response. 
  //  PID cuts are set such that the probability for correct electron id is 90% in each of the two 
  //    detectors
  
// first estimate based on parameterization in the B-> single electron analysis
  Double_t value=0;
  Double_t p1 =11.;
  Double_t p2=0.00007;
  Double_t p3=0.007;
  value=p2+p3*(1.-exp(-TMath::Power(p/p1,4.)));

  value/=0.01; // here remove from TPC+TRD the TRD contribution estimated to be 0.01
// Better estimation based on TRD test beam (as presented by Andrea at Munster)
//  if (p<10.) value*=(1.32-0.18*p+0.076*p*p-0.0037*p*p*p)/100.;
//  if (p>10.) value*=(0.48+0.287*p)/100.;
// From Silvia MASCIOCCHI (Darmstadt) at PWG3 AliceWeek October 2007
  if (p<10.) value*=(0.44-0.06*p+0.031*p*p-0.0008*p*p*p)/100.;
  if (p>10.) value*=(-0.67+0.28*p)/100.;

  return value;
}
//----------------------------------------------------------------------------
Commit	Line	Data
27e190fb	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>
0ac6b6e8	32	#include <TString.h>
27e190fb	33
	34	#include "AliBtoJPSItoEle.h"
	35
	36	ClassImp(AliBtoJPSItoEle)
	37
	38	//----------------------------------------------------------------------------
cbd5890c	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	{
27e190fb	52	// Default constructor
27e190fb	53
0ac6b6e8	54
27e190fb	55	fTrkNum[0] = 0;
	56	fTrkNum[1] = 0;
	57
27e190fb	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
27e190fb	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,
cbd5890c	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	{
27e190fb	100	// Constructor
27e190fb	101
27e190fb	102	fTrkNum[0] = trkNum[0];
	103	fTrkNum[1] = trkNum[1];
	104
0ac6b6e8	105	fV1x = v1[0];
	106	fV1y = v1[1];
	107	fV1z = v1[2];
	108	fV2x = v2[0];
	109	fV2y = v2[1];
	110	fV2z = v2[2];
0ac6b6e8	111
27e190fb	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.;
0ac6b6e8	135
27e190fb	136	}
	137	//----------------------------------------------------------------------------
	138	AliBtoJPSItoEle::~AliBtoJPSItoEle() {}
27e190fb	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(kMJPsikMJPsi-2.kMekMe,2.)-4.kMekMekMekMe)/(2.kMJPsi);
245
246	Double_t beta = P()/Energy();
247	Double_t gamma = Energy()/kMJPsi;
248
249	ctsJPsi = (Ql(1)/gamma-betaTMath::Sqrt(pStarpStar+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(dxdx+dydy);
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(kMekMe+PChild(1)PChild(1));
302	energy[0] = TMath::Sqrt(kMekMe+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	//-----------------------------------------------------------------------------
0ac6b6e8	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	//-----------------------------------------------------------------------------
27e190fb	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
27e190fb	433	value/=0.01; // here remove from TPC+TRD the TRD contribution estimated to be 0.01
0ac6b6e8	434	// Better estimation based on TRD test beam (as presented by Andrea at Munster)
	435	// if (p<10.) value=(1.32-0.18p+0.076pp-0.0037pp*p)/100.;
	436	// if (p>10.) value=(0.48+0.287p)/100.;
	437	// From Silvia MASCIOCCHI (Darmstadt) at PWG3 AliceWeek October 2007
	438	if (p<10.) value=(0.44-0.06p+0.031pp-0.0008pp*p)/100.;
	439	if (p>10.) value=(-0.67+0.28p)/100.;
27e190fb	440
	441	return value;
	442	}
	443	//----------------------------------------------------------------------------
	444
	445
	446
	447