93 #if defined(debug_G4BinaryCascade)
94 #define _CheckChargeAndBaryonNumber_(val) CheckChargeAndBaryonNumber(val)
96 #define _CheckChargeAndBaryonNumber_(val)
111 theImR.push_back(theDecay);
114 theImR.push_back(aAb);
117 theImR.push_back(aSc);
123 theCutOnPAbsorb= 0*
MeV;
137 thePrimaryEscape =
true;
146 projectileA=projectileZ=0;
147 currentInitialEnergy=initial_nuclear_mass=0.;
161 ClearAndDestroy(&theTargetList);
162 ClearAndDestroy(&theSecondaryList);
163 ClearAndDestroy(&theCapturedList);
164 delete thePropagator;
165 delete theCollisionMgr;
167 delete theLateParticle;
169 delete theH1Scatterer;
174 outFile <<
"G4BinaryCascade is an intra-nuclear cascade model in which\n"
175 <<
"an incident hadron collides with a nucleon, forming two\n"
176 <<
"final-state particles, one or both of which may be resonances.\n"
177 <<
"The resonances then decay hadronically and the decay products\n"
178 <<
"are then propagated through the nuclear potential along curved\n"
179 <<
"trajectories until they re-interact or leave the nucleus.\n"
180 <<
"This model is valid for incident pions up to 1.5 GeV and\n"
181 <<
"nucleons up to 10 GeV.\n";
185 outFile <<
"G4BinaryCascade propagtes secondaries produced by a high\n"
186 <<
"energy model through the wounded nucleus.\n"
187 <<
"Secondaries are followed after the formation time and if\n"
188 <<
"within the nucleus are propagated through the nuclear\n"
189 <<
"potential along curved trajectories until they interact\n"
190 <<
"with a nucleon, decay, or leave the nucleus.\n"
191 <<
"An interaction of a secondary with a nucleon produces two\n"
192 <<
"final-state particles, one or both of which may be resonances.\n"
193 <<
"Resonances decay hadronically and the decay products\n"
194 <<
"are in turn propagated through the nuclear potential along curved\n"
195 <<
"trajectories until they re-interact or leave the nucleus.\n"
196 <<
"This model is valid for pions up to 1.5 GeV and\n"
197 <<
"nucleons up to about 3.5 GeV.\n";
212 if(getenv(
"BCDEBUG") )
G4cerr <<
" ######### Binary Cascade Reaction starts ######### "<<
G4endl;
217 if(initial4Momentum.
e()-initial4Momentum.
m()<theBCminP &&
231 if(!getenv(
"I_Am_G4BinaryCascade_Developer") )
238 G4cerr <<
"You are trying to use G4BinaryCascade with " <<definition->GetParticleName()<<
" as projectile."<<
G4endl;
239 G4cerr <<
"G4BinaryCascade should not be used for projectiles other than nucleons or pions."<<
G4endl;
240 G4cerr <<
"If you want to continue, please switch on the developer environment: "<<
G4endl;
241 G4cerr <<
"setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<
G4endl;
242 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade - used for unvalid particle type - Fatal");
247 thePrimaryType = definition;
248 thePrimaryEscape =
false;
252 G4int interactionCounter = 0;
260 ClearAndDestroy(products);
274 initialPosition=GetSpherePoint(1.1*radius, initial4Momentum);
275 kt =
new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
279 secondaries->push_back(kt);
286 }
while(! products );
288 if(++interactionCounter>99)
break;
289 }
while(products->size() == 0);
291 if(products->size()>0)
297 G4ReactionProductVector::iterator iter;
299 for(iter = products->begin(); iter != products->end(); ++iter)
303 (*iter)->GetTotalEnergy(),
304 (*iter)->GetMomentum());
312 if(getenv(
"BCDEBUG") )
G4cerr <<
" ######### Binary Cascade Reaction void, return intial state ######### "<<
G4endl;
318 ClearAndDestroy(products);
324 if(getenv(
"BCDEBUG") )
G4cerr <<
" ######### Binary Cascade Reaction ends ######### "<<
G4endl;
334 #ifdef debug_BIC_Propagate
335 G4cout <<
"G4BinaryCascade Propagate starting -------------------------------------------------------" <<
G4endl;
344 ClearAndDestroy(&theCapturedList);
345 ClearAndDestroy(&theSecondaryList);
346 theSecondaryList.clear();
347 ClearAndDestroy(&theFinalState);
348 std::vector<G4KineticTrack *>::iterator iter;
359 #ifdef debug_BIC_GetExcitationEnergy
360 G4cout <<
"ExcitationEnergy0 " << GetExcitationEnergy() <<
G4endl;
365 G4bool success = BuildLateParticleCollisions(secondaries);
368 products=HighEnergyModelFSProducts(products, secondaries);
369 ClearAndDestroy(secondaries);
372 #ifdef debug_G4BinaryCascade
373 G4cout <<
"G4BinaryCascade::Propagate: warning - high energy model failed energy conservation, returning unchanged high energy final state" <<
G4endl;
383 FindCollisions(&theSecondaryList);
386 if(theCollisionMgr->
Entries() == 0 )
390 #ifdef debug_BIC_return
400 G4bool haveProducts =
false;
401 G4int collisionCount=0;
403 while(theCollisionMgr->
Entries() > 0 && currentZ)
416 if(theCollisionMgr->
Entries() > 0)
420 #ifdef debug_BIC_Propagate_Collisions
421 G4cout <<
" NextCollision * , Time, curtime = " << nextCollision <<
" "
436 if (ApplyCollision(nextCollision))
455 products = FillVoidNucleusProducts(products);
456 #ifdef debug_BIC_return
457 G4cout <<
"return @ Z=0 after collision loop "<<
G4endl;
458 PrintKTVector(&theSecondaryList,std::string(
" theSecondaryList"));
459 G4cout <<
"theTargetList size: " << theTargetList.size() <<
G4endl;
460 PrintKTVector(&theTargetList,std::string(
" theTargetList"));
461 PrintKTVector(&theCapturedList,std::string(
" theCapturedList"));
463 G4cout <<
" ExcitE be4 Correct : " <<GetExcitationEnergy() <<
G4endl;
464 G4cout <<
" Mom Transfered to nucleus : " << theMomentumTransfer <<
" " << theMomentumTransfer.
mag() <<
G4endl;
465 PrintKTVector(&theFinalState,std::string(
" FinalState uncorrected"));
466 G4cout <<
"returned products: " << products->size() <<
G4endl;
485 #ifdef debug_BIC_return
492 #ifdef debug_BIC_Propagate
493 G4cout <<
" Momentum transfer to Nucleus " << theMomentumTransfer <<
" " << theMomentumTransfer.
mag() <<
G4endl;
500 if ( theSecondaryList.size() > 0 )
502 #ifdef debug_G4BinaryCascade
503 G4cerr <<
"G4BinaryCascade: Warning, have active particles at end" <<
G4endl;
504 PrintKTVector(&theSecondaryList,
"active particles @ end added to theFinalState");
507 for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
509 theFinalState.push_back(*iter);
511 theSecondaryList.clear();
514 while ( theCollisionMgr->
Entries() > 0 )
516 #ifdef debug_G4BinaryCascade
517 G4cerr <<
" Warning: remove left over collision(s) " <<
G4endl;
522 #ifdef debug_BIC_Propagate_Excitation
524 PrintKTVector(&theSecondaryList,std::string(
" theSecondaryList"));
525 G4cout <<
"theTargetList size: " << theTargetList.size() <<
G4endl;
527 PrintKTVector(&theCapturedList,std::string(
" theCapturedList"));
529 G4cout <<
" ExcitE be4 Correct : " <<GetExcitationEnergy() <<
G4endl;
530 G4cout <<
" Mom Transfered to nucleus : " << theMomentumTransfer <<
" " << theMomentumTransfer.
mag() <<
G4endl;
531 PrintKTVector(&theFinalState,std::string(
" FinalState uncorrected"));
537 G4double ExcitationEnergy=GetExcitationEnergy();
539 #ifdef debug_BIC_Propagate_finals
540 PrintKTVector(&theFinalState,std::string(
" FinalState be4 corr"));
541 G4cout <<
" Excitation Energy prefinal, #collisions:, out, captured "
542 << ExcitationEnergy <<
" "
543 << collisionCount <<
" "
544 << theFinalState.size() <<
" "
545 << theCapturedList.size()<<
G4endl;
548 if (ExcitationEnergy < 0 )
550 G4int maxtry=5, ntry=0;
553 ExcitationEnergy=GetExcitationEnergy();
554 }
while ( ++ntry < maxtry && ExcitationEnergy < 0 );
557 #ifdef debug_BIC_Propagate_finals
558 PrintKTVector(&theFinalState,std::string(
" FinalState corrected"));
559 G4cout <<
" Excitation Energy final, #collisions:, out, captured "
560 << ExcitationEnergy <<
" "
561 << collisionCount <<
" "
562 << theFinalState.size() <<
" "
563 << theCapturedList.size()<<
G4endl;
567 if ( ExcitationEnergy < 0. )
582 ClearAndDestroy(products);
584 #ifdef debug_BIC_return
595 products= ProductsAddFinalState(products, theFinalState);
597 products= ProductsAddPrecompound(products, precompoundProducts);
602 thePrimaryEscape =
true;
604 #ifdef debug_BIC_return
613 G4double G4BinaryCascade::GetExcitationEnergy()
618 #if defined(debug_G4BinaryCascade) || defined(debug_BIC_GetExcitationEnergy)
619 G4int finalA = theTargetList.size()+theCapturedList.size();
620 G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList);
621 if ( (currentA - finalA) != 0 || (currentZ - finalZ) != 0 )
623 G4cerr <<
"G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} "
624 << currentA <<
" " << finalA <<
" "<< currentZ <<
" " << finalZ <<
G4endl;
635 else if (currentZ==0 )
638 else {nucleusMass = GetFinalNucleusMomentum().
mag()
643 #ifdef debug_G4BinaryCascade
644 G4cout <<
"G4BinaryCascade::GetExcitationEnergy(): Warning - invalid nucleus (A,Z)=("
645 << currentA <<
"," << currentZ <<
")" <<
G4endl;
650 #ifdef debug_BIC_GetExcitationEnergy
652 debug.push_back(
"====> current A, Z");
653 debug.push_back(currentZ);
654 debug.push_back(currentA);
655 debug.push_back(
"====> final A, Z");
656 debug.push_back(finalZ);
657 debug.push_back(finalA);
658 debug.push_back(nucleusMass);
659 debug.push_back(GetFinalNucleusMomentum().mag());
665 excitationE = GetFinalNucleusMomentum().
mag() - nucleusMass;
671 #ifdef debug_BIC_GetExcitationEnergy
673 if ( excitationE < 0 )
675 G4cout <<
"negative ExE final Ion mass " <<nucleusMass<<
G4endl;
677 if(finalZ>.5)
G4cout <<
" Final nuclmom/mass " << Nucl_mom <<
" " << Nucl_mom.
mag()
678 <<
" (A,Z)=("<< finalA <<
","<<finalZ <<
")"
679 <<
" mass " << nucleusMass <<
" "
680 <<
" excitE " << excitationE <<
G4endl;
688 initialExc = theInitial4Mom.
mag()-
690 G4cout <<
"GetExcitationEnergy: Initial nucleus A Z " << A <<
" " << Z <<
" " << initialExc <<
G4endl;
707 void G4BinaryCascade::BuildTargetList()
718 ClearAndDestroy(&theTargetList);
746 theTargetList.push_back(kt);
751 #ifdef debug_BIC_BuildTargetList
752 else {
G4cout <<
"nucleon is hit" << nucleon <<
G4endl;}
759 }
else if (currentZ==0 && currentA>=1 )
764 G4cerr <<
"G4BinaryCascade::BuildTargetList(): Fatal Error - invalid nucleus (A,Z)=("
765 << currentA <<
"," << currentZ <<
")" <<
G4endl;
768 currentInitialEnergy= theInitial4Mom.
e() + theProjectile4Momentum.
e();
770 #ifdef debug_BIC_BuildTargetList
771 G4cout <<
"G4BinaryCascade::BuildTargetList(): nucleus (A,Z)=("
772 << currentA <<
"," << currentZ <<
") mass: " << massInNucleus <<
773 ", theInitial4Mom " << theInitial4Mom <<
774 ", currentInitialEnergy " << currentInitialEnergy <<
G4endl;
784 std::vector<G4KineticTrack *>::iterator iter;
787 projectileA=projectileZ=0;
790 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
792 if((*iter)->GetFormationTime() < StartingTime)
793 StartingTime = (*iter)->GetFormationTime();
798 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
802 G4double FormTime = (*iter)->GetFormationTime() - StartingTime;
803 (*iter)->SetFormationTime(FormTime);
806 FindLateParticleCollision(*iter);
807 lateParticles4Momentum += (*iter)->GetTrackingMomentum();
808 lateA += (*iter)->GetDefinition()->GetBaryonNumber();
809 lateZ +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
813 theSecondaryList.push_back(*iter);
815 theProjectile4Momentum += (*iter)->GetTrackingMomentum();
816 projectileA += (*iter)->GetDefinition()->GetBaryonNumber();
817 projectileZ +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
818 #ifdef debug_BIC_Propagate
819 G4cout <<
" Adding initial secondary " << *iter
820 <<
" time" << (*iter)->GetFormationTime()
821 <<
", state " << (*iter)->GetState() <<
G4endl;
830 theProjectile4Momentum += mom;
834 G4double excitation= theProjectile4Momentum.
e() + initial_nuclear_mass - lateParticles4Momentum.e() - massInNucleus;
835 #ifdef debug_BIC_GetExcitationEnergy
836 G4cout <<
"BIC: Proj.e, nucl initial, nucl final, lateParticles"
837 << theProjectile4Momentum <<
", "
838 << initial_nuclear_mass<<
", " << massInNucleus <<
", "
839 << lateParticles4Momentum <<
G4endl;
840 G4cout <<
"BIC: Proj.e / initial excitation: " << theProjectile4Momentum.
e() <<
" / " << excitation <<
G4endl;
842 success = excitation > 0;
843 #ifdef debug_G4BinaryCascade
845 G4cout <<
"G4BinaryCascade::BuildLateParticleCollisions(): Proj.e / initial excitation: " << theProjectile4Momentum.
e() <<
" / " << excitation <<
G4endl;
855 secondaries->clear();
874 fragment = FindFragments();
877 if(fragment->
GetA() >1.5)
885 else if (theExcitationHandler)
887 precompoundProducts=theExcitationHandler->
BreakItUp(*fragment);
892 if (theTargetList.size() + theCapturedList.size() > 1 ) {
896 std::vector<G4KineticTrack *>::iterator i;
897 if ( theTargetList.size() == 1 ) {i=theTargetList.begin();}
898 if ( theCapturedList.size() == 1 ) {i=theCapturedList.begin();}
903 precompoundProducts->push_back(aNew);
911 precompoundProducts = DecayVoidNucleus();
913 return precompoundProducts;
921 if ( (theTargetList.size()+theCapturedList.size()) > 0 )
924 std::vector<G4KineticTrack *>::iterator aNuc;
926 std::vector<G4double> masses;
929 if ( theTargetList.size() != 0)
931 for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++)
933 G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
934 masses.push_back(mass);
939 if ( theCapturedList.size() != 0)
941 for(aNuc = theCapturedList.begin();
942 aNuc != theCapturedList.end(); aNuc++)
944 G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
945 masses.push_back(mass);
956 if ( eCMS < sumMass )
958 eCMS=sumMass + 2*
MeV*masses.size();
963 std::vector<G4LorentzVector*> * momenta=decay.
Decay(eCMS,masses);
964 std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
967 if ( theTargetList.size() != 0)
969 for ( aNuc=theTargetList.begin();
970 (aNuc != theTargetList.end()) && (aMom!=momenta->end());
976 result->push_back(aNew);
982 if ( theCapturedList.size() != 0)
984 for ( aNuc=theCapturedList.begin();
985 (aNuc != theCapturedList.end()) && (aMom!=momenta->end());
989 (*aNuc)->GetDefinition());
992 result->push_back(aNew);
1008 for(i = 0; i< fs.size(); i++)
1015 products->push_back(aNew);
1017 #ifdef debug_BIC_Propagate_finals
1033 if ( precompoundProducts )
1035 std::vector<G4ReactionProduct *>::iterator j;
1036 for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
1041 #ifdef debug_BIC_Propagate_finals
1044 pProduct *= precompoundLorentzboost;
1045 #ifdef debug_BIC_Propagate_finals
1048 pSumPreco += pProduct;
1049 (*j)->SetTotalEnergy(pProduct.e());
1050 (*j)->SetMomentum(pProduct.vect());
1051 (*j)->SetNewlyAdded(
true);
1052 products->push_back(*j);
1056 precompoundProducts->clear();
1057 delete precompoundProducts;
1065 for(std::vector<G4KineticTrack *>::iterator i = secondaries->begin();
1066 i != secondaries->end(); ++i)
1069 for(std::vector<G4BCAction *>::iterator j = theImR.begin();
1070 j!=theImR.end(); j++)
1073 const std::vector<G4CollisionInitialState *> & aCandList
1074 = (*j)->GetCollisions(*i, theTargetList, theCurrentTime);
1075 for(
size_t count=0; count<aCandList.size(); count++)
1087 void G4BinaryCascade::FindDecayCollision(
G4KineticTrack * secondary)
1090 const std::vector<G4CollisionInitialState *> & aCandList
1091 = theDecay->
GetCollisions(secondary, theTargetList, theCurrentTime);
1092 for(
size_t count=0; count<aCandList.size(); count++)
1099 void G4BinaryCascade::FindLateParticleCollision(
G4KineticTrack * secondary)
1104 if (((
G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(secondary,tin,tout))
1109 }
else if ( tout > 0 )
1122 #ifdef debug_BIC_FindCollision
1123 G4cout <<
"FindLateP Particle, 4-mom, times newState "
1126 <<
" times " << tin <<
" " << tout <<
" "
1130 const std::vector<G4CollisionInitialState *> & aCandList
1131 = theLateParticle->
GetCollisions(secondary, theTargetList, theCurrentTime);
1132 for(
size_t count=0; count<aCandList.size(); count++)
1134 #ifdef debug_BIC_FindCollision
1135 G4cout <<
" Adding a late Col : " << aCandList[count] <<
G4endl;
1148 #ifdef debug_BIC_ApplyCollision
1149 G4cerr <<
"G4BinaryCascade::ApplyCollision start"<<
G4endl;
1150 theCollisionMgr->
Print();
1155 G4bool haveTarget=target_collection.size()>0;
1158 #ifdef debug_G4BinaryCascade
1159 G4cout <<
"G4BinaryCasacde::ApplyCollision(): StateError " << primary <<
G4endl;
1160 PrintKTVector(primary,std::string(
"primay- ..."));
1161 PrintKTVector(&target_collection,std::string(
"... targets"));
1164 theCollisionMgr->
Print();
1181 G4int initialBaryon(0);
1182 G4int initialCharge(0);
1190 G4double initial_Efermi=CorrectShortlivedPrimaryForFermi(primary,target_collection);
1196 #ifdef debug_BIC_ApplyCollision
1197 DebugApplyCollisionFail(collision, products);
1203 G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
1204 G4bool decayCollision= (!haveTarget) && products && products->size() > 1;
1208 #ifdef debug_G4BinaryCascade
1209 G4int lateBaryon(0), lateCharge(0);
1212 if ( lateParticleCollision )
1216 #ifdef debug_G4BinaryCascade
1217 lateBaryon = initialBaryon;
1218 lateCharge = initialCharge;
1220 initialBaryon=initialCharge=0;
1227 if (!lateParticleCollision)
1229 if( !products || products->size()==0 || !CheckPauliPrinciple(products) )
1231 #ifdef debug_BIC_ApplyCollision
1232 if (products)
G4cout <<
" ======Failed Pauli =====" <<
G4endl;
1233 G4cerr <<
"G4BinaryCascade::ApplyCollision blocked"<<
G4endl;
1241 if (! CorrectShortlivedFinalsForFermi(products, initial_Efermi)){
1247 #ifdef debug_BIC_ApplyCollision
1248 DebugApplyCollision(collision, products);
1252 if (products) ClearAndDestroy(products);
1253 if ( decayCollision ) FindDecayCollision(primary);
1259 G4int finalBaryon(0);
1260 G4int finalCharge(0);
1262 for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1264 if ( ! lateParticleCollision )
1266 (*i)->SetState(primary->
GetState());
1268 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1269 finalCharge+=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
1272 if (((
G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout) &&
1273 tin < 0 && tout > 0 )
1275 PrintKTVector((*i),
"particle inside marked not-inside");
1276 G4cout <<
"tin tout: " << tin <<
" " << tout <<
G4endl;
1281 if (((
G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout))
1288 else if ( tout > 0 )
1291 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1292 finalCharge+=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
1297 toFinalState.push_back((*i));
1303 toFinalState.push_back((*i));
1308 if(!toFinalState.empty())
1310 theFinalState.insert(theFinalState.end(),
1311 toFinalState.begin(),toFinalState.end());
1312 std::vector<G4KineticTrack *>::iterator iter1, iter2;
1313 for(iter1 = toFinalState.begin(); iter1 != toFinalState.end();
1316 iter2 = std::find(products->begin(), products->end(),
1318 if ( iter2 != products->end() ) products->erase(iter2);
1324 currentA += finalBaryon-initialBaryon;
1325 currentZ += finalCharge-initialCharge;
1329 oldSecondaries.push_back(primary);
1332 #ifdef debug_G4BinaryCascade
1333 if ( (finalBaryon-initialBaryon-lateBaryon) != 0 || (finalCharge-initialCharge-lateCharge) != 0 )
1335 G4cout <<
"G4BinaryCascade: Error in Balancing: " <<
G4endl;
1336 G4cout <<
"initial/final baryon number, initial/final Charge "
1337 << initialBaryon <<
" "<< finalBaryon <<
" "
1338 << initialCharge <<
" "<< finalCharge <<
" "
1340 <<
", with number of products: "<< products->size() <<
G4endl;
1341 G4cout << G4endl<<
"Initial condition are these:"<<
G4endl;
1354 for(
size_t ii=0; ii< oldTarget.size(); ii++)
1356 oldTarget[ii]->Hit();
1359 UpdateTracksAndCollisions(&oldSecondaries, &oldTarget, products);
1369 G4bool G4BinaryCascade::Absorb()
1377 std::vector<G4KineticTrack *>::iterator iter;
1379 for(iter = theSecondaryList.begin();
1380 iter != theSecondaryList.end(); ++iter)
1385 if(absorber.WillBeAbsorbed(*kt))
1387 absorbList.push_back(kt);
1392 if(absorbList.empty())
1396 for(iter = absorbList.begin(); iter != absorbList.end(); ++iter)
1399 if(!absorber.FindAbsorbers(*kt, theTargetList))
1400 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1402 if(!absorber.FindProducts(*kt))
1403 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1409 while(!CheckPauliPrinciple(products))
1414 ClearAndDestroy(products);
1415 if(!absorber.FindProducts(*kt))
1417 "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1423 toRemove.push_back(kt);
1424 toDelete.push_back(kt);
1426 UpdateTracksAndCollisions(&toRemove, absorbers, products);
1427 ClearAndDestroy(absorbers);
1429 ClearAndDestroy(&toDelete);
1442 std::vector<G4KineticTrack *>::iterator i;
1447 G4int particlesAboveCut=0;
1448 G4int particlesBelowCut=0;
1449 if ( verbose )
G4cout <<
" Capture: secondaries " << theSecondaryList.size() <<
G4endl;
1450 for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1467 ++particlesBelowCut;
1475 if (verbose)
G4cout <<
"Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP "
1476 << particlesAboveCut <<
" " << particlesBelowCut <<
" " << capturedEnergy
1480 if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1483 for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1491 captured.push_back(kt);
1493 theCapturedList.push_back(kt);
1497 UpdateTracksAndCollisions(&captured, NULL, NULL);
1512 fermiMom.
Init(A, Z);
1516 G4KineticTrackVector::iterator i;
1523 for(i = products->begin(); i != products->end(); ++i)
1525 definition = (*i)->GetDefinition();
1551 if(mom.
e() < eFermi )
1560 #ifdef debug_BIC_CheckPauli
1563 for(i = products->begin(); i != products->end(); ++i)
1565 definition = (*i)->GetDefinition();
1574 if ( mom.
e()-mom.
mag()+field > 160*
MeV )
1576 G4cout <<
"momentum problem pFermi=" << pFermi
1577 <<
" mom, mom.m " << mom <<
" " << mom.
mag()
1578 <<
" field " << field <<
G4endl;
1589 void G4BinaryCascade::StepParticlesOut()
1596 while( theSecondaryList.size() > 0 )
1601 std::vector<G4KineticTrack *>::iterator i;
1602 for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1610 ((
G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(kt,tdummy,tStep);
1611 #ifdef debug_BIC_StepParticlesOut
1612 G4cout <<
" minTimeStep, tStep Particle " <<minTimeStep <<
" " <<tStep
1617 PrintKTVector(&theSecondaryList, std::string(
" state ERROR....."));
1618 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1621 if(intersect && tStep<minTimeStep && tStep> 0 )
1623 minTimeStep = tStep;
1626 PrintKTVector(&theSecondaryList, std::string(
" state ERROR....."));
1627 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1634 if(theCollisionMgr->
Entries() > 0)
1638 G4cout <<
" NextCollision * , Time= " << nextCollision <<
" "
1639 <<timeToCollision<<
G4endl;
1641 if ( timeToCollision > minTimeStep )
1643 DoTimeStep(minTimeStep);
1647 if (!DoTimeStep(timeToCollision) )
1659 if ( ApplyCollision(nextCollision))
1671 #ifdef debug_G4BinaryCascade
1672 G4cerr <<
"G4BinaryCascade.cc: Warning - aborting looping particle(s)" <<
G4endl;
1673 PrintKTVector(&theSecondaryList,
" looping particles added to theFinalState");
1677 std::vector<G4KineticTrack *>::iterator iter;
1678 for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
1680 theFinalState.push_back(*iter);
1682 theSecondaryList.clear();
1696 #ifdef debug_BIC_StepParticlesOut
1700 if ( counter > 100 && theCollisionMgr->
Entries() == 0)
1702 #ifdef debug_BIC_StepParticlesOut
1703 PrintKTVector(&theSecondaryList,std::string(
"stepping 100 steps"));
1705 FindCollisions(&theSecondaryList);
1720 G4double G4BinaryCascade::CorrectShortlivedPrimaryForFermi(
1729 if ( std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1736 std::vector<G4KineticTrack *>::iterator titer;
1737 for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
1755 for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1757 G4int PDGcode=(*i)->GetDefinition()->GetPDGEncoding();
1759 final_Efermi+=((
G4RKPropagation *)thePropagator)->GetField(PDGcode,(*i)->GetPosition());
1760 if ( std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1762 resonances.push_back(*i);
1765 if ( resonances.size() > 0 )
1767 G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
1768 for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
1772 G4double newEnergy=mom.
e() + delta_Fermi;
1773 G4double newEnergy2= newEnergy*newEnergy;
1775 if ( newEnergy2 < mass2 )
1788 void G4BinaryCascade::CorrectFinalPandE()
1796 #ifdef debug_BIC_CorrectFinalPandE
1800 if ( theFinalState.size() == 0 )
return;
1802 G4KineticTrackVector::iterator i;
1804 if ( pNucleus.
e() == 0 )
return;
1805 #ifdef debug_BIC_CorrectFinalPandE
1810 for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1812 pFinals += (*i)->Get4Momentum();
1814 #ifdef debug_BIC_CorrectFinalPandE
1815 G4cout <<
"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName()
1816 <<
" 4mom " << (*i)->Get4Momentum()<<
G4endl;
1819 #ifdef debug_BIC_CorrectFinalPandE
1820 G4cout <<
"CorrectFinalPandE pN pF: " <<pNucleus <<
" " <<pFinals <<
G4endl;
1826 #ifdef debug_BIC_CorrectFinalPandE
1827 G4cout <<
"CorrectFinalPandE pCM, CMS pCM " << pCM <<
" " <<toCMS*pCM<<
G4endl;
1828 G4cout <<
"CorrectFinal CMS pN pF " <<toCMS*pNucleus <<
" "
1830 <<
" nucleus initial mass : " <<GetFinal4Momentum().
mag()
1831 <<
" massInNucleus m(nucleus) m(finals) std::sqrt(s): " << massInNucleus <<
" " <<pNucleus.
mag()<<
" "
1832 << pFinals.mag() <<
" " << pCM.
mag() <<
G4endl;
1838 G4double m10 = GetIonMass(currentZ,currentA);
1840 if( s0-(m10+m20)*(m10+m20) < 0 )
1842 #ifdef debug_BIC_CorrectFinalPandE
1843 G4cout <<
"G4BinaryCascade::CorrectFinalPandE() : error! " <<
G4endl;
1845 G4cout <<
"not enough mass to correct: mass, A,Z, mass(nucl), mass(finals) "
1846 << std::sqrt(s0-(m10+m20)*(m10+m20)) <<
" "
1847 << currentA <<
" " << currentZ <<
" "
1848 << m10 <<
" " << m20
1852 PrintKTVector(&theFinalState,
" mass problem");
1858 G4double pInCM = std::sqrt((s0-(m10+m20)*(m10+m20))*(s0-(m10-m20)*(m10-m20))/(4.*s0));
1859 #ifdef debug_BIC_CorrectFinalPandE
1860 G4cout <<
" CorrectFinalPandE pInCM new, CURRENT, ratio : " << pInCM
1861 <<
" " << (pFinals).vect().mag()<<
" " << pInCM/(pFinals).vect().mag() <<
G4endl;
1863 if ( pFinals.vect().mag() > pInCM )
1870 for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1873 G4ThreeVector p3(factor*(toCMS*(*i)->Get4Momentum()).vect());
1877 #ifdef debug_BIC_CorrectFinalPandE
1880 (*i)->Set4Momentum(
p);
1882 #ifdef debug_BIC_CorrectFinalPandE
1883 G4cout <<
"CorrectFinalPandE nucleus corrected mass : " << GetFinal4Momentum() <<
" "
1884 <<GetFinal4Momentum().
mag() << G4endl
1885 <<
" CMS pFinals , mag, 3.mag : " << qFinals <<
" " << qFinals.mag() <<
" " << qFinals.vect().mag()<<
G4endl;
1886 G4cerr <<
" -CorrectFinalPandE 5 " << factor <<
G4endl;
1889 #ifdef debug_BIC_CorrectFinalPandE
1890 else {
G4cerr <<
" -CorrectFinalPandE 6 - no correction done" <<
G4endl; }
1896 void G4BinaryCascade::UpdateTracksAndCollisions(
1902 std::vector<G4KineticTrack *>::iterator iter1, iter2;
1907 if(!oldSecondaries->empty())
1909 for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end();
1912 iter2 = std::find(theSecondaryList.begin(), theSecondaryList.end(),
1914 if ( iter2 != theSecondaryList.end() ) theSecondaryList.erase(iter2);
1924 if(oldTarget->size()!=0)
1928 for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1)
1930 iter2 = std::find(theTargetList.begin(), theTargetList.end(),
1932 theTargetList.erase(iter2);
1940 if(!newSecondaries->empty())
1943 for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end();
1946 theSecondaryList.push_back(*iter1);
1949 PrintKTVector(*iter1,
"undefined in FindCollisions");
1955 FindCollisions(newSecondaries);
1970 ktv(out), wanted_state(astate)
1974 if ( (kt)->GetState() == wanted_state ) ktv->push_back(kt);
1985 #ifdef debug_BIC_DoTimeStep
1987 debug.push_back(
"======> DoTimeStep 1");
debug.dump();
1988 G4cerr <<
"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep
1989 <<
" , time="<<theCurrentTime <<
G4endl;
1990 PrintKTVector(&theSecondaryList, std::string(
"DoTimeStep - theSecondaryList"));
1995 std::vector<G4KineticTrack *>::iterator iter;
1998 std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2003 std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2008 #ifdef debug_BIC_DoTimeStep
2014 thePropagator->
Transport(theSecondaryList, dummy, theTimeStep);
2021 #ifdef debug_BIC_DoTimeStep
2022 G4cout <<
"DoTimeStep : theMomentumTransfer = " << theMomentumTransfer <<
G4endl;
2023 PrintKTVector(&theSecondaryList, std::string(
"DoTimeStep - secondaries aft trsprt"));
2029 std::for_each( kt_outside->begin(),kt_outside->end(),
2036 std::for_each( kt_inside->begin(),kt_inside->end(),
2046 kt_gone_in->clear();
2047 std::for_each( kt_outside->begin(),kt_outside->end(),
2050 kt_gone_out->clear();
2051 std::for_each( kt_inside->begin(),kt_inside->end(),
2054 #ifdef debug_BIC_DoTimeStep
2055 PrintKTVector(fail,std::string(
" Failed to go in/out -> miss_nucleus/captured"));
2056 PrintKTVector(kt_gone_in, std::string(
"recreated kt_gone_in"));
2057 PrintKTVector(kt_gone_out, std::string(
"recreated kt_gone_out"));
2063 std::for_each( kt_outside->begin(),kt_outside->end(),
2066 std::for_each( kt_outside->begin(),kt_outside->end(),
2069 #ifdef debug_BIC_DoTimeStep
2070 PrintKTVector(kt_gone_out, std::string(
"append gone_outs to final state.. theFinalState"));
2073 theFinalState.insert(theFinalState.end(),
2074 kt_gone_out->begin(),kt_gone_out->end());
2078 std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2084 if ( theCollisionMgr->
Entries()> 0 )
2086 if (kt_gone_out->size() )
2089 iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary);
2090 if ( iter != kt_gone_out->end() )
2093 #ifdef debug_BIC_DoTimeStep
2094 G4cout <<
" DoTimeStep - WARNING: deleting current collision!" <<
G4endl;
2098 if ( kt_captured->size() )
2101 iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary);
2102 if ( iter != kt_captured->end() )
2105 #ifdef debug_BIC_DoTimeStep
2106 G4cout <<
" DoTimeStep - WARNING: deleting current collision!" <<
G4endl;
2113 UpdateTracksAndCollisions(kt_gone_out,0 ,0);
2116 if ( kt_captured->size() )
2118 theCapturedList.insert(theCapturedList.end(),
2119 kt_captured->begin(),kt_captured->end());
2123 std::vector<G4KineticTrack *>::iterator i_captured;
2124 for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++)
2126 (*i_captured)->Hit();
2129 UpdateTracksAndCollisions(kt_captured, NULL, NULL);
2132 #ifdef debug_G4BinaryCascade
2135 std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2137 if ( currentZ != (GetTotalCharge(theTargetList)
2138 + GetTotalCharge(theCapturedList)
2139 + GetTotalCharge(*kt_inside)) )
2141 G4cout <<
" error-DoTimeStep aft, A, Z: " << currentA <<
" " << currentZ
2142 <<
" sum(tgt,capt,active) "
2143 << GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList) + GetTotalCharge(*kt_inside)
2144 <<
" targets: " << GetTotalCharge(theTargetList)
2145 <<
" captured: " << GetTotalCharge(theCapturedList)
2146 <<
" active: " << GetTotalCharge(*kt_inside)
2158 theCurrentTime += theTimeStep;
2172 std::vector<G4KineticTrack *>::iterator iter;
2177 G4int secondaries_in(0);
2178 G4int secondaryBarions_in(0);
2179 G4int secondaryCharge_in(0);
2182 for ( iter =in->begin(); iter != in->end(); ++iter)
2185 secondaryCharge_in +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
2186 if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 )
2188 secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber();
2192 secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass();
2198 G4double mass_initial= GetIonMass(currentZ,currentA);
2200 currentZ += secondaryCharge_in;
2201 currentA += secondaryBarions_in;
2206 G4double mass_final= GetIonMass(currentZ,currentA);
2208 G4double correction= secondaryMass_in + mass_initial - mass_final;
2209 if (secondaries_in>1)
2210 {correction /= secondaries_in;}
2212 #ifdef debug_BIC_CorrectBarionsOnBoundary
2213 G4cout <<
"CorrectBarionsOnBoundary,currentZ,currentA,"
2214 <<
"secondaryCharge_in,secondaryBarions_in,"
2215 <<
"energy correction,m_secondry,m_nucl_init,m_nucl_final "
2216 << currentZ <<
" "<< currentA <<
" "
2217 << secondaryCharge_in<<
" "<<secondaryBarions_in<<
" "
2218 << correction <<
" "
2219 << secondaryMass_in <<
" "
2220 << mass_initial <<
" "
2221 << mass_final <<
" "
2223 PrintKTVector(in,std::string(
"in be4 correction"));
2226 for ( iter = in->begin(); iter != in->end(); ++iter)
2228 if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
2230 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
2237 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + barrier);
2239 kt_fail->push_back(*iter);
2240 currentZ -=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
2241 currentA -= (*iter)->GetDefinition()->GetBaryonNumber();
2246 #ifdef debug_BIC_CorrectBarionsOnBoundary
2247 G4cout <<
" CorrectBarionsOnBoundary, aft, A, Z, sec-Z,A,m,m_in_nucleus "
2248 << currentA <<
" " << currentZ <<
" "
2249 << secondaryCharge_in <<
" " << secondaryBarions_in <<
" "
2250 << secondaryMass_in <<
" "
2252 PrintKTVector(in,std::string(
"in AFT correction"));
2259 G4int secondaries_out(0);
2260 G4int secondaryBarions_out(0);
2261 G4int secondaryCharge_out(0);
2264 for ( iter =out->begin(); iter != out->end(); ++iter)
2267 secondaryCharge_out +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
2268 if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 )
2270 secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber();
2274 secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass();
2281 G4double mass_initial= GetIonMass(currentZ,currentA);
2282 currentA -=secondaryBarions_out;
2283 currentZ -=secondaryCharge_out;
2292 G4cerr <<
"G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " <<
2293 secondaryBarions_out <<
" " << secondaryCharge_out <<
G4endl;
2294 PrintKTVector(&theTargetList,
"CorrectBarionsOnBoundary Target");
2295 PrintKTVector(&theCapturedList,
"CorrectBarionsOnBoundary Captured");
2296 PrintKTVector(&theSecondaryList,
"CorrectBarionsOnBoundary Secondaries");
2297 G4cerr <<
"G4BinaryCascade - currentA, currentZ " << currentA <<
" " << currentZ <<
G4endl;
2298 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error");
2300 G4double mass_final=GetIonMass(currentZ,currentA);
2301 G4double correction= mass_initial - mass_final - secondaryMass_out;
2303 if (secondaries_out>1) correction /= secondaries_out;
2304 #ifdef debug_BIC_CorrectBarionsOnBoundary
2305 G4cout <<
"DoTimeStep,currentZ,currentA,"
2306 <<
"secondaries_out,"
2307 <<
"secondaryCharge_out,secondaryBarions_out,"
2308 <<
"energy correction,m_secondry,m_nucl_init,m_nucl_final "
2309 <<
" "<< currentZ <<
" "<< currentA <<
" "
2310 << secondaries_out <<
" "
2311 << secondaryCharge_out<<
" "<<secondaryBarions_out<<
" "
2312 << correction <<
" "
2313 << secondaryMass_out <<
" "
2314 << mass_initial <<
" "
2315 << mass_final <<
" "
2317 PrintKTVector(out,std::string(
"out be4 correction"));
2320 for ( iter = out->begin(); iter != out->end(); ++iter)
2322 if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
2324 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
2335 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() - barrier);
2337 kt_fail->push_back(*iter);
2338 currentZ +=
G4lrint((*iter)->GetDefinition()->GetPDGCharge()/
eplus);
2339 currentA += (*iter)->GetDefinition()->GetBaryonNumber();
2341 #ifdef debug_BIC_CorrectBarionsOnBoundary
2344 G4cout <<
"Not correcting outgoing " << *iter <<
" "
2345 << (*iter)->GetDefinition()->GetPDGEncoding() <<
" "
2346 << (*iter)->GetDefinition()->GetParticleName() <<
G4endl;
2347 PrintKTVector(out,std::string(
"outgoing, one not corrected"));
2353 #ifdef debug_BIC_CorrectBarionsOnBoundary
2354 PrintKTVector(out,std::string(
"out AFTER correction"));
2355 G4cout <<
" DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta "
2356 << currentA <<
" "<< currentZ <<
" "
2357 << secondaryCharge_out <<
" "<< secondaryBarions_out <<
" "<<
2358 secondaryMass_out <<
" "
2359 << massInNucleus <<
" "
2372 G4Fragment * G4BinaryCascade::FindFragments()
2376 #ifdef debug_BIC_FindFragments
2377 G4cout <<
"target, captured, secondary: "
2378 << theTargetList.size() <<
" "
2379 << theCapturedList.size()<<
" "
2380 << theSecondaryList.size()
2384 G4int a = theTargetList.size()+theCapturedList.size();
2386 G4KineticTrackVector::iterator i;
2387 for(i = theTargetList.begin(); i != theTargetList.end(); ++i)
2389 if(
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus) == 1 )
2395 G4int zCaptured = 0;
2397 for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2399 CapturedMomentum += (*i)->Get4Momentum();
2400 if(
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus) == 1 )
2406 G4int z = zTarget+zCaptured;
2408 #ifdef debug_G4BinaryCascade
2409 if ( z != (GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList)) )
2411 G4cout <<
" FindFragment Counting error z a " << z <<
" " <<a <<
" "
2412 << GetTotalCharge(theTargetList) <<
" " << GetTotalCharge(theCapturedList)<<
2414 PrintKTVector(&theTargetList, std::string(
"theTargetList"));
2415 PrintKTVector(&theCapturedList, std::string(
"theCapturedList"));
2429 if ( z < 1 )
return 0;
2432 G4int excitons = theCapturedList.size();
2433 #ifdef debug_BIC_FindFragments
2434 G4cout <<
"Fragment: a= " << a <<
" z= " << z <<
" particles= " << excitons
2435 <<
" Charged= " << zCaptured <<
" holes= " << holes
2436 <<
" excitE= " <<GetExcitationEnergy()
2437 <<
" Final4Momentum= " << GetFinalNucleusMomentum() <<
" capturMomentum= " << CapturedMomentum
2458 G4LorentzVector final4Momentum = theInitial4Mom + theProjectile4Momentum;
2460 for(G4KineticTrackVector::iterator i = theFinalState.begin(); i != theFinalState.end(); ++i)
2462 final4Momentum -= (*i)->Get4Momentum();
2463 finals += (*i)->Get4Momentum();
2466 if((final4Momentum.
vect()/final4Momentum.
e()).mag()>1.0 && currentA > 0)
2468 #ifdef debug_BIC_Final4Momentum
2470 G4cerr <<
"G4BinaryCascade::GetFinal4Momentum - Fatal"<<
G4endl;
2471 G4KineticTrackVector::iterator i;
2472 G4cerr <<
"Total initial 4-momentum " << theProjectile4Momentum <<
G4endl;
2473 G4cerr <<
" GetFinal4Momentum: Initial nucleus "<<theInitial4Mom<<
G4endl;
2474 for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
2476 G4cerr <<
" Final state: "<<(*i)->Get4Momentum()<<(*i)->GetDefinition()->GetParticleName()<<
G4endl;
2479 G4cerr<<
" Final4Momentum = "<<final4Momentum <<
" "<<final4Momentum.
m()<<
G4endl;
2480 G4cerr <<
" current A, Z = "<< currentA<<
", "<<currentZ<<
G4endl;
2486 return final4Momentum;
2497 G4KineticTrackVector::iterator i;
2499 for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2501 CapturedMomentum += (*i)->Get4Momentum();
2507 if ( NucleusMomentum.
e() > 0 )
2511 G4ThreeVector boost= (NucleusMomentum.
vect() -CapturedMomentum.vect())/NucleusMomentum.
e();
2512 if(boost.
mag2()>1.0)
2514 # ifdef debug_BIC_FinalNucleusMomentum
2515 G4cerr <<
"G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<
G4endl;
2517 G4cerr <<
"it 01"<<NucleusMomentum<<
" "<<CapturedMomentum<<
" "<<
G4endl;
2524 precompoundLorentzboost.
set( boost );
2525 #ifdef debug_debug_BIC_FinalNucleusMomentum
2526 G4cout <<
"GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<
" "<<CapturedMomentum<<
" "<<
G4endl;
2528 NucleusMomentum *= nucleusBoost;
2529 #ifdef debug_BIC_FinalNucleusMomentum
2530 G4cout <<
"GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<
G4endl;
2533 return NucleusMomentum;
2550 std::vector<G4KineticTrack *>::iterator iter, jter;
2555 while(!done && tryCount++ <200)
2562 secs = theH1Scatterer->
Scatter(*(*secondaries).front(), aTarget);
2563 #ifdef debug_H1_BinaryCascade
2564 PrintKTVector(secs,
" From Scatter");
2566 for(
size_t ss=0; secs && ss<secs->size(); ss++)
2569 if((*secs)[ss]->GetDefinition()->IsShortLived()) done =
true;
2573 ClearAndDestroy(&theFinalState);
2574 ClearAndDestroy(secondaries);
2577 for(
size_t current=0; secs && current<secs->size(); current++)
2579 if((*secs)[current]->GetDefinition()->IsShortLived())
2583 for(jter=dec->begin(); jter != dec->end(); jter++)
2586 secs->push_back(*jter);
2589 delete (*secs)[current];
2596 theFinalState.push_back((*secs)[current]);
2601 #ifdef debug_H1_BinaryCascade
2602 PrintKTVector(&theFinalState,
" FinalState");
2604 for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
2610 products->push_back(aNew);
2611 #ifdef debug_H1_BinaryCascade
2616 G4cout <<
"final shortlived : ";
2619 G4cout <<
"final un stable : ";
2626 theFinalState.clear();
2651 }
while (
sqr(x1) +
sqr(x2) > 1.);
2677 std::vector<G4KineticTrack *>::iterator i;
2678 for(i = ktv->begin(); i != ktv->end(); ++i)
2687 std::vector<G4ReactionProduct *>::iterator i;
2688 for(i = rpv->begin(); i != rpv->end(); ++i)
2697 if (comment.size() > 0 )
G4cout <<
"G4BinaryCascade::PrintKTVector() " << comment << G4endl;
2699 G4cout <<
" vector: " << ktv <<
", number of tracks: " << ktv->size()
2701 std::vector<G4KineticTrack *>::iterator i;
2704 for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count)
2707 G4cout <<
" track n. " << count;
2711 G4cout <<
"G4BinaryCascade::PrintKTVector():No KineticTrackVector given " <<
G4endl;
2715 void G4BinaryCascade::PrintKTVector(
G4KineticTrack * kt, std::string comment)
2718 if (comment.size() > 0 )
G4cout <<
"G4BinaryCascade::PrintKTVector() "<< comment << G4endl;
2731 G4cout <<
"G4BinaryCascade::PrintKTVector(): No Kinetictrack given" <<
G4endl;
2741 if ( Z > 0 && A >= Z )
2745 }
else if ( A > 0 && Z>0 )
2750 }
else if ( A >= 0 && Z<=0 )
2755 }
else if ( A == 0 && std::abs(Z)<2 )
2762 G4cerr <<
"G4BinaryCascade::GetIonMass() - invalid (A,Z) = ("
2763 << A <<
"," << Z <<
")" <<
G4endl;
2764 throw G4HadronicException(__FILE__, __LINE__,
"G4BinaryCascade::GetIonMass() - giving up");
2773 std::vector<G4KineticTrack *>::iterator iter;
2774 decayKTV.
Decay(&theFinalState);
2776 for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
2779 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2781 Esecondaries +=(*iter)->Get4Momentum().e();
2784 products->push_back(aNew);
2788 for(iter = theCapturedList.begin(); iter != theCapturedList.end(); ++iter)
2791 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2793 Esecondaries +=(*iter)->Get4Momentum().e();
2796 products->push_back(aNew);
2800 while(theCollisionMgr->
Entries() > 0)
2807 if ( lates->size() == 1 ) {
2815 products->push_back(aNew);
2825 decayKTV.
Decay(&theSecondaryList);
2827 for(iter = theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
2830 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2832 Esecondaries +=(*iter)->Get4Momentum().e();
2835 products->push_back(aNew);
2839 for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter)
2841 SumMassNucleons += (*iter)->GetDefinition()->GetPDGMass();
2844 G4double Ekinetic=theProjectile4Momentum.e() + initial_nuclear_mass - Esecondaries - SumMassNucleons;
2847 if (theTargetList.size() ) Ekinetic /= theTargetList.size();
2854 for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter)
2858 G4double p=std::sqrt(
sqr(Ekinetic) + 2.*Ekinetic*mass);
2859 aNew->
SetMomentum(p*(*iter)->Get4Momentum().vect().unit());
2863 products->push_back(aNew);
2871 std::vector<G4KineticTrack *>::iterator iter;
2872 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
2875 aNew->
SetMomentum((*iter)->Get4Momentum().vect());
2879 products->push_back(aNew);
2882 if (currentA == 1 && currentZ == 0) {
2884 }
else if (currentA == 1 && currentZ == 1) {
2886 }
else if (currentA == 2 && currentZ == 1) {
2888 }
else if (currentA == 3 && currentZ == 1) {
2890 }
else if (currentA == 3 && currentZ == 2) {
2892 }
else if (currentA == 4 && currentZ == 2) {
2898 if (fragment != 0) {
2904 products->push_back(theNew);
2909 void G4BinaryCascade::PrintWelcomeMessage()
2911 G4cout <<
"Thank you for using G4BinaryCascade. "<<
G4endl;
2921 for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
2923 G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
2924 if (std::abs(PDGcode)==211 || PDGcode==111 ) havePion=
true;
2927 if ( !products || havePion)
2930 <<
", with NO products! " <<
G4endl;
2931 G4cout << G4endl<<
"Initial condition are these:"<<
G4endl;
2947 G4bool G4BinaryCascade::CheckChargeAndBaryonNumber(
G4String where)
2949 static G4int lastdA(0), lastdZ(0);
2956 std::vector<G4KineticTrack *>::iterator i;
2957 G4int CapturedA(0), CapturedZ(0);
2958 G4int secsA(0), secsZ(0);
2959 for ( i=theCapturedList.begin(); i!=theCapturedList.end(); ++i) {
2960 CapturedA += (*i)->GetDefinition()->GetBaryonNumber();
2961 CapturedZ +=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
2964 for ( i=theSecondaryList.begin(); i!=theSecondaryList.end(); ++i) {
2966 secsA += (*i)->GetDefinition()->GetBaryonNumber();
2967 secsZ +=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
2971 for ( i=theFinalState.begin(); i!=theFinalState.end(); ++i) {
2972 fStateA += (*i)->GetDefinition()->GetBaryonNumber();
2973 fStateZ +=
G4lrint((*i)->GetDefinition()->GetPDGCharge()/
eplus);
2976 G4int deltaA= iStateA - secsA - fStateA -currentA - lateA;
2977 G4int deltaZ= iStateZ - secsZ - fStateZ -currentZ - lateZ;
2979 if (deltaA != 0 || deltaZ!=0 ) {
2980 if (deltaA != lastdA || deltaZ != lastdZ ) {
2981 G4cout <<
"baryon/charge imbalance - " << where << G4endl
2982 <<
"deltaA " <<deltaA<<
", iStateA "<<iStateA<<
", CapturedA "<<CapturedA <<
", secsA "<<secsA
2983 <<
", fStateA "<<fStateA <<
", currentA "<<currentA <<
", lateA "<<lateA << G4endl
2984 <<
"deltaZ "<<deltaZ<<
", iStateZ "<<iStateZ<<
", CapturedZ "<<CapturedZ <<
", secsZ "<<secsZ
2985 <<
", fStateZ "<<fStateZ <<
", currentZ "<<currentZ <<
", lateZ "<<lateZ << G4endl<<
G4endl;
2989 }
else { lastdA=lastdZ=0;}
2998 PrintKTVector(collision->
GetPrimary(),std::string(
" Primary particle"));
3000 PrintKTVector(products,std::string(
" Scatterer products"));
3017 <<
" " << initial << G4endl;;
3020 for (
unsigned int it=0; it < ktv.size(); it++)
3033 <<
" " << initial <<
" Excit " << thisExcitation << G4endl;;
3038 G4int product_barions(0);
3041 for (
unsigned int it=0; it < products->size(); it++)
3053 <<
" " <<
final << G4endl;;
3058 G4int finalA = currentA;
3059 G4int finalZ = currentZ;
3062 finalA -= product_barions;
3063 finalZ -= GetTotalCharge(*products);
3068 G4cout <<
" current/final a,z " << currentA <<
" " << currentZ <<
" "<< finalA<<
" "<< finalZ
3069 <<
" delta-mass " << delta<<
G4endl;
3072 G4cout <<
" initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
3073 <<
" " <<
final <<
" "
3075 << currentInitialEnergy -
final - mass_out
3077 currentInitialEnergy-=
final;
3086 G4ReactionProductVector::iterator iter;
3094 for(iter = products->begin(); iter != products->end(); ++iter)
3104 Efinal += (*iter)->GetTotalEnergy();
3105 pFinal += (*iter)->GetMomentum();
3109 G4cout <<
"BIC E/p delta " <<
G4double GetWeightChange() const
G4bool IsParticipant() const
static G4Triton * TritonDefinition()
Hep3Vector boostVector() const
void Update4Momentum(G4double aEnergy)
static G4He3 * He3Definition()
CascadeState GetState() const
const G4LorentzVector & GetMomentum() const
virtual G4int GetCharge()=0
CLHEP::Hep3Vector G4ThreeVector
const G4HadProjectile * GetPrimaryProjectile() const
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
const G4ThreeVector & GetPosition() const
virtual G4bool StartLoop()=0
G4VPreCompoundModel * GetDeExcitation() const
virtual void Transport(G4KineticTrackVector &theActive, const G4KineticTrackVector &theSpectators, G4double theTimeStep)=0
void RemoveCollision(G4CollisionInitialState *collision)
virtual const G4VNuclearDensity * GetNuclearDensity() const =0
G4KineticTrack * GetPrimary(void)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4bool GetPDGStable() const
G4CollisionInitialState * GetNextCollision()
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
virtual G4int GetMassNumber()=0
static G4Proton * ProtonDefinition()
virtual G4ThreeVector GetMomentumTransfer() const =0
virtual const std::vector< G4CollisionInitialState * > & GetCollisions(G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &, G4double theCurrentTime)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
G4int GetPDGEncoding() const
G4double GetActualMass() const
virtual G4KineticTrackVector * Scatter(const G4KineticTrack &trk1, const G4KineticTrack &trk2)
virtual const G4ThreeVector & GetPosition() const
static void ConstructParticle()
#define _CheckChargeAndBaryonNumber_(val)
virtual void PropagateModelDescription(std::ostream &) const
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState) const
virtual G4double CoulombBarrier()=0
const G4String & GetParticleName() const
std::vector< G4LorentzVector * > * Decay(const G4double, const std::vector< G4double > &) const
void RemoveTracksCollisions(G4KineticTrackVector *ktv)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
G4ParticleDefinition * GetDefinition() const
void SetNewlyAdded(const G4bool f)
G4double GetBarrier(G4int encoding)
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
virtual G4double GetOuterRadius()=0
G4KineticTrackVector * GetFinalState()
void SetMinEnergy(G4double anEnergy)
double precision function energy(A, Z)
G4ParticleDefinition * GetDefinition() const
virtual ~G4BinaryCascade()
G4ExcitationHandler * GetExcitationHandler() const
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
CascadeState SetState(const CascadeState new_state)
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
const G4ParticleDefinition * GetDefinition() const
virtual void Init(G4int theA, G4int theZ)=0
G4bool nucleon(G4int ityp)
static G4PionMinus * PionMinusDefinition()
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
const G4LorentzVector & GetMomentum() const
void SetNucleon(G4Nucleon *aN)
G4double GetKineticEnergy() const
void SetTotalEnergy(const G4double en)
G4double GetFermiMomentum(G4double density)
static G4PionPlus * PionPlusDefinition()
static G4Proton * Proton()
HepLorentzRotation & set(double bx, double by, double bz)
virtual void Init(G4V3DNucleus *theNucleus)=0
G4int GetTargetBaryonNumber()
static G4Neutron * Neutron()
void Set4Momentum(const G4LorentzVector &a4Momentum)
void SetNumberOfParticles(G4int value)
G4KineticTrackVector & GetTargetCollection(void)
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *, G4V3DNucleus *)
const G4LorentzVector & Get4Momentum() const
virtual const std::vector< G4CollisionInitialState * > & GetCollisions(G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &, G4double theCurrentTime)
G4BinaryCascade(G4VPreCompoundModel *ptr=0)
void AddCollision(G4double time, G4KineticTrack *proj, G4KineticTrack *target=NULL)
G4double GetKineticEnergy() const
G4HadronicInteraction * FindModel(const G4String &name)
void operator()(G4KineticTrack *&kt) const
G4bool IsShortLived() const
void SetEnergyChange(G4double anEnergy)
virtual void ModelDescription(std::ostream &) const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)=0
void Init(G4int anA, G4int aZ)
void Decay(G4KineticTrackVector *tracks) const
G4double GetPDGMass() const
G4double GetDensity(const G4ThreeVector &aPosition) const
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4LorentzVector & GetTrackingMomentum() const
static G4HadronicInteractionRegistry * Instance()
G4BCAction * GetGenerator()
Hep3Vector orthogonal() const
G4VPreCompoundModel * theDeExcitation
G4double GetField(G4int encoding, G4ThreeVector pos)
void SetMaxEnergy(const G4double anEnergy)
void SetDeExcitation(G4VPreCompoundModel *ptr)
G4HadFinalState theParticleChange
virtual G4double GetMass()=0
virtual G4ParticleDefinition * GetDefinition() const
Hep3Vector cross(const Hep3Vector &) const
G4V3DNucleus * the3DNucleus
virtual G4Nucleon * GetNextNucleon()=0
HepLorentzRotation inverse() const
void SetNumberOfCharged(G4int value)
G4double GetCollisionTime(void)
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double absoluteLevel)
static G4Deuteron * DeuteronDefinition()
static G4Alpha * AlphaDefinition()
static G4Neutron * NeutronDefinition()
void SetMomentumChange(const G4ThreeVector &aV)
SelectFromKTV(G4KineticTrackVector *out, G4KineticTrack::CascadeState astate)
void AddSecondary(G4DynamicParticle *aP)
G4GLOB_DLL std::ostream G4cerr
G4int GetBaryonNumber() const
CLHEP::HepLorentzVector G4LorentzVector