#include <G4INCLXXInterface.hh>
Inheritance diagram for G4INCLXXInterface:
Public Member Functions | |
G4INCLXXInterface (const G4String &name="INCL++ cascade with G4ExcitationHandler") | |
~G4INCLXXInterface () | |
G4int | operator== (G4INCLXXInterface &right) |
G4int | operator!= (G4INCLXXInterface &right) |
G4ReactionProductVector * | Propagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus) |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus) |
void | DeleteModel () |
Interface for INCL++. This interface handles basic hadron bullet particles (protons, neutrons, pions), as well as light ions.
Example usage in case of protons:
G4INCLXXInterface* inclModel = new G4INCLXXInterface; inclModel -> SetMinEnergy(0.0 * MeV); // Set the energy limits inclModel -> SetMaxEnergy(3.0 * GeV); G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess(); G4ProtonInelasticCrossSection* protonInelasticCrossSection = new G4ProtonInelasticCrossSection(); protonInelasticProcess -> RegisterMe(inclModel); protonInelasticProcess -> AddDataSet(protonInelasticCrossSection); particle = G4Proton::Proton(); processManager = particle -> GetProcessManager(); processManager -> AddDiscreteProcess(protonInelasticProcess);
Definition at line 91 of file G4INCLXXInterface.hh.
G4INCLXXInterface::G4INCLXXInterface | ( | const G4String & | name = "INCL++ cascade with G4ExcitationHandler" |
) |
Definition at line 49 of file G4INCLXXInterface.cc.
References G4INCLXXInterfaceStore::EmitWarning().
00049 : 00050 G4VIntraNuclearTransportModel(nam), 00051 theINCLModel(NULL), 00052 theInterfaceStore(G4INCLXXInterfaceStore::GetInstance()), 00053 complainedAboutBackupModel(false) 00054 { 00055 // Use the environment variable G4INCLXX_NO_DE_EXCITATION to disable de-excitation 00056 if(getenv("G4INCLXX_NO_DE_EXCITATION")) { 00057 G4String message = "de-excitation is completely disabled!"; 00058 theInterfaceStore->EmitWarning(message); 00059 theExcitationHandler = 0; 00060 } else { 00061 theExcitationHandler = new G4ExcitationHandler; 00062 } 00063 00064 theBackupModel = new G4BinaryLightIonReaction; 00065 }
G4INCLXXInterface::~G4INCLXXInterface | ( | ) |
G4HadFinalState * G4INCLXXInterface::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | theNucleus | |||
) | [virtual] |
Main method to apply the INCL physics model.
aTrack | the projectile particle | |
theNucleus | target nucleus |
Implements G4HadronicInteraction.
Definition at line 118 of file G4INCLXXInterface.cc.
References G4INCL::EventInfo::A, G4HadFinalState::AddSecondary(), G4HadronicInteraction::ApplyYourself(), G4INCL::EventInfo::ARem, G4ExcitationHandler::BreakItUp(), G4HadFinalState::Clear(), G4INCL::INCL::configToString(), G4INCL::EventInfo::EKin, G4INCL::EventInfo::EKinRem, G4INCLXXInterfaceStore::EmitWarning(), G4INCL::EventInfo::EStarRem, G4cout, G4endl, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4HadProjectile::GetDefinition(), G4INCLXXInterfaceStore::GetDumpInput(), G4INCLXXInterfaceStore::GetINCLModel(), G4INCLXXInterfaceStore::GetInstance(), G4ParticleTable::GetIon(), G4IonTable::GetIonMass(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4INCLXXInterfaceStore::GetMaxProjMassINCL(), G4HadronicInteraction::GetModelName(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4Nucleus::GetZ_asInt(), isAlive, G4INCL::EventInfo::jxRem, G4INCL::EventInfo::jyRem, G4INCL::EventInfo::jzRem, G4INCL::EventInfo::nParticles, G4INCL::EventInfo::nRemnants, G4INCL::INCL::processEvent(), G4INCL::EventInfo::px, G4INCL::EventInfo::pxRem, G4INCL::EventInfo::py, G4INCL::EventInfo::pyRem, G4INCL::EventInfo::pz, G4INCL::EventInfo::pzRem, G4Fragment::SetAngularMomentum(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), stopAndKill, G4INCL::EventInfo::transparent, G4INCL::EventInfo::Z, and G4INCL::EventInfo::ZRem.
00119 { 00120 // For systems heavier than theMaxProjMassINCL, use another model (typically 00121 // BIC) 00122 const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL(); 00123 if(aTrack.GetDefinition()->GetAtomicMass() > theMaxProjMassINCL 00124 && theNucleus.GetA_asInt() > theMaxProjMassINCL) { 00125 if(!complainedAboutBackupModel) { 00126 complainedAboutBackupModel = true; 00127 std::stringstream ss; 00128 ss << "INCL++ refuses to handle reactions between nuclei with A>" 00129 << theMaxProjMassINCL 00130 << ". A backup model (" 00131 << theBackupModel->GetModelName() 00132 << ") will be used instead."; 00133 G4cout << "[INCL++] Warning: " << ss.str() << G4endl; 00134 } 00135 return theBackupModel->ApplyYourself(aTrack, theNucleus); 00136 } 00137 00138 const G4int maxTries = 200; 00139 00140 // Check if inverse kinematics should be used 00141 const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus); 00142 00143 // If we are running in inverse kinematics, redefine aTrack and theNucleus 00144 G4LorentzRotation *toInverseKinematics = NULL; 00145 G4LorentzRotation *toDirectKinematics = NULL; 00146 G4HadProjectile const *aProjectileTrack = &aTrack; 00147 G4Nucleus *theTargetNucleus = &theNucleus; 00148 if(inverseKinematics) { 00149 G4ParticleTable * const theParticleTable = G4ParticleTable::GetParticleTable(); 00150 const G4int oldTargetA = theNucleus.GetA_asInt(); 00151 const G4int oldTargetZ = theNucleus.GetZ_asInt(); 00152 G4ParticleDefinition *oldTargetDef = theParticleTable->GetIon(oldTargetZ, oldTargetA, 0.0); 00153 const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition(); 00154 00155 if(oldProjectileDef != 0 && oldTargetDef != 0) { 00156 const G4int newTargetA = oldProjectileDef->GetAtomicMass(); 00157 const G4int newTargetZ = oldProjectileDef->GetAtomicNumber(); 00158 00159 if(newTargetA > 0 && newTargetZ > 0) { 00160 // This should give us the same energy per nucleon 00161 theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ); 00162 const G4double theProjectileMass = theParticleTable->GetIonTable()->GetIonMass(oldTargetZ, oldTargetA); 00163 toInverseKinematics = new G4LorentzRotation(aTrack.Get4Momentum().boostVector()); 00164 G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theProjectileMass); 00165 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum); 00166 aProjectileTrack = new G4HadProjectile(swappedProjectileParticle); 00167 } else { 00168 G4String message = "badly defined target after swapping. Falling back to normal (non-swapped) mode."; 00169 theInterfaceStore->EmitWarning(message); 00170 toInverseKinematics = new G4LorentzRotation; 00171 } 00172 } else { 00173 G4String message = "oldProjectileDef or oldTargetDef was null"; 00174 theInterfaceStore->EmitWarning(message); 00175 toInverseKinematics = new G4LorentzRotation; 00176 } 00177 } 00178 00179 // INCL assumes the projectile particle is going in the direction of 00180 // the Z-axis. Here we construct proper rotation to convert the 00181 // momentum vectors of the outcoming particles to the original 00182 // coordinate system. 00183 G4LorentzVector projectileMomentum = aProjectileTrack->Get4Momentum(); 00184 00185 // INCL++ assumes that the projectile is going in the direction of 00186 // the z-axis. In principle, if the coordinate system used by G4 00187 // hadronic framework is defined differently we need a rotation to 00188 // transform the INCL++ reaction products to the appropriate 00189 // frame. Please note that it isn't necessary to apply this 00190 // transform to the projectile because when creating the INCL++ 00191 // projectile particle; \see{toINCLKineticEnergy} needs to use only the 00192 // projectile energy (direction is simply assumed to be along z-axis). 00193 G4RotationMatrix toZ; 00194 toZ.rotateZ(-projectileMomentum.phi()); 00195 toZ.rotateY(-projectileMomentum.theta()); 00196 G4RotationMatrix toLabFrame3 = toZ.inverse(); 00197 G4LorentzRotation toLabFrame(toLabFrame3); 00198 // However, it turns out that the projectile given to us by G4 00199 // hadronic framework is already going in the direction of the 00200 // z-axis so this rotation is actually unnecessary. Both toZ and 00201 // toLabFrame turn out to be unit matrices as can be seen by 00202 // uncommenting the folowing two lines: 00203 // G4cout <<"toZ = " << toZ << G4endl; 00204 // G4cout <<"toLabFrame = " << toLabFrame << G4endl; 00205 00206 theResult.Clear(); // Make sure the output data structure is clean. 00207 theResult.SetStatusChange(stopAndKill); 00208 00209 std::list<G4Fragment> remnants; 00210 00211 G4int nTries = 0; 00212 // INCL can generate transparent events. However, this is meaningful 00213 // only in the standalone code. In Geant4 we must "force" INCL to 00214 // produce a valid cascade. 00215 G4bool eventIsOK = false; 00216 do { 00217 const G4INCL::ParticleSpecies theSpecies = toINCLParticleSpecies(*aProjectileTrack); 00218 const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack); 00219 00220 // The INCL model will be created at the first use 00221 theINCLModel = G4INCLXXInterfaceStore::GetInstance()->GetINCLModel(); 00222 00223 if(theInterfaceStore->GetDumpInput()) { 00224 G4cout << theINCLModel->configToString() << G4endl; 00225 } 00226 const G4INCL::EventInfo eventInfo = theINCLModel->processEvent(theSpecies, kineticEnergy, theTargetNucleus->GetA_asInt(), theTargetNucleus->GetZ_asInt()); 00227 // eventIsOK = !eventInfo.transparent && nTries < maxTries; 00228 eventIsOK = !eventInfo.transparent; 00229 if(eventIsOK) { 00230 00231 // If the collision was run in inverse kinematics, prepare to boost back 00232 // all the reaction products 00233 if(inverseKinematics) { 00234 toDirectKinematics = new G4LorentzRotation(toInverseKinematics->inverse()); 00235 } 00236 00237 for(G4int i = 0; i < eventInfo.nParticles; i++) { 00238 G4int A = eventInfo.A[i]; 00239 G4int Z = eventInfo.Z[i]; 00240 // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl; 00241 G4double kinE = eventInfo.EKin[i]; 00242 G4double px = eventInfo.px[i]; 00243 G4double py = eventInfo.py[i]; 00244 G4double pz = eventInfo.pz[i]; 00245 G4DynamicParticle *p = toG4Particle(A, Z , kinE, px, py, pz); 00246 if(p != 0) { 00247 G4LorentzVector momentum = p->Get4Momentum(); 00248 00249 // Apply the toLabFrame rotation 00250 momentum *= toLabFrame; 00251 // Apply the inverse-kinematics boost 00252 if(inverseKinematics) { 00253 momentum *= *toDirectKinematics; 00254 momentum.setVect(-momentum.vect()); 00255 } 00256 00257 // Set the four-momentum of the reaction products 00258 p->Set4Momentum(momentum); 00259 theResult.AddSecondary(p); 00260 00261 } else { 00262 G4String message = "the model produced a particle that couldn't be converted to Geant4 particle."; 00263 theInterfaceStore->EmitWarning(message); 00264 } 00265 } 00266 00267 for(G4int i = 0; i < eventInfo.nRemnants; i++) { 00268 const G4int A = eventInfo.ARem[i]; 00269 const G4int Z = eventInfo.ZRem[i]; 00270 // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl; 00271 const G4double kinE = eventInfo.EKinRem[i]; 00272 const G4double px = eventInfo.pxRem[i]; 00273 const G4double py = eventInfo.pyRem[i]; 00274 const G4double pz = eventInfo.pzRem[i]; 00275 G4ThreeVector spin( 00276 eventInfo.jxRem[i]*hbar_Planck, 00277 eventInfo.jyRem[i]*hbar_Planck, 00278 eventInfo.jzRem[i]*hbar_Planck 00279 ); 00280 const G4double excitationE = eventInfo.EStarRem[i]; 00281 const G4double nuclearMass = G4NucleiProperties::GetNuclearMass(A, Z) + excitationE; 00282 const G4double scaling = remnant4MomentumScaling(nuclearMass, 00283 kinE, 00284 px, py, pz); 00285 G4LorentzVector fourMomentum(scaling * px, scaling * py, scaling * pz, 00286 nuclearMass + kinE); 00287 if(std::abs(scaling - 1.0) > 0.01) { 00288 std::stringstream ss; 00289 ss << "momentum scaling = " << scaling 00290 << "\n Lorentz vector = " << fourMomentum 00291 << ")\n A = " << A << ", Z = " << Z 00292 << "\n E* = " << excitationE << ", nuclearMass = " << nuclearMass 00293 << "\n remnant i=" << i << ", nRemnants=" << eventInfo.nRemnants 00294 << "\n Reaction was: " << aTrack.GetKineticEnergy()/MeV 00295 << "-MeV " << aTrack.GetDefinition()->GetParticleName() << " + " 00296 << G4ParticleTable::GetParticleTable()->GetIon(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0.0)->GetParticleName() 00297 << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics."; 00298 theInterfaceStore->EmitWarning(ss.str()); 00299 } 00300 00301 // Apply the toLabFrame rotation 00302 fourMomentum *= toLabFrame; 00303 spin *= toLabFrame3; 00304 // Apply the inverse-kinematics boost 00305 if(inverseKinematics) { 00306 fourMomentum *= *toDirectKinematics; 00307 fourMomentum.setVect(-fourMomentum.vect()); 00308 } 00309 00310 G4Fragment remnant(A, Z, fourMomentum); 00311 remnant.SetAngularMomentum(spin); 00312 remnants.push_back(remnant); 00313 } 00314 } 00315 nTries++; 00316 } while(!eventIsOK && nTries < maxTries); 00317 00318 // Clean up the objects that we created for the inverse kinematics 00319 if(inverseKinematics) { 00320 delete aProjectileTrack; 00321 delete theTargetNucleus; 00322 delete toInverseKinematics; 00323 delete toDirectKinematics; 00324 } 00325 00326 if(!eventIsOK) { 00327 std::stringstream ss; 00328 ss << "maximum number of tries exceeded for the proposed " 00329 << aTrack.GetKineticEnergy()/MeV << "-MeV " << aTrack.GetDefinition()->GetParticleName() 00330 << " + " << G4ParticleTable::GetParticleTable()->GetIon(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0.0)->GetParticleName() 00331 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics."; 00332 theInterfaceStore->EmitWarning(ss.str()); 00333 theResult.SetStatusChange(isAlive); 00334 theResult.SetEnergyChange(aTrack.GetKineticEnergy()); 00335 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 00336 return &theResult; 00337 } 00338 00339 // De-excitation: 00340 00341 if(theExcitationHandler != 0) { 00342 for(std::list<G4Fragment>::const_iterator i = remnants.begin(); 00343 i != remnants.end(); i++) { 00344 G4ReactionProductVector *deExcitationResult = theExcitationHandler->BreakItUp((*i)); 00345 00346 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin(); 00347 fragment != deExcitationResult->end(); ++fragment) { 00348 G4ParticleDefinition *def = (*fragment)->GetDefinition(); 00349 if(def != 0) { 00350 G4DynamicParticle *theFragment = new G4DynamicParticle(def, (*fragment)->GetMomentum()); 00351 theResult.AddSecondary(theFragment); 00352 } 00353 } 00354 00355 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin(); 00356 fragment != deExcitationResult->end(); ++fragment) { 00357 delete (*fragment); 00358 } 00359 deExcitationResult->clear(); 00360 delete deExcitationResult; 00361 } 00362 } 00363 00364 remnants.clear(); 00365 00366 return &theResult; 00367 }
void G4INCLXXInterface::DeleteModel | ( | ) | [inline] |
G4int G4INCLXXInterface::operator!= | ( | G4INCLXXInterface & | right | ) | [inline] |
G4int G4INCLXXInterface::operator== | ( | G4INCLXXInterface & | right | ) | [inline] |
G4ReactionProductVector * G4INCLXXInterface::Propagate | ( | G4KineticTrackVector * | theSecondaries, | |
G4V3DNucleus * | theNucleus | |||
) | [virtual] |