HEDISGenerator.cxx

Go to the documentation of this file.

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

 Author: Alfonso Garcia <alfonsog \at nikhef.nl>
         NIKHEF

 For the class documentation see the corresponding header file.

*/
//____________________________________________________________________________

#include "Physics/HEDIS/EventGen/HEDISGenerator.h"
#include "Physics/Hadronization/LeptoHadronization.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Numerical/MathUtils.h"
#include "Framework/Conventions/Constants.h"

#include <RVersion.h>
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
#include <TMCParticle.h>
#else
#include <TMCParticle6.h>
#endif

using namespace genie;
using namespace genie::utils::math;

//___________________________________________________________________________
HEDISGenerator::HEDISGenerator() :
HadronicSystemGenerator("genie::HEDISGenerator")
{
  this->Initialize();
}
//___________________________________________________________________________
HEDISGenerator::HEDISGenerator(string config) :
HadronicSystemGenerator("genie::HEDISGenerator", config)
{
  this->Initialize();
}
//___________________________________________________________________________
HEDISGenerator::~HEDISGenerator()
{

}
//____________________________________________________________________________
void HEDISGenerator::Initialize(void) const
{

}
//___________________________________________________________________________
void HEDISGenerator::ProcessEventRecord(GHepRecord * evrec) const
{
// This method generates the final state hadronic system

  //-- Add the target remnant
  this->AddTargetNucleusRemnant(evrec);
  GHepParticle * target = evrec -> TargetNucleus();
  if(target) evrec->Particle(evrec->RemnantNucleusPosition())->SetStatus(kIStFinalStateNuclearRemnant);

  //-- Add the primary lepton
  this->AddPrimaryLepton(evrec);

  //-- Run the hadronization model and get the fragmentation products
  fHadronizationModel->ProcessEventRecord(evrec);

}
//___________________________________________________________________________
void HEDISGenerator::AddPrimaryLepton(GHepRecord * evrec) const
{

  Interaction * interaction = evrec->Summary();

  // Neutrino 4p
  LongLorentzVector p4v( * evrec->Probe()->P4() );
  LOG("HEDISGenerator", pINFO) << "NEUTRINO @ LAB' =>  E = " << p4v.E() << " //  m = " << p4v.M() << " // p = " << p4v.P();
  LOG("HEDISGenerator", pINFO) << "                  dir = " << p4v.Dx() << " , "  << p4v.Dy() << " , "  << p4v.Dz();

  // Look-up selected kinematics & other needed kinematical params
  long double Q2  = interaction->Kine().Q2(true);
  long double y   = interaction->Kine().y(true);
  long double Ev  = p4v.E();
  long double ml  = interaction->FSPrimLepton()->Mass();
  long double ml2 = powl(ml,2);

  // Compute the final state primary lepton energy and momentum components
  // along and perpendicular the neutrino direction
  long double El  = (1-y)*Ev;
  long double plp = El - 0.5*(Q2+ml2)/Ev;                          // p(//)
  long double plt = sqrtl(fmaxl(0.,El*El-plp*plp-ml2)); // p(-|)
  // Randomize transverse components
  RandomGen * rnd = RandomGen::Instance();
  long double phi  = 2 * constants::kPi * rnd->RndLep().Rndm();
  long double pltx = plt * cosl(phi);
  long double plty = plt * sinl(phi);

  // Lepton 4-momentum in the LAB frame
  LongLorentzVector p4llong( pltx, plty, plp, El );
  p4llong.Rotate(p4v);
  LOG("HEDISGenerator", pINFO) << "LEPTON     @ LAB' =>  E = " << p4llong.E() << " //  m = " << p4llong.M() << " // p = " << p4llong.P();
  LOG("HEDISGenerator", pINFO) << "                    dir = " << p4llong.Dx() << " , "  << p4llong.Dy() << " , "  << p4llong.Dz();

  // Translate from long double to double
  TLorentzVector p4l( (double)p4llong.Px(), (double)p4llong.Py(), (double)p4llong.Pz(), (double)p4llong.E() );

  // Add lepton to EventRecord
  int pdgl = interaction->FSPrimLepton()->PdgCode();
  evrec->AddParticle(pdgl, kIStStableFinalState, evrec->ProbePosition(),-1,-1,-1, p4l, *(evrec->Probe()->X4()));
  evrec->Summary()->KinePtr()->SetFSLeptonP4(p4l);

}
//___________________________________________________________________________
void HEDISGenerator::Configure(const Registry & config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void HEDISGenerator::Configure(string config)
{
  Algorithm::Configure(config);
  this->LoadConfig();
}
//____________________________________________________________________________
void HEDISGenerator::LoadConfig(void)
{
  fHadronizationModel = 0;

  //-- Get the requested hadronization model
  fHadronizationModel =
     dynamic_cast<const EventRecordVisitorI *> (this->SubAlg("Hadronizer"));
  assert(fHadronizationModel);

}