iFatras::HadIntProcessorParametric Class Reference

Parametric implementation of nuclear interactions to be used in Fatras. More...

#include <HadIntProcessorParametric.h>

Inheritance diagram for iFatras::HadIntProcessorParametric:

Collaboration diagram for iFatras::HadIntProcessorParametric:

Public Member Functions
	HadIntProcessorParametric (const std::string &, const std::string &, const IInterface *)
	AlgTool constructor for HadIntProcessorParametric.
virtual	~HadIntProcessorParametric ()
	Destructor.
StatusCode	initialize ()
	AlgTool initailize method.
StatusCode	finalize ()
	AlgTool finalize method.
bool	hadronicInteraction (const Amg::Vector3D &position, const Amg::Vector3D &momentum, double p, double E, double charge, const Trk::MaterialProperties &mprop, double pathCorrection, Trk::ParticleHypothesis particle=Trk::pion) const
	interface for processing of the nuclear interactions
bool	recordHadState (double time, double p, const Amg::Vector3D &vertex, const Amg::Vector3D &particleDir, Trk::ParticleHypothesis particle) const
	interface for processing of the presampled nuclear interaction
bool	doHadronicInteraction (double time, const Amg::Vector3D &position, const Amg::Vector3D &momentum, const Trk::Material *emat, Trk::ParticleHypothesis particle, bool processSecondaries) const
ISF::ISFParticleVector	doHadIntOnLayer (const ISF::ISFParticle *parent, double time, const Amg::Vector3D &position, const Amg::Vector3D &momentum, const Trk::Material *emat, Trk::ParticleHypothesis particle) const

Private Member Functions
ISF::ISFParticleVector	getHadState (const ISF::ISFParticle *parent, double time, double p, const Amg::Vector3D &vertex, const Amg::Vector3D &particleDir, Trk::ParticleHypothesis particle) const
	collect secondaries for layer material update

Static Private Member Functions
static double	absorptionLength (const Trk::MaterialProperties *mat, double p, double q, Trk::ParticleHypothesis particle=Trk::pion)
	interface for calculation of absorption length

Private Attributes
double	m_minimumHadOutEnergy
	hadronic interaction setting
double	m_childMomParam
bool	m_cutChain
bool	m_hadIntFromX0
double	m_hadIntProbScale
double	m_minimumHadInitialEnergy
ServiceHandle< ISF::IParticleBroker >	m_particleBroker
	ISF services & Tools.
ServiceHandle< ISF::ITruthSvc >	m_truthRecordSvc
ServiceHandle< IAtRndmGenSvc >	m_rndGenSvc
	Random Generator service.
int	m_processCode
	MCTruth process code for TruthIncidents created by this tool.
CLHEP::HepRandomEngine *	m_randomEngine
	Random engine.
std::string	m_randomEngineName
	Name of the random number stream.
bool	m_validationMode
ToolHandle< IPhysicsValidationTool >	m_validationTool
bool	m_hadIntValidation
std::string	m_hadIntValidationTreeName
	validation tree name - to be acessed by this from root
std::string	m_hadIntValidationTreeDescription
	validation tree description - second argument in TTree
std::string	m_hadIntValidationTreeFolder
	stream/folder to for the TTree to be written out
TTree *	m_hadIntValidationTree
	Root Validation Tree.
float	m_hadIntPointX
	ntuple variable : hadronic interaction point x coordinate
float	m_hadIntPointY
	ntuple variable : hadronic interaction point y coordinate
float	m_hadIntPointR
	ntuple variable : hadronic interaction point r distance
float	m_hadIntPointZ
	ntuple variable : hadronic interaction point z coordinate
int	m_hadIntMotherPdg
	ntuple variable : hadronic interaction mother Pdg
int	m_hadIntMotherBarcode
	ntuple variable : hadronic interaction mother barcode
float	m_hadIntMotherP
	ntuple variable : hadronic interaction mother momentum
float	m_hadIntMotherPt
	ntuple variable : hadronic interaction mother momentum
float	m_hadIntMotherPhi
	ntuple variable : hadronic interaction mother phi
float	m_hadIntMotherEta
	ntuple variable : hadronic interaction photon eta
int	m_hadIntChildren
	nutple variable : hadronic interaction children numbers
float	m_hadIntChildE
	nutple variable : hadronic interaction children total energy
float	m_hadIntChildP [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child Energy
float	m_hadIntChildPcms [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child Energy
int	m_hadIntChildPdg [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child Pdg
float	m_hadIntChildPhi [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child phi
float	m_hadIntChildEta [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child eta
float	m_hadIntChildTh [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child phi
float	m_hadIntChildThc [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child eta
float	m_hadIntChildDeltaPhi [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child delta phi
float	m_hadIntChildDeltaEta [MAXHADINTCHILDREN] {}
	nutple variable : hadronic interaction child delta eta

Detailed Description

Parametric implementation of nuclear interactions to be used in Fatras.

The parameterisation is gathered from Geant4.

Author

Andre.nosp@m.as.S.nosp@m.alzbu.nosp@m.rger.nosp@m.@cern.nosp@m..ch, Carst.nosp@m.en.M.nosp@m.agass.nosp@m.@cer.nosp@m.n.ch

Definition at line 54 of file HadIntProcessorParametric.h.

Constructor & Destructor Documentation

HadIntProcessorParametric()

iFatras::HadIntProcessorParametric::HadIntProcessorParametric	(	const std::string &	t,
		const std::string &	n,
		const IInterface *	p )

AlgTool constructor for HadIntProcessorParametric.

Definition at line 40 of file HadIntProcessorParametric.cxx.

                                                                                                                       :
  base_class(t,n,p),
  m_minimumHadOutEnergy(200.*CLHEP::MeV),
  m_childMomParam(5.),
  m_cutChain(false),
  m_hadIntFromX0(false),
  m_hadIntProbScale(1.0),
  m_minimumHadInitialEnergy(1000.*CLHEP::MeV),
  m_particleBroker("ISF_ParticleBroker", n),
  m_truthRecordSvc("ISF_TruthRecordSvc", n),
  m_rndGenSvc("AtDSFMTGenSvc", n),
  m_processCode(121),
  m_randomEngine(nullptr),
  m_randomEngineName("FatrasRnd"),
  m_validationMode(false),
  m_validationTool(""),
  m_hadIntValidation(false),
  m_hadIntValidationTreeName("FatrasHadronicInteractions"),
  m_hadIntValidationTreeDescription("Validation output from the HadIntProcessorParametric"),
  m_hadIntValidationTreeFolder("/val/FatrasHadronicInteractions"),
  m_hadIntValidationTree(nullptr),
  m_hadIntPointX(0.),
  m_hadIntPointY(0.),
  m_hadIntPointR(0.),
  m_hadIntPointZ(0.),
  m_hadIntMotherPdg(0),
  m_hadIntMotherBarcode(0),
  m_hadIntMotherP(0.),
  m_hadIntMotherPt(0.),
  m_hadIntMotherPhi(0.),
  m_hadIntMotherEta(0.),
  m_hadIntChildren(0),
  m_hadIntChildE(0.)
{
      // property setting
      declareProperty("MinimumHadronicOutEnergy"            , m_minimumHadOutEnergy);
      declareProperty("HadronicChildMomParam"               , m_childMomParam);
      declareProperty("HadronicInteractionFromX0"           , m_hadIntFromX0);
      declareProperty("HadronicInteractionScaleFactor"      , m_hadIntProbScale);
      declareProperty("MinimumHadronicInitialEnergy"        , m_minimumHadInitialEnergy);
      // the steering --------------------------------------------------------------
      declareProperty("ShortenHadIntChain"                  , m_cutChain);
      declareProperty("ValidationMode"                      , m_validationMode);
      declareProperty("PhysicsValidationTool"               , m_validationTool);
      declareProperty("HadronicInteractionValidation"       , m_hadIntValidation);
      // validation mode
      declareProperty("RandomNumberService"                 , m_rndGenSvc          , "Random number generator");
      declareProperty("RandomStreamName"                    , m_randomEngineName   , "Name of the random number stream");
      // MC Truth Properties
      declareProperty("PhysicsProcessCode"                  , m_processCode        , "MCTruth Physics Process Code" );
      // ISF Services and Tools
      declareProperty("ParticleBroker"                      , m_particleBroker     , "ISF ParticleBroker Svc"       );
      declareProperty("TruthRecordSvc"                      , m_truthRecordSvc     , "ISF Particle Truth Svc"       );
}

~HadIntProcessorParametric()

iFatras::HadIntProcessorParametric::~HadIntProcessorParametric ( )

virtualdefault

Destructor.

Member Function Documentation

absorptionLength()

double iFatras::HadIntProcessorParametric::absorptionLength ( const Trk::MaterialProperties * mat, double p, double q, Trk::ParticleHypothesis particle = Trk::pion )

staticprivate

interface for calculation of absorption length

absorption length

Definition at line 228 of file HadIntProcessorParametric.cxx.

                                                                                  {
  double al = mat->l0();
 
  /* // these parametrization comes from comparison with G4 sampler, but give too many interactions
     // not understood yet
     if ( particle == Trk::pion ) al = ( 0.53+0.2*exp(-p/1000.) ) * mat->l0();
     else if ( particle == Trk::kaon ) al = ( 0.1+0.65*exp(-p/5000.) ) * mat->l0();
     else if ( particle == Trk::proton ) al = 0.57 * mat->l0();
     else al =  mat->l0();
  */
 
  if ( particle == Trk::pion || particle == Trk::kaon || particle == Trk::pi0 || particle == Trk::k0) {
    al *= 1./(1.+ exp(-0.5*(p-270.)*(p-270.)/60./60.));
  }
  if ( particle == Trk::proton || particle == Trk::neutron ) al *=0.7;
  if ( particle == Trk::pion || particle == Trk::pi0) al *=0.9;
 
  return al;
}

doHadIntOnLayer()

ISF::ISFParticleVector iFatras::HadIntProcessorParametric::doHadIntOnLayer	(	const ISF::ISFParticle *	parent,
		double	time,
		const Amg::Vector3D &	position,
		const Amg::Vector3D &	momentum,
		const Trk::Material *	emat,
		Trk::ParticleHypothesis	particle ) const

Definition at line 689 of file HadIntProcessorParametric.cxx.

                                                                             {
  // called from McMaterialEffectsUpdator::interact and
  // McMaterialEffectsUpdator::interactLay methods
  return getHadState(parent, time, momentum.mag(), position, momentum.unit(), particle); // Registers TruthIncident internally
 
}

doHadronicInteraction()

bool iFatras::HadIntProcessorParametric::doHadronicInteraction	(	double	time,
		const Amg::Vector3D &	position,
		const Amg::Vector3D &	momentum,
		const Trk::Material *	emat,
		Trk::ParticleHypothesis	particle,
		bool	processSecondaries ) const

Definition at line 659 of file HadIntProcessorParametric.cxx.

                                                                                                                              {
  // Called by McMaterialEffectsUpdator::interact
  // get parent particle
  const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
  // something is seriously wrong if there is no parent particle
  assert(parent);
 
  ISF::ISFParticleVector ispVec=getHadState(parent, time, momentum.mag(), position, momentum.unit(), particle); // Registers TruthIncident internally
 
 
  // having no secondaries does not necessarily mean the interaction did not take place  : TODO : add flag into ::getHadState
  //  if (!ispVec.size()) return false;
 
  // push onto ParticleStack
  if (processSecondaries && !ispVec.empty() ) {
    for (auto *childParticle : ispVec) {
      //Check that the new ISFParticles have a valid TruthBinding
      if (!childParticle->getTruthBinding()) {
        ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
      }
      m_particleBroker->push(childParticle, parent);
    }
  }
 
  return true;
 
}

finalize()

StatusCode iFatras::HadIntProcessorParametric::finalize

(

)

AlgTool finalize method.

Definition at line 181 of file HadIntProcessorParametric.cxx.

{
    ATH_MSG_DEBUG( "finalize() successful" );
    return StatusCode::SUCCESS;
}

getHadState()

ISF::ISFParticleVector iFatras::HadIntProcessorParametric::getHadState ( const ISF::ISFParticle * parent, double time, double p, const Amg::Vector3D & vertex, const Amg::Vector3D & particleDir, Trk::ParticleHypothesis particle ) const

private

collect secondaries for layer material update

Definition at line 250 of file HadIntProcessorParametric.cxx.

{
  ISF::ISFParticleVector chDef(0);
 
  // sampling of hadronic interaction
  double m = Trk::ParticleMasses::mass[particle];
  double E = sqrt(p*p + m*m);
 
  // get the maximum multiplicity
  double multiplicity_max = 0.25 * E/1000. + 18.;
 
  // multiplicity distribution
  double randx , randy, arg = 0.;
 
  double p1 = 0.;
  double p2 = 0.;
 
  if (E > 15000.) {
    p1 = 8.69;
    p2 = 2.34;
  } else {
    p1 = 6.77;
    p2 = 2.02;
  }
 
  for (;;) {
    randx = 30. * CLHEP::RandFlat::shoot(m_randomEngine);
    randy = 1.  * CLHEP::RandFlat::shoot(m_randomEngine);
    arg = exp(-0.5*( (randx-p1)/p2 + exp(-(randx-p1)/p2) ) );
    if (randy < arg && randx>3 && randx<multiplicity_max) break;
  }
 
  randx *= (1.2-0.4*exp(-0.001*p));     // trying to adjust
 
  int Npart = (int)randx;
 
  // protection against Npart < 3
  if (Npart < 3) {
    ATH_MSG_VERBOSE( "[ had ] Number of particles smaller than 3, parameterisation not valid."
             << " Doing Nothing");
    return chDef;
  }
 
  ATH_MSG_VERBOSE( "[ had ] interaction of " << HepMC::barcode(parent)
           << " with " << Npart << " outgoing particles " );
 
  // record the interaction
 
  // ------ now create the new hadrons ------
  ATH_MSG_DEBUG(  "[ had ] create hadronic shower for particle with PDG ID "
          << parent->pdgCode() << " and barcode "
          << HepMC::barcode(parent) );
 
 
  // create the genParticles
 
  // validation if needed
  if (m_hadIntValidationTree){
    m_hadIntPointX        = vertex.x();
    m_hadIntPointY        = vertex.y();
    m_hadIntPointZ        = vertex.z();
    m_hadIntPointR        = vertex.perp();
    m_hadIntMotherPdg     = parent->pdgCode();
    m_hadIntMotherBarcode = HepMC::barcode(parent); // FIXME barcode-based
    m_hadIntMotherP       = p;
    m_hadIntMotherPt      = p*particleDir.perp();
    m_hadIntMotherPhi     = particleDir.phi();
    m_hadIntMotherEta     = particleDir.eta();
    // reset the children
    m_hadIntChildren      = 0;
  }
 
  ATH_MSG_VERBOSE( "[ had ] incoming particle energy | mass | momentum "
           << E << " | " << m << " | " << p << " | " );
 
  if (m_cutChain && ( HepMC::generations(parent) > 0 || (HepMC::barcode(parent) ==HepMC::UNDEFINED_ID && HepMC::uniqueID(parent) == HepMC::UNDEFINED_ID) ) ) {
    if (m_hadIntValidationTree) m_hadIntValidationTree->Fill();
    ATH_MSG_VERBOSE( "[ had ] interaction initiated by a secondary particle, no children saved " );
    return chDef;
  }
 
  int gen_part = 0;
 
  // new sampling: sample particle type and energy in the CMS frame of outgoing particles
  // creation of shower particles
  double chargedist = 0.;
  std::vector<double> charge(Npart);
  std::vector<Trk::ParticleHypothesis> childType(Npart);
  std::vector<double> newm(Npart);
  std::vector<int> pdgid(Npart);
 
  // children type sampling  : simplified
  //double pif = 0.19;
  //double nef = 0.20;
  //double prf = 0.20;
 
  // sample heavy particles (alpha) but don't save
  double pif = 0.10;
  double nef = 0.30;
  double prf = 0.30;
 
  if ( particle == Trk::pion || particle == Trk::kaon || particle == Trk::pi0 || particle == Trk::k0 ) {
      pif = 0.15;
      nef = 0.25;
      prf = 0.25;
    }
  if ( particle == Trk::proton ) {
    pif = 0.06;
    nef = 0.25;
    prf = 0.35;
  }
  if ( particle == Trk::neutron ) {
    pif = 0.03;
    nef = 0.35;
    prf = 0.17;
  }
 
  for (int i=0; i<Npart; i++) {
    chargedist  = CLHEP::RandFlat::shoot(m_randomEngine);
    if (chargedist<pif) {
      charge[i]=0.;
      childType[i]=Trk::pi0;
      newm[i]=Trk::ParticleMasses::mass[Trk::pi0]; // MeV
      pdgid[i]=111;
      continue;
    }
    if ( chargedist<2*pif) {
      charge[i]=1.;
      childType[i]=Trk::pion;
      newm[i]=Trk::ParticleMasses::mass[Trk::pion]; // MeV
      pdgid[i]=211;
      continue;
    }
    if (chargedist<3*pif) {
      charge[i]=-1.;
      childType[i]=Trk::pion;
      newm[i]=Trk::ParticleMasses::mass[Trk::pion]; // MeV
      pdgid[i]=-211;
      continue;
    }
    if (chargedist<3*pif+nef) {
      charge[i]=0.;
      childType[i]=Trk::neutron;
      newm[i]=ParticleConstants::neutronMassInMeV; // MeV
      pdgid[i]=2112; // neutron
      continue;
    }
    if (chargedist<3*pif+nef+prf) {
      charge[i]=1.;
      childType[i]=Trk::proton;
      newm[i]=Trk::ParticleMasses::mass[Trk::proton]; // MeV
      pdgid[i]=2212;
      continue;
    }
    charge[i]=2.;
    childType[i]=Trk::proton;
    newm[i]=4000.;
    pdgid[i]=20000;
  }
 
  // move the incoming particle type forward
  if ( childType[0] != particle ) {
    for (int i=1; i<Npart; i++) {
      if (childType[i]==particle) {
        childType[i]=childType[0];
        childType[0]=particle;
        double cho = charge[i];
        charge[i]=charge[0];
        charge[0]=parent ? parent->charge() : cho;
    newm[i]=Trk::ParticleMasses::mass[childType[i]]; // MeV
    newm[0]=Trk::ParticleMasses::mass[childType[0]]; // MeV
        break;
      }
    }
  }
 
  /*
  // sample momentum of daughters in CMS frame of outgoing particles  [ exp(-par/p) ]
  std::vector<double> mom;
  mom.clear();mom.reserve(Npart);
  double mom_n = 0.;
  for (int _npart=0; _npart<Npart; _npart++) {
    rand1  = CLHEP::RandFlat::shoot(m_randomEngine);
    mom_n = -log(rand1)/m_childMomParam * p;
    int ipos = _npart;
    while ( ipos>0 && mom_n>mom[ipos-1]) ipos--;
    mom.insert(mom.begin()+ipos,mom_n);
  }
 
  // check if configuration acceptable - if not, resample hardest mom
  double momR = 0.;
  for (int i=1; i<Npart; i++) momR += mom[i];
  if (momR < mom[0]) mom[0] = mom[1]+rand1*(momR-mom[1]);
  */
 
  std::vector<double> mom(Npart);
  std::vector<double> th(Npart);
  std::vector<double> ph(Npart);
 
  // sample first particle energy fraction and random momentum direction
  double eps = 2./Npart;
  double rnd  = CLHEP::RandFlat::shoot(m_randomEngine);
  mom[0] = 0.5*pow(eps,rnd);
  th[0]  = acos( 2*CLHEP::RandFlat::shoot(m_randomEngine)-1.);
  ph[0]  = 2*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
 
  // toss particles around in a way which preserves the total momentum (0.,0.,0.) at this point
  // TODO shoot first particle along the impact direction preferentially
 
  Amg::Vector3D ptemp(mom[0]*sin(th[0])*cos(ph[0]),mom[0]*sin(th[0])*sin(ph[0]),mom[0]*cos(th[0]));
  double ptot = mom[0];
 
  double theta = 0.; double phi = 0.;
  for (int i=1; i<Npart-2; i++) {
    eps = 1./(Npart-i);
    //mom[i] = pow(eps,CLHEP::RandFlat::shoot(m_randomEngine))*(1-ptot);
    mom[i] = ( eps + CLHEP::RandFlat::shoot(m_randomEngine)*(1-eps))*(1-ptot);
    if (ptemp.mag()<1-ptot) {
      while ( fabs(ptemp.mag()-mom[i])>1-ptot-mom[i] ){
    mom[i] =  ( eps + CLHEP::RandFlat::shoot(m_randomEngine)*(1-eps))*(1-ptot);
      }
    }
    // max p remaining
    double p_rem=1-ptot-mom[i];
    double cthmax = fmin(1.,(-ptemp.mag()*ptemp.mag()-mom[i]*mom[i]+p_rem*p_rem)/2/ptemp.mag()/mom[i]);
    //if (cthmax<-1.) std::cout <<"problem in theta sampling:p_rem:ptot:pcurr:"<<p_rem<<","<<ptemp.mag()<<","<<mom[i]<< std::endl;
    double rnd  = CLHEP::RandFlat::shoot(m_randomEngine);
    theta = acos( (cthmax+1.)*rnd-1.);
    phi = 2*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
    HepGeom::Vector3D<double> test(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
    HepGeom::Vector3D<double> dnewHep = HepGeom::RotateZ3D(ptemp.phi())*HepGeom::RotateY3D(ptemp.theta())*test;
    Amg::Vector3D dnew( dnewHep.x(), dnewHep.y(), dnewHep.z() );
    th[i]=dnew.theta();
    ph[i]=dnew.phi();
    ptemp += mom[i]*dnew;
    ptot += mom[i];
  }
 
  eps = 0.5;
  mom[Npart-2] = pow(eps,CLHEP::RandFlat::shoot(m_randomEngine))*(1-ptot);
  mom[Npart-1] = 1-ptot-mom[Npart-2];
 
  if (ptemp.mag()<1-ptot) {
    while (mom[Npart-1]+mom[Npart-2]<ptemp.mag()) {
      mom[Npart-2] = pow(eps,CLHEP::RandFlat::shoot(m_randomEngine))*(1-ptot);
      mom[Npart-1] = 1-ptot-mom[Npart-2];
    }
  }
  if (ptemp.mag()<fabs(mom[Npart-1]-mom[Npart-2]) ) {
    double diff = ptemp.mag()*CLHEP::RandFlat::shoot(m_randomEngine);
    double sum = mom[Npart-1]-mom[Npart-2];
    mom[Npart-2]=0.5*(sum+diff);
    mom[Npart-1]=0.5*(sum-diff);
  }
  double cth =(-ptemp.mag()*ptemp.mag()-mom[Npart-2]*mom[Npart-2]+mom[Npart-1]*mom[Npart-1])/2/ptemp.mag()/mom[Npart-2];
  if (fabs(cth)>1.) cth = (cth>0.) ? 1. : -1.;
 
  theta = acos(cth);
  phi = 2*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
  HepGeom::Vector3D<double> test(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
  HepGeom::Vector3D<double> dnewHep = HepGeom::RotateZ3D(ptemp.phi())*HepGeom::RotateY3D(ptemp.theta())*test;
  Amg::Vector3D dnew( dnewHep.x(), dnewHep.y(), dnewHep.z() );
 
  th[Npart-2]=dnew.theta();
  ph[Npart-2]=dnew.phi();
  ptemp += mom[Npart-2]*dnew;
  Amg::Vector3D dlast = -ptemp;
  th[Npart-1]=dlast.theta();
  ph[Npart-1]=dlast.phi();
 
  // particle sampled, rotate, boost and save final state
  //CLHEP::HepLorentzVector bv(p*particleDir.unit().x(),p*particleDir.unit().y(),p*particleDir.unit().z(),s_currentGenParticle->momentum().e()+mtot);
  double etot = 0.;
  for (int i=0;i<Npart; i++) etot += sqrt(mom[i]*mom[i]+newm[i]*newm[i]);
  double summ = 0.;
  for (int i=0;i<Npart; i++) summ += newm[i];
 
  // std::cout <<"hadronic interaction: current energy, expected :"<< etot <<","<< sqrt(summ*summ+2*summ*p+m*m)<< std::endl;
  // rescale (roughly) to the expected energy
  float scale = sqrt(summ*summ+2*summ*p+m*m)/etot;
  etot = 0.;
  for (int i=0;i<Npart; i++) {
    mom[i] *= scale;
    etot += sqrt(mom[i]*mom[i]+newm[i]*newm[i]);
  }
 
 
  CLHEP::HepLorentzVector bv(p*particleDir.unit().x(),p*particleDir.unit().y(),p*particleDir.unit().z(),sqrt(etot*etot+p*p));
  std::vector<CLHEP::HepLorentzVector> childBoost(Npart);
 
  //std::cout <<"boost vector:"<<p<<","<<bv.e()<<","<<bv.m()<<std::endl;
  //std::cout <<"etot, mother E,m:"<<etot<<","<<E<<","<<m<<std::endl;
 
  Amg::Vector3D in(0.,0.,0.);
  Amg::Vector3D fin(0.,0.,0.);
 
  for (int i=0; i<Npart; i++) {
    Amg::Vector3D dirCms(sin(th[i])*cos(ph[i]),sin(th[i])*sin(ph[i]),cos(th[i]));
    //Amg::Vector3D rotDirCms=HepGeom::RotateZ3D(particleDir.phi())*HepGeom::RotateY3D(particleDir.theta())*dirCms;
    Amg::Vector3D childP = mom[i]*dirCms;
    in += childP;
    CLHEP::HepLorentzVector newp(childP.x(),childP.y(),childP.z(),sqrt(mom[i]*mom[i]+newm[i]*newm[i]));
    CLHEP::HepLorentzVector childPB = newp.boost(bv.boostVector());
    childBoost[i]=childPB;
    fin += Amg::Vector3D(childPB.x(),childPB.y(),childPB.z());
  }
 
  double eout = 0.;
 
  // child particle vector for TruthIncident
  //  Reserve space for as many paricles as created due to hadr. int.
  //  However, the number of child particles for TruthIncident might be
  //  smaller due to (momentum) cuts
  ISF::ISFParticleVector           children(Npart);
  ISF::ISFParticleVector::iterator childrenIt = children.begin();
  unsigned short                numChildren = 0;
 
  for (int i=0; i<Npart; i++) {
    if (pdgid[i]<10000) {
      Amg::Vector3D childP = Amg::Vector3D(childBoost[i].x(),childBoost[i].y(),childBoost[i].z());
      Amg::Vector3D chP = Amg::Vector3D(sin(th[i])*cos(ph[i]),sin(th[i])*sin(ph[i]),cos(th[i]));
 
      eout += childBoost[i].e();
 
      // validation if needed
      if (m_hadIntValidationTree && m_hadIntChildren < MAXHADINTCHILDREN){
    m_hadIntChildPdg[m_hadIntChildren]      = pdgid[i];
    m_hadIntChildP[m_hadIntChildren]        = childP.mag();
    m_hadIntChildPcms[m_hadIntChildren]     = mom[i];
    m_hadIntChildTh[m_hadIntChildren]        = childP.unit().dot(particleDir);
    m_hadIntChildThc[m_hadIntChildren]       =chP.dot(particleDir);
    m_hadIntChildPhi[m_hadIntChildren]      = childP.phi();
    m_hadIntChildEta[m_hadIntChildren]      = childP.eta();
    double deltaPhi = m_hadIntMotherPhi - m_hadIntChildPhi[m_hadIntChildren];
    // rescale the deltaPhi
    deltaPhi -= deltaPhi > M_PI ? M_PI : 0.;
    deltaPhi += deltaPhi < -M_PI ? M_PI : 0.;
    m_hadIntChildDeltaPhi[m_hadIntChildren] = deltaPhi;
    m_hadIntChildDeltaEta[m_hadIntChildren] = m_hadIntMotherEta - m_hadIntChildEta[m_hadIntChildren];
    ++m_hadIntChildren;
      }
 
      if (childP.mag()> m_minimumHadOutEnergy) {
    // get the new particle
    double mass = Trk::ParticleMasses::mass[ childType[i] ];
 
    // create the particle
        const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
        const int id = HepMC::UNDEFINED_ID;
    ISF::ISFParticle *child = new ISF::ISFParticle ( vertex,
                             childP,
                             mass,
                             charge[i],
                             pdgid[i],
                             status,
                             time,
                             *parent,
                                                         id
                                                         );
    // in the validation mode, add process info
    if (m_validationMode) {
      ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
      validInfo->setProcess(m_processCode);
      if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
      else validInfo->setGeneration(1);     // assume parent is a primary track
      child->setUserInformation(validInfo);
    }
    // record child for TruthIncident
    *childrenIt = child;
    ++childrenIt; numChildren++;
      }
 
      gen_part++;
    }
  } // particle loop
 
  children.resize(numChildren);
  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
                   children,
                   m_processCode,
                   parent->nextGeoID(),
                   ISF::fKillsPrimary );
  m_truthRecordSvc->registerTruthIncident( truth);
  // At this point we need to update the properties of the
  // ISFParticles produced in the interaction
  truth.updateChildParticleProperties();
 
  // save info for validation
  if (m_validationMode && m_validationTool.isEnabled()) {
    Amg::Vector3D* nMom = nullptr;
    m_validationTool->saveISFVertexInfo(m_processCode,vertex,*parent,p*particleDir,nMom,children);
    delete nMom;
  }
 
 
  m_hadIntChildE = eout;
 
  if (m_hadIntValidationTree) m_hadIntValidationTree->Fill();
 
  ATH_MSG_VERBOSE( "[ had ] it was kinematically possible to create " << gen_part << " shower particles " );
 
  return children;
 
}

hadronicInteraction()

bool iFatras::HadIntProcessorParametric::hadronicInteraction	(	const Amg::Vector3D &	position,
		const Amg::Vector3D &	momentum,
		double	p,
		double	E,
		double	charge,
		const Trk::MaterialProperties &	mprop,
		double	pathCorrection,
		Trk::ParticleHypothesis	particle = Trk::pion ) const

interface for processing of the nuclear interactions

Definition at line 188 of file HadIntProcessorParametric.cxx.

{
  const Trk::MaterialProperties* extMprop = dynamic_cast<const Trk::MaterialProperties*>(&mprop);
  double prob = 0.;
 
  // m_hadIntProbScale is used later, not here
  if (extMprop && !m_hadIntFromX0) {
 
    double al = absorptionLength(extMprop, p, charge, particle);  // in mm
 
    if (al>0.) prob = exp( (-1.)*pathCorrection*extMprop->thickness()/al );
    else       prob = exp( (-1.)*pathCorrection*extMprop->thicknessInL0());
 
  } else {
    ATH_MSG_VERBOSE( "  [ had ] Nuclear Interaction length not available, using appox. from X0" );
    // using approximation lambda = 0.37 * Z * X0 instead -- giving a warning message
    prob = exp( (-1.)*pathCorrection*mprop.thicknessInX0() / ( 0.37 * mprop.averageZ())  );
  }
 
  if (msgLvl(MSG::VERBOSE))
    ATH_MSG_VERBOSE( "  [ had ] Probability of hadronic interaction = "
    << (1. - prob) * 100. * m_hadIntProbScale << " %." );
 
  // apply a global scalor of the probability
  // (1. - prob) is generally O(0.01), so this is the right way to scale it
  // TODO fix time info (if needed)
  if (CLHEP::RandFlat::shoot(m_randomEngine) < (1. - prob) * m_hadIntProbScale ) return recordHadState(0.,p,
                                                       position,
                                                       momentum.unit(),
                                                       particle); // Registers TruthIncident internally
 
  // no hadronic interactions were computed
  return false;
}

initialize()

StatusCode iFatras::HadIntProcessorParametric::initialize

(

)

AlgTool initailize method.

Definition at line 101 of file HadIntProcessorParametric.cxx.

{
 
    ATH_MSG_DEBUG( "initialize()" );
 
    // get the random generator serice
     if (m_rndGenSvc.retrieve().isFailure()){
          ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
          return StatusCode::FAILURE;
     } else
          ATH_MSG_VERBOSE( "Successfully retrieved " << m_rndGenSvc );
 
    //Get own engine with own seeds:
    m_randomEngine = m_rndGenSvc->GetEngine(m_randomEngineName);
    if (!m_randomEngine) {
        ATH_MSG_FATAL( "Could not get random engine '" << m_randomEngineName << "'" );
        return StatusCode::FAILURE;
    }
 
    // ISF Services
    if (m_particleBroker.retrieve().isFailure()){
        ATH_MSG_FATAL( "Could not retrieve " << m_particleBroker );
        return StatusCode::FAILURE;
    }
    if (m_truthRecordSvc.retrieve().isFailure()){
        ATH_MSG_FATAL( "Could not retrieve " << m_truthRecordSvc );
        return StatusCode::FAILURE;
    }
 
    // the validation setup -------------------------------- PART 3: Hadronic Interaction ---------------------
    ATH_CHECK( m_validationTool.retrieve( DisableTool{ m_validationTool.empty() || !m_validationMode } ) );
 
    if (m_hadIntValidation){
 
      SmartIF<ITHistSvc> tHistSvc{Gaudi::svcLocator()->service("THistSvc")};
      if (!tHistSvc) {
        ATH_MSG_ERROR( "initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
      }
      ATH_MSG_VERBOSE( "Booking hadronic interaction validation TTree ... " );
 
      // create the new Tree
      m_hadIntValidationTree = new TTree(m_hadIntValidationTreeName.c_str(), m_hadIntValidationTreeDescription.c_str());
 
      // counter for boundary surfaces
      m_hadIntValidationTree->Branch("HadIntPointX"    ,  &m_hadIntPointX, "hintX/F");
      m_hadIntValidationTree->Branch("HadIntPointY"    ,  &m_hadIntPointY, "hintY/F");
      m_hadIntValidationTree->Branch("HadIntPointR"    ,  &m_hadIntPointR, "hintR/F");
      m_hadIntValidationTree->Branch("HadIntPointZ"    ,  &m_hadIntPointZ, "hintZ/F");
      m_hadIntValidationTree->Branch("HadIntMotherPdg"    ,  &m_hadIntMotherPdg,       "hintMotherPdg/I");
      m_hadIntValidationTree->Branch("HadIntMotherBarcode"    ,  &m_hadIntMotherBarcode, "hintMotherBarcode/I");
      m_hadIntValidationTree->Branch("HadIntMotherP" ,  &m_hadIntMotherP,    "hintMotherP/F");
      m_hadIntValidationTree->Branch("HadIntMotherPt" ,  &m_hadIntMotherPt,    "hintMotherPt/F");
      m_hadIntValidationTree->Branch("HadIntMotherPhi"    ,  &m_hadIntMotherPhi,       "hintMotherPhi/F");
      m_hadIntValidationTree->Branch("HadIntMotherEta"    ,  &m_hadIntMotherEta,       "hintMohterEta/F");
 
      m_hadIntValidationTree->Branch("HadIntChildren"     ,  &m_hadIntChildren,        "hintcs/I");
      m_hadIntValidationTree->Branch("HadIntChildE  "     ,  &m_hadIntChildE,          "hintce/F");
      m_hadIntValidationTree->Branch("HadIntChildPdg"     ,  m_hadIntChildPdg,         "hintChildPdg[hintcs]/I");
      m_hadIntValidationTree->Branch("HadIntChildP"  ,  m_hadIntChildP,      "hintChildP[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildPcms"  ,  m_hadIntChildPcms,      "hintChildPcms[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildTh"  ,  m_hadIntChildTh,      "hintChildTh[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildThc"  ,  m_hadIntChildThc,      "hintChildThc[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildPhi"     ,  m_hadIntChildPhi,         "hintChildPhi[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildEta"     ,  m_hadIntChildEta,         "hintChildEta[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildDeltaPhi",  m_hadIntChildDeltaPhi,    "hintChildPhi[hintcs]/F");
      m_hadIntValidationTree->Branch("HadIntChildDeltaEta",  m_hadIntChildDeltaEta,    "hintChildEta[hintcs]/F");
 
 
      if ((tHistSvc->regTree(m_hadIntValidationTreeFolder, m_hadIntValidationTree)).isFailure()) {
        ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
        delete m_hadIntValidationTree; m_hadIntValidationTree = nullptr;
      } else {
        ATH_MSG_INFO( "TTree for Hadronic Interactions validation booked." );
      }
    } // ------------- end of validation mode -----------------------------------------------------------------
    ATH_MSG_DEBUG( "finalize() successful" );
    return StatusCode::SUCCESS;
}

recordHadState()

bool iFatras::HadIntProcessorParametric::recordHadState	(	double	time,
		double	p,
		const Amg::Vector3D &	vertex,
		const Amg::Vector3D &	particleDir,
		Trk::ParticleHypothesis	particle ) const

interface for processing of the presampled nuclear interaction

Definition at line 700 of file HadIntProcessorParametric.cxx.

                                                                   {
  // get parent particle
  const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
  // something is seriously wrong if there is no parent particle
  assert(parent);
 
  ISF::ISFParticleVector ispVec=getHadState(parent, time, p, vertex, particleDir, particle); // Registers TruthIncident internally
 
  // having no secondaries does not necessarily mean the interaction did not take place : TODO : add flag into ::getHadState
  //  if (!ispVec.size()) return false;
 
  // push onto ParticleStack
  if (!ispVec.empty() ) {
    for (auto *childParticle : ispVec) {
      //Check that the new ISFParticles have a valid TruthBinding
      if (!childParticle->getTruthBinding()) {
        ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
      }
      m_particleBroker->push(childParticle, parent);
    }
  }
 
  return true;
}

Member Data Documentation

m_childMomParam

double iFatras::HadIntProcessorParametric::m_childMomParam

private

Definition at line 102 of file HadIntProcessorParametric.h.

m_cutChain

bool iFatras::HadIntProcessorParametric::m_cutChain

private

Definition at line 103 of file HadIntProcessorParametric.h.

m_hadIntChildDeltaEta

float iFatras::HadIntProcessorParametric::m_hadIntChildDeltaEta[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child delta eta

Definition at line 155 of file HadIntProcessorParametric.h.

{};   

m_hadIntChildDeltaPhi

float iFatras::HadIntProcessorParametric::m_hadIntChildDeltaPhi[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child delta phi

Definition at line 154 of file HadIntProcessorParametric.h.

{};   

m_hadIntChildE

float iFatras::HadIntProcessorParametric::m_hadIntChildE

mutableprivate

nutple variable : hadronic interaction children total energy

Definition at line 146 of file HadIntProcessorParametric.h.

m_hadIntChildEta

float iFatras::HadIntProcessorParametric::m_hadIntChildEta[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child eta

Definition at line 151 of file HadIntProcessorParametric.h.

{};        

m_hadIntChildP

float iFatras::HadIntProcessorParametric::m_hadIntChildP[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child Energy

Definition at line 147 of file HadIntProcessorParametric.h.

{};     

m_hadIntChildPcms

float iFatras::HadIntProcessorParametric::m_hadIntChildPcms[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child Energy

Definition at line 148 of file HadIntProcessorParametric.h.

{};     

m_hadIntChildPdg

int iFatras::HadIntProcessorParametric::m_hadIntChildPdg[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child Pdg

Definition at line 149 of file HadIntProcessorParametric.h.

{};        

m_hadIntChildPhi

float iFatras::HadIntProcessorParametric::m_hadIntChildPhi[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child phi

Definition at line 150 of file HadIntProcessorParametric.h.

{};        

m_hadIntChildren

int iFatras::HadIntProcessorParametric::m_hadIntChildren

mutableprivate

nutple variable : hadronic interaction children numbers

Definition at line 145 of file HadIntProcessorParametric.h.

m_hadIntChildTh

float iFatras::HadIntProcessorParametric::m_hadIntChildTh[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child phi

Definition at line 152 of file HadIntProcessorParametric.h.

{};        

m_hadIntChildThc

float iFatras::HadIntProcessorParametric::m_hadIntChildThc[MAXHADINTCHILDREN] {}

mutableprivate

nutple variable : hadronic interaction child eta

Definition at line 153 of file HadIntProcessorParametric.h.

{};        

m_hadIntFromX0

bool iFatras::HadIntProcessorParametric::m_hadIntFromX0

private

Definition at line 104 of file HadIntProcessorParametric.h.

m_hadIntMotherBarcode

int iFatras::HadIntProcessorParametric::m_hadIntMotherBarcode

mutableprivate

ntuple variable : hadronic interaction mother barcode

Definition at line 139 of file HadIntProcessorParametric.h.

m_hadIntMotherEta

float iFatras::HadIntProcessorParametric::m_hadIntMotherEta

mutableprivate

ntuple variable : hadronic interaction photon eta

Definition at line 143 of file HadIntProcessorParametric.h.

m_hadIntMotherP

float iFatras::HadIntProcessorParametric::m_hadIntMotherP

mutableprivate

ntuple variable : hadronic interaction mother momentum

Definition at line 140 of file HadIntProcessorParametric.h.

m_hadIntMotherPdg

int iFatras::HadIntProcessorParametric::m_hadIntMotherPdg

mutableprivate

ntuple variable : hadronic interaction mother Pdg

Definition at line 138 of file HadIntProcessorParametric.h.

m_hadIntMotherPhi

float iFatras::HadIntProcessorParametric::m_hadIntMotherPhi

mutableprivate

ntuple variable : hadronic interaction mother phi

Definition at line 142 of file HadIntProcessorParametric.h.

m_hadIntMotherPt

float iFatras::HadIntProcessorParametric::m_hadIntMotherPt

mutableprivate

ntuple variable : hadronic interaction mother momentum

Definition at line 141 of file HadIntProcessorParametric.h.

m_hadIntPointR

float iFatras::HadIntProcessorParametric::m_hadIntPointR

mutableprivate

ntuple variable : hadronic interaction point r distance

Definition at line 136 of file HadIntProcessorParametric.h.

m_hadIntPointX

float iFatras::HadIntProcessorParametric::m_hadIntPointX

mutableprivate

ntuple variable : hadronic interaction point x coordinate

Definition at line 134 of file HadIntProcessorParametric.h.

m_hadIntPointY

float iFatras::HadIntProcessorParametric::m_hadIntPointY

mutableprivate

ntuple variable : hadronic interaction point y coordinate

Definition at line 135 of file HadIntProcessorParametric.h.

m_hadIntPointZ

float iFatras::HadIntProcessorParametric::m_hadIntPointZ

mutableprivate

ntuple variable : hadronic interaction point z coordinate

Definition at line 137 of file HadIntProcessorParametric.h.

m_hadIntProbScale

double iFatras::HadIntProcessorParametric::m_hadIntProbScale

private

Definition at line 105 of file HadIntProcessorParametric.h.

m_hadIntValidation

bool iFatras::HadIntProcessorParametric::m_hadIntValidation

private

Definition at line 126 of file HadIntProcessorParametric.h.

m_hadIntValidationTree

TTree* iFatras::HadIntProcessorParametric::m_hadIntValidationTree

private

Root Validation Tree.

Definition at line 131 of file HadIntProcessorParametric.h.

m_hadIntValidationTreeDescription

std::string iFatras::HadIntProcessorParametric::m_hadIntValidationTreeDescription

private

validation tree description - second argument in TTree

Definition at line 128 of file HadIntProcessorParametric.h.

m_hadIntValidationTreeFolder

std::string iFatras::HadIntProcessorParametric::m_hadIntValidationTreeFolder

private

stream/folder to for the TTree to be written out

Definition at line 129 of file HadIntProcessorParametric.h.

m_hadIntValidationTreeName

std::string iFatras::HadIntProcessorParametric::m_hadIntValidationTreeName

private

validation tree name - to be acessed by this from root

Definition at line 127 of file HadIntProcessorParametric.h.

m_minimumHadInitialEnergy

double iFatras::HadIntProcessorParametric::m_minimumHadInitialEnergy

private

Definition at line 106 of file HadIntProcessorParametric.h.

m_minimumHadOutEnergy

double iFatras::HadIntProcessorParametric::m_minimumHadOutEnergy

private

hadronic interaction setting

Definition at line 101 of file HadIntProcessorParametric.h.

m_particleBroker

ServiceHandle<ISF::IParticleBroker> iFatras::HadIntProcessorParametric::m_particleBroker

private

ISF services & Tools.

Definition at line 110 of file HadIntProcessorParametric.h.

m_processCode

int iFatras::HadIntProcessorParametric::m_processCode

private

MCTruth process code for TruthIncidents created by this tool.

Definition at line 116 of file HadIntProcessorParametric.h.

m_randomEngine

CLHEP::HepRandomEngine* iFatras::HadIntProcessorParametric::m_randomEngine

private

Random engine.

Definition at line 119 of file HadIntProcessorParametric.h.

m_randomEngineName

std::string iFatras::HadIntProcessorParametric::m_randomEngineName

private

Name of the random number stream.

Definition at line 120 of file HadIntProcessorParametric.h.

m_rndGenSvc

ServiceHandle<IAtRndmGenSvc> iFatras::HadIntProcessorParametric::m_rndGenSvc

private

Random Generator service.

Definition at line 114 of file HadIntProcessorParametric.h.

m_truthRecordSvc

ServiceHandle<ISF::ITruthSvc> iFatras::HadIntProcessorParametric::m_truthRecordSvc

private

Definition at line 111 of file HadIntProcessorParametric.h.

m_validationMode

bool iFatras::HadIntProcessorParametric::m_validationMode

private

Definition at line 123 of file HadIntProcessorParametric.h.

m_validationTool

ToolHandle<IPhysicsValidationTool> iFatras::HadIntProcessorParametric::m_validationTool

private

Definition at line 124 of file HadIntProcessorParametric.h.

---

The documentation for this class was generated from the following files:

• HadIntProcessorParametric.h
• HadIntProcessorParametric.cxx