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

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

Definition at line 251 of file gCalcHEDISDiffXsec.cxx.

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

                                                {

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

}

int main	(	int	argc,
		char **	argv
	)

Definition at line 98 of file gCalcHEDISDiffXsec.cxx.

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

                                 {

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

}

void PrintSyntax

(

void

)

vector< double > ReadListFromPath

(

string

path

)

Definition at line 225 of file gCalcHEDISDiffXsec.cxx.

References infile, LOG, and pFATAL.

Referenced by main().

                                             {

  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;

}

Variable Documentation

const double epsilon = 1e-5

Definition at line 76 of file gCalcHEDISDiffXsec.cxx.

Referenced by CheckChargeConservation(), CheckEnergyMomentumConservation(), genie::utils::math::CholeskyDecomposition(), genie::Intranuke::GenerateVertex(), genie::Intranuke2018::GenerateVertex(), osetUtils::isEqual(), genie::BardinIMDRadCorPXSec::Li2(), genie::SPPEventGenerator::Vertex::operator==(), genie::KPhaseSpace::Q2Lim_W_SPP(), genie::KPhaseSpace::Q2Lim_W_SPP_iso(), genie::SmithMonizUtils::Q2QES_SM_lim(), genie::SmithMonizUtils::QEL_EnuMin_SM(), genie::DISHadronicSystemGenerator::SimulateFormationZone(), genie::utils::intranuke::StepParticle(), genie::utils::intranuke2018::StepParticle(), genie::KPhaseSpace::WLim_SPP(), genie::KPhaseSpace::WLim_SPP_iso(), and genie::RosenbluthPXSec::XSec().

int fNu = -1

Definition at line 84 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

string fOutFileName = ""

Definition at line 89 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

string fPathElist = ""

Definition at line 88 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

string fPathXlist = ""

Definition at line 86 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

string fPathYlist = ""

Definition at line 87 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

bool fSaveAll = false

Definition at line 90 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().

int fTgt = -1

Definition at line 85 of file gCalcHEDISDiffXsec.cxx.

Referenced by DecodeCommandLine(), and main().