Geant4: G4ChargeExchange Class Reference

G4ChargeExchange::G4ChargeExchange ( )

Definition at line 48 of file G4ChargeExchange.cc.

References G4Alpha::Alpha(), G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4Deuteron::Deuteron(), G4He3::He3(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4Neutron::Neutron(), G4OmegaMinus::OmegaMinus(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), G4HadronicInteraction::SetMaxEnergy(), G4HadronicInteraction::SetMinEnergy(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4Triton::Triton(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().

00048                                    : G4HadronicInteraction("Charge Exchange")
00049 {
00050   SetMinEnergy( 0.0*GeV );
00051   SetMaxEnergy( 100.*TeV );
00052 
00053   lowEnergyRecoilLimit = 100.*keV;  
00054   lowestEnergyLimit    = 1.*MeV;  
00055 
00056   theProton   = G4Proton::Proton();
00057   theNeutron  = G4Neutron::Neutron();
00058   theAProton  = G4AntiProton::AntiProton();
00059   theANeutron = G4AntiNeutron::AntiNeutron();
00060   thePiPlus   = G4PionPlus::PionPlus();
00061   thePiMinus  = G4PionMinus::PionMinus();
00062   thePiZero   = G4PionZero::PionZero();
00063   theKPlus    = G4KaonPlus::KaonPlus();
00064   theKMinus   = G4KaonMinus::KaonMinus();
00065   theK0S      = G4KaonZeroShort::KaonZeroShort();
00066   theK0L      = G4KaonZeroLong::KaonZeroLong();
00067   theL        = G4Lambda::Lambda();
00068   theAntiL    = G4AntiLambda::AntiLambda();
00069   theSPlus    = G4SigmaPlus::SigmaPlus();
00070   theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
00071   theSMinus   = G4SigmaMinus::SigmaMinus();
00072   theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
00073   theS0       = G4SigmaZero::SigmaZero();
00074   theAS0      = G4AntiSigmaZero::AntiSigmaZero();
00075   theXiMinus  = G4XiMinus::XiMinus();
00076   theXi0      = G4XiZero::XiZero();
00077   theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00078   theAXi0     = G4AntiXiZero::AntiXiZero();
00079   theOmega    = G4OmegaMinus::OmegaMinus();
00080   theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
00081   theD        = G4Deuteron::Deuteron();
00082   theT        = G4Triton::Triton();
00083   theA        = G4Alpha::Alpha();
00084   theHe3      = G4He3::He3();
00085 }

G4ChargeExchange::~G4ChargeExchange ( ) [virtual]

Definition at line 87 of file G4ChargeExchange.cc.

00088 {}

G4HadFinalState * G4ChargeExchange::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	targetNucleus
	)			`[virtual]`

Implements G4HadronicInteraction.

Definition at line 90 of file G4ChargeExchange.cc.

References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4HadProjectile::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), SampleT(), G4DynamicParticle::SetDefinition(), G4HadFinalState::SetEnergyChange(), G4DynamicParticle::SetKineticEnergy(), G4HadFinalState::SetLocalEnergyDeposit(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HadronicInteraction::theParticleChange, and G4HadronicInteraction::verboseLevel.

00092 {
00093   theParticleChange.Clear();
00094   const G4HadProjectile* aParticle = &aTrack;
00095   G4double ekin = aParticle->GetKineticEnergy();
00096 
00097   G4int A = targetNucleus.GetA_asInt();
00098   G4int Z = targetNucleus.GetZ_asInt();
00099 
00100   if(ekin <= lowestEnergyLimit || A < 3) {
00101     theParticleChange.SetEnergyChange(ekin);
00102     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00103     return &theParticleChange;
00104   }
00105 
00106   G4double plab = aParticle->GetTotalMomentum();
00107 
00108   if (verboseLevel > 1)
00109     G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
00110            << plab/GeV << " GeV/c "
00111            << " ekin(MeV) = " << ekin/MeV << "  "
00112            << aParticle->GetDefinition()->GetParticleName() << G4endl;
00113 
00114   // Scattered particle referred to axis of incident particle
00115   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00116 
00117   G4int N = A - Z;
00118   G4int projPDG = theParticle->GetPDGEncoding();
00119   if (verboseLevel > 1)
00120     G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
00121            << " PDGcode= " << projPDG << " on nucleus Z= " << Z
00122            << " A= " << A << " N= " << N
00123            << G4endl;
00124 
00125   G4ParticleDefinition * theDef = 0;
00126 
00127   G4double mass2 = G4NucleiProperties::GetNuclearMass((G4double)A, (G4double)Z);
00128   G4LorentzVector lv1 = aParticle->Get4Momentum();
00129   G4LorentzVector lv0(0.0,0.0,0.0,mass2);
00130 
00131   G4LorentzVector lv  = lv0 + lv1;
00132   G4ThreeVector bst = lv.boostVector();
00133   lv1.boost(-bst);
00134   lv0.boost(-bst);
00135 
00136   // Sample final particles
00137   G4bool theHyperon = false;
00138   G4ParticleDefinition* theRecoil = 0;
00139   G4ParticleDefinition* theSecondary = 0;
00140 
00141   if(theParticle == theProton) {
00142     theSecondary = theNeutron;
00143     Z++;
00144   } else if(theParticle == theNeutron) {
00145     theSecondary = theProton;
00146     Z--;
00147   } else if(theParticle == thePiPlus) {
00148     theSecondary = thePiZero;
00149     Z++;
00150   } else if(theParticle == thePiMinus) {
00151     theSecondary = thePiZero;
00152     Z--;
00153   } else if(theParticle == theKPlus) {
00154     if(G4UniformRand()<0.5) theSecondary = theK0S;
00155     else  theSecondary = theK0L;
00156     Z++;
00157   } else if(theParticle == theKMinus) {
00158     if(G4UniformRand()<0.5) theSecondary = theK0S;
00159     else  theSecondary = theK0L;
00160     Z--;
00161   } else if(theParticle == theK0S || theParticle == theK0L) {
00162     if(G4UniformRand()*A < G4double(Z)) {
00163       theSecondary = theKPlus;
00164       Z--;
00165     } else {
00166       theSecondary = theKMinus;
00167       Z++;
00168     }
00169   } else if(theParticle == theANeutron) {
00170     theSecondary = theAProton;
00171     Z++;
00172   } else if(theParticle == theAProton) {
00173     theSecondary = theANeutron;
00174     Z--;
00175   } else if(theParticle == theL) {
00176     G4double x = G4UniformRand();
00177     if(G4UniformRand()*A < G4double(Z)) {
00178       if(x < 0.2) {
00179         theSecondary = theS0;
00180       } else if (x < 0.4) {
00181         theSecondary = theSPlus;
00182         Z--;
00183       } else if (x < 0.6) {
00184         theSecondary = theProton;
00185         theRecoil = theL;
00186         theHyperon = true;
00187         A--;
00188       } else if (x < 0.8) {
00189         theSecondary = theProton;
00190         theRecoil = theS0;
00191         theHyperon = true;
00192         A--;
00193       } else {
00194         theSecondary = theNeutron;
00195         theRecoil = theSPlus;
00196         theHyperon = true;
00197         A--;
00198       }
00199     } else {
00200       if(x < 0.2) {
00201         theSecondary = theS0;
00202       } else if (x < 0.4) {
00203         theSecondary = theSMinus;
00204         Z++;
00205       } else if (x < 0.6) {
00206         theSecondary = theNeutron;
00207         theRecoil = theL;
00208         A--;
00209         theHyperon = true;
00210       } else if (x < 0.8) {
00211         theSecondary = theNeutron;
00212         theRecoil = theS0;
00213         theHyperon = true;
00214         A--;
00215       } else {
00216         theSecondary = theProton;
00217         theRecoil = theSMinus;
00218         theHyperon = true;
00219         A--;
00220       }
00221     }
00222   }
00223 
00224   if (Z == 1 && A == 2) theDef = theD;
00225   else if (Z == 1 && A == 3) theDef = theT;
00226   else if (Z == 2 && A == 3) theDef = theHe3;
00227   else if (Z == 2 && A == 4) theDef = theA;
00228   else {
00229     theDef = 
00230       G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
00231   }
00232   if(!theSecondary) { return &theParticleChange; }
00233 
00234   G4double m11 = theSecondary->GetPDGMass();
00235   G4double m21 = theDef->GetPDGMass();
00236   if(theRecoil)  { m21 += theRecoil->GetPDGMass(); }
00237   else           { theRecoil = theDef; }
00238 
00239   G4double etot = lv0.e() + lv1.e();
00240 
00241   // kinematiacally impossible
00242   if(etot < m11 + m21) {
00243     theParticleChange.SetEnergyChange(ekin);
00244     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00245     return &theParticleChange;
00246   }
00247 
00248   G4ThreeVector p1 = lv1.vect();
00249   G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
00250   //  G4double e2 = etot - e1;
00251   G4double ptot = std::sqrt(e1*e1 - m11*m11);
00252 
00253   G4double tmax = 4.0*ptot*ptot;
00254   G4double g2 = GeV*GeV; 
00255 
00256   G4double t = g2*SampleT(tmax/g2, A);
00257 
00258   if(verboseLevel>1)
00259     G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
00260            << " ptot= " << ptot << G4endl;
00261 
00262   // Sampling in CM system
00263   G4double phi  = G4UniformRand()*twopi;
00264   G4double cost = 1. - 2.0*t/tmax;
00265   if(std::abs(cost) > 1.0) cost = 1.0;
00266   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
00267 
00268   //if (verboseLevel > 1)
00269   //  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
00270 
00271   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
00272   v1 *= ptot;
00273   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
00274   G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
00275 
00276   nlv0.boost(bst);
00277   nlv1.boost(bst);
00278 
00279   theParticleChange.SetStatusChange(stopAndKill);
00280   theParticleChange.SetEnergyChange(0.0);
00281   G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
00282   theParticleChange.AddSecondary(aSec);
00283 
00284   G4double erec =  nlv0.e() - m21;
00285 
00286   //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl;  
00287 
00288   if(theHyperon) {
00289     theParticleChange.SetLocalEnergyDeposit(erec);
00290     aSec = new G4DynamicParticle();
00291     aSec->SetDefinition(theRecoil);
00292     aSec->SetKineticEnergy(0.0);
00293   } else if(erec > lowEnergyRecoilLimit) {
00294     aSec = new G4DynamicParticle(theRecoil, nlv0);
00295     theParticleChange.AddSecondary(aSec);
00296   } else {
00297     if(erec < 0.0) erec = 0.0;
00298     theParticleChange.SetLocalEnergyDeposit(erec);
00299   }
00300   return &theParticleChange;
00301 }

G4double G4ChargeExchange::SampleT	(	G4double	p,
		G4double	A
	)

Definition at line 303 of file G4ChargeExchange.cc.

References G4UniformRand.

Referenced by ApplyYourself().

00304 {
00305   G4double aa, bb, cc, dd;
00306   if (A <= 62.) {
00307     aa = std::pow(A, 1.63);
00308     bb = 14.5*std::pow(A, 0.66);
00309     cc = 1.4*std::pow(A, 0.33);
00310     dd = 10.;
00311   } else {
00312     aa = std::pow(A, 1.33);
00313     bb = 60.*std::pow(A, 0.33);
00314     cc = 0.4*std::pow(A, 0.40);
00315     dd = 10.;
00316   }
00317   G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;
00318   G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;
00319   
00320   G4double t;
00321   G4double y = bb;
00322   if(G4UniformRand()*(x1 + x2) < x2) y = dd;
00323 
00324   do {t = -std::log(G4UniformRand())/y;} while (t > tmax);
00325 
00326   return t;
00327 }

void G4ChargeExchange::SetLowestEnergyLimit ( G4double value ) [inline]

Definition at line 112 of file G4ChargeExchange.hh.

00113 {
00114   lowestEnergyLimit = value;
00115 }

void G4ChargeExchange::SetRecoilKinEnergyLimit ( G4double value ) [inline]

Definition at line 107 of file G4ChargeExchange.hh.

00108 {
00109   lowEnergyRecoilLimit = value;
00110 }


Public Member Functions
	G4ChargeExchange ()
virtual	~G4ChargeExchange ()
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void	SetLowestEnergyLimit (G4double value)
void	SetRecoilKinEnergyLimit (G4double value)
G4double	SampleT (G4double p, G4double A)

G4ChargeExchange Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation