PhysInteractionSelector.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*
 Copyright (c) 2003-2020, The GENIE Collaboration
 For the full text of the license visit http://copyright.genie-mc.org

 Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
 University of Liverpool & STFC Rutherford Appleton Laboratory 
*/
//____________________________________________________________________________

#include <vector>
#include <sstream>
#include <cstdlib>
#include <iomanip>

#include <TMath.h>
#include <TLorentzVector.h>

#include "Framework/EventGen/XSecAlgorithmI.h"
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/PhysInteractionSelector.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/EventGeneratorI.h"
#include "Framework/EventGen/InteractionList.h"
#include "Framework/EventGen/InteractionGeneratorMap.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Numerical/Spline.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/PrintUtils.h"

using std::vector;
using std::endl;
using std::setw;
using std::setprecision;
using std::setfill;
using std::ostringstream;
using namespace genie;
//using namespace genie::units;

//___________________________________________________________________________
PhysInteractionSelector::PhysInteractionSelector() :
InteractionSelectorI("genie::PhysInteractionSelector")
{

}
//___________________________________________________________________________
PhysInteractionSelector::PhysInteractionSelector(string config) :
InteractionSelectorI("genie::PhysInteractionSelector", config)
{

}
//___________________________________________________________________________
PhysInteractionSelector::~PhysInteractionSelector()
{

}
//___________________________________________________________________________
EventRecord * PhysInteractionSelector::SelectInteraction
     (const InteractionGeneratorMap * igmap, const TLorentzVector & p4) const
{
  if(!igmap) {
     LOG("IntSel", pERROR)
        << "\n*** NULL InteractionGeneratorMap! Can't select interaction";
     return 0;
  }
  if(igmap->size() <= 0) {
     LOG("IntSel", pERROR)
       << "\n*** Empty InteractionGeneratorMap! Can't select interaction";
     return 0;
  }

  // Get the list of spline objects
  // Should have been constructed at the job initialization
  XSecSplineList * xssl = 0;
  if (fUseSplines) xssl = XSecSplineList::Instance();

  const InteractionList & ilst = igmap->GetInteractionList();
  vector<double> xseclist(ilst.size());

  string istate = ilst[0]->InitState().AsString();
  ostringstream msg;
  msg << "Selecting an interaction for the given initial state = "
      << istate << " at E = " << p4.E() << " GeV";

  LOG("IntSel", pNOTICE)
             << utils::print::PrintFramedMesg(msg.str(), 0, '=');
  LOG("IntSel", pNOTICE)
     << "Computing xsecs for all relevant modeled interactions:";

  unsigned int i=0;
  InteractionList::const_iterator intliter = ilst.begin();

  ostringstream xsec_table_printout;

  xsec_table_printout
      << " |"  << setfill('-') << setw(112) << "|" << endl
      << " | " << setfill(' ') << setw(80) << "interaction"
      << " | cross-section (1E-38*cm^2) |" << endl
      << " |"  << setfill('-') << setw(112) << "|" << endl;

  for( ; intliter != ilst.end(); ++intliter) {

     Interaction * interaction = new Interaction(**intliter);
     interaction->InitStatePtr()->SetProbeP4(p4);

     SLOG("IntSel", pDEBUG)
           << "Computing xsec for: \n  " << interaction->AsString();

     // get the cross section for this interaction
     const XSecAlgorithmI * xsec_alg =
               igmap->FindGenerator(interaction)->CrossSectionAlg();
     assert(xsec_alg);

     double xsec = 0; // cross section for this interaction

     bool spline_computed = xssl->SplineExists(xsec_alg, interaction);
     bool eval = fUseSplines && spline_computed;
     if (eval) {
           const InitialState & init = interaction->InitState();
           const ProcessInfo & proc  = interaction->ProcInfo();
           double E = init.ProbeE(kRfLab);
           if(TMath::IsNaN(E)) {
                 BLOG("IntSel", pFATAL) << *interaction;
                 BLOG("IntSel", pFATAL) << "E = " << E;
                 abort();
           }
           const Spline * spl = xssl->GetSpline(xsec_alg,interaction);
           if(spl->ClosestKnotValueIsZero(E,"-")) xsec = 0;
           else xsec = spl->Evaluate(E);
     } else {
           xsec = xsec_alg->Integral(interaction);
     }
     TMath::Max(0., xsec);
/*
     LOG("IntSel", pNOTICE)
       << interaction->AsString()
       << " --> xsec " << (eval ? "[**interp**]" : "[**calc**]")
       << " = " << xsec/genie::units::cm2 << " cm^2";
*/
     xsec_table_printout
           << " | " << setfill(' ') << setw(80) << interaction->AsString()
           << " | " << setfill(' ') << setw(26) << xsec/(1E-38*genie::units::cm2)
           << " | " << endl;

     xseclist[i++] = xsec;
     delete interaction;

  } // loop over interaction that can be generated

  xsec_table_printout
      << " |"  << setfill('-') << setw(112) << "|" << endl;

  LOG("IntSel", pNOTICE)
    << "\n" << xsec_table_printout.str();

  // select an interaction

  LOG("IntSel", pINFO)
            << "Selecting an entry from the Interaction List";
  double xsec_sum  = 0;
  for(unsigned int iint = 0; iint < xseclist.size(); iint++) {
     xsec_sum       += xseclist[iint];
     xseclist[iint]  = xsec_sum;

     SLOG("IntSel", pINFO)
             << "Sum{xsec}(0->" << iint << ") = " << xsec_sum;
  }
  RandomGen * rnd = RandomGen::Instance();
  double R = xsec_sum * rnd->RndISel().Rndm();

  LOG("IntSel", pINFO)
      << "Generating Rndm (0. -> max = " << xsec_sum << ") = " << R;

  for(unsigned int iint = 0; iint < xseclist.size(); iint++) {

     SLOG("IntSel", pDEBUG)
               << "Sum{xsec}(0->" << iint <<") = " << xseclist[iint];

     if( R < xseclist[iint] ) {
       Interaction * selected_interaction = new Interaction (*ilst[iint]);
       selected_interaction->InitStatePtr()->SetProbeP4(p4);

       // set the cross section for the selected interaction (just extract it
       // from the array of summed xsecs rather than recomputing it)
       double xsec_pedestal = (iint > 0) ? xseclist[iint-1] : 0.;
       double xsec = xseclist[iint] - xsec_pedestal;
       assert(xsec>0);

       LOG("IntSel", pNOTICE)
         << "Selected interaction: " << selected_interaction->AsString();

       // bootstrap the event record
       EventRecord * evrec = new EventRecord;
       evrec->AttachSummary(selected_interaction);
       evrec->SetXSec(xsec);

       return evrec;
     }
  }
  LOG("IntSel", pERROR) << "Could not select interaction";
  return 0;
}
//___________________________________________________________________________
void PhysInteractionSelector::Configure(const Registry & config)
{
  Algorithm::Configure(config);
  this->LoadConfigData();
}
//____________________________________________________________________________
void PhysInteractionSelector::Configure(string param_set)
{
  Algorithm::Configure(param_set);
  this->LoadConfigData();
}
//____________________________________________________________________________
void PhysInteractionSelector::LoadConfigData(void)
{
  //check whether the user prefers the cross sections to be calculated or
  //evaluated from a spline object constructed at the job initialization
  fUseSplines = false ;
  GetParam( "UseStoredXSecs", fUseSplines ) ;

}
//___________________________________________________________________________