Implementation of nuclear interactions of hadrons in the tracker layers (based on FTF model of Geant4). More...

Inheritance diagram for fastsim::NuclearInteractionFTF:

Public Member Functions
void	interact (fastsim::Particle &particle, const SimplifiedGeometry &layer, std::vector< std::unique_ptr< fastsim::Particle > > &secondaries, const RandomEngineAndDistribution &random) override
	Perform the interaction. More...

	NuclearInteractionFTF (const std::string &name, const edm::ParameterSet &cfg)
	Constructor. More...

	~NuclearInteractionFTF () override
	Default destructor. More...

Public Member Functions inherited from fastsim::InteractionModel
const std::string	getName ()
	Return (unique) name of this interaction. More...

	InteractionModel (std::string name)
	Constructor. More...

virtual void	registerProducts (edm::ProducerBase &producer) const
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...

virtual void	storeProducts (edm::Event &iEvent)
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...

virtual	~InteractionModel ()
	Default destructor. More...

Private Attributes
const double	corrfactors_el [numHadrons][npoints]
	elastic interaction length corrections per particle and energy More...

const double	corrfactors_inel [numHadrons][npoints]
	inelastic interaction length corrections per particle and energy More...

G4LorentzVector	curr4Mom

int	currIdx
	Index of particle in vector of Geant4 particles. More...

const G4ParticleDefinition *	currParticle

G4Track *	currTrack

G4Step *	dummyStep

size_t	index
	Index for energy of particle (vectors parametrized as a function of energy) More...

double	intLengthElastic
	Elastic interaction length. More...

double	intLengthInelastic
	Inelastic interaction length. More...

const double	nuclElLength [numHadrons]
	elastic interaction length in Silicon at 1 TeV per particle More...

const double	nuclInelLength [numHadrons]
	inelastic interaction length in Silicon at 1 TeV per particle More...

G4Nucleus	targetNucleus

G4CascadeInterface *	theBertiniCascade

double	theBertiniLimit
	Bertini cascade for low-energy hadrons. More...

G4ThreeVector	theBoost

G4GeneratorPrecompoundInterface *	theCascade

G4DiffuseElastic *	theDiffuseElastic

double	theDistCut
	Cut on deltaR for ClosestChargedDaughter algorithm (FastSim tracking) More...

double	theEnergyLimit
	Minimum energy of interacting particle. More...

G4TheoFSGenerator *	theHadronicModel

G4LundStringFragmentation *	theLund

G4HadProjectile	theProjectile

G4ExcitedStringDecay *	theStringDecay

G4FTFModel *	theStringModel

G4PhysicsLogVector *	vect

G4ThreeVector	vectProj

Static Private Attributes
static const int	npoints = 21
	Number of steps to cover range of particle energies. More...

static const int	numHadrons = 30
	Number of G4 hadrons. More...

static const G4ParticleDefinition *	theG4Hadron [numHadrons] = { 0 }
	The Geant4 particles. More...

static int	theId [numHadrons] = {0}
	The pdgIDs of the Geant4 particles. More...

Detailed Description

Implementation of nuclear interactions of hadrons in the tracker layers (based on FTF model of Geant4).

Computes the probability for hadrons to interact with a nucleon of the tracker material (inelastically) and then sample nuclear interaction using FTF model of Geant4 Also, another implementation of nuclear interactions can be used that is based on FullSim events (NuclearInteraction).

Definition at line 104 of file NuclearInteractionFTF.cc.

Constructor & Destructor Documentation

fastsim::NuclearInteractionFTF::NuclearInteractionFTF	(	const std::string &	name,
		const edm::ParameterSet &	cfg
	)

Constructor.

Definition at line 246 of file NuclearInteractionFTF.cc.

References currIdx, currParticle, currTrack, dummyStep, edm::ParameterSet::getParameter(), GeV, mps_fire::i, index, fastsim::initializeOnce, intLengthElastic, intLengthInelastic, BPhysicsValidation_cfi::Lambda, MeV, npoints, numHadrons, phase2TrackerDigitizer_cfi::SigmaZero, targetNucleus, theBertiniCascade, theBertiniLimit, theDiffuseElastic, theDistCut, theEnergyLimit, theG4Hadron, theHadronicModel, theId, theLund, theStringDecay, theStringModel, and vect.

     : fastsim::InteractionModel(name)
     , curr4Mom(0.,0.,0.,0.)
     , vectProj(0.,0.,1.)
     , theBoost(0.,0.,0.)
 {
     // Angular distance of particle before/after interaction. If it too large, a daughter is created instead
     theDistCut = cfg.getParameter<double>("distCut");
     // Set energy limits for processes
     theBertiniLimit = cfg.getParameter<double>("bertiniLimit") * CLHEP::GeV;
     theEnergyLimit = cfg.getParameter<double>("energyLimit") * CLHEP::GeV;
 
     // FTF model
     theHadronicModel = new G4TheoFSGenerator("FTF");
     theStringModel = new G4FTFModel();
     G4GeneratorPrecompoundInterface* cascade 
       = new G4GeneratorPrecompoundInterface(new CMSDummyDeexcitation());
     theLund = new G4LundStringFragmentation();
     theStringDecay = new G4ExcitedStringDecay(theLund);
     theStringModel->SetFragmentationModel(theStringDecay);
 
     theHadronicModel->SetTransport(cascade);
     theHadronicModel->SetHighEnergyGenerator(theStringModel);
     theHadronicModel->SetMinEnergy(theEnergyLimit);
 
     // Bertini Cascade 
     theBertiniCascade = new G4CascadeInterface();
 
     theDiffuseElastic = new G4DiffuseElastic();
 
     // Geant4 particles and cross sections
     std::call_once(initializeOnce, [this] () {
         theG4Hadron[0] = G4Proton::Proton();
         theG4Hadron[1] = G4Neutron::Neutron();
         theG4Hadron[2] = G4PionPlus::PionPlus();
         theG4Hadron[3] = G4PionMinus::PionMinus();
         theG4Hadron[4] = G4KaonPlus::KaonPlus();
         theG4Hadron[5] = G4KaonMinus::KaonMinus();
         theG4Hadron[6] = G4KaonZeroLong::KaonZeroLong();
         theG4Hadron[7] = G4KaonZeroShort::KaonZeroShort();
         theG4Hadron[8] = G4KaonZero::KaonZero();
         theG4Hadron[9] = G4AntiKaonZero::AntiKaonZero();
         theG4Hadron[10]= G4Lambda::Lambda();
         theG4Hadron[11]= G4OmegaMinus::OmegaMinus();
         theG4Hadron[12]= G4SigmaMinus::SigmaMinus();
         theG4Hadron[13]= G4SigmaPlus::SigmaPlus();
         theG4Hadron[14]= G4SigmaZero::SigmaZero();
         theG4Hadron[15]= G4XiMinus::XiMinus();
         theG4Hadron[16]= G4XiZero::XiZero();
         theG4Hadron[17]= G4AntiProton::AntiProton();
         theG4Hadron[18]= G4AntiNeutron::AntiNeutron();
         theG4Hadron[19]= G4AntiLambda::AntiLambda();
         theG4Hadron[20]= G4AntiOmegaMinus::AntiOmegaMinus();
         theG4Hadron[21]= G4AntiSigmaMinus::AntiSigmaMinus();
         theG4Hadron[22]= G4AntiSigmaPlus::AntiSigmaPlus();
         theG4Hadron[23]= G4AntiSigmaZero::AntiSigmaZero();
         theG4Hadron[24]= G4AntiXiMinus::AntiXiMinus();
         theG4Hadron[25]= G4AntiXiZero::AntiXiZero();
         theG4Hadron[26]= G4AntiAlpha::AntiAlpha();
         theG4Hadron[27]= G4AntiDeuteron::AntiDeuteron();
         theG4Hadron[28]= G4AntiTriton::AntiTriton();
         theG4Hadron[29]= G4AntiHe3::AntiHe3();
 
         // other Geant4 particles
         G4ParticleDefinition* ion = G4GenericIon::GenericIon();
         ion->SetProcessManager(new G4ProcessManager(ion));
         G4DecayPhysics decays;
         decays.ConstructParticle();  
         G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
         partTable->SetVerboseLevel(0);
         partTable->SetReadiness();
 
         for(int i=0; i<numHadrons; ++i) {
             theId[i] = theG4Hadron[i]->GetPDGEncoding();
         }
     }); 
 
     // local objects
     vect = new G4PhysicsLogVector(npoints-1,100*MeV,TeV);
     intLengthElastic = intLengthInelastic = 0.0;
     currIdx = 0;
     index = 0;
     currTrack = nullptr;
     currParticle = theG4Hadron[0];
 
     // fill projectile particle definitions
     dummyStep = new G4Step();
     dummyStep->SetPreStepPoint(new G4StepPoint());
 
     // target is always Silicon
     targetNucleus.SetParameters(28, 14);
 }

fastsim::NuclearInteractionFTF::~NuclearInteractionFTF ( )

override

Default destructor.

Definition at line 339 of file NuclearInteractionFTF.cc.

References theLund, theStringDecay, theStringModel, and vect.

                                                      {
   delete theStringDecay;
   delete theStringModel;
   delete theLund;
   delete vect;
 }

Member Function Documentation

void fastsim::NuclearInteractionFTF::interact	(	fastsim::Particle &	particle,
		const SimplifiedGeometry &	layer,
		std::vector< std::unique_ptr< fastsim::Particle > > &	secondaries,
		const RandomEngineAndDistribution &	random
	)

overridevirtual

Perform the interaction.

Parameters

particle	The particle that interacts with the matter.
layer	The detector layer that interacts with the particle.
secondaries	Particles that are produced in the interaction (if any).
random	The Random Engine.

Implements fastsim::InteractionModel.

Definition at line 347 of file NuclearInteractionFTF.cc.

References funct::abs(), fastsim::Particle::charge(), corrfactors_el, corrfactors_inel, funct::cos(), curr4Mom, currIdx, currParticle, currTrack, DEFINE_EDM_PLUGIN, dir, dummyStep, MillePedeFileConverter_cfg::e, RandomEngineAndDistribution::flatShoot(), fastsim::Particle::getMotherDeltaR(), fastsim::Particle::getMotherPdgId(), fastsim::SimplifiedGeometry::getNuclearInteractionThicknessFactor(), fastsim::SimplifiedGeometry::getThickness(), GeV, index, intLengthElastic, intLengthInelastic, ResonanceBuilder::mass, fastsim::Particle::momentum(), npoints, nuclElLength, nuclInelLength, numHadrons, fastsim::Particle::pdgId(), fastsim::Particle::position(), mps_fire::result, fastsim::Particle::simTrackIndex(), funct::sin(), mathSSE::sqrt(), targetNucleus, theBertiniCascade, theBertiniLimit, theBoost, theDiffuseElastic, theDistCut, theEnergyLimit, theG4Hadron, theHadronicModel, theId, theProjectile, vect, and vectProj.

 {
     int thePid = particle.pdgId();
     //
     // no valid PDGid
     //
     if(abs(thePid) <= 100 || abs(thePid) >= 1000000)
     {
     return;
     }
 
     // particle lost all its enrgy in previous interaction
     if (particle.momentum().E() < 1E-10)
     {
     return;
     }
 
     double radLengths = layer.getThickness(particle.position(),particle.momentum());
     // TEC layers have some fudge factor for nuclear interactions
     radLengths *= layer.getNuclearInteractionThicknessFactor();
     //
     // no material
     //
     if(radLengths < 1E-10)
     {
     return;
     }
 
     // get the G4 hadron
     if(thePid != theId[currIdx]) {
         currParticle = nullptr;
         currIdx = 0;
         for(; currIdx<numHadrons; ++currIdx) {
             if(theId[currIdx] == thePid) {
                 currParticle = theG4Hadron[currIdx];
                 // neutral kaons
                 if(7 == currIdx || 8 == currIdx) {
                   currParticle = theG4Hadron[9];
                   if(random.flatShoot() > 0.5) { currParticle = theG4Hadron[10]; }
                 }   
                 break;
             }
         } 
     }
 
     // no valid G4 hadron found
     if(!currParticle){
     return;
     }
 
     // fill projectile for Geant4
     double mass = currParticle->GetPDGMass();
     double ekin = CLHEP::GeV * particle.momentum().e() - mass;
 
     // check interaction length
     intLengthElastic   = nuclElLength[currIdx];
     intLengthInelastic = nuclInelLength[currIdx];
     if(0.0 == intLengthInelastic){
     return;
     }
 
     // apply corrections
     if(ekin <= vect->Energy(0)){
         intLengthElastic   *= corrfactors_el[currIdx][0];
         intLengthInelastic *= corrfactors_inel[currIdx][0];
     }else if(ekin < vect->Energy(npoints-1)){
         index = vect->FindBin(ekin, index);
         double e1 = vect->Energy(index);
         double e2 = vect->Energy(index+1);
         intLengthElastic   *= ((corrfactors_el[currIdx][index] * (e2 - ekin) + 
               corrfactors_el[currIdx][index+1] * (ekin - e1)) / (e2 - e1));
         intLengthInelastic *= ((corrfactors_inel[currIdx][index] * (e2 - ekin) + 
               corrfactors_inel[currIdx][index+1] * (ekin - e1)) / (e2 - e1));
     }
 
     double currInteractionLength = -G4Log(random.flatShoot()) *
                 intLengthElastic * intLengthInelastic / (intLengthElastic + intLengthInelastic); 
 
     // Check position of nuclear interaction
     if(currInteractionLength > radLengths)
     {
     return; 
     }
 
     // fill projectile for Geant4
     double px = particle.momentum().px();
     double py = particle.momentum().py();
     double pz = particle.momentum().pz();
     double ptot = sqrt(px*px + py*py + pz*pz);
     double norm = 1./ptot;
     G4ThreeVector dir(px*norm, py*norm, pz*norm);
 
     G4DynamicParticle* dynParticle = new G4DynamicParticle(theG4Hadron[currIdx],dir,ekin);
     currTrack = new G4Track(dynParticle, 0.0, vectProj);
     currTrack->SetStep(dummyStep);
 
     theProjectile.Initialise(*currTrack); 
     delete currTrack;
 
     G4HadFinalState* result;
 
     // elastic interaction
     if(random.flatShoot() > intLengthElastic / (intLengthElastic + intLengthInelastic)){
         result = theDiffuseElastic->ApplyYourself(theProjectile, targetNucleus);
         G4ThreeVector dirnew = result->GetMomentumChange().unit();
         double cost = (dir*dirnew);
         double sint = std::sqrt((1. - cost) * (1. + cost));
 
         // This vector is already in units of GeV (FastSim standard even if it is a G4Vector)
         curr4Mom.set(ptot*dirnew.x(),ptot*dirnew.y(),ptot*dirnew.z(),particle.momentum().e());
 
         // Always create a daughter if the kink is large enough
         if(sint > theDistCut){ 
             // Create secondary
             secondaries.emplace_back(new fastsim::Particle(thePid
                 ,particle.position()
                 ,XYZTLorentzVector(curr4Mom.px(), curr4Mom.py(), curr4Mom.pz(), curr4Mom.e())));
 
             // ClosestChargedDaughter thing for FastSim tracking
             if(particle.charge() != 0){
                 secondaries.back()->setMotherDeltaR(particle.momentum());
                 secondaries.back()->setMotherPdgId(thePid);
                 secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
             }
 
             // The particle is destroyed
             particle.momentum().SetXYZT(0., 0., 0., 0.);
         }else{
             // ... or just rotate particle
             particle.momentum().SetXYZT(curr4Mom.px(), curr4Mom.py(), curr4Mom.pz(), curr4Mom.e());
         }
 
     // inelastic interaction
     }else{
         // Bertini cascade for low-energy hadrons (except light anti-nuclei)
         // FTFP is applied above energy limit and for all anti-hyperons and anti-ions 
         if(ekin <= theBertiniLimit && currIdx < 17){ 
             result = theBertiniCascade->ApplyYourself(theProjectile, targetNucleus);
         }else{
             result = theHadronicModel->ApplyYourself(theProjectile, targetNucleus);
         }
 
         if(result){
             int nsec = result->GetNumberOfSecondaries();
             if(nsec > 0){
                 result->SetTrafoToLab(theProjectile.GetTrafoToLab());
 
                 // Generate angle
                 double phi = random.flatShoot()*CLHEP::twopi;
 
                 // rotate and store secondaries
                 for (int j=0; j<nsec; ++j) {
                     const G4DynamicParticle* dp = result->GetSecondary(j)->GetParticle();
                     int daughterID = dp->GetParticleDefinition()->GetPDGEncoding();
 
                     // rotate around primary direction
                     curr4Mom = dp->Get4Momentum();
                     curr4Mom.rotate(phi, vectProj);
                     curr4Mom *= result->GetTrafoToLab();
 
                     // prompt 2-gamma decay for pi0, eta, eta'p
                     // don't have a charge so the closest daughter thing does not have to be done
                     if(111 == daughterID || 221 == daughterID || 331 == daughterID){
                         theBoost = curr4Mom.boostVector();
                         double e = 0.5*dp->GetParticleDefinition()->GetPDGMass();
                         double fi  = random.flatShoot()*CLHEP::twopi; 
                         double cth = 2.*random.flatShoot() - 1.0;
                         double sth = sqrt((1.0 - cth)*(1.0 + cth)); 
                         G4LorentzVector lv(e*sth*cos(fi),e*sth*sin(fi),e*cth,e);
                         lv.boost(theBoost);
 
                         // create secondaries
                         secondaries.emplace_back(new fastsim::Particle(22
                                 ,particle.position()
                                 ,XYZTLorentzVector(lv.px()/CLHEP::GeV, lv.py()/CLHEP::GeV, lv.pz()/CLHEP::GeV, lv.e()/CLHEP::GeV)));
 
                         curr4Mom -= lv;
                         if(curr4Mom.e() > theEnergyLimit) { 
                             secondaries.emplace_back(new fastsim::Particle(22
                                 ,particle.position()
                                 ,XYZTLorentzVector(curr4Mom.px()/CLHEP::GeV, curr4Mom.py()/CLHEP::GeV, curr4Mom.pz()/CLHEP::GeV, curr4Mom.e()/CLHEP::GeV)));
                         }
 
                         // The mother particle is destroyed
                         particle.momentum().SetXYZT(0., 0., 0., 0.);
                     }else{
                         // Copied from the original code, however I am not sure if all that is correct!
                         // Energy of particles increases during the interaction!
                         // std::cout << particle.momentum().e() << "(" << dynParticle->GetParticleDefinition()->GetPDGMass()/CLHEP::GeV << ") -> " << curr4Mom.e()/CLHEP::GeV << "(" << dp->GetParticleDefinition()->GetPDGMass()/CLHEP::GeV << ")" << std::endl;
 
                         if(curr4Mom.e() > theEnergyLimit + dp->GetParticleDefinition()->GetPDGMass()){
                             // Create secondary
                             secondaries.emplace_back(new fastsim::Particle(daughterID
                                 ,particle.position()
                                 ,XYZTLorentzVector(curr4Mom.px()/CLHEP::GeV, curr4Mom.py()/CLHEP::GeV, curr4Mom.pz()/CLHEP::GeV, curr4Mom.e()/CLHEP::GeV)));
 
                             // ClosestChargedDaughter thing for FastSim tracking
                             if(particle.charge() != 0 && std::abs(particle.charge()-dp->GetParticleDefinition()->GetPDGCharge()) < 1E-10){
                                 secondaries.back()->setMotherDeltaR(particle.momentum());
                                 secondaries.back()->setMotherPdgId(particle.getMotherDeltaR() == -1 ? particle.pdgId() : particle.getMotherPdgId());
                                 secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
                             }
 
                             // The mother particle is destroyed
                             particle.momentum().SetXYZT(0., 0., 0., 0.);
                         }
                     }
                 }
             }
 
             // Clear the final state
             result->Clear();
         }
     }
 }

Member Data Documentation

const double fastsim::NuclearInteractionFTF::corrfactors_el[numHadrons][npoints]

private

elastic interaction length corrections per particle and energy

Definition at line 195 of file NuclearInteractionFTF.cc.

Referenced by interact().

const double fastsim::NuclearInteractionFTF::corrfactors_inel[numHadrons][npoints]

private

inelastic interaction length corrections per particle and energy

Definition at line 161 of file NuclearInteractionFTF.cc.

Referenced by interact().

G4LorentzVector fastsim::NuclearInteractionFTF::curr4Mom

private

Definition at line 145 of file NuclearInteractionFTF.cc.

Referenced by interact().

int fastsim::NuclearInteractionFTF::currIdx

private

Index of particle in vector of Geant4 particles.

Definition at line 157 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

const G4ParticleDefinition* fastsim::NuclearInteractionFTF::currParticle

private

Definition at line 141 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

G4Track* fastsim::NuclearInteractionFTF::currTrack

private

Definition at line 140 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

G4Step* fastsim::NuclearInteractionFTF::dummyStep

private

Definition at line 139 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

size_t fastsim::NuclearInteractionFTF::index

private

Index for energy of particle (vectors parametrized as a function of energy)

Definition at line 158 of file NuclearInteractionFTF.cc.

Referenced by BeautifulSoup.PageElement::insert(), interact(), and NuclearInteractionFTF().

double fastsim::NuclearInteractionFTF::intLengthElastic

private

Elastic interaction length.

Definition at line 154 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

double fastsim::NuclearInteractionFTF::intLengthInelastic

private

Inelastic interaction length.

Definition at line 155 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

const int fastsim::NuclearInteractionFTF::npoints = 21

staticprivate

Number of steps to cover range of particle energies.

Definition at line 124 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

const double fastsim::NuclearInteractionFTF::nuclElLength[numHadrons]

private

Initial value:

= {
        9.248, 9.451, 11.545, 11.671, 32.081, 34.373, 34.373, 32.081, 0, 0, 15.739, 20.348, 15.739, 15.739, 0, 20.349, 0, 9.7514, 12.864, 15.836, 20.516, 15.836, 15.744, 0, 20.517, 20.44, 4.129, 6.0904, 4.5204, 4.5204
        }

elastic interaction length in Silicon at 1 TeV per particle

Definition at line 234 of file NuclearInteractionFTF.cc.

Referenced by interact().

const double fastsim::NuclearInteractionFTF::nuclInelLength[numHadrons]

private

Initial value:

= {
        4.5606, 4.4916, 5.7511, 5.7856, 6.797, 6.8373, 6.8171, 6.8171, 0, 0, 4.6926, 4.6926, 4.6926, 4.6926, 0, 4.6926, 4.6926, 4.3171, 4.3171, 4.6926, 4.6926, 4.6926, 4.6926, 0, 4.6926, 4.6926, 2.509, 2.9048, 2.5479, 2.5479
        }

inelastic interaction length in Silicon at 1 TeV per particle

Definition at line 229 of file NuclearInteractionFTF.cc.

Referenced by interact().

const int fastsim::NuclearInteractionFTF::numHadrons = 30

staticprivate

Number of G4 hadrons.

Definition at line 123 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().

G4Nucleus fastsim::NuclearInteractionFTF::targetNucleus

private

Definition at line 143 of file NuclearInteractionFTF.cc.

Referenced by interact(), and NuclearInteractionFTF().