KovalenkoQELCharmPXSec.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*
 Copyright (c) 2003-2024, The GENIE Collaboration
 For the full text of the license visit http://copyright.genie-mc.org

 Costas Andreopoulos <c.andreopoulos \at cern.ch>
 University of Liverpool
*/
//____________________________________________________________________________

#include <TMath.h>
#include <Math/Integrator.h>

#include "Framework/Algorithm/AlgConfigPool.h"
#include "Physics/XSectionIntegration/XSecIntegratorI.h"
#include "Physics/Charm/XSection/KovalenkoQELCharmPXSec.h"
#include "Framework/Conventions/GBuild.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Conventions/RefFrame.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Messenger/Messenger.h"
/////#include "Numerical/IntegratorI.h"
#include "Physics/PartonDistributions/PDF.h"
#include "Physics/PartonDistributions/PDFModelI.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/KineUtils.h"
#include "Framework/Numerical/GSLUtils.h"

using namespace genie;
using namespace genie::constants;

//____________________________________________________________________________
KovalenkoQELCharmPXSec::KovalenkoQELCharmPXSec() :
XSecAlgorithmI("genie::KovalenkoQELCharmPXSec")
{

}
//____________________________________________________________________________
KovalenkoQELCharmPXSec::KovalenkoQELCharmPXSec(string config) :
XSecAlgorithmI("genie::KovalenkoQELCharmPXSec", config)
{

}
//____________________________________________________________________________
KovalenkoQELCharmPXSec::~KovalenkoQELCharmPXSec()
{

}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::XSec(
                  const Interaction * interaction, KinePhaseSpace_t kps) const
{
  if(! this -> ValidProcess    (interaction) ) return 0.;
  if(! this -> ValidKinematics (interaction) ) return 0.;

  //----- get kinematics & init state - compute auxiliary vars
  const Kinematics &   kinematics  = interaction->Kine();
  const InitialState & init_state  = interaction->InitState();
  const Target &       target      = init_state.Tgt();

  //neutrino energy & momentum transfer
  double E   = init_state.ProbeE(kRfHitNucRest);
  double E2  = E * E;
  double Q2  = kinematics.Q2();

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("QELCharmXSec", pDEBUG) << "E = " << E << ", Q2 = " << Q2;
#endif

  //resonance mass & nucleon mass
  double MR    = this->MRes  (interaction);
  double MR2   = TMath::Power(MR,2);
  double Mnuc  = target.HitNucMass();
  double Mnuc2 = TMath::Power(Mnuc,2);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("QELCharmXSec", pDEBUG) << "M(RES) = " << MR;
#endif

  //----- Calculate the differential cross section dxsec/dQ^2
  double Gf        = kGF2 / (2*kPi);
  double vR        = (MR2 - Mnuc2 + Q2) / (2*Mnuc);
  double xiR       = this->xiBar(Q2, Mnuc, vR);
  double vR2       = vR*vR;
  double vR_E      = vR/E;
  double Q2_4E2    = Q2/(4*E2);
  double Q2_2MExiR = Q2/(2*Mnuc*E*xiR);
  double Z         = this->ZR(interaction);
  double D         = this->DR(interaction);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("QELCharmXSec", pDEBUG)
     << "Z = " << Z << ", D = " << D << ". xiR = " << xiR << ", vR = " << vR;
#endif

  double xsec = Gf*Z*D * (1 - vR_E + Q2_4E2 + Q2_2MExiR) *
                                     TMath::Sqrt(vR2 + Q2) / (vR*xiR);

  //----- The algorithm computes dxsec/dQ2
  //      Check whether variable tranformation is needed
  if(kps!=kPSQ2fE) {
    double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
    xsec *= J;
  }

  //----- If requested return the free nucleon xsec even for input nuclear tgt
  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;

  //----- Nuclear cross section (simple scaling here)
  int nuc   = target.HitNucPdg();
  int NNucl = (pdg::IsProton(nuc)) ? target.Z() : target.N();
  xsec *= NNucl;

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("QELCharmXSec", pINFO)
     << "dsigma/dQ2(E=" << E << ", Q2=" << Q2 << ") = "
                             << xsec / (1E-40*units::cm2) << " x 1E-40 cm^2";
#endif

  return xsec;
}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::ZR(const Interaction * interaction) const
{
  const XclsTag &      xcls       = interaction->ExclTag();
  const InitialState & init_state = interaction->InitState();

  int pdgc = xcls.CharmHadronPdg();
  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );

  if      ( pdgc == kPdgLambdaPc && isN ) return fScLambdaP;
  else if ( pdgc == kPdgSigmaPc  && isN ) return fScSigmaP;
  else if ( pdgc == kPdgSigmaPPc && isP ) return fScSigmaPP;
  else                                    abort();
}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::DR(const Interaction * interaction) const
{
  const InitialState & init_state = interaction -> InitState();

  // Compute PDFs
  PDF pdfs;
  pdfs.SetModel(fPDFModel);   // <-- attach algorithm

  // Compute integration area = [xi_bar_plus, xi_bar_minus]
  const Kinematics & kinematics = interaction->Kine();

  double Q2     = kinematics.Q2();
  double Mnuc   = init_state.Tgt().HitNucMass();
  double Mnuc2  = TMath::Power(Mnuc,2);
  double MR     = this->MRes(interaction);
  double DeltaR = this->ResDM(interaction);

  double vR_minus  = ( TMath::Power(MR-DeltaR,2) - Mnuc2 + Q2 ) / (2*Mnuc);
  double vR_plus   = ( TMath::Power(MR+DeltaR,2) - Mnuc2 + Q2 ) / (2*Mnuc);

  LOG("QELCharmXSec", pDEBUG)
    << "vR = [plus: " << vR_plus << ", minus: " << vR_minus << "]";

  double xi_bar_minus  = 0.999;
  double xi_bar_plus  = this->xiBar(Q2, Mnuc, vR_plus);

  LOG("QELCharmXSec", pDEBUG)
    << "Integration limits = [" << xi_bar_plus << ", " << xi_bar_minus << "]";

  int pdgc = init_state.Tgt().HitNucPdg();

  ROOT::Math::IBaseFunctionOneDim * integrand = new
          utils::gsl::wrap::KovQELCharmIntegrand(&pdfs,Q2,pdgc);
  ROOT::Math::IntegrationOneDim::Type ig_type =
          utils::gsl::Integration1DimTypeFromString("adaptive");

  double abstol   = 1;    // We mostly care about relative tolerance
  double reltol   = 1E-4;
  int    nmaxeval = 100000;
  ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval);
  double D = ig.Integral(xi_bar_plus, xi_bar_minus);

  delete integrand;

  return D;
}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::xiBar(double Q2, double Mnuc, double v) const
{
  double Mo2 = fMo*fMo;
  double v2  = v *v;
  double xi  = (Q2/Mnuc) / (v + TMath::Sqrt(v2+Q2));
  double xib = xi * ( 1 + (1 + Mo2/(Q2+Mo2))*Mo2/Q2 );
  return xib;
}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::ResDM(const Interaction * interaction) const
{
// Resonance Delta M obeys the constraint DM(R+/-) <= |M(R+/-) - M(R)|
// where M(R-) <= M(R) <= M(R+) are the masses of the neighboring
// resonances R+, R-.
// Get the values from the algorithm conf. registry, and if they do not exist
// set them to default values (Eq.(20) in Sov.J.Nucl.Phys.52:934 (1990)

  const XclsTag & xcls = interaction->ExclTag();

  int pdgc = xcls.CharmHadronPdg();

  bool isLambda = (pdgc == kPdgLambdaPc);
  bool isSigma  = (pdgc == kPdgSigmaPc || pdgc == kPdgSigmaPPc);

  if      ( isLambda ) return fResDMLambda;
  else if ( isSigma  ) return fResDMSigma;
  else                 abort();

  return 0;
}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::MRes(const Interaction * interaction) const
{
  const XclsTag & xcls = interaction->ExclTag();

  int pdgc  = xcls.CharmHadronPdg();
  double MR = PDGLibrary::Instance()->Find(pdgc)->Mass();
  return MR;
}
//____________________________________________________________________________
double KovalenkoQELCharmPXSec::Integral(const Interaction * interaction) const
{
  double xsec = fXSecIntegrator->Integrate(this,interaction);
  return xsec;
}
//____________________________________________________________________________
bool KovalenkoQELCharmPXSec::ValidProcess(
                                        const Interaction * interaction) const
{
  // Make sure we are dealing with one of the following channels:
  // v + n --> mu- + Lambda_{c}^{+} (2285)
  // v + n --> mu- + Sigma_{c}^{+} (2455)
  // v + p --> mu- + Sigma_{c}^{++} (2455)

  if(interaction->TestBit(kISkipProcessChk)) return true;

  const XclsTag &      xcls       = interaction->ExclTag();
  const InitialState & init_state = interaction->InitState();
  const ProcessInfo &  proc_info  = interaction->ProcInfo();

  bool is_exclusive_charm = (xcls.IsCharmEvent() && !xcls.IsInclusiveCharm());
  if(!is_exclusive_charm) return false;

  if(!proc_info.IsQuasiElastic()) return false;
  if(!proc_info.IsWeak())         return false;

  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );

  int pdgc = xcls.CharmHadronPdg();

  bool can_handle = (
     (pdgc == kPdgLambdaPc && isN) || /* v + n -> l + #Lambda_{c}^{+} */
     (pdgc == kPdgSigmaPc  && isN) || /* v + n -> l + #Sigma_{c}^{+}  */
     (pdgc == kPdgSigmaPPc && isP)    /* v + p -> l + #Sigma_{c}^{++} */
  );
  return can_handle;
}
//____________________________________________________________________________
bool KovalenkoQELCharmPXSec::ValidKinematics(
                                        const Interaction * interaction) const
{
  if(interaction->TestBit(kISkipKinematicChk)) return true;

  const InitialState & init_state  = interaction->InitState();
  double E = init_state.ProbeE(kRfHitNucRest);

  //resonance, final state primary lepton & nucleon mass
  double MR    = this -> MRes  (interaction);
  double ml    = interaction->FSPrimLepton()->Mass();
  double Mnuc  = init_state.Tgt().HitNucP4Ptr()->M();
  double Mnuc2 = TMath::Power(Mnuc,2);

  //resonance threshold
  double ER = ( TMath::Power(MR+ml,2) - Mnuc2 ) / (2*Mnuc);

  if(E <= ER) return false;

  return true;
}
//____________________________________________________________________________
void KovalenkoQELCharmPXSec::Configure(const Registry & config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void KovalenkoQELCharmPXSec::Configure(string param_set)
{
  Algorithm::Configure(param_set);
  this->LoadConfig();
}
//____________________________________________________________________________
void KovalenkoQELCharmPXSec::LoadConfig(void)
{
  fPDFModel   = 0;
  ///fIntegrator = 0;

  // Get config values or set defaults
  GetParamDef( "Scale-LambdaP", fScLambdaP,  0.8 * 0.0102 ) ;
  GetParamDef( "Scale-SigmaP",  fScSigmaP ,  0.8 * 0.0028 ) ;
  GetParamDef( "Scale-SigmaPP", fScSigmaPP,  0.8 * 0.0159 ) ;
  GetParamDef( "Res-DeltaM-Lambda", fResDMLambda,  0.56 ) ;      //GeV
  GetParamDef( "Res-DeltaM-Sigma",  fResDMSigma,   0.20 ) ;      //GeV
  GetParamDef( "Mo",                fMo,           sqrt(0.1) );  //GeV

  // get PDF model and integrator

  fPDFModel = dynamic_cast<const PDFModelI *>(this->SubAlg("PDF-Set"));
  assert(fPDFModel);

  // load XSec Integrator
  fXSecIntegrator =
      dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
  assert(fXSecIntegrator);

  // load numerical integrator for integrand in diff x-section calc.
//  fIntegrator = dynamic_cast<const IntegratorI *>(this->SubAlg("Integrator"));
//  assert(fIntegrator);
}
//____________________________________________________________________________
// Auxiliary scalar function for internal integration
//____________________________________________________________________________
utils::gsl::wrap::KovQELCharmIntegrand::KovQELCharmIntegrand(
           PDF * pdf, double Q2, int nucleon_pdgc) :
ROOT::Math::IBaseFunctionOneDim()
{
  fPDF  = pdf;
  fQ2   = TMath::Max(0.3, Q2);
  fPdgC = nucleon_pdgc;
}
//____________________________________________________________________________
utils::gsl::wrap::KovQELCharmIntegrand::~KovQELCharmIntegrand()
{

}
//____________________________________________________________________________
unsigned int utils::gsl::wrap::KovQELCharmIntegrand::NDim(void) const
{
  return 1;
}
//____________________________________________________________________________
double utils::gsl::wrap::KovQELCharmIntegrand::DoEval(double xin) const
{
  if(xin<0 || xin>1) return 0.;

  fPDF->Calculate(xin, fQ2);
  bool isP = pdg::IsProton(fPdgC);
  double f = (isP) ? fPDF->DownValence() : fPDF->UpValence();

  LOG("QELCharmXSec", pDEBUG)
        << "f(xin = " << xin << ", Q2 = " << fQ2 << ") = " << f;

  return f;
}
//____________________________________________________________________________
ROOT::Math::IBaseFunctionOneDim *
  utils::gsl::wrap::KovQELCharmIntegrand::Clone(void) const
{
  return new utils::gsl::wrap::KovQELCharmIntegrand(fPDF, fQ2, fPdgC);
}
//____________________________________________________________________________