99 #define _CheckChargeAndBaryonNumber_(val) CheckChargeAndBaryonNumber(val)
102 #define _CheckChargeAndBaryonNumber_(val)
106 #define _DebugEpConservation(val) DebugEpConservation(val)
109 #define _DebugEpConservation(val)
112 #ifdef G4MULTITHREADED
169 #ifdef G4MULTITHREADED
174 #ifdef G4MULTITHREADED
204 outFile <<
"G4BinaryCascade is an intra-nuclear cascade model in which\n"
205 <<
"an incident hadron collides with a nucleon, forming two\n"
206 <<
"final-state particles, one or both of which may be resonances.\n"
207 <<
"The resonances then decay hadronically and the decay products\n"
208 <<
"are then propagated through the nuclear potential along curved\n"
209 <<
"trajectories until they re-interact or leave the nucleus.\n"
210 <<
"This model is valid for incident pions up to 1.5 GeV and\n"
211 <<
"nucleons up to 10 GeV.\n"
212 <<
"The remaining excited nucleus is handed on to ";
220 outFile <<
"G4ExcitationHandler";
225 outFile <<
"void.\n";
231 outFile <<
"G4BinaryCascade propagtes secondaries produced by a high\n"
232 <<
"energy model through the wounded nucleus.\n"
233 <<
"Secondaries are followed after the formation time and if\n"
234 <<
"within the nucleus are propagated through the nuclear\n"
235 <<
"potential along curved trajectories until they interact\n"
236 <<
"with a nucleon, decay, or leave the nucleus.\n"
237 <<
"An interaction of a secondary with a nucleon produces two\n"
238 <<
"final-state particles, one or both of which may be resonances.\n"
239 <<
"Resonances decay hadronically and the decay products\n"
240 <<
"are in turn propagated through the nuclear potential along curved\n"
241 <<
"trajectories until they re-interact or leave the nucleus.\n"
242 <<
"This model is valid for pions up to 1.5 GeV and\n"
243 <<
"nucleons up to about 3.5 GeV.\n"
244 <<
"The remaining excited nucleus is handed on to ";
252 outFile <<
"G4ExcitationHandler";
257 outFile <<
"void.\n";
274 if(std::getenv(
"BCDEBUG") )
G4cerr <<
" ######### Binary Cascade Reaction starts ######### "<<
G4endl;
279 if(initial4Momentum.
e()-initial4Momentum.
m()<
theBCminP &&
293 if(!std::getenv(
"I_Am_G4BinaryCascade_Developer") )
300 G4cerr <<
"You are trying to use G4BinaryCascade with " <<definition->GetParticleName()<<
" as projectile."<<
G4endl;
301 G4cerr <<
"G4BinaryCascade should not be used for projectiles other than nucleons or pions."<<
G4endl;
302 G4cerr <<
"If you want to continue, please switch on the developer environment: "<<
G4endl;
303 G4cerr <<
"setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<
G4endl;
304 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade - used for unvalid particle type - Fatal");
316 G4int interactionCounter = 0,collisionLoopMaxCount;
333 collisionLoopMaxCount = 200;
339 kt =
new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
343 secondaries->push_back(kt);
351 }
while(! products && --collisionLoopMaxCount>0);
353 if(++interactionCounter>99)
break;
355 }
while(products && products->size() == 0);
357 if(products && products->size()>0)
363 G4ReactionProductVector::iterator iter;
365 for(iter = products->begin(); iter != products->end(); ++iter)
369 (*iter)->GetTotalEnergy(),
370 (*iter)->GetMomentum());
373 if(time < 0.0) { time = 0.0; }
374 aNew.SetTime(timePrimary + time);
375 aNew.SetCreatorModelType((*iter)->GetCreatorModel());
383 if(std::getenv(
"BCDEBUG") )
G4cerr <<
" ######### Binary Cascade Reaction void, return intial state ######### "<<
G4endl;
398 if(std::getenv(
"BCDEBUG") )
G4cerr <<
" ######### Binary Cascade Reaction ends ######### "<<
G4endl;
408 #ifdef debug_BIC_Propagate
409 G4cout <<
"G4BinaryCascade Propagate starting -------------------------------------------------------" <<
G4endl;
423 std::vector<G4KineticTrack *>::iterator iter;
434 #ifdef debug_BIC_GetExcitationEnergy
447 #ifdef debug_G4BinaryCascade
448 G4cout <<
"G4BinaryCascade::Propagate: warning - high energy model failed energy conservation, returning unchanged high energy final state" <<
G4endl;
466 #ifdef debug_BIC_return
476 G4bool haveProducts =
false;
477 G4int collisionCount=0;
478 G4int collisionLoopMaxCount=1000000;
494 #ifdef debug_BIC_Propagate_Collisions
495 G4cout <<
" NextCollision * , Time, curtime = " << nextCollision <<
" "
537 #ifdef debug_BIC_return
538 G4cout <<
"return @ Z=0 after collision loop "<<
G4endl;
547 G4cout <<
"returned products: " << products->size() <<
G4endl;
568 #ifdef debug_BIC_return
575 #ifdef debug_BIC_Propagate
586 #ifdef debug_G4BinaryCascade
587 G4cerr <<
"G4BinaryCascade: Warning, have active particles at end" <<
G4endl;
600 #ifdef debug_G4BinaryCascade
601 G4cerr <<
" Warning: remove left over collision(s) " <<
G4endl;
606 #ifdef debug_BIC_Propagate_Excitation
623 #ifdef debug_BIC_Propagate_finals
625 G4cout <<
" Excitation Energy prefinal, #collisions:, out, captured "
626 << ExcitationEnergy <<
" "
627 << collisionCount <<
" "
632 if (ExcitationEnergy < 0 )
634 G4int maxtry=5, ntry=0;
638 }
while ( ++ntry < maxtry && ExcitationEnergy < 0 );
642 #ifdef debug_BIC_Propagate_finals
644 G4cout <<
" Excitation Energy final, #collisions:, out, captured "
645 << ExcitationEnergy <<
" "
646 << collisionCount <<
" "
652 if ( ExcitationEnergy < 0. )
654 #ifdef debug_G4BinaryCascade
655 G4cerr <<
"G4BinaryCascade-Warning: negative excitation energy ";
664 #ifdef debug_BIC_return
665 G4cout <<
" negative Excitation E return empty products " << products <<
G4endl;
688 #ifdef debug_BIC_return
702 #if defined(debug_G4BinaryCascade) || defined(debug_BIC_GetExcitationEnergy)
707 G4cerr <<
"G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} "
727 #ifdef debug_G4BinaryCascade
728 G4cout <<
"G4BinaryCascade::GetExcitationEnergy(): Warning - invalid nucleus (A,Z)=("
734 #ifdef debug_BIC_GetExcitationEnergy
755 #ifdef debug_BIC_GetExcitationEnergy
757 if ( excitationE < 0 )
759 G4cout <<
"negative ExE final Ion mass " <<nucleusMass<<
G4endl;
761 if(finalZ>.5)
G4cout <<
" Final nuclmom/mass " << Nucl_mom <<
" " << Nucl_mom.
mag()
762 <<
" (A,Z)=("<< finalA <<
","<<finalZ <<
")"
763 <<
" mass " << nucleusMass <<
" "
764 <<
" excitE " << excitationE <<
G4endl;
773 G4cout <<
"GetExcitationEnergy: Initial nucleus A Z " << A <<
" " << Z <<
" " << initialExc <<
G4endl;
834 #ifdef debug_BIC_BuildTargetList
835 else {
G4cout <<
"nucleon is hit" << nucleon <<
G4endl;}
847 G4cerr <<
"G4BinaryCascade::BuildTargetList(): Fatal Error - invalid nucleus (A,Z)=("
853 #ifdef debug_BIC_BuildTargetList
854 G4cout <<
"G4BinaryCascade::BuildTargetList(): nucleus (A,Z)=("
867 std::vector<G4KineticTrack *>::iterator iter;
873 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
875 if((*iter)->GetFormationTime() < StartingTime)
876 StartingTime = (*iter)->GetFormationTime();
881 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
885 G4double FormTime = (*iter)->GetFormationTime() - StartingTime;
886 (*iter)->SetFormationTime(FormTime);
890 lateParticles4Momentum += (*iter)->GetTrackingMomentum();
891 lateA += (*iter)->GetDefinition()->GetBaryonNumber();
899 projectileA += (*iter)->GetDefinition()->GetBaryonNumber();
901 #ifdef debug_BIC_Propagate
902 G4cout <<
" Adding initial secondary " << *iter
903 <<
" time" << (*iter)->GetFormationTime()
904 <<
", state " << (*iter)->GetState() <<
G4endl;
918 #ifdef debug_BIC_GetExcitationEnergy
919 G4cout <<
"BIC: Proj.e, nucl initial, nucl final, lateParticles"
922 << lateParticles4Momentum <<
G4endl;
925 success = excitation > 0;
926 #ifdef debug_G4BinaryCascade
938 secondaries->clear();
981 std::vector<G4KineticTrack *>::iterator i;
989 precompoundProducts->push_back(aNew);
999 #ifdef debug_BIC_DeexcitationProducts
1003 if ( precompoundProducts )
1005 std::vector<G4ReactionProduct *>::iterator j;
1006 for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
1009 Preco_momentum += pProduct;
1012 G4cout <<
"finalNuclMom / sum preco products" << fragment_momentum <<
" / " << Preco_momentum
1013 <<
" delta E "<< fragment_momentum.
e() - Preco_momentum.
e() <<
G4endl;
1017 return precompoundProducts;
1028 std::vector<G4KineticTrack *>::iterator aNuc;
1030 std::vector<G4double> masses;
1038 masses.push_back(mass);
1049 masses.push_back(mass);
1060 if ( eCMS < sumMass )
1062 eCMS=sumMass + 2*
MeV*masses.size();
1067 std::vector<G4LorentzVector*> * momenta=decay.
Decay(eCMS,masses);
1068 std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
1081 result->push_back(aNew);
1094 (*aNuc)->GetDefinition());
1098 result->push_back(aNew);
1114 #ifdef debug_BIC_Propagate_finals
1117 for(i = 0; i< fs.size(); i++)
1125 products->push_back(aNew);
1127 #ifdef debug_BIC_Propagate_finals
1137 #ifdef debug_BIC_Propagate_finals
1138 G4cout <<
" Final state momentum " << mom_fs <<
G4endl;
1149 if ( precompoundProducts )
1151 std::vector<G4ReactionProduct *>::iterator j;
1152 for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
1157 #ifdef debug_BIC_Propagate_finals
1158 G4cout <<
"BIC: pProduct be4 boost " <<pProduct <<
G4endl;
1161 #ifdef debug_BIC_Propagate_finals
1162 G4cout <<
"BIC: pProduct aft boost " <<pProduct <<
G4endl;
1164 pSumPreco += pProduct;
1165 (*j)->SetTotalEnergy(pProduct.e());
1166 (*j)->SetMomentum(pProduct.vect());
1167 (*j)->SetNewlyAdded(
true);
1168 products->push_back(*j);
1172 precompoundProducts->clear();
1173 delete precompoundProducts;
1181 for(std::vector<G4KineticTrack *>::iterator i = secondaries->begin();
1182 i != secondaries->end(); ++i)
1184 for(std::vector<G4BCAction *>::iterator j =
theImR.begin();
1188 const std::vector<G4CollisionInitialState *> & aCandList
1190 for(
size_t count=0; count<aCandList.size(); count++)
1205 const std::vector<G4CollisionInitialState *> & aCandList
1207 for(
size_t count=0; count<aCandList.size(); count++)
1224 }
else if ( tout > 0 )
1237 #ifdef debug_BIC_FindCollision
1238 G4cout <<
"FindLateP Particle, 4-mom, times newState "
1241 <<
" times " << tin <<
" " << tout <<
" "
1245 const std::vector<G4CollisionInitialState *> & aCandList
1247 for(
size_t count=0; count<aCandList.size(); count++)
1249 #ifdef debug_BIC_FindCollision
1250 G4cout <<
" Adding a late Col : " << aCandList[count] <<
G4endl;
1263 #ifdef debug_BIC_ApplyCollision
1264 G4cerr <<
"G4BinaryCascade::ApplyCollision start"<<
G4endl;
1270 G4bool haveTarget=target_collection.size()>0;
1273 #ifdef debug_G4BinaryCascade
1274 G4cout <<
"G4BinaryCasacde::ApplyCollision(): StateError " << primary <<
G4endl;
1276 PrintKTVector(&target_collection,std::string(
"... targets"));
1296 G4int initialBaryon(0);
1297 G4int initialCharge(0);
1311 #ifdef debug_BIC_ApplyCollision
1318 G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
1319 G4bool decayCollision= (!haveTarget) && products && products->size() > 1;
1323 #ifdef debug_G4BinaryCascade
1324 G4int lateBaryon(0), lateCharge(0);
1327 if ( lateParticleCollision )
1331 #ifdef debug_G4BinaryCascade
1332 lateBaryon = initialBaryon;
1333 lateCharge = initialCharge;
1335 initialBaryon=initialCharge=0;
1342 if (!lateParticleCollision)
1346 #ifdef debug_BIC_ApplyCollision
1347 if (products)
G4cout <<
" ======Failed Pauli =====" <<
G4endl;
1348 G4cerr <<
"G4BinaryCascade::ApplyCollision blocked"<<
G4endl;
1362 #ifdef debug_BIC_ApplyCollision
1374 G4int finalBaryon(0);
1375 G4int finalCharge(0);
1377 for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1379 if ( ! lateParticleCollision )
1381 (*i)->SetState(primary->
GetState());
1383 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1384 finalCharge+=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
1388 tin < 0 && tout > 0 )
1391 G4cout <<
"tin tout: " << tin <<
" " << tout <<
G4endl;
1403 else if ( tout > 0 )
1406 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1407 finalCharge+=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
1412 toFinalState.push_back((*i));
1418 toFinalState.push_back((*i));
1423 if(!toFinalState.empty())
1426 toFinalState.begin(),toFinalState.end());
1427 std::vector<G4KineticTrack *>::iterator iter1, iter2;
1428 for(iter1 = toFinalState.begin(); iter1 != toFinalState.end();
1431 iter2 = std::find(products->begin(), products->end(),
1433 if ( iter2 != products->end() ) products->erase(iter2);
1439 currentA += finalBaryon-initialBaryon;
1440 currentZ += finalCharge-initialCharge;
1444 oldSecondaries.push_back(primary);
1447 #ifdef debug_G4BinaryCascade
1448 if ( (finalBaryon-initialBaryon-lateBaryon) != 0 || (finalCharge-initialCharge-lateCharge) != 0 )
1450 G4cout <<
"G4BinaryCascade: Error in Balancing: " <<
G4endl;
1451 G4cout <<
"initial/final baryon number, initial/final Charge "
1452 << initialBaryon <<
" "<< finalBaryon <<
" "
1453 << initialCharge <<
" "<< finalCharge <<
" "
1455 <<
", with number of products: "<< products->size() <<
G4endl;
1456 G4cout << G4endl<<
"Initial condition are these:"<<
G4endl;
1469 for(
size_t ii=0; ii< oldTarget.size(); ii++)
1471 oldTarget[ii]->Hit();
1492 std::vector<G4KineticTrack *>::iterator iter;
1502 absorbList.push_back(kt);
1507 if(absorbList.empty())
1511 for(iter = absorbList.begin(); iter != absorbList.end(); ++iter)
1515 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1518 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1521 G4int maxLoopCount = 1000;
1527 "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1529 if ( --maxLoopCount < 0 )
throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1534 toRemove.push_back(kt);
1535 toDelete.push_back(kt);
1553 std::vector<G4KineticTrack *>::iterator i;
1558 G4int particlesAboveCut=0;
1559 G4int particlesBelowCut=0;
1578 ++particlesBelowCut;
1586 if (verbose)
G4cout <<
"Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP "
1587 << particlesAboveCut <<
" " << particlesBelowCut <<
" " << capturedEnergy
1591 if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*
theCutOnP)
1602 captured.push_back(kt);
1623 fermiMom.
Init(A, Z);
1627 G4KineticTrackVector::iterator i;
1634 for(i = products->begin(); i != products->end(); ++i)
1636 definition = (*i)->GetDefinition();
1662 if(mom.
e() < eFermi )
1671 #ifdef debug_BIC_CheckPauli
1674 for(i = products->begin(); i != products->end(); ++i)
1676 definition = (*i)->GetDefinition();
1685 if ( mom.
e()-mom.
mag()+field > 160*
MeV )
1687 G4cout <<
"momentum problem pFermi=" << pFermi
1688 <<
" mom, mom.m " << mom <<
" " << mom.
mag()
1689 <<
" field " << field <<
G4endl;
1713 std::vector<G4KineticTrack *>::iterator i;
1723 #ifdef debug_BIC_StepParticlesOut
1724 G4cout <<
" minTimeStep, tStep Particle " <<minTimeStep <<
" " <<tStep
1730 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1733 if(intersect && tStep<minTimeStep && tStep> 0 )
1735 minTimeStep = tStep;
1739 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1752 if ( timeToCollision > minTimeStep )
1782 #ifdef debug_G4BinaryCascade
1783 G4cerr <<
"G4BinaryCascade.cc: Warning - aborting looping particle(s)" <<
G4endl;
1788 std::vector<G4KineticTrack *>::iterator iter;
1807 #ifdef debug_BIC_StepParticlesOut
1813 #ifdef debug_BIC_StepParticlesOut
1824 #ifdef debug_BIC_return
1825 G4cout <<
"return @ Z=0 after collision loop "<<
G4endl;
1858 if (
std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1865 std::vector<G4KineticTrack *>::iterator titer;
1866 for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
1884 for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1886 G4int PDGcode=(*i)->GetDefinition()->GetPDGEncoding();
1889 if (
std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1891 resonances.push_back(*i);
1894 if ( resonances.size() > 0 )
1896 G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
1897 for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
1901 G4double newEnergy=mom.
e() + delta_Fermi;
1902 G4double newEnergy2= newEnergy*newEnergy;
1904 if ( newEnergy2 < mass2 )
1926 #ifdef debug_BIC_CorrectFinalPandE
1932 G4KineticTrackVector::iterator i;
1934 if ( pNucleus.
e() == 0 )
return;
1935 #ifdef debug_BIC_CorrectFinalPandE
1942 pFinals += (*i)->Get4Momentum();
1944 #ifdef debug_BIC_CorrectFinalPandE
1945 G4cout <<
"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName()
1946 <<
" 4mom " << (*i)->Get4Momentum()<<
G4endl;
1949 #ifdef debug_BIC_CorrectFinalPandE
1950 G4cout <<
"CorrectFinalPandE pN pF: " <<pNucleus <<
" " <<pFinals <<
G4endl;
1956 #ifdef debug_BIC_CorrectFinalPandE
1957 G4cout <<
"CorrectFinalPandE pCM, CMS pCM " << pCM <<
" " <<toCMS*pCM<<
G4endl;
1958 G4cout <<
"CorrectFinal CMS pN pF " <<toCMS*pNucleus <<
" "
1961 <<
" massInNucleus m(nucleus) m(finals) std::sqrt(s): " <<
massInNucleus <<
" " <<pNucleus.
mag()<<
" "
1962 << pFinals.mag() <<
" " << pCM.
mag() <<
G4endl;
1970 if( s0-(m10+m20)*(m10+m20) < 0 )
1972 #ifdef debug_BIC_CorrectFinalPandE
1973 G4cout <<
"G4BinaryCascade::CorrectFinalPandE() : error! " <<
G4endl;
1975 G4cout <<
"not enough mass to correct: mass^2, A,Z, mass(nucl), mass(finals) "
1976 << (s0-(m10+m20)*(m10+m20)) <<
" "
1978 << m10 <<
" " << m20
1988 G4double pInCM = std::sqrt((s0-(m10+m20)*(m10+m20))*(s0-(m10-m20)*(m10-m20))/(4.*s0));
1989 #ifdef debug_BIC_CorrectFinalPandE
1990 G4cout <<
" CorrectFinalPandE pInCM new, CURRENT, ratio : " << pInCM
1991 <<
" " << (pFinals).vect().mag()<<
" " << pInCM/(pFinals).vect().mag() <<
G4endl;
1993 if ( pFinals.vect().mag() > pInCM )
2003 G4ThreeVector p3(factor*(toCMS*(*i)->Get4Momentum()).vect());
2007 #ifdef debug_BIC_CorrectFinalPandE
2010 (*i)->Set4Momentum(p);
2012 #ifdef debug_BIC_CorrectFinalPandE
2015 <<
" CMS pFinals , mag, 3.mag : " << qFinals <<
" " << qFinals.mag() <<
" " << qFinals.vect().mag()<<
G4endl;
2016 G4cerr <<
" -CorrectFinalPandE 5 " << factor <<
G4endl;
2019 #ifdef debug_BIC_CorrectFinalPandE
2020 else {
G4cerr <<
" -CorrectFinalPandE 6 - no correction done" <<
G4endl; }
2032 std::vector<G4KineticTrack *>::iterator iter1, iter2;
2037 if(!oldSecondaries->empty())
2039 for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end();
2054 if(oldTarget->size()!=0)
2058 for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1)
2070 if(!newSecondaries->empty())
2073 for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end();
2115 #ifdef debug_BIC_DoTimeStep
2118 G4cerr <<
"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep
2125 std::vector<G4KineticTrack *>::iterator iter;
2138 #ifdef debug_BIC_DoTimeStep
2151 #ifdef debug_BIC_DoTimeStep
2161 std::for_each( kt_outside->begin(),kt_outside->end(),
2168 std::for_each( kt_inside->begin(),kt_inside->end(),
2178 kt_gone_in->clear();
2179 std::for_each( kt_outside->begin(),kt_outside->end(),
2182 kt_gone_out->clear();
2183 std::for_each( kt_inside->begin(),kt_inside->end(),
2186 #ifdef debug_BIC_DoTimeStep
2187 PrintKTVector(fail,std::string(
" Failed to go in/out -> miss_nucleus/captured"));
2188 PrintKTVector(kt_gone_in, std::string(
"recreated kt_gone_in"));
2189 PrintKTVector(kt_gone_out, std::string(
"recreated kt_gone_out"));
2195 std::for_each( kt_outside->begin(),kt_outside->end(),
2198 std::for_each( kt_outside->begin(),kt_outside->end(),
2201 #ifdef debug_BIC_DoTimeStep
2202 PrintKTVector(kt_gone_out, std::string(
"append gone_outs to final state.. theFinalState"));
2206 kt_gone_out->begin(),kt_gone_out->end());
2218 if (kt_gone_out->size() )
2221 iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary);
2222 if ( iter != kt_gone_out->end() )
2225 #ifdef debug_BIC_DoTimeStep
2226 G4cout <<
" DoTimeStep - WARNING: deleting current collision!" <<
G4endl;
2230 if ( kt_captured->size() )
2233 iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary);
2234 if ( iter != kt_captured->end() )
2237 #ifdef debug_BIC_DoTimeStep
2238 G4cout <<
" DoTimeStep - WARNING: deleting current collision!" <<
G4endl;
2248 if ( kt_captured->size() )
2251 kt_captured->begin(),kt_captured->end());
2255 std::vector<G4KineticTrack *>::iterator i_captured;
2256 for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++)
2258 (*i_captured)->Hit();
2264 #ifdef debug_G4BinaryCascade
2274 <<
" sum(tgt,capt,active) "
2304 std::vector<G4KineticTrack *>::iterator iter;
2309 G4int secondaries_in(0);
2310 G4int secondaryBarions_in(0);
2311 G4int secondaryCharge_in(0);
2314 for ( iter =in->begin(); iter != in->end(); ++iter)
2317 secondaryCharge_in +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
2318 if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 )
2320 secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber();
2324 secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass();
2340 G4double correction= secondaryMass_in + mass_initial - mass_final;
2341 if (secondaries_in>1)
2342 {correction /= secondaries_in;}
2344 #ifdef debug_BIC_CorrectBarionsOnBoundary
2345 G4cout <<
"CorrectBarionsOnBoundary,currentZ,currentA,"
2346 <<
"secondaryCharge_in,secondaryBarions_in,"
2347 <<
"energy correction,m_secondry,m_nucl_init,m_nucl_final "
2349 << secondaryCharge_in<<
" "<<secondaryBarions_in<<
" "
2350 << correction <<
" "
2351 << secondaryMass_in <<
" "
2352 << mass_initial <<
" "
2353 << mass_final <<
" "
2357 for ( iter = in->begin(); iter != in->end(); ++iter)
2359 if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
2361 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
2368 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + barrier);
2370 kt_fail->push_back(*iter);
2372 currentA -= (*iter)->GetDefinition()->GetBaryonNumber();
2377 #ifdef debug_BIC_CorrectBarionsOnBoundary
2378 G4cout <<
" CorrectBarionsOnBoundary, aft, Z, A, sec-Z,A,m,m_in_nucleus "
2380 << secondaryCharge_in <<
" " << secondaryBarions_in <<
" "
2381 << secondaryMass_in <<
" "
2390 G4int secondaries_out(0);
2391 G4int secondaryBarions_out(0);
2392 G4int secondaryCharge_out(0);
2395 for ( iter =out->begin(); iter != out->end(); ++iter)
2398 secondaryCharge_out +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
2399 if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 )
2401 secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber();
2405 secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass();
2423 G4cerr <<
"G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " <<
2424 secondaryBarions_out <<
" " << secondaryCharge_out <<
G4endl;
2429 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error");
2432 G4double correction= mass_initial - mass_final - secondaryMass_out;
2435 if (secondaries_out>1) correction /= secondaries_out;
2436 #ifdef debug_BIC_CorrectBarionsOnBoundary
2437 G4cout <<
"DoTimeStep,(current Z,A),"
2438 <<
"(secondaries out,Charge,Barions),"
2439 <<
"* energy correction,(m_secondry,m_nucl_init,m_nucl_final) "
2441 << secondaries_out <<
","
2442 << secondaryCharge_out<<
","<<secondaryBarions_out<<
") * "
2443 << correction <<
" ("
2444 << secondaryMass_out <<
", "
2445 << mass_initial <<
", "
2446 << mass_final <<
")"
2451 for ( iter = out->begin(); iter != out->end(); ++iter)
2453 if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
2455 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
2466 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() - barrier);
2468 kt_fail->push_back(*iter);
2470 currentA += (*iter)->GetDefinition()->GetBaryonNumber();
2472 #ifdef debug_BIC_CorrectBarionsOnBoundary
2475 G4cout <<
"Not correcting outgoing " << *iter <<
" "
2476 << (*iter)->GetDefinition()->GetPDGEncoding() <<
" "
2477 << (*iter)->GetDefinition()->GetParticleName() <<
G4endl;
2478 PrintKTVector(out,std::string(
"outgoing, one not corrected"));
2484 #ifdef debug_BIC_CorrectBarionsOnBoundary
2486 G4cout <<
" DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta "
2488 << secondaryCharge_out <<
" "<< secondaryBarions_out <<
" "<<
2489 secondaryMass_out <<
" "
2507 #ifdef debug_BIC_FindFragments
2508 G4cout <<
"target, captured, secondary: "
2517 G4KineticTrackVector::iterator i;
2520 if(
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus) == 1 )
2526 G4int zCaptured = 0;
2530 CapturedMomentum += (*i)->Get4Momentum();
2531 if(
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus) == 1 )
2537 G4int z = zTarget+zCaptured;
2539 #ifdef debug_G4BinaryCascade
2542 G4cout <<
" FindFragment Counting error z a " << z <<
" " <<a <<
" "
2560 if ( z < 1 )
return 0;
2564 #ifdef debug_BIC_FindFragments
2565 G4cout <<
"Fragment: a= " << a <<
" z= " << z <<
" particles= " << excitons
2566 <<
" Charged= " << zCaptured <<
" holes= " << holes
2593 final4Momentum -= (*i)->Get4Momentum();
2594 finals += (*i)->Get4Momentum();
2597 if(final4Momentum.
e()> 0 && (final4Momentum.
vect()/final4Momentum.
e()).mag()>1.0 &&
currentA > 0)
2599 #ifdef debug_BIC_Final4Momentum
2601 G4cerr <<
"G4BinaryCascade::GetFinal4Momentum - Fatal"<<
G4endl;
2602 G4KineticTrackVector::iterator i;
2603 G4cerr <<
"Total initial 4-momentum " << theProjectile4Momentum <<
G4endl;
2607 G4cerr <<
" Final state: "<<(*i)->Get4Momentum()<<(*i)->GetDefinition()->GetParticleName()<<
G4endl;
2610 G4cerr<<
" Final4Momentum = "<<final4Momentum <<
" "<<final4Momentum.
m()<<
G4endl;
2617 return final4Momentum;
2628 G4KineticTrackVector::iterator i;
2632 CapturedMomentum += (*i)->Get4Momentum();
2638 if ( NucleusMomentum.
e() > 0 )
2643 if(boost.
mag2()>1.0)
2645 # ifdef debug_BIC_FinalNucleusMomentum
2646 G4cerr <<
"G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<
G4endl;
2648 G4cerr <<
"it 01"<<NucleusMomentum<<
" "<<CapturedMomentum<<
" "<<
G4endl;
2656 #ifdef debug_debug_BIC_FinalNucleusMomentum
2657 G4cout <<
"GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<
" "<<CapturedMomentum<<
" "<<
G4endl;
2659 NucleusMomentum *= nucleusBoost;
2660 #ifdef debug_BIC_FinalNucleusMomentum
2661 G4cout <<
"GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<
G4endl;
2664 return NucleusMomentum;
2681 std::vector<G4KineticTrack *>::iterator iter, jter;
2686 while(!done && tryCount++ <200)
2694 #ifdef debug_H1_BinaryCascade
2697 for(
size_t ss=0; secs && ss<secs->size(); ss++)
2700 if((*secs)[ss]->GetDefinition()->IsShortLived()) done =
true;
2708 for(
size_t current=0; secs && current<secs->size(); current++)
2710 if((*secs)[current]->GetDefinition()->IsShortLived())
2714 for(jter=dec->begin(); jter != dec->end(); jter++)
2717 secs->push_back(*jter);
2720 delete (*secs)[current];
2732 #ifdef debug_H1_BinaryCascade
2742 products->push_back(aNew);
2743 #ifdef debug_H1_BinaryCascade
2748 G4cout <<
"final shortlived : ";
2751 G4cout <<
"final un stable : ";
2783 }
while (
sqr(x1) +
sqr(x2) > 1.);
2809 std::vector<G4KineticTrack *>::iterator i;
2810 for(i = ktv->begin(); i != ktv->end(); ++i)
2819 std::vector<G4ReactionProduct *>::iterator i;
2820 for(i = rpv->begin(); i != rpv->end(); ++i)
2829 if (comment.size() > 0 )
G4cout <<
"G4BinaryCascade::PrintKTVector() " << comment <<
G4endl;
2831 G4cout <<
" vector: " << ktv <<
", number of tracks: " << ktv->size()
2833 std::vector<G4KineticTrack *>::iterator i;
2836 for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count)
2839 G4cout <<
" track n. " << count;
2843 G4cout <<
"G4BinaryCascade::PrintKTVector():No KineticTrackVector given " <<
G4endl;
2850 if (comment.size() > 0 )
G4cout <<
"G4BinaryCascade::PrintKTVector() "<< comment <<
G4endl;
2863 G4cout <<
"G4BinaryCascade::PrintKTVector(): No Kinetictrack given" <<
G4endl;
2873 if ( Z > 0 && A >= Z )
2877 }
else if ( A > 0 && Z>0 )
2882 }
else if ( A >= 0 && Z<=0 )
2887 }
else if ( A == 0 )
2894 G4cerr <<
"G4BinaryCascade::GetIonMass() - invalid (A,Z) = ("
2895 << A <<
"," << Z <<
")" <<
G4endl;
2896 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::GetIonMass() - giving up");
2907 std::vector<G4KineticTrack *>::iterator iter;
2908 std::vector<G4ReactionProduct *>::iterator rpiter;
2914 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2917 Esecondaries +=(*iter)->Get4Momentum().e();
2918 psecondaries +=(*iter)->Get4Momentum();
2921 products->push_back(aNew);
2933 if ( lates->size() == 1 ) {
2945 products->push_back(aNew);
2967 (*iter)->Update4Momentum((*iter)->Get4Momentum().vect()+transferCorrection);
2968 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2971 Esecondaries +=(*iter)->Get4Momentum().e();
2972 psecondaries +=(*iter)->Get4Momentum();
2974 products->push_back(aNew);
2981 (*iter)->Update4Momentum((*iter)->Get4Momentum().vect()+transferCorrection);
2982 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2985 Esecondaries +=(*iter)->Get4Momentum().e();
2986 psecondaries +=(*iter)->Get4Momentum();
2988 products->push_back(aNew);
2995 SumMassNucleons += (*iter)->GetDefinition()->GetPDGMass();
2996 pNucleons += (*iter)->Get4Momentum();
3000 #ifdef debug_BIC_FillVoidnucleus
3002 psecondaries - pNucleons;
3012 for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
3013 TotalEkin+=(*rpiter)->GetKineticEnergy();
3017 correction=1. + (Ekinetic-Ekineticrdm)/TotalEkin;
3019 #ifdef debug_G4BinaryCascade
3021 G4cout <<
"BLIC::FillVoidNucleus() fail correction, Ekinetic, TotalEkin " << Ekinetic <<
""<< TotalEkin <<
G4endl;
3025 for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
3026 (*rpiter)->SetKineticEnergy((*rpiter)->GetKineticEnergy()*correction);
3027 (*rpiter)->SetMomentum((*rpiter)->GetTotalMomentum() * (*rpiter)->GetMomentum().unit());
3031 Ekinetic=Ekineticrdm*correction;
3044 products->push_back(aNew);
3049 for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
3050 psecondaries +=
G4LorentzVector((*rpiter)->GetMomentum(),(*rpiter)->GetTotalEnergy() );
3062 SumMom=initial4Mom.
vect()-SumMom;
3065 std::vector<G4ReactionProduct *>::reverse_iterator
reverse;
3066 while ( SumMom.
mag() > 0.1*
MeV && loopcount++ < 10)
3068 G4int index=products->size();
3069 for (reverse=products->rbegin(); reverse!=products->rend(); ++
reverse, --index){
3070 SumMom=initial4Mom.
vect();
3071 for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
3072 SumMom-=(*rpiter)->GetMomentum();
3075 G4double p=((*reverse)->GetMomentum()).mag();
3076 (*reverse)->SetMomentum( p*(((*reverse)->GetMomentum()+SumMom).unit()));
3087 std::vector<G4KineticTrack *>::iterator iter;
3088 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
3091 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
3098 products->push_back(aNew);
3117 if (fragment != 0) {
3125 products->push_back(theNew);
3132 G4cout <<
"Thank you for using G4BinaryCascade. "<<
G4endl;
3142 for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
3144 G4int PDGcode=
std::abs((*i)->GetDefinition()->GetPDGEncoding());
3145 if (
std::abs(PDGcode)==211 || PDGcode==111 ) havePion=
true;
3148 if ( !products || havePion)
3151 G4cout <<
" Collision " << collision <<
", type: "<<
typeid(action).
name()
3152 <<
", with NO products! " <<
G4endl;
3153 G4cout << G4endl<<
"Initial condition are these:"<<
G4endl;
3171 static G4int lastdA(0), lastdZ(0);
3178 std::vector<G4KineticTrack *>::iterator i;
3179 G4int CapturedA(0), CapturedZ(0);
3180 G4int secsA(0), secsZ(0);
3182 CapturedA += (*i)->GetDefinition()->GetBaryonNumber();
3183 CapturedZ +=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
3188 secsA += (*i)->GetDefinition()->GetBaryonNumber();
3189 secsZ +=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
3194 fStateA += (*i)->GetDefinition()->GetBaryonNumber();
3195 fStateZ +=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
3201 #ifdef debugCheckChargeAndBaryonNumberverbose
3202 G4cout << where <<
" A: iState= "<< iStateA<<
", secs= "<< secsA<<
", fState= "<< fStateA<<
", current= "<<
currentA<<
", late= " <<lateA <<
G4endl;
3203 G4cout << where <<
" Z: iState= "<< iStateZ<<
", secs= "<< secsZ<<
", fState= "<< fStateZ<<
", current= "<<
currentZ<<
", late= " <<lateZ <<
G4endl;
3206 if (deltaA != 0 || deltaZ!=0 ) {
3207 if (deltaA != lastdA || deltaZ != lastdZ ) {
3208 G4cout <<
"baryon/charge imbalance - " << where << G4endl
3209 <<
"deltaA " <<deltaA<<
", iStateA "<<iStateA<<
", CapturedA "<<CapturedA <<
", secsA "<<secsA
3210 <<
", fStateA "<<fStateA <<
", currentA "<<
currentA <<
", lateA "<<lateA << G4endl
3211 <<
"deltaZ "<<deltaZ<<
", iStateZ "<<iStateZ<<
", CapturedZ "<<CapturedZ <<
", secsZ "<<secsZ
3212 <<
", fStateZ "<<fStateZ <<
", currentZ "<<
currentZ <<
", lateZ "<<lateZ << G4endl<<
G4endl;
3216 }
else { lastdA=lastdZ=0;}
3244 <<
" " << initial <<
G4endl;;
3247 for (
unsigned int it=0;
it < ktv.size();
it++)
3260 <<
" " << initial <<
" Excit " << thisExcitation <<
G4endl;;
3265 G4int product_barions(0);
3268 for (
unsigned int it=0;
it < products->size();
it++)
3280 <<
" " <<
final <<
G4endl;;
3289 finalA -= product_barions;
3294 <<
" delta-mass " << delta<<
G4endl;
3298 <<
" " <<
final <<
" "
3311 G4ReactionProductVector::iterator iter;
3319 for(iter = products->begin(); iter != products->end(); ++iter)
3322 G4cout <<
" Secondary E - Ekin / p " <<
3323 (*iter)->GetDefinition()->GetParticleName() <<
" " <<
3324 (*iter)->GetTotalEnergy() <<
" - " <<
3325 (*iter)->GetKineticEnergy()<<
" / " <<
3326 (*iter)->GetMomentum().x() <<
" " <<
3327 (*iter)->GetMomentum().y() <<
" " <<
3328 (*iter)->GetMomentum().z() <<
G4endl;
3329 Efinal += (*iter)->GetTotalEnergy();
3330 pFinal += (*iter)->GetMomentum();
3334 G4cout <<
"BIC E/p delta " <<
3352 std::vector<G4KineticTrack *>::iterator ktiter;
3356 G4cout <<
" Secondary E - Ekin / p " <<
3357 (*ktiter)->GetDefinition()->GetParticleName() <<
" " <<
3358 (*ktiter)->Get4Momentum().e() <<
" - " <<
3359 (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() <<
" / " <<
3360 (*ktiter)->Get4Momentum().vect() <<
G4endl;
3361 psecs += (*ktiter)->Get4Momentum();
3367 G4cout <<
" Target E - Ekin / p " <<
3368 (*ktiter)->GetDefinition()->GetParticleName() <<
" " <<
3369 (*ktiter)->Get4Momentum().e() <<
" - " <<
3370 (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() <<
" / " <<
3371 (*ktiter)->Get4Momentum().vect() <<
G4endl;
3372 ptgts += (*ktiter)->Get4Momentum();
3378 G4cout <<
" Captured E - Ekin / p " <<
3379 (*ktiter)->GetDefinition()->GetParticleName() <<
" " <<
3380 (*ktiter)->Get4Momentum().e() <<
" - " <<
3381 (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() <<
" / " <<
3382 (*ktiter)->Get4Momentum().vect() <<
G4endl;
3383 pcpts += (*ktiter)->Get4Momentum();
3389 G4cout <<
" Finals E - Ekin / p " <<
3390 (*ktiter)->GetDefinition()->GetParticleName() <<
" " <<
3391 (*ktiter)->Get4Momentum().e() <<
" - " <<
3392 (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() <<
" / " <<
3393 (*ktiter)->Get4Momentum().vect() <<
G4endl;
3394 pfins += (*ktiter)->Get4Momentum();
3397 G4cout <<
" Secondaries " << psecs <<
", Targets " << ptgts << G4endl
3398 <<
" Captured " << pcpts <<
", Finals " << pfins << G4endl