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
00030
00031
00032
00033
00034 #include "G4DynamicParticle.hh"
00035 #include "G4TheoFSGenerator.hh"
00036 #include "G4ReactionProductVector.hh"
00037 #include "G4ReactionProduct.hh"
00038 #include "G4IonTable.hh"
00039
00040 G4TheoFSGenerator::G4TheoFSGenerator(const G4String& name)
00041 : G4HadronicInteraction(name)
00042 , theTransport(0), theHighEnergyGenerator(0)
00043 , theQuasielastic(0), theProjectileDiffraction(0)
00044 {
00045 theParticleChange = new G4HadFinalState;
00046 }
00047
00048 G4TheoFSGenerator::~G4TheoFSGenerator()
00049 {
00050 delete theParticleChange;
00051 }
00052
00053 void G4TheoFSGenerator::ModelDescription(std::ostream& outFile) const
00054 {
00055 outFile << GetModelName() <<" consists of a " << theHighEnergyGenerator->GetModelName()
00056 << " string model and of "
00057 << ".\n"
00058 << "The string model simulates the interaction of\n"
00059 << "an incident hadron with a nucleus, forming \n"
00060 << "excited strings, decays these strings into hadrons,\n"
00061 << "and leaves an excited nucleus.\n"
00062 << "The string model:\n";
00063 theHighEnergyGenerator->ModelDescription(outFile);
00064
00065 }
00066
00067
00068 G4HadFinalState * G4TheoFSGenerator::ApplyYourself(const G4HadProjectile & thePrimary, G4Nucleus &theNucleus)
00069 {
00070
00071 theParticleChange->Clear();
00072 theParticleChange->SetStatusChange(stopAndKill);
00073
00074
00075
00076
00077 G4DynamicParticle aTemp(const_cast<G4ParticleDefinition *>(thePrimary.GetDefinition()),
00078 thePrimary.Get4Momentum().vect());
00079 const G4DynamicParticle * aPart = &aTemp;
00080
00081 if ( theQuasielastic ) {
00082
00083 if ( theQuasielastic->GetFraction(theNucleus, *aPart) > G4UniformRand() )
00084 {
00085
00086 G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, *aPart);
00087
00088 if (result)
00089 {
00090 for(unsigned int i=0; i<result->size(); i++)
00091 {
00092 G4DynamicParticle * aNew =
00093 new G4DynamicParticle(result->operator[](i)->GetDefinition(),
00094 result->operator[](i)->Get4Momentum().e(),
00095 result->operator[](i)->Get4Momentum().vect());
00096 theParticleChange->AddSecondary(aNew);
00097 delete result->operator[](i);
00098 }
00099 delete result;
00100
00101 } else
00102 {
00103 theParticleChange->SetStatusChange(isAlive);
00104 theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy());
00105 theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit());
00106
00107 }
00108 return theParticleChange;
00109 }
00110 }
00111 if ( theProjectileDiffraction) {
00112
00113 if ( theProjectileDiffraction->GetFraction(theNucleus, *aPart) > G4UniformRand() )
00114
00115
00116
00117 {
00118
00119 G4KineticTrackVector *result= theProjectileDiffraction->Scatter(theNucleus, *aPart);
00120
00121 if (result)
00122 {
00123 for(unsigned int i=0; i<result->size(); i++)
00124 {
00125 G4DynamicParticle * aNew =
00126 new G4DynamicParticle(result->operator[](i)->GetDefinition(),
00127 result->operator[](i)->Get4Momentum().e(),
00128 result->operator[](i)->Get4Momentum().vect());
00129 theParticleChange->AddSecondary(aNew);
00130 delete result->operator[](i);
00131 }
00132 delete result;
00133
00134 } else
00135 {
00136 theParticleChange->SetStatusChange(isAlive);
00137 theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy());
00138 theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit());
00139
00140 }
00141 return theParticleChange;
00142 }
00143 }
00144 G4KineticTrackVector * theInitialResult =
00145 theHighEnergyGenerator->Scatter(theNucleus, *aPart);
00146
00147
00148 #ifdef DEBUG_initial_result
00149 G4double E_out(0);
00150 G4IonTable * ionTable=G4ParticleTable::GetParticleTable()->GetIonTable();
00151 std::vector<G4KineticTrack *>::iterator ir_iter;
00152 for(ir_iter=theInitialResult->begin(); ir_iter!=theInitialResult->end(); ir_iter++)
00153 {
00154
00155 E_out += (*ir_iter)->Get4Momentum().e();
00156 }
00157 G4double init_mass= ionTable->GetIonMass(theNucleus.GetZ_asInt(),theNucleus.GetA_asInt());
00158 G4double init_E=aPart->Get4Momentum().e();
00159
00160 const std::vector<G4Nucleon> & thy = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
00161 G4int resZ(0),resA(0);
00162 G4double delta_m(0);
00163 for(size_t them=0; them<thy.size(); them++)
00164 {
00165 if(thy[them].AreYouHit()) {
00166 ++resA;
00167 if ( thy[them].GetDefinition() == G4Proton::Proton() ) {
00168 ++resZ;
00169 delta_m +=G4Proton::Proton()->GetPDGMass();
00170 } else {
00171 delta_m +=G4Neutron::Neutron()->GetPDGMass();
00172 }
00173 }
00174 }
00175 G4double final_mass(0);
00176 if ( theNucleus.GetA_asInt() ) {
00177 final_mass=ionTable->GetIonMass(theNucleus.GetZ_asInt()-resZ,theNucleus.GetA_asInt()- resA);
00178 }
00179 G4double E_excit=init_mass + init_E - final_mass - E_out;
00180 G4cout << " Corrected delta mass " << init_mass - final_mass - delta_m << G4endl;
00181 G4cout << "initial E, mass = " << init_E << ", " << init_mass << G4endl;
00182 G4cout << " final E, mass = " << E_out <<", " << final_mass << " excitation_E " << E_excit << G4endl;
00183 #endif
00184
00185 G4ReactionProductVector * theTransportResult = NULL;
00186 G4int hitCount = 0;
00187 const std::vector<G4Nucleon>& they = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
00188 for(size_t them=0; them<they.size(); them++)
00189 {
00190 if(they[them].AreYouHit()) hitCount ++;
00191 }
00192 if(hitCount != theHighEnergyGenerator->GetWoundedNucleus()->GetMassNumber() )
00193 {
00194 theTransport->SetPrimaryProjectile(thePrimary);
00195 theTransportResult =
00196 theTransport->Propagate(theInitialResult, theHighEnergyGenerator->GetWoundedNucleus());
00197 if ( !theTransportResult ) {
00198 G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
00199 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
00200 }
00201 }
00202 else
00203 {
00204 theTransportResult = theDecay.Propagate(theInitialResult, theHighEnergyGenerator->GetWoundedNucleus());
00205 if ( !theTransportResult ) {
00206 G4cout << "G4TheoFSGenerator: null ptr from decay propagate " << G4endl;
00207 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from decay propagate");
00208 }
00209 }
00210
00211
00212 unsigned int i;
00213 for(i=0; i<theTransportResult->size(); i++)
00214 {
00215 G4DynamicParticle * aNew =
00216 new G4DynamicParticle(theTransportResult->operator[](i)->GetDefinition(),
00217 theTransportResult->operator[](i)->GetTotalEnergy(),
00218 theTransportResult->operator[](i)->GetMomentum());
00219
00220 theParticleChange->AddSecondary(aNew);
00221 delete theTransportResult->operator[](i);
00222 }
00223
00224
00225 delete theTransportResult;
00226
00227
00228 return theParticleChange;
00229 }
00230
00231 std::pair<G4double, G4double> G4TheoFSGenerator::GetEnergyMomentumCheckLevels() const
00232 {
00233 if ( theHighEnergyGenerator ) {
00234 return theHighEnergyGenerator->GetEnergyMomentumCheckLevels();
00235 } else {
00236 return std::pair<G4double, G4double>(DBL_MAX, DBL_MAX);
00237 }
00238 }