gCalcHEDISDiffXsec.cxx File Reference

#include "TMath.h"
#include "TSystem.h"
#include "TTree.h"
#include "TFile.h"
#include "Framework/EventGen/XSecAlgorithmI.h"
#include "Framework/EventGen/InteractionList.h"
#include "Framework/EventGen/EventGeneratorI.h"
#include "Framework/EventGen/GEVGDriver.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/KineUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Interaction/InitialState.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Conventions/Units.h"
#include <fstream>
#include <iostream>

Include dependency graph for gCalcHEDISDiffXsec.cxx:

Go to the source code of this file.

Functions
void	PrintSyntax (void)
void	DecodeCommandLine (int argc, char *argv[])
vector< double >	ReadListFromPath (string path)
int	main (int argc, char **argv)

Variables
const double	epsilon = 1e-5
int	fNu = -1
int	fTgt = -1
string	fPathXlist = ""
string	fPathYlist = ""
string	fPathElist = ""
string	fOutFileName = ""
bool	fSaveAll = false

Function Documentation

DecodeCommandLine()

void DecodeCommandLine	(	int	argc,
		char *	argv[] )

Definition at line 251 of file gCalcHEDISDiffXsec.cxx.

                                                {
 
  // Common run options. Set defaults and read.
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);
 
  // Parse run options for this app
  CmdLnArgParser parser(argc,argv);
 
  if( parser.OptionExists('p') ){ 
    fNu = parser.ArgAsInt('p');
    LOG("gcalchedisdiffxsec", pINFO) << "Probe = " << fNu;
  }          
  else {
    LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input neutrino type!";
    PrintSyntax();
    exit(1);
  }
 
  if( parser.OptionExists('t') ){ 
    fTgt = parser.ArgAsInt('t');
    LOG("gcalchedisdiffxsec", pINFO) << "Target = " << fTgt;
  }          
  else {
    LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input target type!";
    PrintSyntax();
    exit(1);
  }
  
  if( parser.OptionExists('x') ){ 
    fPathXlist = parser.ArgAsString('x');
    LOG("gcalchedisdiffxsec", pINFO) << "PathXlist = " << fPathXlist;
  }          
  else {
    LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input path to Xlist!";
    PrintSyntax();
    exit(1);
  }
 
  if( parser.OptionExists('y') ){ 
    fPathYlist = parser.ArgAsString('y');
    LOG("gcalchedisdiffxsec", pINFO) << "PathYlist = " << fPathYlist;
  }          
  else {
    LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input path to Ylist!";
    PrintSyntax();
    exit(1);
  }
 
  if( parser.OptionExists('e') ){ 
    fPathElist = parser.ArgAsString('e');
    LOG("gcalchedisdiffxsec", pINFO) << "PathElist = " << fPathElist;
  }          
  else {
    LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input path to Elist!";
    PrintSyntax();
    exit(1);
  }
 
  if( parser.OptionExists('s') ) {
    fSaveAll = true;
    LOG("gcalchedisdiffxsec", pINFO) << "SaveAll = " << fSaveAll;
  }          
 
  if( parser.OptionExists('o') ){ 
    fOutFileName = parser.ArgAsString('o');
    LOG("gcalchedisdiffxsec", pINFO) << "OutFileName = " << fOutFileName;
  }          
  else {
    LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified output file name!";
    PrintSyntax();
    exit(1);
  }
 
}

References genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsString(), fNu, fOutFileName, fPathElist, fPathXlist, fPathYlist, fSaveAll, fTgt, genie::RunOpt::Instance(), LOG, genie::CmdLnArgParser::OptionExists(), pFATAL, pINFO, PrintSyntax(), and genie::RunOpt::ReadFromCommandLine().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv )

Definition at line 98 of file gCalcHEDISDiffXsec.cxx.

                                 {
 
  DecodeCommandLine(argc,argv); 
 
  RunOpt::Instance()->BuildTune();
  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
 
  string fChannel = RunOpt::Instance()->EventGeneratorList();
 
  GEVGDriver evg_driver;
  InitialState init_state(fTgt, fNu);
  evg_driver.SetEventGeneratorList(fChannel);
  evg_driver.Configure(init_state);
 
  vector<double> ve = ReadListFromPath(fPathElist);
  vector<double> vx = ReadListFromPath(fPathXlist);
  vector<double> vy = ReadListFromPath(fPathYlist);
 
  double widthx = (TMath::Log10(vx[1])-TMath::Log10(vx[0]));
 
  LOG("gcalchedisdiffxsec", pINFO) << widthx;
 
  int quark;
  double ei;
  double bx;
  double by;
  double dxsec;
  
  TString treename = Form("diffxsec_nu%d_tgt%d_%s",fNu,fTgt,fChannel.c_str());
 
  LOG("gcalchedisdiffxsec", pINFO) << treename;
 
  TTree * tree = new TTree(treename,treename);
  tree->Branch( "Quark",    &quark, "Quark/I"    );
  tree->Branch( "Ei",       &ei,    "Ei/D"       );
  tree->Branch( "Bx",       &bx,    "Bx/D"       );
  tree->Branch( "By",       &by,    "By/D"       );
  tree->Branch( "DiffXsec", &dxsec, "DiffXsec/D" );
 
  const InteractionList * intlst   = evg_driver.Interactions();
 
  InteractionList::const_iterator intliter;
  for(intliter = intlst->begin(); intliter != intlst->end(); ++intliter) {
 
    const Interaction * interaction = *intliter;
 
    quark = 10000 * interaction->InitState().Tgt().HitQrkPdg();
    if (!interaction->InitState().Tgt().HitSeaQrk()) {
      if ( quark>0 ) quark += 100;
      else           quark -= 100;
    }
    quark += 1 * interaction->ExclTag().FinalQuarkPdg();
 
    LOG("gcalchedisdiffxsec", pINFO) << "Current interaction: " << interaction->AsString();
 
    const XSecAlgorithmI * xsec_alg = evg_driver.FindGenerator(interaction)->CrossSectionAlg();
 
    for ( unsigned i=0; i<ve.size(); i++ ) {
 
      ei = ve[i];
 
      LOG("gcalchedisdiffxsec", pINFO) << "Energy: " << ei << " [GeV]";
 
      TLorentzVector p4(0,0,ei,ei);
      interaction->InitStatePtr()->SetProbeP4(p4);
 
      for ( unsigned j=0; j<vy.size(); j++ ) {
 
        double z = vy[j]; // z = (eo-em)/(ei-em)
 
        by = (1-z)*(ei-10.)/ei; // y = 1 - eo/ei;
 
        LOG("gcalchedisdiffxsec", pDEBUG) << "  z: " << z << "  y: " << by;
 
        if      ( by==1. ) by -= 1e-4; 
        else if ( by==0. ) {
          by = 5./ei;
          double by_prev = (1-vy[j-1])*(ei-10.)/ei;
          LOG("gcalchedisdiffxsec", pDEBUG) << "  by: " << by << "  by_prev: " << by_prev << "  " << vy[j-1];
          if (by>by_prev) by = 1e-4;
        }
        
        if ( by<0 || by>1 ) continue;
        
        interaction->KinePtr()->Sety(by);
        
        if (fSaveAll) {
          for ( unsigned k=0; k<vx.size(); k++ ) {
            bx = vx[k];
            interaction->KinePtr()->Setx(bx);
            utils::kinematics::UpdateWQ2FromXY(interaction);
            dxsec = xsec_alg->XSec(interaction, kPSxyfE) / units::cm2;
            LOG("gcalchedisdiffxsec", pDEBUG) << "x: " << bx << "  y: " << by << " -> d2sigmadxdy[E=" << ei << "GeV] = " << dxsec << " cm2";
            tree->Fill();
          }
        }
        else {
          dxsec = 0.;
          for ( unsigned k=0; k<vx.size(); k++ ) {
            bx = vx[k];
            interaction->KinePtr()->Setx(bx);
            utils::kinematics::UpdateWQ2FromXY(interaction);
            dxsec += xsec_alg->XSec(interaction, kPSxyfE) * bx;
          }
          dxsec *=  widthx*TMath::Log(10) / units::cm2;
          LOG("gcalchedisdiffxsec", pDEBUG) << "  y: " << by << " -> d2sigmady[E=" << ei << "GeV] = " << dxsec << " cm2";
          tree->Fill();
        }
        
      }
    
    }
 
  }
 
  TFile * outfile = new TFile(fOutFileName.c_str(),"RECREATE");
  tree->Write(tree->GetName());
  outfile->Close();
 
}

References genie::Interaction::AsString(), genie::RunOpt::BuildTune(), genie::units::cm2, genie::GEVGDriver::Configure(), genie::EventGeneratorI::CrossSectionAlg(), DecodeCommandLine(), e, genie::RunOpt::EventGeneratorList(), genie::Interaction::ExclTag(), genie::XclsTag::FinalQuarkPdg(), genie::GEVGDriver::FindGenerator(), fNu, fOutFileName, fPathElist, fPathXlist, fPathYlist, fSaveAll, fTgt, genie::Target::HitQrkPdg(), genie::Target::HitSeaQrk(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::RunOpt::Instance(), genie::GEVGDriver::Interactions(), genie::Interaction::KinePtr(), genie::kPSxyfE, LOG, genie::utils::app_init::MesgThresholds(), pDEBUG, pINFO, ReadListFromPath(), genie::GEVGDriver::SetEventGeneratorList(), genie::InitialState::SetProbeP4(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::InitialState::Tgt(), genie::utils::kinematics::UpdateWQ2FromXY(), and genie::XSecAlgorithmI::XSec().

PrintSyntax()

void PrintSyntax

(

void

)

Definition at line 327 of file gCalcHEDISDiffXsec.cxx.

{
  LOG("gcalchedisdiffxsec", pNOTICE)
      << "\n\n" << "Syntax:" << "\n"
      << "   gcalchedisdiffxsec -p nu -t tgt -o root_file\n"
      << "            -x pathxlist\n"
      << "            -y pathylist\n"
      << "            -e pathelist\n"
      << "            --tune genie_tune\n"
      << "            --event-generator-list list_name\n"
      << "            [--message-thresholds xml_file]\n";
}

References LOG, and pNOTICE.

Referenced by DecodeCommandLine().

ReadListFromPath()

vector< double > ReadListFromPath

(

string

path

)

Definition at line 225 of file gCalcHEDISDiffXsec.cxx.

                                             {
 
  vector<double> list;
 
  std::ifstream infile(path.c_str());
 
  //check to see that the file was opened correctly:
  if (!infile.is_open()) {
    LOG("gcalchedisdiffxsec", pFATAL) << "There was a problem opening the input file!";
    exit(1);//exit or do additional error checking
  }
 
  double val = 0.;
  while (infile >> val) list.push_back(val);
 
  return list;
 
}

References infile, LOG, and pFATAL.

Referenced by main().

Variable Documentation

epsilon

const double epsilon = 1e-5

Definition at line 76 of file gCalcHEDISDiffXsec.cxx.

Referenced by CheckChargeConservation(), CheckEnergyMomentumConservation(), genie::utils::math::CholeskyDecomposition(), genie::Intranuke2018::GenerateVertex(), genie::Intranuke2025::GenerateVertex(), genie::Intranuke::GenerateVertex(), osetUtils::isEqual(), osetUtils::isEqual(), genie::BardinIMDRadCorPXSec::Li2(), genie::SPPEventGenerator::Vertex::operator==(), genie::DISHadronicSystemGenerator::SimulateFormationZone(), genie::utils::intranuke2018::StepParticle(), genie::utils::intranuke2025::StepParticle(), genie::utils::intranuke::StepParticle(), and genie::RosenbluthPXSec::XSec().

fNu

int fNu = -1

Definition at line 84 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

fOutFileName

string fOutFileName = ""

Definition at line 89 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

fPathElist

string fPathElist = ""

Definition at line 88 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

fPathXlist

string fPathXlist = ""

Definition at line 86 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

fPathYlist

string fPathYlist = ""

Definition at line 87 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

fSaveAll

bool fSaveAll = false

Definition at line 90 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

fTgt

int fTgt = -1

Definition at line 85 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().