00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <algorithm>
00027 #include <vector>
00028 #include <cmath>
00029 #include <numeric>
00030
00031 #include "G4BinaryLightIonReaction.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "G4LorentzVector.hh"
00035 #include "G4LorentzRotation.hh"
00036 #include "G4ReactionProductVector.hh"
00037 #include "G4ping.hh"
00038 #include "G4Delete.hh"
00039 #include "G4Neutron.hh"
00040 #include "G4VNuclearDensity.hh"
00041 #include "G4FermiMomentum.hh"
00042 #include "G4HadTmpUtil.hh"
00043 #include "G4PreCompoundModel.hh"
00044 #include "G4HadronicInteractionRegistry.hh"
00045
00046
00047
00048
00049 G4BinaryLightIonReaction::G4BinaryLightIonReaction(G4VPreCompoundModel* ptr)
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 }
00067
00068 G4BinaryLightIonReaction::~G4BinaryLightIonReaction()
00069 {}
00070
00071 void G4BinaryLightIonReaction::ModelDescription(std::ostream& outFile) const
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 }
00079
00080
00081 struct ReactionProduct4Mom
00082 {
00083 G4LorentzVector operator()(G4LorentzVector a,G4ReactionProduct* b) {return a + G4LorentzVector(b->GetMomentum(), b->GetTotalEnergy() );}
00084 };
00085
00086 G4HadFinalState *G4BinaryLightIonReaction::
00087 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus )
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
00100 G4LorentzRotation toBreit(mom.boostVector());
00101
00102 G4bool swapped=SetLighterAsProjectile(mom, toBreit);
00103
00104 G4ReactionProductVector * result = 0;
00105 G4ReactionProductVector * cascaders=0;
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 if( (mom.t()-mom.mag())/pA < 50*MeV )
00117 {
00118
00119
00120 cascaders=FuseNucleiAndPrompound(mom);
00121 }
00122 else
00123 {
00124 result=Interact(mom,toBreit);
00125
00126 if(! result )
00127 {
00128 {
00129
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
00153 G4double theStatisticalExEnergy = GetProjectileExcitation();
00154
00155
00156 pInitialState = mom;
00157
00158 pInitialState.setT(pInitialState.getT() +
00159 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA));
00160
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
00172 std::vector<G4ReactionProduct *>::iterator iter;
00173
00174
00175
00176
00177
00178
00179
00180 delete result;
00181 result=0;
00182 G4LorentzVector momentum(pInitialState-pFinalState);
00183 G4int loopcount(0);
00184
00185 while (std::abs(momentum.e()-pspectators.e()) > 10*MeV)
00186 {
00187 G4LorentzVector pCorrect(pInitialState-pspectators);
00188
00189
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
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
00256 G4LorentzRotation toZ;
00257 toZ.rotateZ(-1*mom.phi());
00258 toZ.rotateY(-1*mom.theta());
00259 G4LorentzRotation toLab(toZ.inverse());
00260
00261
00262
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
00284 theResult.AddSecondary(aNew);
00285 Etot += tmp.e();
00286
00287
00288
00289
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 }
00306
00307
00308
00309
00310 G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector(
00311 G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom)
00312
00313 {
00314 const int nAttemptScale = 2500;
00315 const double ErrLimit = 1.E-6;
00316 if (Output->empty())
00317 return TRUE;
00318 G4LorentzVector SumMom(0,0,0,0);
00319 G4double SumMass = 0;
00320 G4double TotalCollisionMass = TotalCollisionMom.m();
00321 size_t i = 0;
00322
00323 for(i = 0; i < Output->size(); i++)
00324 {
00325 SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
00326 SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
00327 }
00328
00329
00330 if (SumMass > TotalCollisionMass) return FALSE;
00331 SumMass = SumMom.m2();
00332 if (SumMass < 0) return FALSE;
00333 SumMass = std::sqrt(SumMass);
00334
00335
00336 G4ThreeVector Beta = -SumMom.boostVector();
00337
00338
00339 for(i = 0; i < Output->size(); i++)
00340 {
00341 G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
00342 mom *= Beta;
00343 (*Output)[i]->SetMomentum(mom.vect());
00344 (*Output)[i]->SetTotalEnergy(mom.e());
00345 }
00346
00347
00348
00349 G4double Scale = 0,OldScale=0;
00350 G4double factor = 1.;
00351 G4int cAttempt = 0;
00352 G4double Sum = 0;
00353 G4bool success = false;
00354 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
00355 {
00356 Sum = 0;
00357 for(i = 0; i < Output->size(); i++)
00358 {
00359 G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
00360 HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
00361 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
00362 HadronMom.setE(E);
00363 (*Output)[i]->SetMomentum(HadronMom.vect());
00364 (*Output)[i]->SetTotalEnergy(HadronMom.e());
00365 Sum += E;
00366 }
00367 OldScale=Scale;
00368 Scale = TotalCollisionMass/Sum - 1;
00369
00370 if (std::abs(Scale) <= ErrLimit
00371 || OldScale == Scale)
00372 {
00373 if (debug_G4BinaryLightIonReactionResults) G4cout << "E/p corrector: " << cAttempt << G4endl;
00374 success = true;
00375 break;
00376 }
00377 if ( cAttempt > 10 )
00378 {
00379
00380 factor=std::max(1.,std::log(std::abs(OldScale/(OldScale-Scale))));
00381
00382 }
00383 }
00384
00385 if( (!success) && debug_G4BinaryLightIonReactionResults)
00386 {
00387 G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl;
00388 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
00389 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
00390 }
00391
00392
00393 Beta = TotalCollisionMom.boostVector();
00394
00395 for(i = 0; i < Output->size(); i++)
00396 {
00397 G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
00398 mom *= Beta;
00399 (*Output)[i]->SetMomentum(mom.vect());
00400 (*Output)[i]->SetTotalEnergy(mom.e());
00401 }
00402 return TRUE;
00403 }
00404 G4bool G4BinaryLightIonReaction::SetLighterAsProjectile(G4LorentzVector & mom,const G4LorentzRotation & toBreit)
00405 {
00406 G4bool swapped = false;
00407 if(tA<pA)
00408 {
00409 swapped = true;
00410 G4int tmp(0);
00411 tmp = tA; tA=pA; pA=tmp;
00412 tmp = tZ; tZ=pZ; pZ=tmp;
00413 G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
00414 G4LorentzVector it(m1, G4ThreeVector(0,0,0));
00415 mom = toBreit*it;
00416 }
00417 return swapped;
00418 }
00419 G4ReactionProductVector * G4BinaryLightIonReaction::FuseNucleiAndPrompound(const G4LorentzVector & mom)
00420 {
00421 G4Fragment aPreFrag;
00422 aPreFrag.SetA(pA+tA);
00423 aPreFrag.SetZ(pZ+tZ);
00424 aPreFrag.SetNumberOfParticles(pA);
00425 aPreFrag.SetNumberOfCharged(pZ);
00426 aPreFrag.SetNumberOfHoles(0);
00427 G4ThreeVector plop(0.,0., mom.vect().mag());
00428 G4double m_nucl=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA);
00429 G4LorentzVector aL(mom.t()+m_nucl, plop);
00430 aPreFrag.SetMomentum(aL);
00431 G4ParticleDefinition * preFragDef;
00432 preFragDef = G4ParticleTable::GetParticleTable()
00433 ->FindIon(pZ+tZ,pA+tA,0,pZ+tZ);
00434 aPreFrag.SetParticleDefinition(preFragDef);
00435
00436
00437
00438 G4ReactionProductVector * cascaders = theProjectileFragmentation->DeExcite(aPreFrag);
00439
00440 for(size_t count = 0; count<cascaders->size(); count++)
00441 {
00442 cascaders->operator[](count)->SetNewlyAdded(true);
00443
00444 }
00445
00446 return cascaders;
00447 }
00448 G4ReactionProductVector * G4BinaryLightIonReaction::Interact(G4LorentzVector & mom, const G4LorentzRotation & toBreit)
00449 {
00450 G4ReactionProductVector * result = 0;
00451 G4double projectileMass(0);
00452 G4LorentzVector it;
00453
00454 G4int tryCount(0);
00455 do
00456 {
00457 ++tryCount;
00458 projectile3dNucleus = new G4Fancy3DNucleus;
00459 projectile3dNucleus->Init(pA, pZ);
00460 projectile3dNucleus->CenterNucleons();
00461 projectileMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
00462 projectile3dNucleus->GetCharge(),projectile3dNucleus->GetMassNumber());
00463 it=toBreit * G4LorentzVector(projectileMass,G4ThreeVector(0,0,0));
00464
00465 target3dNucleus = new G4Fancy3DNucleus;
00466 target3dNucleus->Init(tA, tZ);
00467 G4double impactMax = target3dNucleus->GetOuterRadius()+projectile3dNucleus->GetOuterRadius();
00468
00469 G4double aX=(2.*G4UniformRand()-1.)*impactMax;
00470 G4double aY=(2.*G4UniformRand()-1.)*impactMax;
00471 G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi);
00472
00473 G4KineticTrackVector * initalState = new G4KineticTrackVector;
00474 projectile3dNucleus->StartLoop();
00475 G4Nucleon * aNuc;
00476 G4LorentzVector tmpV(0,0,0,0);
00477 G4LorentzVector nucleonMom(1./pA*mom);
00478 nucleonMom.setZ(nucleonMom.vect().mag());
00479 nucleonMom.setX(0);
00480 nucleonMom.setY(0);
00481 theFermi.Init(pA,pZ);
00482 while( (aNuc=projectile3dNucleus->GetNextNucleon()) )
00483 {
00484 G4LorentzVector p4 = aNuc->GetMomentum();
00485 tmpV+=p4;
00486 G4ThreeVector nucleonPosition(aNuc->GetPosition());
00487 G4double density=(projectile3dNucleus->GetNuclearDensity())->GetDensity(nucleonPosition);
00488 nucleonPosition += pos;
00489 G4KineticTrack * it1 = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom );
00490 it1->SetState(G4KineticTrack::outside);
00491 G4double pfermi= theFermi.GetFermiMomentum(density);
00492 G4double mass = aNuc->GetDefinition()->GetPDGMass();
00493 G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass;
00494 it1->SetProjectilePotential(-Efermi);
00495 initalState->push_back(it1);
00496 }
00497
00498 result=theModel->Propagate(initalState, target3dNucleus);
00499 if( result && result->size()==0)
00500 {
00501 delete result;
00502 result=0;
00503 }
00504 if ( ! result )
00505 {
00506 delete target3dNucleus;
00507 delete projectile3dNucleus;
00508 }
00509
00510
00511
00512
00513 } while (! result && tryCount< 150);
00514 return result;
00515 }
00516 G4double G4BinaryLightIonReaction::GetProjectileExcitation()
00517 {
00518 spectatorA=spectatorZ=0;
00519
00520 G4Nucleon * aNuc;
00521
00522
00523
00524
00525
00526
00527 G4double theStatisticalExEnergy = 0;
00528 projectile3dNucleus->StartLoop();
00529 while( (aNuc=projectile3dNucleus->GetNextNucleon()) )
00530 {
00531
00532 if(!aNuc->AreYouHit())
00533 {
00534 spectatorA++;
00535 spectatorZ+=G4lrint(aNuc->GetDefinition()->GetPDGCharge()/eplus);
00536 }
00537 else
00538 {
00539 G4ThreeVector aPosition(aNuc->GetPosition());
00540 G4double localDensity = projectile3dNucleus->GetNuclearDensity()->GetDensity(aPosition);
00541 G4double localPfermi = theFermi.GetFermiMomentum(localDensity);
00542 G4double nucMass = aNuc->GetDefinition()->GetPDGMass();
00543 G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
00544 G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag());
00545 theStatisticalExEnergy += deltaE;
00546 }
00547 }
00548 return theStatisticalExEnergy;
00549 }
00550
00551 G4LorentzVector G4BinaryLightIonReaction::SortResult(G4ReactionProductVector * result, G4ReactionProductVector * spectators,G4ReactionProductVector * cascaders)
00552 {
00553 unsigned int i(0);
00554
00555 G4LorentzVector pspectators(0,0,0,0);
00556 pFinalState=G4LorentzVector(0,0,0,0);
00557 for(i=0; i<result->size(); i++)
00558 {
00559 if( (*result)[i]->GetNewlyAdded() )
00560 {
00561 pFinalState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
00562 cascaders->push_back((*result)[i]);
00563 }
00564 else {
00565
00566 pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
00567 spectators->push_back((*result)[i]);
00568
00569
00570 }
00571
00572
00573
00574
00575
00576 }
00577
00578 return pspectators;
00579 }
00580
00581 void G4BinaryLightIonReaction::DeExciteSpectatorNucleus(G4ReactionProductVector * spectators, G4ReactionProductVector * cascaders,
00582 G4double theStatisticalExEnergy, G4LorentzVector & pSpectators)
00583 {
00584
00585 G4ReactionProductVector * proFrag = 0;
00586 G4LorentzVector pFragment(0.,0.,0.,0.);
00587
00588 G4LorentzRotation boost_fragments;
00589
00590
00591
00592 G4LorentzVector pFragments(0,0,0,0);
00593
00594 if(spectatorZ>0 && spectatorA>1)
00595 {
00596
00597 G4Fragment aProRes;
00598 aProRes.SetA(spectatorA);
00599 aProRes.SetZ(spectatorZ);
00600 aProRes.SetNumberOfParticles(0);
00601 aProRes.SetNumberOfCharged(0);
00602 aProRes.SetNumberOfHoles(pA-spectatorA);
00603 G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(spectatorZ,spectatorA);
00604 pFragment=G4LorentzVector(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) );
00605 aProRes.SetMomentum(pFragment);
00606 G4ParticleDefinition * resDef;
00607 resDef = G4ParticleTable::GetParticleTable()->FindIon(spectatorZ,spectatorA,0,spectatorZ);
00608 aProRes.SetParticleDefinition(resDef);
00609
00610 proFrag = theHandler->BreakItUp(aProRes);
00611
00612 boost_fragments = G4LorentzRotation(pSpectators.boostVector());
00613
00614
00615
00616
00617
00618
00619 G4ReactionProductVector::iterator ispectator;
00620 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
00621 {
00622 delete *ispectator;
00623 }
00624 }
00625 else if(spectatorA!=0)
00626 {
00627 G4ReactionProductVector::iterator ispectator;
00628 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
00629 {
00630 (*ispectator)->SetNewlyAdded(true);
00631 cascaders->push_back(*ispectator);
00632 pFragments+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy());
00633
00634
00635
00636
00637 }
00638 }
00639
00640 delete spectators;
00641
00642
00643 G4ReactionProductVector::iterator ii;
00644 if(proFrag)
00645 {
00646 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
00647 {
00648 (*ii)->SetNewlyAdded(true);
00649 G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy());
00650 tmp *= boost_fragments;
00651 (*ii)->SetMomentum(tmp.vect());
00652 (*ii)->SetTotalEnergy(tmp.e());
00653
00654 pFragments += tmp;
00655 }
00656 }
00657
00658
00659
00660
00661
00662 G4LorentzVector pCas=pInitialState - pFragments;
00663
00664
00665 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas);
00666 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
00667 {
00668 G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl;
00669 }
00670
00671
00672 if(proFrag)
00673 {
00674 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
00675 {
00676 cascaders->push_back(*ii);
00677 }
00678 delete proFrag;
00679 }
00680
00681 if ( ! EnergyIsCorrect )
00682 {
00683
00684 if (! EnergyAndMomentumCorrector(cascaders,pInitialState))
00685 {
00686 if(debug_G4BinaryLightIonReactionResults)
00687 G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl;
00688 }
00689 }
00690
00691 }
00692