SpectralFunc1d.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*
 Copyright (c) 2003-2025, 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
 
 For the class documentation see the corresponding header file.
 
 Important revisions after version 2.0.0 :
 
*/
//____________________________________________________________________________
 
#include <sstream>
 
#include <TSystem.h>
 
#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/Messenger/Messenger.h"
#include "Physics/NuclearState/SpectralFunc1d.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Numerical/Spline.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Physics/NuclearState/NuclearUtils.h"
 
using std::ostringstream;
 
using namespace genie;
using namespace genie::constants;
using namespace genie::controls;
using namespace genie::utils;
 
//____________________________________________________________________________
SpectralFunc1d::SpectralFunc1d() :
NuclearModelI("genie::SpectralFunc1d")
{
 
}
//____________________________________________________________________________
SpectralFunc1d::SpectralFunc1d(string config) :
NuclearModelI("genie::SpectralFunc1d", config)
{
 
}
//____________________________________________________________________________
SpectralFunc1d::~SpectralFunc1d()
{
  this->CleanUp();
}
//____________________________________________________________________________
bool SpectralFunc1d::GenerateNucleon(const Target & target) const
{
  RandomGen * rnd = RandomGen::Instance();
  int Z = target.Z();
 
  map<int, double>::const_iterator  dbl_it;
  map<int, Spline*>::const_iterator spl_it;
 
  // Select fermi momentum from the integrated (over removal energies) s/f.
  //
  spl_it = fSFk.find(Z);
  dbl_it = fMaxProb.find(Z);
  if(spl_it == fSFk.end() || dbl_it == fMaxProb.end()) {
    fCurrRemovalEnergy = 0.;
    fCurrMomentum.SetXYZ(0.,0.,0.);
    return false;
  }
 
  double prob_max = dbl_it->second;
  LOG("SpectralFunc1", pDEBUG) << "Max probability = " << prob_max;
 
  double p = 0;
  unsigned int niter = 0;
  while(1) {
    if(niter > kRjMaxIterations) {
       LOG("SpectralFunc1", pWARN)
        << "Couldn't generate a hit nucleon after " << niter << " iterations";
       return false;
    }
    niter++;
 
    if(fUseRFGMomentumCutoff) p = fPCutOff * rnd->RndGen().Rndm();
    else p = rnd->RndGen().Rndm();
 
    double prob  = spl_it->second->Evaluate(p);
    double probg = prob_max * rnd->RndGen().Rndm();
    LOG("SpectralFunc1", pDEBUG) << "Trying p = " << p << " / prob = " << prob;
 
    bool accept = (probg<prob);
    if(accept) break;
  }
 
  LOG("SpectralFunc1", pINFO) << "|p,nucleon| = " << p;
 
  double costheta = -1. + 2. * rnd->RndGen().Rndm();
  double sintheta = TMath::Sqrt(1.-costheta*costheta);
  double fi       = 2 * kPi * rnd->RndGen().Rndm();
  double cosfi    = TMath::Cos(fi);
  double sinfi    = TMath::Sin(fi);
 
  double px = p*sintheta*cosfi;
  double py = p*sintheta*sinfi;
  double pz = p*costheta;
 
  fCurrMomentum.SetXYZ(px,py,pz);
 
  // Set removal energy
  // Do it either in the same way as in the FG model or by using the average
  // removal energy for the seleced pF as calculated from the s/f itself
  //
  if(fUseRFGRemovalE) {
    dbl_it = fNucRmvE.find(Z);
    if(dbl_it != fNucRmvE.end()) fCurrRemovalEnergy = dbl_it->second;
    else fCurrRemovalEnergy = nuclear::BindEnergyPerNucleon(target);
  } else {
    spl_it = fSFw.find(Z);
    if(spl_it==fSFw.end()) {
       fCurrRemovalEnergy = 0.;
       fCurrMomentum.SetXYZ(0.,0.,0.);
       return false;
    } else fCurrRemovalEnergy = spl_it->second->Evaluate(p);
  }
 
  return true;
}
//____________________________________________________________________________
double SpectralFunc1d::Prob(
                  double p, double /*w*/, const Target & target) const
{
  if(fUseRFGMomentumCutoff) { if(p > fPCutOff) return 0; }
 
  int Z = target.Z();
  map<int, Spline*>::const_iterator spl_it = fSFk.find(Z);
 
  if(spl_it == fSFk.end()) return 0;
  else return TMath::Max(0.,spl_it->second->Evaluate(p));
}
//____________________________________________________________________________
void SpectralFunc1d::Configure(const Registry & config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void SpectralFunc1d::Configure(string param_set)
{
  Algorithm::Configure(param_set);
  this->LoadConfig();
}
//____________________________________________________________________________
void SpectralFunc1d::LoadConfig(void)
{
  
  NuclearModelI::LoadConfig() ; 
 
  LOG("SpectralFunc1", pDEBUG) << "Loading coonfiguration for SpectralFunc1d";
 
  this->CleanUp();
 
  // Load spectral function data.
  // Hopefully analytical expressions will be available soon.
  // Currently I have spectral functions for C12 and Fe56 only.
  //
  string data_dir =
      string(gSystem->Getenv("GENIE")) + string("/data/evgen/nucl/spectral_functions/");
 
  string c12_sf1dk_file  = data_dir + "benhar-sf1dk-12c.data";
  string fe56_sf1dk_file = data_dir + "benhar-sf1dk-56fe.data";
  string c12_sf1dw_file  = data_dir + "benhar-sf1dw-12c.data";
  string fe56_sf1dw_file = data_dir + "benhar-sf1dw-56fe.data";
 
  Spline * spl = 0;
 
  spl = new Spline(c12_sf1dk_file);
  fSFk.insert(map<int, Spline*>::value_type(6,spl));
  spl = new Spline(fe56_sf1dk_file);
  fSFk.insert(map<int, Spline*>::value_type(26,spl));
 
  spl = new Spline(c12_sf1dw_file);
  fSFw.insert(map<int, Spline*>::value_type(6,spl));
  spl = new Spline(fe56_sf1dw_file);
  fSFw.insert(map<int, Spline*>::value_type(26,spl));
 
  // scan for max.
  map<int, Spline*>::const_iterator spliter;
  int    n    = 100;
  double pmin =  0.000;
  double dp   =  0.010;
  for(spliter = fSFk.begin(); spliter != fSFk.end(); ++spliter) {
    double prob_max = 0;
    int Z = spliter->first;
    spl = spliter->second;
    for(int i=0; i<n; i++) {
       double p = pmin + i*dp;
       prob_max = TMath::Max(prob_max, spl->Evaluate(p));
    }
    fMaxProb.insert(map<int,double>::value_type(Z,prob_max));
  }
 
  // Check whether to use the same removal energies as in the FG model or
  // to use the average removal energy for the selected fermi momentum
  // (computed from the spectral function itself)
  GetParam( "UseRFGRemovalE", fUseRFGRemovalE ) ;
 
  // Check whether to use the same momentum cutoff as in the FG model
  GetParam("UseRFGMomentumCutoff", fUseRFGMomentumCutoff ) ;
 
  //Get the momentum cutoff
  GetParam( "RFG-MomentumCutOff", fPCutOff ) ;
 
  // Removal energies as used in the FG model
  // Load removal energy for specific nuclei from either the algorithm's
  // configuration file or the UserPhysicsOptions file.
  // If none is used use Wapstra's semi-empirical formula.
  for(int Z=1; Z<140; Z++) {
    for(int A=Z; A<3*Z; A++) {
      ostringstream key;
      int pdgc = pdg::IonPdgCode(A,Z);
      key   << "RFG-NucRemovalE@Pdg="     << pdgc;
      RgKey rgkey   = key.str();
      double eb ;
      if ( GetParam( rgkey, eb, false ) ) {
        eb = TMath::Max(eb, 0.);
        LOG("BodekRitchie", pNOTICE)
          << "Nucleus: " << pdgc << " -> using Eb =  " << eb << " GeV";
        fNucRmvE.insert(map<int,double>::value_type(Z,eb));
      }
    }
  }
}
//____________________________________________________________________________
void SpectralFunc1d::CleanUp(void)
{
  map<int, Spline*>::iterator spliter;
 
  for(spliter = fSFk.begin(); spliter != fSFk.end(); ++spliter) {
    Spline * spl = spliter->second;
    if(spl) delete spl;
  }
  for(spliter = fSFw.begin(); spliter != fSFw.end(); ++spliter) {
    Spline * spl = spliter->second;
    if(spl) delete spl;
  }
  fSFk.clear();
  fSFw.clear();
  fNucRmvE.clear();
  fMaxProb.clear();
}
//____________________________________________________________________________