SRCNuclearRecoil.cxx

Go to the documentation of this file.

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

 Author: Afroditi Papadopoulou <apapadop \at mit.edu>
         Massachusetts Institute of Technology - October 04, 2019

 @ October 4, 2019 - Afroditi Papadopoulou (AP)
   Created this new module that controls the addition of the recoil nucleon in the event record 
   and extracts its kinematics 
*/
//____________________________________________________________________________

#include <cstdlib>

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

#include "Framework/Algorithm/AlgFactory.h"
#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Units.h"
#include "Physics/NuclearState/SRCNuclearRecoil.h"

#include "Physics/NuclearState/NuclearModel.h"
#include "Physics/NuclearState/NuclearModelI.h"
#include "Framework/EventGen/EVGThreadException.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepStatus.h"
#include "Framework/GHEP/GHepFlags.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Physics/NuclearState/FermiMomentumTablePool.h"
#include "Physics/NuclearState/FermiMomentumTable.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/Utils/KineUtils.h"
#include "Physics/NuclearState/NuclearUtils.h"

using namespace genie;
using namespace genie::constants;

//___________________________________________________________________________
SRCNuclearRecoil::SRCNuclearRecoil() :
SecondNucleonEmissionI("genie::SRCNuclearRecoil")
{

}
//___________________________________________________________________________
SRCNuclearRecoil::SRCNuclearRecoil(string config) :
  SecondNucleonEmissionI("genie::SRCNuclearRecoil", config )
{

}

//___________________________________________________________________________

SRCNuclearRecoil::~SRCNuclearRecoil()
{

}

//___________________________________________________________________________

void SRCNuclearRecoil::ProcessEventRecord(GHepRecord * evrec) const
{

  const Interaction *  interaction = evrec       -> Summary();
  const InitialState & init_state  = interaction -> InitState();
  const Target       & tgt         = init_state.Tgt();

  // do nothing for non-nuclear targets
  if(! tgt.IsNucleus()) return;

  // access the hit nucleon and target nucleus at the GHEP record
  GHepParticle * nucleon = evrec->HitNucleon();
  GHepParticle * nucleus = evrec->TargetNucleus();
  assert(nucleon);
  assert(nucleus);

  // Set this to either a proton or neutron to eject a secondary particle
  int eject_nucleon_pdg = this->SRCRecoilPDG( *nucleon, tgt );

  // Ejection of secondary particle
  if (eject_nucleon_pdg != 0) { EmitSecondNucleon(evrec,eject_nucleon_pdg); }

}

//___________________________________________________________________________

int SRCNuclearRecoil::SRCRecoilPDG( const GHepParticle & nucleon, const Target & tgt) const {

      int eject_nucleon_pdg = 0;

      // const int nucleus_pdgc = tgt->Pdg();
      const int nucleon_pdgc = nucleon.Pdg();

      double pN2 = TMath::Power(nucleon.P4()->Rho(),2.); // (momentum of struck nucleon)^2

      double kF = fNuclModel -> LocalFermiMomentum( tgt, 
                                                    nucleon_pdgc, 
                                                    nucleon.X4()->Vect().Mag() );

      if (TMath::Sqrt(pN2) > kF) {
        double Pp = (nucleon_pdgc == kPdgProton) ? fPPPairPercentage : fPNPairPercentage;
        RandomGen * rnd = RandomGen::Instance();
        double prob = rnd->RndGen().Rndm();
        eject_nucleon_pdg = (prob > Pp) ? kPdgNeutron : kPdgProton;
      }
      
      return eject_nucleon_pdg;
}
//___________________________________________________________________________
void SRCNuclearRecoil::Configure(const Registry & config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void SRCNuclearRecoil::Configure(string config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void SRCNuclearRecoil::LoadConfig(void)
{

  SecondNucleonEmissionI::LoadConfig() ;

  this->GetParamDef("PNPairPercentage",       fPNPairPercentage,    0.95);

  if (fPNPairPercentage < 0. || fPNPairPercentage > 1.) { 

        LOG("SRCNuclearRecoil", pFATAL)
        << "PNPairPercentage either less than 0 or greater than 1: Exiting" ;

        exit(78); 
  }

  fPPPairPercentage = 1. - fPNPairPercentage;

}
//____________________________________________________________________________