genie::HEDISStrucFunc Class Reference

#include <HEDISStrucFunc.h>

Collaboration diagram for genie::HEDISStrucFunc:

Classes
struct	Cleaner
class	HEDISStrucFuncTable

Public Types
enum	StrucFuncType {   kMHTUndefined = 0 , kSFT1 , kSFT2 , kSFT3 ,   kSFnumber }
typedef enum genie::HEDISStrucFunc::StrucFuncType	HEDISStrucFuncType_t

Public Member Functions
SF_xQ2	EvalQrkSFLO (const Interaction *in, double x, double Q2)
SF_xQ2	EvalNucSFLO (const Interaction *in, double x, double Q2)
SF_xQ2	EvalNucSFNLO (const Interaction *in, double x, double Q2)

Static Public Member Functions
static HEDISStrucFunc *	Instance (SF_info sfinfo)

Private Member Functions
	HEDISStrucFunc (SF_info sfinfo)
	HEDISStrucFunc (const HEDISStrucFunc &)
	~HEDISStrucFunc ()
void	CreateQrkSF (const Interaction *in, string sfFile)
void	CreateNucSF (const Interaction *in, string sfFile)
string	QrkSFName (const Interaction *in)
string	NucSFName (const Interaction *in)
int	QrkSFCode (const Interaction *in)
int	NucSFCode (const Interaction *in)

Private Attributes
map< int, HEDISStrucFuncTable >	fQrkSFLOTables
map< int, HEDISStrucFuncTable >	fNucSFLOTables
map< int, HEDISStrucFuncTable >	fNucSFNLOTables
SF_info	fSF
vector< double >	sf_x_array
vector< double >	sf_q2_array

Static Private Attributes
static HEDISStrucFunc *	fgInstance = 0

Friends
struct	Cleaner

Detailed Description

Definition at line 191 of file HEDISStrucFunc.h.

Member Typedef Documentation

HEDISStrucFuncType_t

typedef enum genie::HEDISStrucFunc::StrucFuncType genie::HEDISStrucFunc::HEDISStrucFuncType_t

Member Enumeration Documentation

StrucFuncType

enum genie::HEDISStrucFunc::StrucFuncType

Enumerator
kMHTUndefined
kSFT1
kSFT2
kSFT3
kSFnumber

Definition at line 199 of file HEDISStrucFunc.h.

                                 {  
        kMHTUndefined = 0,
        kSFT1, 
        kSFT2, 
        kSFT3, 
        kSFnumber, 
      } HEDISStrucFuncType_t;

Constructor & Destructor Documentation

HEDISStrucFunc() [1/2]

HEDISStrucFunc::HEDISStrucFunc ( SF_info sfinfo )

private

Definition at line 47 of file HEDISStrucFunc.cxx.

{
 
  fSF = sfinfo;
 
  string basedir = "";
  if ( gSystem->Getenv("HEDIS_SF_DATA_PATH")==NULL ) basedir = string(gSystem->Getenv("GENIE")) + "/data/evgen/hedis-sf";
  else                                               basedir = string(gSystem->Getenv("HEDIS_SF_DATA_PATH"));
  LOG("HEDISStrucFunc", pWARN) << "Base directory: " << basedir;
 
  if ( gSystem->AccessPathName( basedir.c_str(), kReadPermission ) ) {
      LOG("HEDISStrucFunc", pFATAL) << "Base directory doesnt exist or you dont have read permission.";
      LOG("HEDISStrucFunc", pFATAL) << "Remember!!!";
      LOG("HEDISStrucFunc", pFATAL) << "Path to base directory is defined with the enviroment variable HEDIS_SF_DATA_PATH.";
      LOG("HEDISStrucFunc", pFATAL) << "If not defined, default location is $GENIE/data/evgen/hedis-sf";
      assert(0);
  }
 
  // Name of the directory where SF tables are (or will be) saved. 
  string SFname = basedir + "/" + RunOpt::Instance()->Tune()->Name();
 
  // Check that the directory where SF tables are stored exists
  LOG("HEDISStrucFunc", pWARN) << "SF are (or will be) in following directory: " << SFname;
  if ( gSystem->AccessPathName( SFname.c_str(), kReadPermission ) ) {
    LOG("HEDISStrucFunc", pWARN) << "Bad news! Directory doesnt exists.";
    LOG("HEDISStrucFunc", pWARN) << "HEDIS package requires computation of SF";        
    LOG("HEDISStrucFunc", pWARN) << "This will be SLOW!!!!!";        
    if ( gSystem->mkdir(SFname.c_str())==0 ) {
      LOG("HEDISStrucFunc", pINFO) << "Creating Metafile with input information";
      std::ofstream meta_stream((SFname+"/Inputs.txt").c_str());
      meta_stream << fSF;
      meta_stream.close();      
    }
    else {
      LOG("HEDISStrucFunc", pFATAL) << "You dont have write permission in the following directory:";
      LOG("HEDISStrucFunc", pFATAL) << SFname;
      LOG("HEDISStrucFunc", pFATAL) << "Remember!!!";
      LOG("HEDISStrucFunc", pFATAL) << "Path to base directory is defined with the enviroment variable HEDIS_SF_DATA_PATH.";
      LOG("HEDISStrucFunc", pFATAL) << "If not defined, default location is $GENIE/data/evgen/hedis-sf";
      assert(0);
    }
  }
  else {
    LOG("HEDISStrucFunc", pWARN) << "Good news! Directory already exists.";
    LOG("HEDISStrucFunc", pWARN) << "Precomputed SF stored in that directory will be used";
    LOG("HEDISStrucFunc", pINFO) << "Comparing info from Metafile and selected Tune (CommonParam.xml)";
    SF_info cm;
    std::ifstream meta_stream((SFname+"/Inputs.txt").c_str(), std::ios::in);
    meta_stream >> cm;
    meta_stream.close();
    if (cm==fSF) {
      LOG("HEDISStrucFunc", pINFO) << "Info from MetaFile and Tune match";
    }
    else {
      LOG("HEDISStrucFunc", pFATAL) << "Info from MetaFile and Tune doesnt match";
      LOG("HEDISStrucFunc", pFATAL) << "MetaFile Path : " << SFname << "/Inputs.txt";        
      LOG("HEDISStrucFunc", pFATAL) << "From Tune : ";        
      LOG("HEDISStrucFunc", pFATAL) << cm;        
      LOG("HEDISStrucFunc", pFATAL) << "From MetaFile : ";        
      LOG("HEDISStrucFunc", pFATAL) << fSF;        
      assert(0);
    }
  }
 
  fSF.Sin2ThW = fSF.Rho==0. ? fSF.Sin2ThW : 1.-fSF.MassW*fSF.MassW/fSF.MassZ/fSF.MassZ/(1+fSF.Rho);
  LOG("HEDISStrucFunc", pINFO) << "Weingber angle computation:";
  LOG("HEDISStrucFunc", pINFO) << "Rho     = " << fSF.Rho;
  LOG("HEDISStrucFunc", pINFO) << "Sin2ThW = " << fSF.Sin2ThW;
 
  // initialising lhapdf
#ifdef __GENIE_LHAPDF6_ENABLED__
  LOG("HEDISStrucFunc", pINFO) << "Initialising LHAPDF6...";
  pdf = LHAPDF::mkPDF(fSF.LHAPDFset, fSF.LHAPDFmember);
  xPDFmin  = pdf->xMin();
  Q2PDFmin = pdf->q2Min();
  Q2PDFmax = pdf->q2Max();
  LOG("HEDISStrucFunc", pINFO) << "PDF info:";
  LOG("HEDISStrucFunc", pINFO) << "OrderQCD = " << pdf->orderQCD();
  LOG("HEDISStrucFunc", pINFO) << "FlavorScheme = " << pdf->info().get_entry("FlavorScheme");
  LOG("HEDISStrucFunc", pINFO) << "NumFlavors = " << pdf->info().get_entry("NumFlavors");
  LOG("HEDISStrucFunc", pINFO) << "Xmin = " << xPDFmin << "  Xmax = " << pdf->xMax() << "  Q2min = " << Q2PDFmin << "  Q2max = " << Q2PDFmax;
  LOG("HEDISStrucFunc", pINFO) << "MZ = " << pdf->info().get_entry("MZ");
  for (int i=1; i<7; i++) {
    mPDFQrk[i] = pdf->quarkMass(i);
    LOG("HEDISStrucFunc", pINFO) << "M" << i << " = " << mPDFQrk[i];
  }
#endif
#ifdef __GENIE_LHAPDF5_ENABLED__
  LOG("HEDISStrucFunc", pINFO) << "Initialising LHAPDF5...";
  LHAPDF::initPDFByName(fSF.LHAPDFset, LHAPDF::LHGRID, fSF.LHAPDFmember);
  xPDFmin  = LHAPDF::getXmin(0);
  Q2PDFmin = LHAPDF::getQ2min(0);
  Q2PDFmax = LHAPDF::getQ2max(0);
  LOG("HEDISStrucFunc", pINFO) << "PDF info:";
  LOG("HEDISStrucFunc", pINFO) << "Xmin = " << xPDFmin << "  Xmax = " << LHAPDF::getXmax(0) << "  Q2min = " << Q2PDFmin << "  Q2max = " << Q2PDFmax;
  for (int i=1; i<7; i++) {
    mPDFQrk[i] = LHAPDF::getQMass(i);
    LOG("HEDISStrucFunc", pINFO) << "M" << i << " = " << mPDFQrk[i];
  }
#endif
 
  // checking if LHAPDF boundaries matches boundaries defined by user
  if ( fSF.XGridMin < xPDFmin ) {
    LOG("HEDISStrucFunc", pWARN) << "Lower boundary in X is smaller than input PDF";
    LOG("HEDISStrucFunc", pWARN) << "xPDFmin = "  << xPDFmin;
    LOG("HEDISStrucFunc", pWARN) << "xGridMin = " << fSF.XGridMin;
    LOG("HEDISStrucFunc", pWARN) << "DANGER: An extrapolation will be done!!!!!";
  }
  else if ( fSF.Q2GridMin < Q2PDFmin ) {
    LOG("HEDISStrucFunc", pWARN) << "Lower boundary in Q2 is smaller than input PDF";
    LOG("HEDISStrucFunc", pWARN) << "Q2PDFmin = "  << Q2PDFmin;
    LOG("HEDISStrucFunc", pWARN) << "Q2GridMin = " << fSF.Q2GridMin;
    LOG("HEDISStrucFunc", pWARN) << "DANGER: An extrapolation will be done!!!!!";
  }
  else if ( fSF.Q2GridMax > Q2PDFmax ) {
    LOG("HEDISStrucFunc", pWARN) << "Upper boundary in Q2 is bigger than input PDF";
    LOG("HEDISStrucFunc", pWARN) << "Q2PDFmax = "  << Q2PDFmax;
    LOG("HEDISStrucFunc", pWARN) << "Q2GridMax = " << fSF.Q2GridMax;
    LOG("HEDISStrucFunc", pWARN) << "DANGER: An extrapolation will be done!!!!!";
  }
 
  // define arrays to fill from data files
  double dlogq2 = TMath::Abs( TMath::Log10(fSF.Q2GridMin)-TMath::Log10(fSF.Q2GridMax) ) / fSF.NGridQ2;
  double dlogx  = TMath::Abs( TMath::Log10(fSF.XGridMin)-TMath::Log10(1.) ) / fSF.NGridX;
  LOG("HEDISStrucFunc", pINFO) << "Grid x,Q2 :" << fSF.NGridX << " , " << fSF.NGridQ2;
  for ( double logq2 = TMath::Log10(fSF.Q2GridMin); logq2<TMath::Log10(fSF.Q2GridMax); logq2+= dlogq2 ) {
    double q2 = TMath::Power( 10, logq2 + 0.5*dlogq2 );
    if (q2>fSF.Q2GridMax) continue;
    sf_q2_array.push_back(q2);
    LOG("HEDISStrucFunc", pDEBUG) << "q2: " << sf_q2_array.back();
  }
  for ( double logx = TMath::Log10(fSF.XGridMin); logx<TMath::Log10(1.); logx+= dlogx ) {
    double x = TMath::Power( 10, logx + 0.5*dlogx );
    if ( x>1. ) continue;
    sf_x_array.push_back(x);
    LOG("HEDISStrucFunc", pDEBUG) << "x: " << sf_x_array.back();
  }
 
  // Change to variables that are suitable for BLI2DNonUnifGrid
  int nx = sf_q2_array.size();
  int ny = sf_x_array.size();
  double x[nx];
  double y[ny];
  double z[nx*ny];
  for (int i=0; i<nx; i++) x[i] = sf_q2_array[i];
  for (int j=0; j<ny; j++) y[j] = sf_x_array[j];
  
  vector<InteractionType_t> inttype;
  inttype.push_back(kIntWeakCC);
  inttype.push_back(kIntWeakNC);
  vector<InitialState> init_state;
  init_state.push_back(InitialState(1, 2, kPdgNuE));
  init_state.push_back(InitialState(1, 2, kPdgAntiNuE));
 
  HEDISInteractionListGenerator * helist = new HEDISInteractionListGenerator();
  InteractionList * ilist = helist->CreateHEDISlist(init_state,inttype);
 
  // Load structure functions for each quark at LO
  for(InteractionList::iterator in=ilist->begin(); in!=ilist->end(); ++in) {
 
    string sfFile = SFname + "/QrkSF_LO_" + QrkSFName(*in) + ".dat";
    // Make sure data files are available
    LOG("HEDISStrucFunc", pINFO) << "Checking if file " << sfFile << " exists...";        
    if ( gSystem->AccessPathName( sfFile.c_str()) ) {
      LOG("HEDISStrucFunc", pWARN) << "File doesnt exist. SF table will be computed.";        
      CreateQrkSF( *in, sfFile );
    }
    else if ( atoi(gSystem->GetFromPipe(("wc -w "+sfFile+" | awk '{print $1}'").c_str()))!=kSFT3*nx*ny ) {
      LOG("HEDISStrucFunc", pWARN) << "File does not contain all the need points. SF table will be recomputed.";        
      gSystem->Exec(("rm "+sfFile).c_str());
      CreateQrkSF( *in, sfFile );
    }
    std::ifstream sf_stream(sfFile.c_str(), std::ios::in);
    // Loop over F1,F2,F3
    for(int sf = 1; sf < kSFnumber; ++sf) {
      // Loop over x/Q2 bins
      for ( int ij=0; ij<nx*ny; ij++ ) sf_stream >> z[ij];
      // Create SF tables with BLI2DNonUnifGrid using x,Q2 binning
      fQrkSFLOTables[QrkSFCode(*in)].Table[(HEDISStrucFuncType_t)sf] = new genie::BLI2DNonUnifGrid( nx, ny, x, y, z );
    }
  }
 
  if (fSF.IsNLO) {
#ifdef __GENIE_APFEL_ENABLED__
    // initialising APFEL framework
    LOG("HEDISStrucFunc", pINFO) << "Initialising APFEL..." ; 
    APFEL::SetPDFSet(fSF.LHAPDFset);
    APFEL::SetReplica(fSF.LHAPDFmember);
    if (fSF.Scheme=="BGR") {
      APFEL::SetMassScheme("FONLL-B");
      APFEL::SetPoleMasses(mPDFQrk[4],mPDFQrk[5],mPDFQrk[6]);
    } 
    else if (fSF.Scheme=="CSMS") {
      APFEL::SetMassScheme("ZM-VFNS");
      APFEL::SetPoleMasses(mPDFQrk[4],mPDFQrk[5],mPDFQrk[5]+0.1);
    }
    else if (fSF.Scheme=="GGHR") {
      APFEL::SetMassScheme("FFNS5");
      APFEL::SetPoleMasses(mPDFQrk[4],mPDFQrk[5],mPDFQrk[6]);
      APFEL::SetMaxFlavourPDFs(5);
      APFEL::SetMaxFlavourAlpha(5);
    }
    else {
      LOG("HEDISStrucFunc", pERROR) << "Mass Scheme is not set properly";
      assert(0);
    }
    APFEL::SetQLimits(TMath::Sqrt(fSF.Q2GridMin),TMath::Sqrt(fSF.Q2GridMax));
    APFEL::SetMaxFlavourPDFs(6);
    APFEL::SetMaxFlavourAlpha(6);
    APFEL::SetNumberOfGrids(3);
    APFEL::SetGridParameters(1,90,3,fSF.XGridMin);
    APFEL::SetGridParameters(2,50,5,1e-1);
    APFEL::SetGridParameters(3,40,5,8e-1);
    APFEL::SetPerturbativeOrder(1);
    APFEL::SetAlphaQCDRef(pdf->alphasQ(fSF.MassZ),fSF.MassZ);    
    APFEL::SetProtonMass(kProtonMass);
    APFEL::SetWMass(fSF.MassW);
    APFEL::SetZMass(fSF.MassZ);
    APFEL::SetSin2ThetaW(fSF.Sin2ThW);
    APFEL::SetCKM(fSF.Vud, fSF.Vus, fSF.Vub,
                  fSF.Vcd, fSF.Vcs, fSF.Vcb,
                  fSF.Vtd, fSF.Vts, fSF.Vtb);
#endif
 
    // Load structure functions for each quark at LO
    int nch = -1;
    for(InteractionList::iterator in=ilist->begin(); in!=ilist->end(); ++in) {
 
      if ( nch==NucSFCode(*in) ) continue;
      nch = NucSFCode(*in);
      
      string sfFile = SFname + "/NucSF_NLO_" + NucSFName(*in) + ".dat";
      // Make sure data files are available
      LOG("HEDISStrucFunc", pINFO) << "Checking if file " << sfFile << " exists...";        
      if ( gSystem->AccessPathName( sfFile.c_str()) ) {
#ifdef __GENIE_APFEL_ENABLED__
        LOG("HEDISStrucFunc", pWARN) << "File doesnt exist. SF table will be computed.";        
        CreateNucSF( *in, sfFile );
#else
        LOG("HEDISStrucFunc", pERROR) << "File doesnt exist. APFEL is needed for NLO SF";        
        assert(0);
#endif
      }
      else if ( atoi(gSystem->GetFromPipe(("wc -w "+sfFile+" | awk '{print $1}'").c_str()))!=kSFT3*nx*ny ) {
#ifdef __GENIE_APFEL_ENABLED__
        LOG("HEDISStrucFunc", pWARN) << "File does not contain all the need points. SF table will be recomputed.";        
        gSystem->Exec(("rm "+sfFile).c_str());
        CreateNucSF( *in, sfFile );
#else
        LOG("HEDISStrucFunc", pERROR) << "File does not contain all the need points. APFEL is needed for NLO SF";        
        assert(0);
#endif
      }
      std::ifstream sf_stream(sfFile.c_str(), std::ios::in);
      // Loop over F1,F2,F3
      for(int sf = 1; sf < kSFnumber; ++sf) {
        // Loop over x/Q2 bins
        for ( int ij=0; ij<nx*ny; ij++ ) sf_stream >> z[ij];
        // Create SF tables with BLI2DNonUnifGrid using x,Q2 binning
        fNucSFNLOTables[nch].Table[(HEDISStrucFuncType_t)sf] = new genie::BLI2DNonUnifGrid( nx, ny, x, y, z );
      }        
 
      //compute structure functions for each nucleon at LO using quark grids
      LOG("HEDISStrucFunc", pDEBUG) << "Creating LO " << sfFile;              
      vector <int> qcodes;
      for(InteractionList::iterator in2=ilist->begin(); in2!=ilist->end(); ++in2) {
        if (NucSFCode(*in2)==nch) qcodes.push_back(QrkSFCode(*in2));
      }
      // Loop over F1,F2,F3
      for(int sf = 1; sf < kSFnumber; ++sf) {
        int ij = 0;
        // Loop over Q2 bins
        for (int i=0; i<nx; i++) {
          // Loop over x bins
          for (int j=0; j<ny; j++) {
            double sum = 0.;
            // NucSF = sum_qrks QrkSF
            for(vector<int>::const_iterator iq=qcodes.begin(); iq!=qcodes.end(); ++iq) 
              sum += fQrkSFLOTables[*iq].Table[(HEDISStrucFuncType_t)sf]->Evaluate(x[i],y[j]);
            z[ij] = sum;
            ij++;
          }
        }
        // Create SF tables with BLI2DNonUnifGrid using x,Q2 binning
        fNucSFLOTables[nch].Table[(HEDISStrucFuncType_t)sf] = new genie::BLI2DNonUnifGrid( nx, ny, x, y, z );
      }
    }
  }
  
  fgInstance = 0;
 
}

References genie::HEDISInteractionListGenerator::CreateHEDISlist(), CreateNucSF(), CreateQrkSF(), e, fgInstance, fNucSFLOTables, fNucSFNLOTables, fQrkSFLOTables, fSF, genie::RunOpt::Instance(), genie::kIntWeakCC, genie::kIntWeakNC, genie::kPdgAntiNuE, genie::kPdgNuE, genie::constants::kProtonMass, kSFnumber, kSFT3, LOG, mPDFQrk, genie::TuneId::Name(), NucSFCode(), NucSFName(), pDEBUG, pERROR, pFATAL, pINFO, pWARN, Q2PDFmax, Q2PDFmin, QrkSFCode(), QrkSFName(), sf_q2_array, sf_x_array, genie::RunOpt::Tune(), and xPDFmin.

Referenced by HEDISStrucFunc(), and Instance().

HEDISStrucFunc() [2/2]

genie::HEDISStrucFunc::HEDISStrucFunc ( const HEDISStrucFunc & )

private

References HEDISStrucFunc().

~HEDISStrucFunc()

HEDISStrucFunc::~HEDISStrucFunc ( )

private

Definition at line 340 of file HEDISStrucFunc.cxx.

{
 
}

Member Function Documentation

CreateNucSF()

void genie::HEDISStrucFunc::CreateNucSF ( const Interaction * in, string sfFile )

private

Referenced by HEDISStrucFunc().

CreateQrkSF()

void HEDISStrucFunc::CreateQrkSF ( const Interaction * in, string sfFile )

private

Definition at line 358 of file HEDISStrucFunc.cxx.

{
 
  // variables used to tag the SF for particular channel
  bool iscc        = in->ProcInfo().IsWeakCC();
  bool isnu        = pdg::IsNeutrino(in->InitState().ProbePdg());
  bool ispr        = pdg::IsProton(in->InitState().Tgt().HitNucPdg());
  bool sea_iq      = in->InitState().Tgt().HitSeaQrk();
  int pdg_iq       = in->InitState().Tgt().HitQrkPdg();
  int pdg_fq       = in->ExclTag().FinalQuarkPdg();
  double mass_nucl = in->InitState().Tgt().HitNucMass();
 
  // up and down quark swicth depending on proton or neutron interaction  
  int qrkd = 0;
  int qrku = 0;
  if ( ispr ) { qrkd = 1 ; qrku = 2; }
  else        { qrkd = 2 ; qrku = 1; }
 
  // variables associated to the PDF and coupling of the quarks
  int qpdf1 = -999;                                     // number used to compute PDF
  int qpdf2 = -999;                                     // auxiliary number used to compute PDF for valence quarks
  double Cp2 = -999;                                    // couping for F1,F2
  double Cp3 = -999;                                    // couping for F3
  double sign3 = isnu ? +1. : -1.;  // sign change for nu/nubar in F3
  if ( iscc ) {
    if ( isnu ) {
      if      ( pdg_iq== 1 && !sea_iq && pdg_fq== 2 ) { qpdf1 =  qrkd; qpdf2 = -qrkd; Cp2 = 2*fSF.Vud*fSF.Vud; Cp3 =  2*fSF.Vud*fSF.Vud; }
      else if ( pdg_iq== 1 && !sea_iq && pdg_fq== 4 ) { qpdf1 =  qrkd; qpdf2 = -qrkd; Cp2 = 2*fSF.Vcd*fSF.Vcd; Cp3 =  2*fSF.Vcd*fSF.Vcd; }
      else if ( pdg_iq== 1 && !sea_iq && pdg_fq== 6 ) { qpdf1 =  qrkd; qpdf2 = -qrkd; Cp2 = 2*fSF.Vtd*fSF.Vtd; Cp3 =  2*fSF.Vtd*fSF.Vtd; }
      else if ( pdg_iq== 1 &&  sea_iq && pdg_fq== 2 ) { qpdf1 = -qrkd;                Cp2 = 2*fSF.Vud*fSF.Vud; Cp3 =  2*fSF.Vud*fSF.Vud; }
      else if ( pdg_iq== 1 &&  sea_iq && pdg_fq== 4 ) { qpdf1 = -qrkd;                Cp2 = 2*fSF.Vcd*fSF.Vcd; Cp3 =  2*fSF.Vcd*fSF.Vcd; }
      else if ( pdg_iq== 1 &&  sea_iq && pdg_fq== 6 ) { qpdf1 = -qrkd;                Cp2 = 2*fSF.Vtd*fSF.Vtd; Cp3 =  2*fSF.Vtd*fSF.Vtd; }
      else if ( pdg_iq== 3 &&  sea_iq && pdg_fq== 2 ) { qpdf1 =  3;                   Cp2 = 2*fSF.Vus*fSF.Vus; Cp3 =  2*fSF.Vus*fSF.Vus; }
      else if ( pdg_iq== 3 &&  sea_iq && pdg_fq== 4 ) { qpdf1 =  3;                   Cp2 = 2*fSF.Vcs*fSF.Vcs; Cp3 =  2*fSF.Vcs*fSF.Vcs; }
      else if ( pdg_iq== 3 &&  sea_iq && pdg_fq== 6 ) { qpdf1 =  3;                   Cp2 = 2*fSF.Vts*fSF.Vts; Cp3 =  2*fSF.Vts*fSF.Vts; }
      else if ( pdg_iq== 5 &&  sea_iq && pdg_fq== 2 ) { qpdf1 =  5;                   Cp2 = 2*fSF.Vub*fSF.Vub; Cp3 =  2*fSF.Vub*fSF.Vub; }
      else if ( pdg_iq== 5 &&  sea_iq && pdg_fq== 4 ) { qpdf1 =  5;                   Cp2 = 2*fSF.Vcb*fSF.Vcb; Cp3 =  2*fSF.Vcb*fSF.Vcb; }
      else if ( pdg_iq== 5 &&  sea_iq && pdg_fq== 6 ) { qpdf1 =  5;                   Cp2 = 2*fSF.Vtb*fSF.Vtb; Cp3 =  2*fSF.Vtb*fSF.Vtb; }
      else if ( pdg_iq==-2 &&  sea_iq && pdg_fq==-1 ) { qpdf1 = -qrku;                Cp2 = 2*fSF.Vud*fSF.Vud; Cp3 = -2*fSF.Vud*fSF.Vud; }
      else if ( pdg_iq==-2 &&  sea_iq && pdg_fq==-3 ) { qpdf1 = -qrku;                Cp2 = 2*fSF.Vus*fSF.Vus; Cp3 = -2*fSF.Vus*fSF.Vus; }
      else if ( pdg_iq==-2 &&  sea_iq && pdg_fq==-5 ) { qpdf1 = -qrku;                Cp2 = 2*fSF.Vub*fSF.Vub; Cp3 = -2*fSF.Vub*fSF.Vub; }
      else if ( pdg_iq==-4 &&  sea_iq && pdg_fq==-1 ) { qpdf1 = -4;                   Cp2 = 2*fSF.Vcd*fSF.Vcd; Cp3 = -2*fSF.Vcd*fSF.Vcd; }
      else if ( pdg_iq==-4 &&  sea_iq && pdg_fq==-3 ) { qpdf1 = -4;                   Cp2 = 2*fSF.Vcs*fSF.Vcs; Cp3 = -2*fSF.Vcs*fSF.Vcs; }
      else if ( pdg_iq==-4 &&  sea_iq && pdg_fq==-5 ) { qpdf1 = -4;                   Cp2 = 2*fSF.Vcb*fSF.Vcb; Cp3 = -2*fSF.Vcb*fSF.Vcb; }
    }
    else {
      if      ( pdg_iq== 2 && !sea_iq && pdg_fq== 1 ) { qpdf1 =  qrku; qpdf2 = -qrku; Cp2 = 2*fSF.Vud*fSF.Vud; Cp3 =  2*fSF.Vud*fSF.Vud; }
      else if ( pdg_iq== 2 && !sea_iq && pdg_fq== 3 ) { qpdf1 =  qrku; qpdf2 = -qrku; Cp2 = 2*fSF.Vus*fSF.Vus; Cp3 =  2*fSF.Vus*fSF.Vus; }
      else if ( pdg_iq== 2 && !sea_iq && pdg_fq== 5 ) { qpdf1 =  qrku; qpdf2 = -qrku; Cp2 = 2*fSF.Vub*fSF.Vub; Cp3 =  2*fSF.Vub*fSF.Vub; }
      else if ( pdg_iq== 2 &&  sea_iq && pdg_fq== 1 ) { qpdf1 = -qrku;                Cp2 = 2*fSF.Vud*fSF.Vud; Cp3 =  2*fSF.Vud*fSF.Vud; }
      else if ( pdg_iq== 2 &&  sea_iq && pdg_fq== 3 ) { qpdf1 = -qrku;                Cp2 = 2*fSF.Vus*fSF.Vus; Cp3 =  2*fSF.Vus*fSF.Vus; }
      else if ( pdg_iq== 2 &&  sea_iq && pdg_fq== 5 ) { qpdf1 = -qrku;                Cp2 = 2*fSF.Vub*fSF.Vub; Cp3 =  2*fSF.Vub*fSF.Vub; }
      else if ( pdg_iq== 4 &&  sea_iq && pdg_fq== 1 ) { qpdf1 =  4;                   Cp2 = 2*fSF.Vcd*fSF.Vcd; Cp3 =  2*fSF.Vcd*fSF.Vcd; }
      else if ( pdg_iq== 4 &&  sea_iq && pdg_fq== 3 ) { qpdf1 =  4;                   Cp2 = 2*fSF.Vcs*fSF.Vcs; Cp3 =  2*fSF.Vcs*fSF.Vcs; }
      else if ( pdg_iq== 4 &&  sea_iq && pdg_fq== 5 ) { qpdf1 =  4;                   Cp2 = 2*fSF.Vcb*fSF.Vcb; Cp3 =  2*fSF.Vcb*fSF.Vcb; }
      else if ( pdg_iq==-1 &&  sea_iq && pdg_fq==-2 ) { qpdf1 = -qrkd;                Cp2 = 2*fSF.Vud*fSF.Vud; Cp3 = -2*fSF.Vud*fSF.Vud; }
      else if ( pdg_iq==-1 &&  sea_iq && pdg_fq==-4 ) { qpdf1 = -qrkd;                Cp2 = 2*fSF.Vcd*fSF.Vcd; Cp3 = -2*fSF.Vcd*fSF.Vcd; }
      else if ( pdg_iq==-1 &&  sea_iq && pdg_fq==-6 ) { qpdf1 = -qrkd;                Cp2 = 2*fSF.Vtd*fSF.Vtd; Cp3 = -2*fSF.Vtd*fSF.Vtd; }
      else if ( pdg_iq==-3 &&  sea_iq && pdg_fq==-2 ) { qpdf1 = -3;                   Cp2 = 2*fSF.Vus*fSF.Vus; Cp3 = -2*fSF.Vus*fSF.Vus; }
      else if ( pdg_iq==-3 &&  sea_iq && pdg_fq==-4 ) { qpdf1 = -3;                   Cp2 = 2*fSF.Vcs*fSF.Vcs; Cp3 = -2*fSF.Vcs*fSF.Vcs; }
      else if ( pdg_iq==-3 &&  sea_iq && pdg_fq==-6 ) { qpdf1 = -3;                   Cp2 = 2*fSF.Vts*fSF.Vts; Cp3 = -2*fSF.Vts*fSF.Vts; }
      else if ( pdg_iq==-5 &&  sea_iq && pdg_fq==-2 ) { qpdf1 = -5;                   Cp2 = 2*fSF.Vub*fSF.Vub; Cp3 = -2*fSF.Vub*fSF.Vub; }
      else if ( pdg_iq==-5 &&  sea_iq && pdg_fq==-4 ) { qpdf1 = -5;                   Cp2 = 2*fSF.Vcb*fSF.Vcb; Cp3 = -2*fSF.Vcb*fSF.Vcb; }
      else if ( pdg_iq==-5 &&  sea_iq && pdg_fq==-6 ) { qpdf1 = -5;                   Cp2 = 2*fSF.Vtb*fSF.Vtb; Cp3 = -2*fSF.Vtb*fSF.Vtb; }
    }
  }
  else {           
    double c2u = TMath::Power( 1./2. - 4./3.*fSF.Sin2ThW,2) + 1./4.;
    double c2d = TMath::Power(-1./2. + 2./3.*fSF.Sin2ThW,2) + 1./4.;
    double c3u = 1./2. - 4./3.*fSF.Sin2ThW;
    double c3d = 1./2. - 2./3.*fSF.Sin2ThW;
    if      ( pdg_iq== 1 && !sea_iq && pdg_fq== 1 ) { qpdf1 =  qrkd; qpdf2 = -qrkd; Cp2 = c2d; Cp3 =  c3d; }
    else if ( pdg_iq== 2 && !sea_iq && pdg_fq== 2 ) { qpdf1 =  qrku; qpdf2 = -qrku; Cp2 = c2u; Cp3 =  c3u; }
    else if ( pdg_iq== 1 &&  sea_iq && pdg_fq== 1 ) { qpdf1 = -qrkd;                Cp2 = c2d; Cp3 =  c3d; }
    else if ( pdg_iq== 2 &&  sea_iq && pdg_fq== 2 ) { qpdf1 = -qrku;                Cp2 = c2u; Cp3 =  c3u; }
    else if ( pdg_iq== 3 &&  sea_iq && pdg_fq== 3 ) { qpdf1 =  3;                   Cp2 = c2d; Cp3 =  c3d; }
    else if ( pdg_iq== 4 &&  sea_iq && pdg_fq== 4 ) { qpdf1 =  4;                   Cp2 = c2u; Cp3 =  c3u; }
    else if ( pdg_iq== 5 &&  sea_iq && pdg_fq== 5 ) { qpdf1 =  5;                   Cp2 = c2d; Cp3 =  c3d; }
    else if ( pdg_iq==-1 &&  sea_iq && pdg_fq==-1 ) { qpdf1 = -qrkd;                Cp2 = c2d; Cp3 = -c3d; }
    else if ( pdg_iq==-2 &&  sea_iq && pdg_fq==-2 ) { qpdf1 = -qrku;                Cp2 = c2u; Cp3 = -c3u; }
    else if ( pdg_iq==-3 &&  sea_iq && pdg_fq==-3 ) { qpdf1 = -3;                   Cp2 = c2d; Cp3 = -c3d; }
    else if ( pdg_iq==-4 &&  sea_iq && pdg_fq==-4 ) { qpdf1 = -4;                   Cp2 = c2u; Cp3 = -c3u; }
    else if ( pdg_iq==-5 &&  sea_iq && pdg_fq==-5 ) { qpdf1 = -5;                   Cp2 = c2d; Cp3 = -c3d; }
  }   
 
  // open file in which SF will be stored
  std::ofstream sf_stream(sfFile.c_str());
 
  // loop over 3 different SF: F1,F2,F3
  for(int sf = 1; sf < 4; sf++) {
    for ( unsigned int i=0; i<sf_q2_array.size(); i++ ) {
      double Q2 = sf_q2_array[i];
      for ( unsigned int j=0; j<sf_x_array.size(); j++ ) {
        double x = sf_x_array[j];
 
        double z = x; // this variable is introduce in case you want to apply scaling
 
        // W threshold
        if      (fSF.QrkThrs==1) { 
            if ( Q2*(1/z-1)+mass_nucl*mass_nucl <= TMath::Power(mass_nucl+mPDFQrk[TMath::Abs(pdg_fq)],2) ) { sf_stream << 0. << "  "; continue; }
        } 
        // W threshold and slow rescaling
        else if (fSF.QrkThrs==2) {
            if ( Q2*(1/z-1)+mass_nucl*mass_nucl <= TMath::Power(mass_nucl+mPDFQrk[TMath::Abs(pdg_fq)],2) ) { sf_stream << 0. << "  "; continue; }
            z *= 1+mPDFQrk[TMath::Abs(pdg_fq)]*mPDFQrk[TMath::Abs(pdg_fq)]/Q2;
        }
        // Slow rescaling
        else if (fSF.QrkThrs==3) {
            z *= 1+mPDFQrk[TMath::Abs(pdg_fq)]*mPDFQrk[TMath::Abs(pdg_fq)]/Q2;
        }
 
        // Fill x,Q2 used to extract PDF. If values outside boundaries then freeze them.
        double xPDF = TMath::Max( z, xPDFmin );
        double Q2PDF = TMath::Max( Q2, Q2PDFmin );
        Q2PDF = TMath::Min( Q2PDF, Q2PDFmax  );
 
        // Extract PDF requiring then to be higher than zero
#ifdef __GENIE_LHAPDF6_ENABLED__
        double fPDF = fmax( pdf->xfxQ2(qpdf1, xPDF, Q2PDF)/z , 0.);
        if (qpdf2!= -999) fPDF -= fmax( pdf->xfxQ2(qpdf2, xPDF, Q2PDF)/z , 0.);
#endif
#ifdef __GENIE_LHAPDF5_ENABLED__
        double fPDF = fmax( LHAPDF::xfx(xPDF, TMath::Sqrt(Q2PDF), qpdf1)/z , 0.);
        if (qpdf2!= -999) fPDF -= fmax( LHAPDF::xfx(xPDF, TMath::Sqrt(Q2PDF), qpdf2)/z , 0.);
#endif
 
        // Compute SF
        double tmp = -999;
        if      ( sf==1 ) tmp = fPDF*Cp2/2;
        else if ( sf==2 ) tmp = fPDF*Cp2*z;
        else if ( sf==3 ) tmp = fPDF*Cp3*sign3;
 
        // Save SF for particular x and Q2 in file
        LOG("HEDISStrucFunc", pDEBUG) << "QrkSFLO" << sf << "[x=" << x << "," << Q2 << "] = " << tmp;
        sf_stream << tmp << "  ";
        
      }
    }
  }
 
  // Close file in which SF are stored
  sf_stream.close();
 
}

References genie::Interaction::ExclTag(), genie::XclsTag::FinalQuarkPdg(), fSF, genie::Target::HitNucMass(), genie::Target::HitNucPdg(), genie::Target::HitQrkPdg(), genie::Target::HitSeaQrk(), genie::Interaction::InitState(), genie::pdg::IsNeutrino(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), LOG, mPDFQrk, pDEBUG, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), Q2PDFmax, Q2PDFmin, sf_q2_array, sf_x_array, genie::InitialState::Tgt(), and xPDFmin.

Referenced by HEDISStrucFunc().

EvalNucSFLO()

SF_xQ2 HEDISStrucFunc::EvalNucSFLO	(	const Interaction *	in,
		double	x,
		double	Q2 )

Definition at line 621 of file HEDISStrucFunc.cxx.

{
  int code = NucSFCode(in);
  SF_xQ2 sf;
  sf.F1 = fNucSFLOTables[code].Table[kSFT1]->Evaluate(Q2,x);
  sf.F2 = fNucSFLOTables[code].Table[kSFT2]->Evaluate(Q2,x);
  sf.F3 = fNucSFLOTables[code].Table[kSFT3]->Evaluate(Q2,x);
  return sf;
}

References genie::SF_xQ2::F1, genie::SF_xQ2::F2, genie::SF_xQ2::F3, fNucSFLOTables, kSFT1, kSFT2, kSFT3, and NucSFCode().

Referenced by genie::HEDISPXSec::XSec().

EvalNucSFNLO()

SF_xQ2 HEDISStrucFunc::EvalNucSFNLO	(	const Interaction *	in,
		double	x,
		double	Q2 )

Definition at line 631 of file HEDISStrucFunc.cxx.

{
  int code = NucSFCode(in);
  SF_xQ2 sf;
  sf.F1 = fNucSFNLOTables[code].Table[kSFT1]->Evaluate(Q2,x);
  sf.F2 = fNucSFNLOTables[code].Table[kSFT2]->Evaluate(Q2,x);
  sf.F3 = fNucSFNLOTables[code].Table[kSFT3]->Evaluate(Q2,x);
  return sf;
}

References genie::SF_xQ2::F1, genie::SF_xQ2::F2, genie::SF_xQ2::F3, fNucSFNLOTables, kSFT1, kSFT2, kSFT3, and NucSFCode().

Referenced by genie::HEDISPXSec::XSec().

EvalQrkSFLO()

SF_xQ2 HEDISStrucFunc::EvalQrkSFLO	(	const Interaction *	in,
		double	x,
		double	Q2 )

Definition at line 611 of file HEDISStrucFunc.cxx.

{
  int code = QrkSFCode(in);
  SF_xQ2 sf;
  sf.F1 = fQrkSFLOTables[code].Table[kSFT1]->Evaluate(Q2,x);
  sf.F2 = fQrkSFLOTables[code].Table[kSFT2]->Evaluate(Q2,x);
  sf.F3 = fQrkSFLOTables[code].Table[kSFT3]->Evaluate(Q2,x);
  return sf;
}

References genie::SF_xQ2::F1, genie::SF_xQ2::F2, genie::SF_xQ2::F3, fQrkSFLOTables, kSFT1, kSFT2, kSFT3, and QrkSFCode().

Referenced by genie::HEDISPXSec::XSec().

Instance()

HEDISStrucFunc * HEDISStrucFunc::Instance ( SF_info sfinfo )

static

Definition at line 345 of file HEDISStrucFunc.cxx.

{
  if(fgInstance == 0) {
    LOG("HEDISStrucFunc", pINFO) << "Late initialization";
    static HEDISStrucFunc::Cleaner cleaner;
    cleaner.DummyMethodAndSilentCompiler();
    fgInstance = new HEDISStrucFunc(sfinfo);
  }  
  return fgInstance;
}

References genie::HEDISStrucFunc::Cleaner::DummyMethodAndSilentCompiler(), fgInstance, HEDISStrucFunc(), LOG, and pINFO.

Referenced by genie::HEDISPXSec::XSec().

NucSFCode()

int HEDISStrucFunc::NucSFCode ( const Interaction * in )

private

Definition at line 603 of file HEDISStrucFunc.cxx.

{
  int code = 100*pdg::IsNeutrino(in->InitState().ProbePdg());
  code    += 10*in->ProcInfo().IsWeakCC();
  code    += 1*pdg::IsProton(in->InitState().Tgt().HitNucPdg());
  return code;
}

References genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsNeutrino(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by EvalNucSFLO(), EvalNucSFNLO(), and HEDISStrucFunc().

NucSFName()

string HEDISStrucFunc::NucSFName ( const Interaction * in )

private

Definition at line 584 of file HEDISStrucFunc.cxx.

{
  string sin = pdg::IsNeutrino(in->InitState().ProbePdg()) ? "nu_" : "nubar_";
  sin += in->ProcInfo().IsWeakCC() ? "cc_" : "nc_";
  sin += pdg::IsProton(in->InitState().Tgt().HitNucPdg()) ? "p" : "n";
  return sin;
}

References genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsNeutrino(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by HEDISStrucFunc().

QrkSFCode()

int HEDISStrucFunc::QrkSFCode ( const Interaction * in )

private

Definition at line 592 of file HEDISStrucFunc.cxx.

{
  int code = 10000000*pdg::IsNeutrino(in->InitState().ProbePdg());
  code    += 1000000*in->ProcInfo().IsWeakCC();
  code    += 100000*pdg::IsProton(in->InitState().Tgt().HitNucPdg());
  code    += 10000*in->InitState().Tgt().HitSeaQrk();
  code    += 100*(6+in->InitState().Tgt().HitQrkPdg());
  code    += 1*(6+in->ExclTag().FinalQuarkPdg());
  return code;
}

References genie::Interaction::ExclTag(), genie::Target::HitNucPdg(), genie::Target::HitSeaQrk(), genie::Interaction::InitState(), genie::pdg::IsNeutrino(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by EvalQrkSFLO(), and HEDISStrucFunc().

QrkSFName()

string HEDISStrucFunc::QrkSFName ( const Interaction * in )

private

Definition at line 573 of file HEDISStrucFunc.cxx.

{
  string sin = pdg::IsNeutrino(in->InitState().ProbePdg()) ? "nu_" : "nubar_";
  sin += in->ProcInfo().IsWeakCC() ? "cc_" : "nc_";
  sin += pdg::IsProton(in->InitState().Tgt().HitNucPdg()) ? "p_" : "n_";
  sin += "iq"+to_string(in->InitState().Tgt().HitQrkPdg());
  sin += in->InitState().Tgt().HitSeaQrk() ? "sea_" : "val_";
  sin += "fq"+to_string(in->ExclTag().FinalQuarkPdg());
  return sin;
}

References genie::Interaction::ExclTag(), genie::XclsTag::FinalQuarkPdg(), genie::Target::HitNucPdg(), genie::Target::HitQrkPdg(), genie::Target::HitSeaQrk(), genie::Interaction::InitState(), genie::pdg::IsNeutrino(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by HEDISStrucFunc().

Friends And Related Symbol Documentation

Cleaner

friend struct Cleaner

friend

Definition at line 265 of file HEDISStrucFunc.h.

Member Data Documentation

fgInstance

HEDISStrucFunc * HEDISStrucFunc::fgInstance = 0

staticprivate

Definition at line 244 of file HEDISStrucFunc.h.

Referenced by HEDISStrucFunc(), Instance(), and genie::HEDISStrucFunc::Cleaner::~Cleaner().

fNucSFLOTables

map<int, HEDISStrucFuncTable> genie::HEDISStrucFunc::fNucSFLOTables

private

Definition at line 248 of file HEDISStrucFunc.h.

Referenced by EvalNucSFLO(), and HEDISStrucFunc().

fNucSFNLOTables

map<int, HEDISStrucFuncTable> genie::HEDISStrucFunc::fNucSFNLOTables

private

Definition at line 249 of file HEDISStrucFunc.h.

Referenced by EvalNucSFNLO(), and HEDISStrucFunc().

fQrkSFLOTables

map<int, HEDISStrucFuncTable> genie::HEDISStrucFunc::fQrkSFLOTables

private

Definition at line 247 of file HEDISStrucFunc.h.

Referenced by EvalQrkSFLO(), and HEDISStrucFunc().

fSF

SF_info genie::HEDISStrucFunc::fSF

private

Definition at line 251 of file HEDISStrucFunc.h.

Referenced by CreateQrkSF(), and HEDISStrucFunc().

sf_q2_array

vector<double> genie::HEDISStrucFunc::sf_q2_array

private

Definition at line 253 of file HEDISStrucFunc.h.

Referenced by CreateQrkSF(), and HEDISStrucFunc().

sf_x_array

vector<double> genie::HEDISStrucFunc::sf_x_array

private

Definition at line 252 of file HEDISStrucFunc.h.

Referenced by CreateQrkSF(), and HEDISStrucFunc().

---

The documentation for this class was generated from the following files:

• HEDISStrucFunc.h
• HEDISStrucFunc.cxx