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::ProducesCollector) 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] = {nullptr}
	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 101 of file NuclearInteractionFTF.cc.

Constructor & Destructor Documentation

◆ NuclearInteractionFTF()

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

Constructor.

Definition at line 298 of file NuclearInteractionFTF.cc.

References looper::cfg, currIdx, currParticle, currTrack, bJpsiMuMuTrigSettings_cff::decays, dummyStep, mps_fire::i, index, fastsim::initializeOnce, intLengthElastic, intLengthInelastic, BPhysicsValidation_cfi::Lambda, 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, []() {
     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);
 }

◆ ~NuclearInteractionFTF()

fastsim::NuclearInteractionFTF::~NuclearInteractionFTF ( )

override

Default destructor.

Definition at line 386 of file NuclearInteractionFTF.cc.

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

Member Function Documentation

◆ interact()

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 393 of file NuclearInteractionFTF.cc.

References funct::abs(), fastsim::Particle::charge(), corrfactors_el, corrfactors_inel, funct::cos(), DeadROC_duringRun::dir, Calorimetry_cff::dp, MillePedeFileConverter_cfg::e, StorageManager_cfg::e1, EcalCondDBWriter_cfi::Energy, RandomEngineAndDistribution::flatShoot(), fastsim::Particle::getMotherDeltaR(), fastsim::Particle::getMotherPdgId(), dqmiolumiharvest::j, EgHLTOffHistBins_cfi::mass, fastsim::Particle::momentum(), npoints, nuclElLength, nuclInelLength, numHadrons, fastsim::Particle::pdgId(), fastsim::Particle::position(), multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, mps_fire::result, fastsim::Particle::simTrackIndex(), funct::sin(), and mathSSE::sqrt().

                                                                                          {
   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

◆ corrfactors_el

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

private

elastic interaction length corrections per particle and energy

Definition at line 220 of file NuclearInteractionFTF.cc.

◆ corrfactors_inel

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

private

inelastic interaction length corrections per particle and energy

Definition at line 160 of file NuclearInteractionFTF.cc.

◆ curr4Mom

G4LorentzVector fastsim::NuclearInteractionFTF::curr4Mom

private

Definition at line 144 of file NuclearInteractionFTF.cc.

◆ currIdx

int fastsim::NuclearInteractionFTF::currIdx

private

Index of particle in vector of Geant4 particles.

Definition at line 156 of file NuclearInteractionFTF.cc.

Public Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ NuclearInteractionFTF()

◆ ~NuclearInteractionFTF()

Member Function Documentation

◆ interact()

Member Data Documentation

◆ corrfactors_el

◆ corrfactors_inel

◆ curr4Mom

◆ currIdx

◆ currParticle

◆ currTrack

◆ dummyStep

◆ index

◆ intLengthElastic

◆ intLengthInelastic

◆ npoints

◆ nuclElLength

◆ nuclInelLength

◆ numHadrons

◆ targetNucleus

◆ theBertiniCascade

◆ theBertiniLimit

◆ theBoost

◆ theCascade

◆ theDiffuseElastic

◆ theDistCut

◆ theEnergyLimit

◆ theG4Hadron

◆ theHadronicModel

◆ theId

◆ theLund

◆ theProjectile

◆ theStringDecay

◆ theStringModel

◆ vect

◆ vectProj