#include <NuclearInteractionFTFSimulator.h>

Inheritance diagram for NuclearInteractionFTFSimulator:

Public Member Functions
	NuclearInteractionFTFSimulator (unsigned int distAlgo, double distCut, double elimit, double eth)
	Constructor. More...

	~NuclearInteractionFTFSimulator () override
	Default Destructor. More...

Public Member Functions inherited from MaterialEffectsSimulator
RHEP_const_iter	beginDaughters () const
	Returns const iterator to the beginning of the daughters list. More...

int	closestDaughterId ()
	The id of the closest charged daughter (filled for nuclear interactions only) More...

double	eMass () const
	Electron mass in GeV/c2. More...

RHEP_const_iter	endDaughters () const
	Returns const iterator to the end of the daughters list. More...

double	excitE () const
	Mean excitation energy (in GeV) More...

	MaterialEffectsSimulator (double A=28.0855, double Z=14.0000, double density=2.329, double radLen=9.360)

unsigned	nDaughters () const
	Returns the number of daughters. More...

XYZVector	orthogonal (const XYZVector &) const
	A vector orthogonal to another one (because it's not in XYZTLorentzVector) More...

double	radLenIncm () const
	One radiation length in cm. More...

double	rho () const
	Density in g/cm3. More...

virtual void	save ()
	Used by NuclearInteractionSimulator to save last sampled event. More...

void	setNormalVector (const GlobalVector &normal)
	Sets the vector normal to the surface traversed. More...

double	theA () const
	A. More...

double	theZ () const
	Z. More...

void	updateState (ParticlePropagator &myTrack, double radlen, RandomEngineAndDistribution const *)
	Compute the material effect (calls the sub class) More...

virtual	~MaterialEffectsSimulator ()

Private Member Functions
void	compute (ParticlePropagator &Particle, RandomEngineAndDistribution const *) override
	Generate a nuclear interaction according to the probability that it happens. More...

double	distanceToPrimary (const RawParticle &Particle, const RawParticle &aDaughter) const

void	saveDaughter (ParticlePropagator &Particle, const G4LorentzVector &lv, int pdgid)

Private Attributes
G4LorentzVector	curr4Mom

int	currIdx

const G4ParticleDefinition *	currParticle

G4Track *	currTrack

double	distMin

G4Step *	dummyStep

size_t	index

double	intLengthElastic

double	intLengthInelastic

G4Nucleus	targetNucleus

G4CascadeInterface *	theBertiniCascade

double	theBertiniLimit

G4ThreeVector	theBoost

G4GeneratorPrecompoundInterface *	theCascade

G4DiffuseElastic *	theDiffuseElastic

unsigned int	theDistAlgo

double	theDistCut

double	theEnergyLimit

G4TheoFSGenerator *	theHadronicModel

G4LundStringFragmentation *	theLund

G4HadProjectile	theProjectile

G4ExcitedStringDecay *	theStringDecay

G4FTFModel *	theStringModel

G4PhysicsLogVector *	vect

G4ThreeVector	vectProj

Static Private Attributes
static const G4ParticleDefinition *	theG4Hadron [numHadrons] = { 0 }

static int	theId [numHadrons] = {0}

Additional Inherited Members
Public Types inherited from MaterialEffectsSimulator
typedef std::vector< RawParticle >::const_iterator	RHEP_const_iter

Protected Attributes inherited from MaterialEffectsSimulator
std::vector< RawParticle >	_theUpdatedState

double	A

double	density

double	radLen

double	radLengths

int	theClosestChargedDaughterId

GlobalVector	theNormalVector

double	Z

Detailed Description

Definition at line 40 of file NuclearInteractionFTFSimulator.h.

Constructor & Destructor Documentation

NuclearInteractionFTFSimulator::NuclearInteractionFTFSimulator	(	unsigned int	distAlgo,
		double	distCut,
		double	elimit,
		double	eth
	)

Constructor.

Definition at line 153 of file NuclearInteractionFTFSimulator.cc.

References currIdx, currParticle, currTrack, dummyStep, mps_fire::i, index, initializeOnce, intLengthElastic, intLengthInelastic, BPhysicsValidation_cfi::Lambda, MeV, npoints, numHadrons, phase2TrackerDigitizer_cfi::SigmaZero, targetNucleus, theBertiniCascade, theDiffuseElastic, theEnergyLimit, theG4Hadron, theHadronicModel, theId, theLund, theStringDecay, theStringModel, and vect.

                                                                     :
   curr4Mom(0.,0.,0.,0.),
   vectProj(0.,0.,1.),
   theBoost(0.,0.,0.),
   theBertiniLimit(elimit),
   theEnergyLimit(eth),
   theDistCut(distCut),
   distMin(1E99),
   theDistAlgo(distAlgo)
 {
   // 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);
 }

NuclearInteractionFTFSimulator::~NuclearInteractionFTFSimulator ( )

override

Default Destructor.

Definition at line 245 of file NuclearInteractionFTFSimulator.cc.

References theLund, theStringDecay, theStringModel, and vect.

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

Member Function Documentation

void NuclearInteractionFTFSimulator::compute	(	ParticlePropagator &	Particle,
		RandomEngineAndDistribution const *	random
	)

overrideprivatevirtual

Generate a nuclear interaction according to the probability that it happens.

Implements MaterialEffectsSimulator.

Definition at line 253 of file NuclearInteractionFTFSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, corrfactors_el, corrfactors_inel, funct::cos(), curr4Mom, currIdx, currParticle, currTrack, dir, distMin, reco::dp, dummyStep, MillePedeFileConverter_cfg::e, reco::e1, reco::e2, RandomEngineAndDistribution::flatShoot(), GeV, index, intLengthElastic, intLengthInelastic, ResonanceBuilder::mass, RawParticle::momentum(), npoints, nuclElLength, nuclInelLength, numHadrons, phi, RawParticle::pid(), MaterialEffectsSimulator::radLengths, mps_fire::result, saveDaughter(), funct::sin(), mathSSE::sqrt(), targetNucleus, theBertiniCascade, theBertiniLimit, theBoost, MaterialEffectsSimulator::theClosestChargedDaughterId, theDiffuseElastic, theDistCut, theEnergyLimit, theG4Hadron, theHadronicModel, theId, theProjectile, vect, and vectProj.

 {
   //std::cout << "#### Primary " << Particle.pid() << " E(GeV)= " 
   //        << Particle.momentum().e() << std::endl;
 
   int thePid = Particle.pid(); 
   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;
       }
     }
   }
   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));
   }
   /*
   std::cout << " Primary " <<  currParticle->GetParticleName() 
         << "  E(GeV)= " << e*fact << std::endl;
   */
 
   double currInteractionLength = -G4Log(random->flatShoot())*intLengthElastic*intLengthInelastic
     /(intLengthElastic + intLengthInelastic); 
   /*
   std::cout << "*NuclearInteractionFTFSimulator::compute: R(X0)= " << radLengths
         << " Rnuc(X0)= " << theNuclIntLength[currIdx] << "  IntLength(X0)= " 
             << currInteractionLength << std::endl;
   */
   // 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);
   /*
   std::cout << " Primary " <<  currParticle->GetParticleName() 
         << "  E(GeV)= " << e*fact << "  P(GeV/c)= (" 
         << px << " " << py << " " << pz << ")" << std::endl;
   */
 
   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 + intLengthInelastic) > intLengthElastic) {
 
     result = theDiffuseElastic->ApplyYourself(theProjectile, targetNucleus);
     G4ThreeVector dirnew = result->GetMomentumChange().unit();
     double cost = (dir*dirnew);
     double sint = std::sqrt((1. - cost)*(1. + cost));
 
     curr4Mom.set(ptot*dirnew.x(),ptot*dirnew.y(),ptot*dirnew.z(),Particle.momentum().e());
 
     // Always create a daughter if the kink is large engough 
     if (sint > theDistCut) { 
       saveDaughter(Particle, curr4Mom, thePid); 
     } else {
       Particle.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(0 < nsec) {
 
     result->SetTrafoToLab(theProjectile.GetTrafoToLab());
     _theUpdatedState.clear();
 
     //std::cout << "   " << nsec << " secondaries" << std::endl;
     // Generate angle
     double phi = random->flatShoot()*CLHEP::twopi;
     theClosestChargedDaughterId = -1;
     distMin = 1e99;
 
     // rotate and store secondaries
     for (int j=0; j<nsec; ++j) {
 
       const G4DynamicParticle* dp = result->GetSecondary(j)->GetParticle();
       int thePid = dp->GetParticleDefinition()->GetPDGEncoding();
 
       // rotate around primary direction
       curr4Mom = dp->Get4Momentum();
       curr4Mom.rotate(phi, vectProj);
       curr4Mom *= result->GetTrafoToLab();
       /*
         std::cout << j << ". " << dp->GetParticleDefinition()->GetParticleName() 
         << "  " << thePid
         << "  " << curr4Mom*fact << std::endl;
       */
       // prompt 2-gamma decay for pi0, eta, eta'p
       if(111 == thePid || 221 == thePid || 331 == thePid) {
         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);
         saveDaughter(Particle, lv, 22); 
         curr4Mom -= lv;
         if(curr4Mom.e() > theEnergyLimit) { 
           saveDaughter(Particle, curr4Mom, 22); 
         } 
       } else {
         if(curr4Mom.e() > theEnergyLimit + dp->GetParticleDefinition()->GetPDGMass()) { 
           saveDaughter(Particle, curr4Mom, thePid); 
         }
       }
     }
     result->Clear();
       }
     }
   }
 }

double NuclearInteractionFTFSimulator::distanceToPrimary	(	const RawParticle &	Particle,
		const RawParticle &	aDaughter
	)		const

private

Definition at line 437 of file NuclearInteractionFTFSimulator.cc.

References RawParticle::charge(), SoftLeptonByDistance_cfi::distance, and theDistAlgo.

Referenced by saveDaughter().

 {
   double distance = 2E99;
   // Compute the distance only for charged primaries
   if ( fabs(Particle.charge()) > 1E-6 ) { 
 
     // The secondary must have the same charge
     double chargeDiff = fabs(aDaughter.charge()-Particle.charge());
     if ( fabs(chargeDiff) < 1E-6 ) {
 
       // Here are two distance definitions * to be tuned *
       switch ( theDistAlgo ) { 
     
       case 1:
     // sin(theta12)
     distance = (aDaughter.Vect().Unit().Cross(Particle.Vect().Unit())).R();
     break;
     
       case 2: 
     // sin(theta12) * p1/p2
     distance = (aDaughter.Vect().Cross(Particle.Vect())).R()
       /aDaughter.Vect().Mag2();
     break;
     
       default:
     // Should not happen
     break;  
       }
     }
   } 
   return distance;
 }

void NuclearInteractionFTFSimulator::saveDaughter	(	ParticlePropagator &	Particle,
		const G4LorentzVector &	lv,
		int	pdgid
	)

private

Definition at line 418 of file NuclearInteractionFTFSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, SoftLeptonByDistance_cfi::distance, distanceToPrimary(), distMin, fact, training_settings::idx, MaterialEffectsSimulator::theClosestChargedDaughterId, and theDistCut.

Referenced by compute().

 {
   unsigned int idx = _theUpdatedState.size();   
   _theUpdatedState.push_back(Particle);
   _theUpdatedState[idx].SetXYZT(lv.px()*fact,lv.py()*fact,lv.pz()*fact,lv.e()*fact);
   _theUpdatedState[idx].setID(pdgid);
 
   // Store the closest daughter index (for later tracking purposes, so charged particles only) 
   double distance = distanceToPrimary(Particle,_theUpdatedState[idx]);
   // Find the closest daughter, if closer than a given upper limit.
   if ( distance < distMin && distance < theDistCut ) {
     distMin = distance;
     theClosestChargedDaughterId = idx;
   }
   // std::cout << _theUpdatedState[idx] << std::endl;
 }