#include <NuclearInteractionFTFSimulator.h>

Inheritance diagram for NuclearInteractionFTFSimulator:

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

	~NuclearInteractionFTFSimulator ()
	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 *)
	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

int	numHadrons

G4Nucleus	targetNucleus

G4ThreeVector	theBoost

G4GeneratorPrecompoundInterface *	theCascade

unsigned int	theDistAlgo

double	theDistCut

double	theEnergyLimit

std::vector< const G4ParticleDefinition * >	theG4Hadron

G4TheoFSGenerator *	theHadronicModel

std::vector< int >	theId

G4LundStringFragmentation *	theLund

std::vector< double >	theNuclIntLength

G4HadProjectile	theProjectile

G4ExcitedStringDecay *	theStringDecay

G4FTFModel *	theStringModel

G4ThreeVector	vectProj

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 36 of file NuclearInteractionFTFSimulator.h.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 63 of file NuclearInteractionFTFSimulator.cc.

References currIdx, currParticle, currTrack, dummyStep, GeV, i, numHadrons, targetNucleus, theBoost, theEnergyLimit, theG4Hadron, theHadronicModel, theId, theLund, theNuclIntLength, theStringDecay, theStringModel, and vectProj.

                   :
   theDistCut(distCut),
   distMin(1E99),
   theDistAlgo(distAlgo)
 {
   theEnergyLimit = 1*CLHEP::GeV;
   currIdx = 0;
 
   // 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);
 
   // Geant4 particles
   numHadrons = 30;
   theG4Hadron.resize(numHadrons,0);
   theG4Hadron[0] = G4Proton::Proton();
   theG4Hadron[1] = G4Neutron::Neutron();
   theG4Hadron[2] = G4PionPlus::PionPlus();
   theG4Hadron[3] = G4PionMinus::PionMinus();
   theG4Hadron[4] = G4AntiProton::AntiProton();
   theG4Hadron[5] = G4KaonPlus::KaonPlus();
   theG4Hadron[6] = G4KaonMinus::KaonMinus();
   theG4Hadron[7] = G4KaonZeroLong::KaonZeroLong();
   theG4Hadron[8] = G4KaonZeroShort::KaonZeroShort();
   theG4Hadron[9] = G4KaonZero::KaonZero();
   theG4Hadron[10]= G4AntiKaonZero::AntiKaonZero();
   theG4Hadron[11]= G4Lambda::Lambda();
   theG4Hadron[12]= G4OmegaMinus::OmegaMinus();
   theG4Hadron[13]= G4SigmaMinus::SigmaMinus();
   theG4Hadron[14]= G4SigmaPlus::SigmaPlus();
   theG4Hadron[15]= G4SigmaZero::SigmaZero();
   theG4Hadron[16]= G4XiMinus::XiMinus();
   theG4Hadron[17]= G4XiZero::XiZero();
   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();
 
   G4GenericIon::GenericIon();
   G4DecayPhysics decays;
   decays.ConstructParticle();  
   G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
   partTable->SetReadiness();
 
   // interaction length in units of radiation length 
   // computed for 5 GeV projectile energy, default value for K0_L
   theNuclIntLength.resize(numHadrons,6.46);
   theNuclIntLength[0] = 4.528;
   theNuclIntLength[1] = 4.524;
   theNuclIntLength[2] = 4.493;
   theNuclIntLength[3] = 4.493;
   theNuclIntLength[4] = 3.593;
   theNuclIntLength[5] = 7.154;
   theNuclIntLength[6] = 5.889;
   theNuclIntLength[11]= 4.986;
   theNuclIntLength[12]= 4.983;
   theNuclIntLength[13]= 4.986;
   theNuclIntLength[14]= 4.986;
   theNuclIntLength[15]= 4.986;
   theNuclIntLength[16]= 4.986;
   theNuclIntLength[17]= 4.986;
   theNuclIntLength[18]= 3.597;
   theNuclIntLength[19]= 3.608;
   theNuclIntLength[20]= 3.639;
   theNuclIntLength[21]= 3.613;
   theNuclIntLength[22]= 3.613;
   theNuclIntLength[23]= 3.613;
   theNuclIntLength[24]= 3.62;
   theNuclIntLength[25]= 3.62;
   theNuclIntLength[26]= 1.971;
   theNuclIntLength[27]= 2.301;
   theNuclIntLength[28]= 1.997;
   theNuclIntLength[29]= 1.997;
 
   // list of PDG codes
   theId.resize(numHadrons,0);
 
   // local objects
   currIdx = 0;
   currTrack = 0;
   currParticle = theG4Hadron[0];
   vectProj.set(0.0,0.0,1.0);  
   theBoost.set(0.0,0.0,1.0);  
 
   // fill projectile particle definitions
   dummyStep = new G4Step();
   dummyStep->SetPreStepPoint(new G4StepPoint());
   for(int i=0; i<numHadrons; ++i) {
     theId[i] = theG4Hadron[i]->GetPDGEncoding();
   }
 
   // target is always Silicon
   targetNucleus.SetParameters(28, 14);
 }

NuclearInteractionFTFSimulator::~NuclearInteractionFTFSimulator ( )

Default Destructor.

Definition at line 177 of file NuclearInteractionFTFSimulator.cc.

References theLund, theStringDecay, and theStringModel.

                                                                 {
 
   delete theStringDecay;
   delete theStringModel;
   delete theLund;
 }

Member Function Documentation

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

privatevirtual

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

Implements MaterialEffectsSimulator.

Definition at line 184 of file NuclearInteractionFTFSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, funct::cos(), curr4Mom, currIdx, currParticle, currTrack, dir, distMin, dummyStep, alignCSCRings::e, RandomEngineAndDistribution::flatShoot(), GeV, j, RawParticle::momentum(), numHadrons, phi, RawParticle::pid(), MaterialEffectsSimulator::radLengths, query::result, saveDaughter(), funct::sin(), mathSSE::sqrt(), targetNucleus, theBoost, MaterialEffectsSimulator::theClosestChargedDaughterId, theEnergyLimit, theG4Hadron, theHadronicModel, theId, theNuclIntLength, theProjectile, and vectProj.

 {
   //std::cout << "#### Primary " << Particle.pid() << " E(GeV)= " 
   //        << Particle.momentum().e() << std::endl;
 
   int thePid = Particle.pid(); 
   if(thePid != theId[currIdx]) {
     currParticle = 0;
     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 e = CLHEP::GeV*Particle.momentum().e();
   double mass = currParticle->GetPDGMass();
   /*
   std::cout << " Primary " <<  currParticle->GetParticleName() 
         << "  E(GeV)= " << e*fact << std::endl;
   */
   if(e <= theEnergyLimit + mass) { return; }
 
   double  currInteractionLength = -G4Log(random->flatShoot())*theNuclIntLength[currIdx]; 
   /*
   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 norm = 1.0/sqrt(px*px + py*py + pz*pz);
   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,e-mass);
   currTrack = new G4Track(dynParticle, 0.0, vectProj);
   currTrack->SetStep(dummyStep);
 
   theProjectile.Initialise(*currTrack); 
   delete currTrack;
 
   G4HadFinalState* 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;
       saveDaughter(Particle, curr4Mom, 22); 
     } else {
       saveDaughter(Particle, curr4Mom, thePid); 
     }
       }
     }
   }
 }

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

private

Definition at line 315 of file NuclearInteractionFTFSimulator.cc.

References RawParticle::charge(), HLT_25ns14e33_v1_cff::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 296 of file NuclearInteractionFTFSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, HLT_25ns14e33_v1_cff::distance, distanceToPrimary(), distMin, fact, customizeTrackingMonitorSeedNumber::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;
 }