PhotonCOHPXSec.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

 Alfonso Garcia <aagarciasoto \at km3net.de>
 IFIC & Harvard University
*/
//____________________________________________________________________________

#include "Physics/HELepton/XSection/PhotonCOHPXSec.h"
#include "Physics/XSectionIntegration/XSecIntegratorI.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/RefFrame.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/KineUtils.h"

using namespace genie;
using namespace genie::constants;

const double a  = 0.523 * (units::fermi);
const double r0 = 1.126 * (units::fermi);

//____________________________________________________________________________
PhotonCOHPXSec::PhotonCOHPXSec() :
XSecAlgorithmI("genie::PhotonCOHPXSec")
{
  born = new Born();
}
//____________________________________________________________________________
PhotonCOHPXSec::PhotonCOHPXSec(string config) :
XSecAlgorithmI("genie::PhotonCOHPXSec", config)
{

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

}
//____________________________________________________________________________
double PhotonCOHPXSec::XSec(
   const Interaction * interaction, KinePhaseSpace_t kps) const
{

  if(! this -> ValidProcess    (interaction) ) return 0.;

  const InitialState & init_state = interaction -> InitState();
  const Kinematics &   kinematics = interaction -> Kine();

  int probepdg = init_state.ProbePdg();
  double E     = init_state.ProbeE(kRfLab);

  double mlout = 0;
  if      (pdg::IsNuE  (TMath::Abs(probepdg))) mlout = kElectronMass;
  else if (pdg::IsNuMu (TMath::Abs(probepdg))) mlout = kMuonMass;
  else if (pdg::IsNuTau(TMath::Abs(probepdg))) mlout = kTauMass;
  double mlout2 = mlout*mlout;

  int A        = init_state.Tgt().A(); 
  int Z        = init_state.Tgt().Z(); 
  double Mtgt = Z*kProtonMass + (A-Z)*kNeutronMass;

  double n1 = kinematics.GetKV(kKVn1);
  double n2 = kinematics.GetKV(kKVn2);
  double n3 = kinematics.GetKV(kKVn3);

  double mL = mlout+kMw;
  
  double Delta = TMath::Power(2.*E*Mtgt-mL*mL,2)-4.*TMath::Power(Mtgt*mL,2);
  if (Delta<0) return 0.;

  double s12_min = E/(2.*E+Mtgt)*(mL*mL+2.*E*Mtgt-TMath::Sqrt(Delta));
  double s12_max = E/(2.*E+Mtgt)*(mL*mL+2.*E*Mtgt+TMath::Sqrt(Delta));
  double s12 = TMath::Exp( TMath::Log(s12_min)+(TMath::Log(s12_max)-TMath::Log(s12_min))*n2);

  double Q2_min = TMath::Power(s12,2)*Mtgt/2./E/(2.*E*Mtgt-s12);
  double Q2_max = s12 - mL*mL;
  double Q2 = TMath::Exp( TMath::Log(Q2_min) + (TMath::Log(Q2_max)-TMath::Log(Q2_min))*n3 );

  double s = s12 - Q2;
  double s13 = s12/2.*((1.+kMw2/s-mlout2/s)-TMath::Sqrt(born->Lambda(1.,kMw2/s,mlout2/s))*n1);

  double ME_T = born->PXSecPhoton_T(s12,s13,Q2,mlout2) * (1.-s12/2./E/Mtgt-TMath::Power(s12/2./E,2)/Q2);
  double ME_L = born->PXSecPhoton_L(s12,s13,Q2,mlout2) * TMath::Power(1.-s12/4./E/Mtgt,2);

  double ME   = ME_T+ME_L;
  double dps2 = 1/32./kPi2/s12 * TMath::Sqrt( born->Lambda(1.,kMw2/s,mlout2/s) ) * (TMath::Log(Q2_max)-TMath::Log(Q2_min)) * (TMath::Log(s12_max)-TMath::Log(s12_min));    
  double dP   = born->GetReAlpha()*TMath::Power(Z,2)*F2_Q( TMath::Sqrt(Q2), r0*TMath::Power(A,1./3.) );

  double xsec = ME * dps2 * dP;
   
  if(kps!=kPSn1n2n3fE) {
      LOG("PhotonCOHPXSec", pWARN)
          << "Doesn't support transformation from "
          << KinePhaseSpace::AsString(kPSn1n2n3fE) << " to "
          << KinePhaseSpace::AsString(kps);
      xsec = 0;
  }

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

  LOG("PhotonCOHPXSec", pINFO) << "dxsec/dn1dn2dn3 (E= " << E << ", n1= " << n1 << ", n2=" << n2 << ", n3=" << n3 << ") = " << xsec;

  return xsec;
}
//____________________________________________________________________________
double PhotonCOHPXSec::F2_Q(double Q, double r0) const
{
  // Analytic Woods-Saxon, A.3 of https://arxiv.org/pdf/1807.10973.pdf
  double FF = 0.0;
  double coth = 1./TMath::TanH(kPi*Q*a);
  FF = 3.*kPi*a/(TMath::Power(r0,2)+TMath::Power(kPi*a,2)) / (Q*r0* TMath::SinH(kPi*Q*a));
  FF *= (kPi*a*coth*TMath::Sin(Q*r0) - r0 * TMath::Cos(Q*r0));        
  return TMath::Power(FF,2);
}
//____________________________________________________________________________
double PhotonCOHPXSec::Integral(const Interaction * interaction) const
{
  double xsec = fXSecIntegrator->Integrate(this,interaction);
  return xsec;
}
//____________________________________________________________________________
bool PhotonCOHPXSec::ValidProcess(const Interaction* interaction) const
{
  if(interaction->TestBit(kISkipProcessChk)) return true;

  const ProcessInfo & proc_info  = interaction->ProcInfo();
  if(!proc_info.IsPhotonCoherent()) return false;

  const InitialState & init_state = interaction -> InitState();
  if(!pdg::IsLepton(init_state.ProbePdg())) return false;

  if(init_state.Tgt().HitNucIsSet()) return false;
  
  return true;
}
//____________________________________________________________________________


//____________________________________________________________________________
void PhotonCOHPXSec::Configure(const Registry & config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void PhotonCOHPXSec::Configure(string config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void PhotonCOHPXSec::LoadConfig(void)
{

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

}