G4INCLXXInterface Class Reference

INCL++ intra-nuclear cascade with G4ExcitationHandler for de-excitation. More...

#include <G4INCLXXInterface.hh>

Inheritance diagram for G4INCLXXInterface:

G4VIntraNuclearTransportModel G4HadronicInteraction

Public Member Functions

 G4INCLXXInterface (const G4String &name="INCL++ cascade with G4ExcitationHandler")
 ~G4INCLXXInterface ()
G4int operator== (G4INCLXXInterface &right)
G4int operator!= (G4INCLXXInterface &right)
G4ReactionProductVectorPropagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
void DeleteModel ()

Detailed Description

INCL++ intra-nuclear cascade with G4ExcitationHandler for de-excitation.

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);
The same setup procedure is needed for neutron, pion and generic-ion inelastic processes as well.

Definition at line 91 of file G4INCLXXInterface.hh.


Constructor & Destructor Documentation

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 (  ) 

Definition at line 67 of file G4INCLXXInterface.cc.

00068 {
00069   delete theBackupModel;
00070   delete theExcitationHandler;
00071 }


Member Function Documentation

G4HadFinalState * G4INCLXXInterface::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus theNucleus 
) [virtual]

Main method to apply the INCL physics model.

Parameters:
aTrack the projectile particle
theNucleus target nucleus
Returns:
the output of the INCL physics model

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]

Definition at line 114 of file G4INCLXXInterface.hh.

00114                      {
00115     delete theINCLModel;
00116     theINCLModel = NULL;
00117   }

G4int G4INCLXXInterface::operator!= ( G4INCLXXInterface right  )  [inline]

Definition at line 100 of file G4INCLXXInterface.hh.

00100                                              {
00101     return (this != &right);
00102   }

G4int G4INCLXXInterface::operator== ( G4INCLXXInterface right  )  [inline]

Definition at line 96 of file G4INCLXXInterface.hh.

00096                                              {
00097     return (this == &right);
00098   }

G4ReactionProductVector * G4INCLXXInterface::Propagate ( G4KineticTrackVector theSecondaries,
G4V3DNucleus theNucleus 
) [virtual]

Implements G4VIntraNuclearTransportModel.

Definition at line 369 of file G4INCLXXInterface.cc.

00369                                                                                             {
00370   return 0;
00371 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:15 2013 for Geant4 by  doxygen 1.4.7