BMPTCalc.cxx

Go to the documentation of this file.

#ifndef bmptcalc_cxx
#define bmptcalc_cxx
#include "MessageService/MsgService.h"
#include "MCReweight/BMPTCalc.h"
#include "MCReweight/BMPTEvent.h"
#include "TMath.h"

//****************************
BMPTCalc::BMPTCalc()
{
}

//****************************
BMPTCalc::BMPTCalc(BMPTConfig *config)
{
  fConfig = *config;
}

//****************************
BMPTCalc::~BMPTCalc()
{

}

//****************************
double BMPTCalc::YLD(BMPTEvent *ev,BMPTConfig *config)
{

  //Return yield of secondaries (tertiaries, etc.) produced as a
  //function of the secondary's energy, transverse momentum, and
  //distance along target. (Eqn. 12)
  //    E         Secondary energy (GeV)
  //    PID       Particle type (GEANT3 convention)
  //    PX        Secondary X momentum (GeV/c)
  //    PY        Secondary Y momentum (GeV/c)
  //    PZ        Secondary Z momentum (GeV/c)
  //    Z         Distance along target (cm)

  BMPTConfig *conf = &fConfig;
  if(config!=0) conf = config;

  double px = ev->px;
  double py = ev->py;
  double pz = ev->pz;

  double AvNo = 6.02e23;  //Avogadro's number
  double M = ev->Mass();
  double P = pz*sqrt((px/pz)*(px/pz) + (py/pz)*(py/pz) + 1.0);
  double E = P*sqrt((M/P)*(M/P) + 1.0);

  //Work in log's to extend dynamic range
  double DYDZ =
    TMath::Log(AvNo) + TMath::Log(conf->GetTARGET_RHO()) + 
    TMath::Log(conf->GetLAMBDA_P()) + 3.0*TMath::Log(P) + 
    TMath::Log(this->XSec(ev,config)) + 
    TMath::Log(this->FZ(ev,config)) - 
    TMath::Log(100.0) - TMath::Log(conf->GetA_TARGET()) - TMath::Log(E);
  return TMath::Exp(DYDZ);

}

//****************************
double BMPTCalc::FZ(BMPTEvent *ev,BMPTConfig *config)
{

  BMPTConfig *conf = &fConfig;
  if(config!=0) conf = config;

  double fz = 
    ( TMath::Exp(-ev->z/conf->GetLAMBDA_P())/conf->GetLAMBDA_P() ) *
    ( ( 1.0 + this->AH(ev,config)*(ev->z/conf->GetLAMBDA_P())) ) *
    ( TMath::Exp( -this->ZPRIME(ev)/conf->GetLAMBDA_S() ) );
  return fz;
  
}

//****************************
double BMPTCalc::AH(BMPTEvent *ev,BMPTConfig *config)
{
  
  //Factor to account for tertiary production See Eqn. 17 and Table 3
  BMPTConfig *conf = &fConfig;
  if(config!=0) conf = config;

  double x_lab = sqrt(ev->Momentum()*ev->Momentum() + 
                      ev->Mass()*ev->Mass())/conf->GetBEAM_P();

  //Pions:
  if(ev->pid==8 || ev->pid==9) {
    double ATertiaryH = 0.80;
    double BTertiaryH = 7.30;
    return ATertiaryH * TMath::Power(1.0 - x_lab,BTertiaryH);
  }
  
  //Kaons:
  if(ev->pid==10 || ev->pid==11 || 
     ev->pid==12 || ev->pid==16){
    double ATertiaryH = 1.56;
    double BTertiaryH = 10.1;
    return ATertiaryH * TMath::Power(1.0 - x_lab,BTertiaryH);
  }

  return 1.0;

}

//****************************
double BMPTCalc::ZPRIME(BMPTEvent *ev)
{
  //assume (0,0,-0.296) is front, centre of target
  double z0 = -0.296; //LE beam
  double tZ = 0.9;  //length of target (m)
  double tX = 6.4e-3;  //width of target (m)
  double tY = 20.0e-3; //height of target (m)
  double x1 = ev->x/100.;  //x,y,z are in cm
  double y1 = ev->y/100.;
  double z1 = ev->z/100.;
  double x2 = x1;
  double y2 = y1;
  double z2 = z1;
  double dxdz = ev->px/ev->pz;
  double dydz = ev->py/ev->pz;
  double a1 = (tX/2. - x1)/dxdz;
  double a2 = (-tX/2. - x1)/dxdz;
  double a = a1; if(a2>a1) a = a2; 
  double b1 = (tY/2. - y1)/dydz;
  double b2 = (-tY/2. - y1)/dydz;
  double b = b1; if(b2>b1) b = b2;
  double c = tZ - z0 - z1;
  if(c<a&&c<b) z2 += c;
  else if(a<c&&a<b) z2 += a;
  else if(b<a&&b<c) z2 += b;
  x2 += dxdz*(z2-z1);
  y2 += dydz*(z2-z1);
  return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
}

//****************************
double BMPTCalc::XSec(BMPTEvent *ev,BMPTConfig *config)
{
  
  //Return a weight (E * d3sigma/dp3 ) for a particle produced in a
  //hadron-nucleus interaction according to the BMPT (hep-ph/0101163)
  //parameterization
  BMPTConfig *conf = &fConfig;
  if(config!=0) conf = config;

  double XF    = ev->pz/conf->GetBEAM_P();
  double XR    = XF;
  double PT    = sqrt(ev->px*ev->px + ev->py*ev->py);

  double CPA   = conf->GetCPA(ev->pid);
  double CPB   = conf->GetCPB(ev->pid);
  double SMA   = conf->GetSMA(ev->pid);
  double SMB   = conf->GetSMB(ev->pid);
  double ALPHA = conf->GetALPHA(ev->pid);
  double BETA  = conf->GetBETA(ev->pid);
  double GAMMA = conf->GetGAMMA(ev->pid);
  double DELTA = conf->GetDELTA(ev->pid);
  double R0    = conf->GetR0(ev->pid);
  double R1    = conf->GetR1(ev->pid);

  double A_PRIME = SMA/TMath::Power(XR,GAMMA);
  double B_PRIME = SMA*SMA/(2.0*TMath::Power(XR,DELTA));

  //Calculate the production cross-section for pi/K/p/pbar
  //Work in log's first to accomodate large dynamic range
  double ED3SDP3 = 0;

  if (ev->pid==8  || ev->pid==9  ||  // pi+ / pi-
      ev->pid==10 || ev->pid==11 ||  // K+/K-/K0S/K0L
      ev->pid==12 || ev->pid==16 ) {
    ED3SDP3 =
      (  TMath::Log(CPA) ) +
      (  ALPHA*TMath::Log(1.0 - XR) ) +
      (  TMath::Log(1.0 + CPB*XR) ) +
      ( -BETA*TMath::Log(XR) ) +
      (  TMath::Log( 1.0 + A_PRIME*PT + B_PRIME*PT*PT ) ) +
      ( -A_PRIME*PT );
  }
  else if (ev->pid==14) {      // proton
    ED3SDP3 =
      ( TMath::Log(CPA) ) +
      ( TMath::Log( 1.0 + CPB*XR) ) +
      ( SMB * PT * PT * TMath::Log(1.0 - XR) ) +
      ( TMath::Log( 1.0 + SMA*PT + (0.5*SMA*SMA)*PT*PT) ) +
      ( -SMA*PT );
  }
  else if (ev->pid==15) {       // anti-proton
    ED3SDP3 =
      ( TMath::Log(CPA) ) +
      ( ALPHA * TMath::Log(1.0 - XR ) ) +
      ( -BETA*TMath::Log(XR) ) +
      ( TMath::Log(1.0 + SMA*PT + (0.5*SMA*SMA)*PT*PT ) ) +
      ( -SMA*PT );
  }

  //Calculate the ratios of negatives to positives
  double R = 1.0;
  if(ev->pid==9) R = R0 * TMath::Power((1.0 + XR),R1);
  else if(ev->pid==10 || ev->pid==12 || 
          ev->pid==16) R = R0 * TMath::Power((1.0 - XR),R1);
    
  //Calculate production off of Be
  ED3SDP3 = TMath::Exp(ED3SDP3)/R;
  
  //K0's are a special case
  if(ev->pid==10 || ev->pid==16)
    ED3SDP3 = ED3SDP3 * ( 0.25*(R + 3.0) );

  //Scale to different target nuclei  
  ALPHA = (0.74-0.55*XF+0.26*XF*XF)*(0.98+0.21*PT*PT);
  ED3SDP3 = ED3SDP3*TMath::Power((conf->GetA_TARGET()/conf->GetA_BE()),ALPHA);

  return ED3SDP3;

}

//****************************
double BMPTCalc::Reweight(BMPTEvent *ev,BMPTConfig *rwtConfig,
                          BMPTConfig *config)
{  
  if(!rwtConfig) return 1;
  if(!ev)        return 1;

  double x0 = XSec(ev,config);
  double x1 = XSec(ev,rwtConfig);

  if(x0!=0) return x1/x0;
  return 1;
}

#endif

---