#include <G4BinaryLightIonReaction.hh>
Inheritance diagram for G4BinaryLightIonReaction:
Public Member Functions | |
G4BinaryLightIonReaction (G4VPreCompoundModel *ptr=0) | |
virtual | ~G4BinaryLightIonReaction () |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus) |
void | SetPrecompound (G4VPreCompoundModel *ptr) |
void | SetDeExcitation (G4ExcitationHandler *ptr) |
virtual void | ModelDescription (std::ostream &) const |
Definition at line 34 of file G4BinaryLightIonReaction.hh.
G4BinaryLightIonReaction::G4BinaryLightIonReaction | ( | G4VPreCompoundModel * | ptr = 0 |
) |
Definition at line 49 of file G4BinaryLightIonReaction.cc.
References G4HadronicInteractionRegistry::FindModel(), G4VPreCompoundModel::GetExcitationHandler(), and G4HadronicInteractionRegistry::Instance().
00050 : G4HadronicInteraction("Binary Light Ion Cascade"), 00051 theProjectileFragmentation(ptr), 00052 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0), 00053 projectile3dNucleus(0),target3dNucleus(0) 00054 { 00055 if(!ptr) { 00056 G4HadronicInteraction* p = 00057 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO"); 00058 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p); 00059 if(!pre) { pre = new G4PreCompoundModel(); } 00060 theProjectileFragmentation = pre; 00061 } 00062 theModel = new G4BinaryCascade(theProjectileFragmentation); 00063 theHandler = theProjectileFragmentation->GetExcitationHandler(); 00064 00065 debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0; 00066 }
G4BinaryLightIonReaction::~G4BinaryLightIonReaction | ( | ) | [virtual] |
G4HadFinalState * G4BinaryLightIonReaction::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | theNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 87 of file G4BinaryLightIonReaction.cc.
References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4cerr, G4cout, G4endl, G4lrint(), G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4IonTable::GetIonMass(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4HadFinalState::GetNumberOfSecondaries(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4HadProjectile::GetTotalEnergy(), G4Nucleus::GetZ_asInt(), isAlive, G4DynamicParticle::Set4Momentum(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), and stopAndKill.
00088 { 00089 static G4int eventcounter=0; 00090 eventcounter++; 00091 if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number starts ######### "<<eventcounter<<G4endl; 00092 G4ping debug("debug_G4BinaryLightIonReaction"); 00093 pA=aTrack.GetDefinition()->GetBaryonNumber(); 00094 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus); 00095 tA=targetNucleus.GetA_asInt(); 00096 tZ=targetNucleus.GetZ_asInt(); 00097 00098 G4LorentzVector mom(aTrack.Get4Momentum()); 00099 //G4cout << "proj mom : " << mom << G4endl; 00100 G4LorentzRotation toBreit(mom.boostVector()); 00101 00102 G4bool swapped=SetLighterAsProjectile(mom, toBreit); 00103 //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl; 00104 G4ReactionProductVector * result = 0; 00105 G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector; 00106 // G4double m_nucl(0); // to check energy balance 00107 00108 00109 // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA); 00110 // G4cout << "Entering the decision point " 00111 // << (mom.t()-mom.mag())/pA << " " 00112 // << pA<<" "<< pZ<<" " 00113 // << tA<<" "<< tZ<<G4endl 00114 // << " "<<mom.t()-mom.mag()<<" " 00115 // << mom.t()- m1<<G4endl; 00116 if( (mom.t()-mom.mag())/pA < 50*MeV ) 00117 { 00118 // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl; 00119 // m_nucl = mom.mag(); 00120 cascaders=FuseNucleiAndPrompound(mom); 00121 } 00122 else 00123 { 00124 result=Interact(mom,toBreit); 00125 00126 if(! result ) 00127 { 00128 { 00129 // abort!! 00130 00131 G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl; 00132 G4cerr << " Primary " << aTrack.GetDefinition() 00133 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() 00134 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") " 00135 << ", kinetic energy " << aTrack.GetKineticEnergy() 00136 << G4endl; 00137 G4cerr << " Target nucleus (A,Z)=(" 00138 << (swapped?pA:tA) << "," 00139 << (swapped?pZ:tZ) << ")" << G4endl; 00140 G4cerr << " if frequent, please submit above information as bug report" 00141 << G4endl << G4endl; 00142 00143 theResult.Clear(); 00144 theResult.SetStatusChange(isAlive); 00145 theResult.SetEnergyChange(aTrack.GetKineticEnergy()); 00146 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 00147 return &theResult; 00148 00149 } 00150 } 00151 00152 // Calculate excitation energy, 00153 G4double theStatisticalExEnergy = GetProjectileExcitation(); 00154 00155 00156 pInitialState = mom; 00157 //G4cout << "pInitialState from aTrack : " << pInitialState; 00158 pInitialState.setT(pInitialState.getT() + 00159 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA)); 00160 //G4cout << "target nucleus added : " << pInitialState << G4endl; 00161 00162 delete target3dNucleus;target3dNucleus=0; 00163 delete projectile3dNucleus;projectile3dNucleus=0; 00164 00165 G4ReactionProductVector * spectators= new G4ReactionProductVector; 00166 00167 cascaders = new G4ReactionProductVector; 00168 00169 G4LorentzVector pspectators=SortResult(result,spectators,cascaders); 00170 00171 // pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom); 00172 std::vector<G4ReactionProduct *>::iterator iter; 00173 00174 //G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl; 00175 // if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0) 00176 // { 00177 // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" << 00178 // spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl; 00179 // } 00180 delete result; 00181 result=0; 00182 G4LorentzVector momentum(pInitialState-pFinalState); 00183 G4int loopcount(0); 00184 //G4cout << "momentum, pspectators : " << momentum << " / " << pspectators << G4endl; 00185 while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) 00186 { 00187 G4LorentzVector pCorrect(pInitialState-pspectators); 00188 //G4cout << "BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl; 00189 // Correct outgoing casacde particles.... to have momentum of (initial state - spectators) 00190 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect); 00191 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults) 00192 { 00193 G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl; 00194 } 00195 pFinalState=G4LorentzVector(0,0,0,0); 00196 unsigned int i; 00197 for(i=0; i<cascaders->size(); i++) 00198 { 00199 pFinalState += G4LorentzVector( (*cascaders)[i]->GetMomentum(), (*cascaders)[i]->GetTotalEnergy() ); 00200 } 00201 momentum=pInitialState-pFinalState; 00202 if (++loopcount > 10 ) 00203 { 00204 if ( momentum.vect().mag() - momentum.e()> 10*keV ) 00205 { 00206 G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl; 00207 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()"); 00208 } else { 00209 break; 00210 } 00211 00212 } 00213 } 00214 00215 if (spectatorA > 0 ) 00216 { 00217 // check spectator momentum 00218 if ( momentum.vect().mag() - momentum.e()> 10*keV ) 00219 { 00220 00221 G4ReactionProductVector::iterator ispectator; 00222 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) 00223 { 00224 delete *ispectator; 00225 } 00226 delete spectators; 00227 00228 G4cout << "G4BinaryLightIonReaction.cc: mom check: " << momentum 00229 << " 3.mag "<< momentum.vect().mag() << G4endl 00230 << " .. pInitialState/pFinalState/spectators " << pInitialState <<" " 00231 << pFinalState << " " << pspectators << G4endl 00232 << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl; 00233 G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl; 00234 G4cout << " Primary " << aTrack.GetDefinition() 00235 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() 00236 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") " 00237 << ", kinetic energy " << aTrack.GetKineticEnergy() 00238 << G4endl; 00239 G4cout << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt() 00240 << "," << targetNucleus.GetZ_asInt() << ")" << G4endl; 00241 G4cout << " if frequent, please submit above information as bug report" 00242 << G4endl << G4endl; 00243 00244 theResult.Clear(); 00245 theResult.SetStatusChange(isAlive); 00246 theResult.SetEnergyChange(aTrack.GetKineticEnergy()); 00247 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 00248 return &theResult; 00249 } 00250 00251 00252 DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum); 00253 } 00254 } 00255 // Rotate to lab 00256 G4LorentzRotation toZ; 00257 toZ.rotateZ(-1*mom.phi()); 00258 toZ.rotateY(-1*mom.theta()); 00259 G4LorentzRotation toLab(toZ.inverse()); 00260 00261 // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped. 00262 // theResult.Clear(); 00263 theResult.Clear(); 00264 theResult.SetStatusChange(stopAndKill); 00265 G4double Etot(0); 00266 size_t i=0; 00267 for(i=0; i<cascaders->size(); i++) 00268 { 00269 if((*cascaders)[i]->GetNewlyAdded()) 00270 { 00271 G4DynamicParticle * aNew = 00272 new G4DynamicParticle((*cascaders)[i]->GetDefinition(), 00273 (*cascaders)[i]->GetTotalEnergy(), 00274 (*cascaders)[i]->GetMomentum() ); 00275 G4LorentzVector tmp = aNew->Get4Momentum(); 00276 if(swapped) 00277 { 00278 tmp*=toBreit.inverse(); 00279 tmp.setVect(-tmp.vect()); 00280 } 00281 tmp *= toLab; 00282 aNew->Set4Momentum(tmp); 00283 //G4cout << "result[" << i << "], 4vect: " << tmp << G4endl; 00284 theResult.AddSecondary(aNew); 00285 Etot += tmp.e(); 00286 // G4cout << "LIBIC: Secondary " << aNew->GetDefinition()->GetParticleName() 00287 // <<" "<< aNew->GetMomentum() 00288 // <<" "<< aNew->GetTotalEnergy() 00289 // << G4endl; 00290 } 00291 delete (*cascaders)[i]; 00292 } 00293 delete cascaders; 00294 00295 #ifdef debug_BLIR_result 00296 G4cout << "Result analysis, secondaries" << theResult.GetNumberOfSecondaries() << G4endl; 00297 G4cout << " Energy conservation initial/primary/nucleus/final/delta(init-final) " 00298 << aTrack.GetTotalEnergy() + m_nucl << aTrack.GetTotalEnergy() << m_nucl <<Etot 00299 << aTrack.GetTotalEnergy() + m_nucl - Etot; 00300 #endif 00301 00302 if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### "<<eventcounter<<G4endl; 00303 00304 return &theResult; 00305 }
void G4BinaryLightIonReaction::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 71 of file G4BinaryLightIonReaction.cc.
00072 { 00073 outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n" 00074 << "using G4BinaryCasacde to model the interaction of a light\n" 00075 << "nucleus with a nucleus.\n" 00076 << "The lighter of the two nuclei is treated like a set of projectiles\n" 00077 << "which are transported simultanously through the heavier nucleus.\n"; 00078 }
void G4BinaryLightIonReaction::SetDeExcitation | ( | G4ExcitationHandler * | ptr | ) | [inline] |
Definition at line 74 of file G4BinaryLightIonReaction.hh.
References G4VPreCompoundModel::SetExcitationHandler().
00075 { 00076 theProjectileFragmentation->SetExcitationHandler(ptr); 00077 theHandler = ptr; 00078 }
void G4BinaryLightIonReaction::SetPrecompound | ( | G4VPreCompoundModel * | ptr | ) | [inline] |
Definition at line 69 of file G4BinaryLightIonReaction.hh.
References G4VPreCompoundModel::GetExcitationHandler().
00070 { 00071 if(ptr) { theProjectileFragmentation = ptr; } 00072 theHandler = theProjectileFragmentation->GetExcitationHandler(); 00073 }