#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

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

Member Enumeration Documentation

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::HEDISStrucFunc ( SF_info sfinfo )

private

Definition at line 47 of file HEDISStrucFunc.cxx.

References genie::units::cm, genie::HEDISInteractionListGenerator::CreateHEDISlist(), e, genie::RunOpt::Instance(), genie::kIntWeakCC, genie::kIntWeakNC, genie::kPdgAntiNuE, genie::kPdgNuE, genie::constants::kProtonMass, LOG, mPDFQrk, genie::TuneId::Name(), pDEBUG, pERROR, pFATAL, pINFO, pWARN, Q2PDFmax, Q2PDFmin, genie::RunOpt::Tune(), and xPDFmin.

 {
 
   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;
 
 }

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

private

HEDISStrucFunc::~HEDISStrucFunc ( )

private

Definition at line 340 of file HEDISStrucFunc.cxx.

 {
 
 }

Member Function Documentation

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

private

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

private

Definition at line 358 of file HEDISStrucFunc.cxx.

References genie::Interaction::ExclTag(), genie::XclsTag::FinalQuarkPdg(), 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(), genie::utils::kinematics::Q2(), Q2PDFmax, Q2PDFmin, genie::InitialState::Tgt(), and xPDFmin.

 {
 
   // 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();
 
 }

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

Definition at line 621 of file HEDISStrucFunc.cxx.

References genie::SF_xQ2::F1, genie::SF_xQ2::F2, and genie::SF_xQ2::F3.

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

 {
   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;
 }

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

Definition at line 631 of file HEDISStrucFunc.cxx.

References genie::SF_xQ2::F1, genie::SF_xQ2::F2, and genie::SF_xQ2::F3.

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

 {
   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;
 }

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

Definition at line 611 of file HEDISStrucFunc.cxx.

References genie::SF_xQ2::F1, genie::SF_xQ2::F2, and genie::SF_xQ2::F3.

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

 {
   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;
 }

HEDISStrucFunc * HEDISStrucFunc::Instance ( SF_info sfinfo )

static

Definition at line 345 of file HEDISStrucFunc.cxx.

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

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

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

int HEDISStrucFunc::NucSFCode ( const Interaction * in )

private

Definition at line 603 of file HEDISStrucFunc.cxx.

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().

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

string HEDISStrucFunc::NucSFName ( const Interaction * in )

private

Definition at line 584 of file HEDISStrucFunc.cxx.

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().

 {
   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;
 }

int HEDISStrucFunc::QrkSFCode ( const Interaction * in )

private

Definition at line 592 of file HEDISStrucFunc.cxx.

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().

 {
   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;
 }

string HEDISStrucFunc::QrkSFName ( const Interaction * in )

private

Definition at line 573 of file HEDISStrucFunc.cxx.

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().

 {
   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;
 }

Friends And Related Function Documentation

friend struct Cleaner

friend

Definition at line 265 of file HEDISStrucFunc.h.

Member Data Documentation

HEDISStrucFunc * HEDISStrucFunc::fgInstance = 0

staticprivate

Definition at line 244 of file HEDISStrucFunc.h.

Referenced by genie::HEDISStrucFunc::Cleaner::~Cleaner().

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

private

Definition at line 248 of file HEDISStrucFunc.h.

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

private

Definition at line 249 of file HEDISStrucFunc.h.

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

private

Definition at line 247 of file HEDISStrucFunc.h.

SF_info genie::HEDISStrucFunc::fSF

private

Definition at line 251 of file HEDISStrucFunc.h.

vector<double> genie::HEDISStrucFunc::sf_q2_array

private

Definition at line 253 of file HEDISStrucFunc.h.

vector<double> genie::HEDISStrucFunc::sf_x_array

private

Definition at line 252 of file HEDISStrucFunc.h.

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

Classes

Public Types

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Friends

Detailed Description

Member Typedef Documentation

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation