00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <iostream>
00030 #include <signal.h>
00031
00032 #include "G4RPGPionSuppression.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "Randomize.hh"
00035 #include "G4HadReentrentException.hh"
00036
00037 G4RPGPionSuppression::G4RPGPionSuppression()
00038 : G4RPGReaction() {}
00039
00040
00041 G4bool G4RPGPionSuppression::
00042 ReactionStage(const G4HadProjectile* ,
00043 G4ReactionProduct& modifiedOriginal,
00044 G4bool& incidentHasChanged,
00045 const G4DynamicParticle* ,
00046 G4ReactionProduct& targetParticle,
00047 G4bool& targetHasChanged,
00048 const G4Nucleus& targetNucleus,
00049 G4ReactionProduct& currentParticle,
00050 G4FastVector<G4ReactionProduct,256>& vec,
00051 G4int& vecLen,
00052 G4bool ,
00053 G4ReactionProduct& )
00054 {
00055
00056
00057
00058
00059 G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00060 G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00061 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00062 G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
00063 2.0*targetMass*etOriginal );
00064 G4double eAvailable = cmEnergy - mOriginal - targetMass;
00065 for (G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;
00066
00067 const G4double atomicWeight = targetNucleus.GetA_asInt();
00068 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00069 const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
00070
00071 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00072 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00073 G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
00074 G4ParticleDefinition *aProton = G4Proton::Proton();
00075 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00076 G4double piMass = aPiPlus->GetPDGMass()/GeV;
00077 G4double nucleonMass = aNeutron->GetPDGMass()/GeV;
00078
00079 const G4bool antiTest =
00080 modifiedOriginal.GetDefinition() != G4AntiProton::AntiProton() &&
00081 modifiedOriginal.GetDefinition() != G4AntiNeutron::AntiNeutron() &&
00082 modifiedOriginal.GetDefinition() != G4AntiLambda::AntiLambda() &&
00083 modifiedOriginal.GetDefinition() != G4AntiSigmaPlus::AntiSigmaPlus() &&
00084 modifiedOriginal.GetDefinition() != G4AntiSigmaMinus::AntiSigmaMinus() &&
00085 modifiedOriginal.GetDefinition() != G4AntiXiZero::AntiXiZero() &&
00086 modifiedOriginal.GetDefinition() != G4AntiXiMinus::AntiXiMinus() &&
00087 modifiedOriginal.GetDefinition() != G4AntiOmegaMinus::AntiOmegaMinus();
00088
00089 if( antiTest && (
00090 currentParticle.GetDefinition() == aPiPlus ||
00091 currentParticle.GetDefinition() == aPiZero ||
00092 currentParticle.GetDefinition() == aPiMinus ) &&
00093 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
00094 ( G4UniformRand() <= atomicWeight/300.0 ) )
00095 {
00096 if (eAvailable > nucleonMass - piMass) {
00097 if( G4UniformRand() > atomicNumber/atomicWeight )
00098 currentParticle.SetDefinitionAndUpdateE( aNeutron );
00099 else
00100 currentParticle.SetDefinitionAndUpdateE( aProton );
00101 incidentHasChanged = true;
00102 }
00103 }
00104
00105 if( antiTest && (
00106 targetParticle.GetDefinition() == aPiPlus ||
00107 targetParticle.GetDefinition() == aPiZero ||
00108 targetParticle.GetDefinition() == aPiMinus ) &&
00109 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
00110 ( G4UniformRand() <= atomicWeight/300.0 ) )
00111 {
00112 if (eAvailable > nucleonMass - piMass) {
00113 if( G4UniformRand() > atomicNumber/atomicWeight )
00114 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00115 else
00116 targetParticle.SetDefinitionAndUpdateE( aProton );
00117 targetHasChanged = true;
00118 }
00119 }
00120
00121 for( G4int i=0; i<vecLen; ++i )
00122 {
00123 if( antiTest && (
00124 vec[i]->GetDefinition() == aPiPlus ||
00125 vec[i]->GetDefinition() == aPiZero ||
00126 vec[i]->GetDefinition() == aPiMinus ) &&
00127 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
00128 ( G4UniformRand() <= atomicWeight/300.0 ) )
00129 {
00130 if (eAvailable > nucleonMass - piMass) {
00131 if( G4UniformRand() > atomicNumber/atomicWeight )
00132 vec[i]->SetDefinitionAndUpdateE( aNeutron );
00133 else
00134 vec[i]->SetDefinitionAndUpdateE( aProton );
00135 }
00136 }
00137 }
00138
00139 return true;
00140 }
00141
00142
00143