Geant4  10.00.p02
G4BinaryLightIonReaction.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 #include <algorithm>
27 #include <vector>
28 #include <cmath>
29 #include <numeric>
30 
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4LorentzVector.hh"
35 #include "G4LorentzRotation.hh"
37 #include "G4ping.hh"
38 #include "G4Delete.hh"
39 #include "G4Neutron.hh"
40 #include "G4VNuclearDensity.hh"
41 #include "G4FermiMomentum.hh"
42 #include "G4HadTmpUtil.hh"
43 #include "G4PreCompoundModel.hh"
45 
46 //#define debug_G4BinaryLightIonReaction
47 //#define debug_BLIR_finalstate
48 
50 : G4HadronicInteraction("Binary Light Ion Cascade"),
51  theProjectileFragmentation(ptr),
52  pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
53  projectile3dNucleus(0),target3dNucleus(0)
54 {
55  if(!ptr) {
58  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
59  if(!pre) { pre = new G4PreCompoundModel(); }
61  }
64 
65  debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0;
66 }
67 
69 {}
70 
72 {
73  outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
74  << "using G4BinaryCasacde to model the interaction of a light\n"
75  << "nucleus with a nucleus.\n"
76  << "The lighter of the two nuclei is treated like a set of projectiles\n"
77  << "which are transported simultanously through the heavier nucleus.\n";
78 }
79 
80 //--------------------------------------------------------------------------------
82 {
84 };
85 
87 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus )
88 {
89  if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction starts ######### " << G4endl;
90  G4ping debug("debug_G4BinaryLightIonReaction");
91  pA=aTrack.GetDefinition()->GetBaryonNumber();
93  tA=targetNucleus.GetA_asInt();
94  tZ=targetNucleus.GetZ_asInt();
95 
96  G4LorentzVector mom(aTrack.Get4Momentum());
97  //G4cout << "proj mom : " << mom << G4endl;
98  G4LorentzRotation toBreit(mom.boostVector());
99 
100  G4bool swapped=SetLighterAsProjectile(mom, toBreit);
101  //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl;
102  G4ReactionProductVector * result = 0;
103  G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector;
104 // G4double m_nucl(0); // to check energy balance
105 
106 
107  // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
108  // G4cout << "Entering the decision point "
109  // << (mom.t()-mom.mag())/pA << " "
110  // << pA<<" "<< pZ<<" "
111  // << tA<<" "<< tZ<<G4endl
112  // << " "<<mom.t()-mom.mag()<<" "
113  // << mom.t()- m1<<G4endl;
114  if( (mom.t()-mom.mag())/pA < 50*MeV )
115  {
116  // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl;
117  // m_nucl = mom.mag();
118  cascaders=FuseNucleiAndPrompound(mom);
119  if( !cascaders )
120  {
121 
122  // abort!! happens for too low energy for nuclei to fuse
123 
124  theResult.Clear();
127  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
128  return &theResult;
129  }
130  }
131  else
132  {
133  result=Interact(mom,toBreit);
134 
135  if(! result )
136  {
137  // abort!!
138 
139  G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl;
140  G4cerr << " Primary " << aTrack.GetDefinition()
141  << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
142  << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
143  << ", kinetic energy " << aTrack.GetKineticEnergy()
144  << G4endl;
145  G4cerr << " Target nucleus (A,Z)=("
146  << (swapped?pA:tA) << ","
147  << (swapped?pZ:tZ) << ")" << G4endl;
148  G4cerr << " if frequent, please submit above information as bug report"
149  << G4endl << G4endl;
150 
151  theResult.Clear();
154  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
155  return &theResult;
156  }
157 
158  // Calculate excitation energy,
159  G4double theStatisticalExEnergy = GetProjectileExcitation();
160 
161 
162  pInitialState = mom;
163  //G4cout << "pInitialState from aTrack : " << pInitialState;
164  pInitialState.setT(pInitialState.getT() +
166  //G4cout << "target nucleus added : " << pInitialState << G4endl;
167 
170 
172 
173  cascaders = new G4ReactionProductVector;
174 
175  G4LorentzVector pspectators=SortResult(result,spectators,cascaders);
176 
177  // pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom);
178  std::vector<G4ReactionProduct *>::iterator iter;
179 
180  //G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl;
181  // if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0)
182  // {
183  // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" <<
184  // spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl;
185  // }
186  delete result;
187  result=0;
189  G4int loopcount(0);
190  //G4cout << "momentum, pspectators : " << momentum << " / " << pspectators << G4endl;
191  while (std::abs(momentum.e()-pspectators.e()) > 10*MeV)
192  {
193  G4LorentzVector pCorrect(pInitialState-pspectators);
194  //G4cout << "BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl;
195  // Correct outgoing casacde particles.... to have momentum of (initial state - spectators)
196  G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
197  if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
198  {
199  G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl;
200  }
201  pFinalState=G4LorentzVector(0,0,0,0);
202  for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
203  {
204  pFinalState += G4LorentzVector( (*iter)->GetMomentum(), (*iter)->GetTotalEnergy() );
205  }
206  momentum=pInitialState-pFinalState;
207  if (++loopcount > 10 )
208  {
209  if ( momentum.vect().mag() - momentum.e()> 10*keV )
210  {
211  G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl;
212  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
213  } else {
214  break;
215  }
216 
217  }
218  }
219 
220  if (spectatorA > 0 )
221  {
222  // check spectator momentum
223  if ( momentum.vect().mag() - momentum.e()> 10*keV )
224  {
225 
226  for (iter=spectators->begin();iter!=spectators->end();iter++)
227  {
228  delete *iter;
229  }
230  delete spectators;
231  for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
232  {
233  delete *iter;
234  }
235  delete cascaders;
236 
237  G4cout << "G4BinaryLightIonReaction.cc: mom check: " << momentum
238  << " 3.mag "<< momentum.vect().mag() << G4endl
239  << " .. pInitialState/pFinalState/spectators " << pInitialState <<" "
240  << pFinalState << " " << pspectators << G4endl
241  << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl;
242  G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl;
243  G4cout << " Primary " << aTrack.GetDefinition()
244  << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
245  << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
246  << ", kinetic energy " << aTrack.GetKineticEnergy()
247  << G4endl;
248  G4cout << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt()
249  << "," << targetNucleus.GetZ_asInt() << ")" << G4endl;
250  G4cout << " if frequent, please submit above information as bug report"
251  << G4endl << G4endl;
252 
253  theResult.Clear();
256  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
257  return &theResult;
258  }
259 
260 
261  DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
262  } else {
263  delete spectators;
264  }
265  }
266  // Rotate to lab
267  G4LorentzRotation toZ;
268  toZ.rotateZ(-1*mom.phi());
269  toZ.rotateY(-1*mom.theta());
270  G4LorentzRotation toLab(toZ.inverse());
271 
272  // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped.
273  // theResult.Clear();
274  theResult.Clear();
276  G4double Etot(0);
277  G4ReactionProductVector::iterator iter;
278  for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
279  {
280  if((*iter)->GetNewlyAdded())
281  {
282  G4DynamicParticle * aNew =
283  new G4DynamicParticle((*iter)->GetDefinition(),
284  (*iter)->GetTotalEnergy(),
285  (*iter)->GetMomentum() );
286  G4LorentzVector tmp = aNew->Get4Momentum();
287  if(swapped)
288  {
289  tmp*=toBreit.inverse();
290  tmp.setVect(-tmp.vect());
291  }
292  tmp *= toLab;
293  aNew->Set4Momentum(tmp);
294  //G4cout << "result[" << i << "], 4vect: " << tmp << G4endl;
295  theResult.AddSecondary(aNew);
296  Etot += tmp.e();
297  // G4cout << "LIBIC: Secondary " << aNew->GetDefinition()->GetParticleName()
298  // <<" "<< aNew->GetMomentum()
299  // <<" "<< aNew->GetTotalEnergy()
300  // << G4endl;
301  }
302  delete *iter;
303  }
304  delete cascaders;
305 
306 #ifdef debug_BLIR_result
307  G4cout << "Result analysis, secondaries" << theResult.GetNumberOfSecondaries() << G4endl;
308  G4cout << " Energy conservation initial/primary/nucleus/final/delta(init-final) "
309  << aTrack.GetTotalEnergy() + m_nucl << aTrack.GetTotalEnergy() << m_nucl <<Etot
310  << aTrack.GetTotalEnergy() + m_nucl - Etot;
311 #endif
312 
313  if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### " << G4endl;
314 
315  return &theResult;
316 }
317 
318 //--------------------------------------------------------------------------------
319 
320 //****************************************************************************
322  G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom)
323 //****************************************************************************
324 {
325  const int nAttemptScale = 2500;
326  const double ErrLimit = 1.E-6;
327  if (Output->empty())
328  return TRUE;
329  G4LorentzVector SumMom(0,0,0,0);
330  G4double SumMass = 0;
331  G4double TotalCollisionMass = TotalCollisionMom.m();
332  size_t i = 0;
333  // Calculate sum hadron 4-momenta and summing hadron mass
334  for(i = 0; i < Output->size(); i++)
335  {
336  SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
337  SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
338  }
339  // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass "
340  // << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl;
341  if (SumMass > TotalCollisionMass) return FALSE;
342  SumMass = SumMom.m2();
343  if (SumMass < 0) return FALSE;
344  SumMass = std::sqrt(SumMass);
345 
346  // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
347  G4ThreeVector Beta = -SumMom.boostVector();
348  //G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl;
349  //--old Output->Boost(Beta);
350  for(i = 0; i < Output->size(); i++)
351  {
352  G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
353  mom *= Beta;
354  (*Output)[i]->SetMomentum(mom.vect());
355  (*Output)[i]->SetTotalEnergy(mom.e());
356  }
357 
358  // Scale total c.m.s. hadron energy (hadron system mass).
359  // It should be equal interaction mass
360  G4double Scale = 0,OldScale=0;
361  G4double factor = 1.;
362  G4int cAttempt = 0;
363  G4double Sum = 0;
364  G4bool success = false;
365  for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
366  {
367  Sum = 0;
368  for(i = 0; i < Output->size(); i++)
369  {
370  G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
371  HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
372  G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
373  HadronMom.setE(E);
374  (*Output)[i]->SetMomentum(HadronMom.vect());
375  (*Output)[i]->SetTotalEnergy(HadronMom.e());
376  Sum += E;
377  }
378  OldScale=Scale;
379  Scale = TotalCollisionMass/Sum - 1;
380  // G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl;
381  if (std::abs(Scale) <= ErrLimit
382  || OldScale == Scale) // protect 'frozen' situation and divide by 0 in calculating new factor below
383  {
384  if (debug_G4BinaryLightIonReactionResults) G4cout << "E/p corrector: " << cAttempt << G4endl;
385  success = true;
386  break;
387  }
388  if ( cAttempt > 10 )
389  {
390  // G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl;
391  factor=std::max(1.,std::log(std::abs(OldScale/(OldScale-Scale))));
392  // G4cout << " ? factor ? " << factor << G4endl;
393  }
394  }
395 
396  if( (!success) && debug_G4BinaryLightIonReactionResults)
397  {
398  G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl;
399  G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
400  G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
401  }
402 
403  // Compute c.m.s. interaction velocity and KTV back boost
404  Beta = TotalCollisionMom.boostVector();
405  //--old Output->Boost(Beta);
406  for(i = 0; i < Output->size(); i++)
407  {
408  G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
409  mom *= Beta;
410  (*Output)[i]->SetMomentum(mom.vect());
411  (*Output)[i]->SetTotalEnergy(mom.e());
412  }
413  return TRUE;
414 }
416 {
417  G4bool swapped = false;
418  if(tA<pA)
419  {
420  swapped = true;
421  G4int tmp(0);
422  tmp = tA; tA=pA; pA=tmp;
423  tmp = tZ; tZ=pZ; pZ=tmp;
425  G4LorentzVector it(m1, G4ThreeVector(0,0,0));
426  mom = toBreit*it;
427  }
428  return swapped;
429 }
431 {
432  // Check if kinematically nuclei can fuse.
435  G4LorentzVector pCompound(mom.e()+mTarget,mom.vect());
436  G4double m2Compound=pCompound.m2();
437  if (m2Compound < sqr(mFused) ) {
438  //G4cout << "G4BLIC: projectile p, mTarget, mFused, mCompound, delta: " <<mom << " " << mTarget << " " << mFused
439  // << " " << sqrt(m2Compound)<< " " << sqrt(m2Compound) - mFused << G4endl;
440  return 0;
441  }
442 
443  G4Fragment aPreFrag;
444  aPreFrag.SetA(pA+tA);
445  aPreFrag.SetZ(pZ+tZ);
446  aPreFrag.SetNumberOfParticles(pA);
447  aPreFrag.SetNumberOfCharged(pZ);
448  aPreFrag.SetNumberOfHoles(0);
449  //GF FIXME: whyusing plop in z direction? this will not conserve momentum?
450  //G4ThreeVector plop(0.,0., mom.vect().mag());
451  //G4LorentzVector aL(mom.t()+mTarget, plop);
452  G4LorentzVector aL(mom.t()+mTarget,mom.vect());
453  aPreFrag.SetMomentum(aL);
454 
455 
456  //G4cout << "Fragment INFO "<< pA+tA <<" "<<pZ+tZ<<" "
457  // << aL <<" "<<G4endl << aPreFrag << G4endl;
459  //G4double tSum = 0;
460  for(size_t count = 0; count<cascaders->size(); count++)
461  {
462  cascaders->operator[](count)->SetNewlyAdded(true);
463  //tSum += cascaders->operator[](count)->GetKineticEnergy();
464  }
465  // G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl;
466  return cascaders;
467 }
469 {
470  G4ReactionProductVector * result = 0;
471  G4double projectileMass(0);
472  G4LorentzVector it;
473 
474  G4int tryCount(0);
475  do
476  {
477  ++tryCount;
483  it=toBreit * G4LorentzVector(projectileMass,G4ThreeVector(0,0,0));
484 
488  // G4cout << "out radius - nucleus - projectile " << target3dNucleus->GetOuterRadius()/fermi << " - " << projectile3dNucleus->GetOuterRadius()/fermi << G4endl;
489  G4double aX=(2.*G4UniformRand()-1.)*impactMax;
490  G4double aY=(2.*G4UniformRand()-1.)*impactMax;
491  G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi);
492 
493  G4KineticTrackVector * initalState = new G4KineticTrackVector;
495  G4Nucleon * aNuc;
496  G4LorentzVector tmpV(0,0,0,0);
497  G4LorentzVector nucleonMom(1./pA*mom);
498  nucleonMom.setZ(nucleonMom.vect().mag());
499  nucleonMom.setX(0);
500  nucleonMom.setY(0);
501  theFermi.Init(pA,pZ);
502  while( (aNuc=projectile3dNucleus->GetNextNucleon()) )
503  {
504  G4LorentzVector p4 = aNuc->GetMomentum();
505  tmpV+=p4;
506  G4ThreeVector nucleonPosition(aNuc->GetPosition());
507  G4double density=(projectile3dNucleus->GetNuclearDensity())->GetDensity(nucleonPosition);
508  nucleonPosition += pos;
509  G4KineticTrack * it1 = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom );
511  G4double pfermi= theFermi.GetFermiMomentum(density);
512  G4double mass = aNuc->GetDefinition()->GetPDGMass();
513  G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass;
514  it1->SetProjectilePotential(-Efermi);
515  initalState->push_back(it1);
516  }
517 
518  result=theModel->Propagate(initalState, target3dNucleus);
519  if( result && result->size()==0)
520  {
521  delete result;
522  result=0;
523  }
524  if ( ! result )
525  {
526  delete target3dNucleus;
527  delete projectile3dNucleus;
528  }
529 
530  // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>());
531  // delete initalState;
532 
533  } while (! result && tryCount< 150);
534  return result;
535 }
537 {
539 
540  G4Nucleon * aNuc;
541  // targetNucleus->StartLoop();
542  // while( (aNuc=targetNucleus->GetNextNucleon()) )
543  // {
544  // G4cout << " tgt Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl;
545  // }
546  // the projectileNucleus excitation energy estimate...
547  G4double theStatisticalExEnergy = 0;
549  while( (aNuc=projectile3dNucleus->GetNextNucleon()) )
550  {
551  // G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl;
552  if(!aNuc->AreYouHit())
553  {
554  spectatorA++;
556  }
557  else
558  {
559  G4ThreeVector aPosition(aNuc->GetPosition());
560  G4double localDensity = projectile3dNucleus->GetNuclearDensity()->GetDensity(aPosition);
561  G4double localPfermi = theFermi.GetFermiMomentum(localDensity);
562  G4double nucMass = aNuc->GetDefinition()->GetPDGMass();
563  G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
564  G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag());
565  theStatisticalExEnergy += deltaE;
566  }
567  }
568  return theStatisticalExEnergy;
569 }
570 
572 {
573  unsigned int i(0);
574  // G4int spectA(0),spectZ(0);
575  G4LorentzVector pspectators(0,0,0,0);
576  pFinalState=G4LorentzVector(0,0,0,0);
577  for(i=0; i<result->size(); i++)
578  {
579  if( (*result)[i]->GetNewlyAdded() )
580  {
581  pFinalState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
582  cascaders->push_back((*result)[i]);
583  }
584  else {
585  // G4cout <<" spectator ... ";
586  pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
587  spectators->push_back((*result)[i]);
588  // spectA++;
589  // spectZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()/eplus);
590  }
591 
592  // G4cout << (*result)[i]<< " "
593  // << (*result)[i]->GetDefinition()->GetParticleName() << " "
594  // << (*result)[i]->GetMomentum()<< " "
595  // << (*result)[i]->GetTotalEnergy() << G4endl;
596  }
597  //G4cout << "pFinalState / pspectators" << pFinalState << " / " << pspectators << G4endl;
598  return pspectators;
599 }
600 
602  G4double theStatisticalExEnergy, G4LorentzVector & pSpectators)
603 {
604  // call precompound model
605  G4ReactionProductVector * proFrag = 0;
606  G4LorentzVector pFragment(0.,0.,0.,0.);
607  // G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
608  G4LorentzRotation boost_fragments;
609  // G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
610  // G4LorentzRotation boost_spectator_mom(-momentum.boostVector());
611  // G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl;
612  G4LorentzVector pFragments(0,0,0,0);
613 
614  if(spectatorZ>0 && spectatorA>1)
615  {
616  // Make the fragment
617  G4Fragment aProRes;
618  aProRes.SetA(spectatorA);
619  aProRes.SetZ(spectatorZ);
620  aProRes.SetNumberOfParticles(0);
621  aProRes.SetNumberOfCharged(0);
622  aProRes.SetNumberOfHoles(pA-spectatorA);
624  pFragment=G4LorentzVector(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) );
625  aProRes.SetMomentum(pFragment);
626 
627  proFrag = theHandler->BreakItUp(aProRes);
628 
629  boost_fragments = G4LorentzRotation(pSpectators.boostVector());
630 
631  // G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE "
632  // << spectatorA <<" "<< spectatorZ <<" "<< mFragment <<" "
633  // << momentum.mag() <<" "<< momentum.mag() - mFragment
634  // << " "<<theStatisticalExEnergy
635  // << " "<< boost_fragments*pFragment<< G4endl;
636  G4ReactionProductVector::iterator ispectator;
637  for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
638  {
639  delete *ispectator;
640  }
641  }
642  else if(spectatorA!=0)
643  {
644  G4ReactionProductVector::iterator ispectator;
645  for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
646  {
647  (*ispectator)->SetNewlyAdded(true);
648  cascaders->push_back(*ispectator);
649  pFragments+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy());
650  // G4cout << "from spectator "
651  // << (*ispectator)->GetDefinition()->GetParticleName() << " "
652  // << (*ispectator)->GetMomentum()<< " "
653  // << (*ispectator)->GetTotalEnergy() << G4endl;
654  }
655  }
656  // / if (spectators)
657  delete spectators;
658 
659  // collect the evaporation part and boost to spectator frame
660  G4ReactionProductVector::iterator ii;
661  if(proFrag)
662  {
663  for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
664  {
665  (*ii)->SetNewlyAdded(true);
666  G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy());
667  tmp *= boost_fragments;
668  (*ii)->SetMomentum(tmp.vect());
669  (*ii)->SetTotalEnergy(tmp.e());
670  // result->push_back(*ii);
671  pFragments += tmp;
672  }
673  }
674 
675  // G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum
676  // <<" "<< pFragments-momentum << G4endl;
677 
678  // correct p/E of Cascade secondaries
679  G4LorentzVector pCas=pInitialState - pFragments;
680 
681  //G4cout <<" Going to correct from " << pFinalState << " to " << pCas << G4endl;
682  G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas);
683  if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
684  {
685  G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl;
686  }
687 
688  // Add deexcitation secondaries
689  if(proFrag)
690  {
691  for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
692  {
693  cascaders->push_back(*ii);
694  }
695  delete proFrag;
696  }
697  //G4cout << "EnergyIsCorrect? " << EnergyIsCorrect << G4endl;
698  if ( ! EnergyIsCorrect )
699  {
700  //G4cout <<" ! EnergyIsCorrect " << pFinalState << " to " << pInitialState << G4endl;
701  if (! EnergyAndMomentumCorrector(cascaders,pInitialState))
702  {
704  G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl;
705  }
706  }
707 
708 }
709 
Definition: G4ping.hh:33
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static const double MeV
Definition: G4SIunits.hh:193
G4BinaryLightIonReaction(G4VPreCompoundModel *ptr=0)
virtual void ModelDescription(std::ostream &) const
const G4LorentzVector & GetMomentum() const
Definition: G4Nucleon.hh:71
CLHEP::Hep3Vector G4ThreeVector
G4Fancy3DNucleus * projectile3dNucleus
std::ofstream outFile
Definition: GammaRayTel.cc:68
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
G4ReactionProductVector * FuseNucleiAndPrompound(const G4LorentzVector &mom)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:355
G4double a
Definition: TRTMaterials.hh:39
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
void SetProjectilePotential(const G4double aPotential)
int G4int
Definition: G4Types.hh:78
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState) const
void SetA(G4double value)
Definition: G4Fragment.hh:314
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4double density
Definition: TRTMaterials.hh:39
G4ExcitationHandler * GetExcitationHandler() const
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
CascadeState SetState(const CascadeState new_state)
const G4ParticleDefinition * GetDefinition() const
bool G4bool
Definition: G4Types.hh:79
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1232
G4VPreCompoundModel * theProjectileFragmentation
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:276
G4double GetKineticEnergy() const
#define FALSE
Definition: globals.hh:52
void Init(G4int theA, G4int theZ)
G4double GetFermiMomentum(G4double density)
G4bool SetLighterAsProjectile(G4LorentzVector &mom, const G4LorentzRotation &toBreit)
#define TRUE
Definition: globals.hh:55
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:364
G4ReactionProductVector * Interact(G4LorentzVector &mom, const G4LorentzRotation &)
void DeExciteSpectatorNucleus(G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders, G4double theStatisticalExEnergy, G4LorentzVector &momentum)
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *, G4V3DNucleus *)
const G4LorentzVector & Get4Momentum() const
G4LorentzVector Get4Momentum() const
G4HadronicInteraction * FindModel(const G4String &name)
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
G4double GetTotalEnergy() const
void Init(G4int anA, G4int aZ)
G4double GetPDGMass() const
G4double GetDensity(const G4ThreeVector &aPosition) const
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define debug
static G4HadronicInteractionRegistry * Instance()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4LorentzVector SortResult(G4ReactionProductVector *result, G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders)
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
virtual G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
void SetZ(G4double value)
Definition: G4Fragment.hh:308
static const double keV
Definition: G4SIunits.hh:195
G4double GetOuterRadius()
G4LorentzVector operator()(G4LorentzVector a, G4ReactionProduct *b)
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:369
T sqr(const T &x)
Definition: templates.hh:145
const G4VNuclearDensity * GetNuclearDensity() const
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
G4Nucleon * GetNextNucleon()
void SetMomentumChange(const G4ThreeVector &aV)
G4bool EnergyAndMomentumCorrector(G4ReactionProductVector *products, G4LorentzVector &TotalCollisionMom)
G4int GetNumberOfSecondaries() const
static const G4double pos
void AddSecondary(G4DynamicParticle *aP)
static const double fermi
Definition: G4SIunits.hh:93
G4GLOB_DLL std::ostream G4cerr
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector
G4ExcitationHandler * theHandler