Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 #include "G4Log.hh"
46 
47 //#define debug_G4BinaryLightIonReaction
48 //#define debug_BLIR_finalstate
49 //#define debug_BLIR_result
50 
52 : G4HadronicInteraction("Binary Light Ion Cascade"),
53  theProjectileFragmentation(ptr),
54  pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
55  projectile3dNucleus(0),target3dNucleus(0)
56 {
57  if(!ptr) {
60  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
61  if(!pre) { pre = new G4PreCompoundModel(); }
62  theProjectileFragmentation = pre;
63  }
64  theModel = new G4BinaryCascade(theProjectileFragmentation);
65  theHandler = theProjectileFragmentation->GetExcitationHandler();
66 
67  debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0;
68 }
69 
71 {}
72 
73 void G4BinaryLightIonReaction::ModelDescription(std::ostream& outFile) const
74 {
75  outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
76  << "using G4BinaryCasacde to model the interaction of a light\n"
77  << "nucleus with a nucleus.\n"
78  << "The lighter of the two nuclei is treated like a set of projectiles\n"
79  << "which are transported simultanously through the heavier nucleus.\n";
80 }
81 
82 //--------------------------------------------------------------------------------
84 {
86 };
87 
89 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus )
90 {
91  if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction starts ######### " << G4endl;
92  G4ping debug("debug_G4BinaryLightIonReaction");
93  pA=aTrack.GetDefinition()->GetBaryonNumber();
94  pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus);
95  tA=targetNucleus.GetA_asInt();
96  tZ=targetNucleus.GetZ_asInt();
97 
98  G4LorentzVector mom(aTrack.Get4Momentum());
99  //G4cout << "proj mom : " << mom << G4endl;
100  G4LorentzRotation toBreit(mom.boostVector());
101 
102  G4bool swapped=SetLighterAsProjectile(mom, toBreit);
103  //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl;
105  G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector;
106 // G4double m_nucl(0); // to check energy balance
107 
108 
109  // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
110  // G4cout << "Entering the decision point "
111  // << (mom.t()-mom.mag())/pA << " "
112  // << pA<<" "<< pZ<<" "
113  // << tA<<" "<< tZ<<G4endl
114  // << " "<<mom.t()-mom.mag()<<" "
115  // << mom.t()- m1<<G4endl;
116  if( (mom.t()-mom.mag())/pA < 50*MeV )
117  {
118  // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl;
119  // m_nucl = mom.mag();
120  cascaders=FuseNucleiAndPrompound(mom);
121  if( !cascaders )
122  {
123 
124  // abort!! happens for too low energy for nuclei to fuse
125 
126  theResult.Clear();
127  theResult.SetStatusChange(isAlive);
128  theResult.SetEnergyChange(aTrack.GetKineticEnergy());
129  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
130  return &theResult;
131  }
132  }
133  else
134  {
135  result=Interact(mom,toBreit);
136 
137  if(! result )
138  {
139  // abort!!
140 
141  G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl;
142  G4cerr << " Primary " << aTrack.GetDefinition()
143  << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
144  << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
145  << ", kinetic energy " << aTrack.GetKineticEnergy()
146  << G4endl;
147  G4cerr << " Target nucleus (A,Z)=("
148  << (swapped?pA:tA) << ","
149  << (swapped?pZ:tZ) << ")" << G4endl;
150  G4cerr << " if frequent, please submit above information as bug report"
151  << G4endl << G4endl;
152 
153  theResult.Clear();
154  theResult.SetStatusChange(isAlive);
155  theResult.SetEnergyChange(aTrack.GetKineticEnergy());
156  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
157  return &theResult;
158  }
159 
160  // Calculate excitation energy,
161  G4double theStatisticalExEnergy = GetProjectileExcitation();
162 
163 
164  pInitialState = mom;
165  //G4cout << "BLIC: pInitialState from aTrack : " << pInitialState;
166  pInitialState.setT(pInitialState.getT() +
168  //G4cout << "BLIC: target nucleus added : " << pInitialState << G4endl;
169 
170  delete target3dNucleus;target3dNucleus=0;
171  delete projectile3dNucleus;projectile3dNucleus=0;
172 
174 
175  cascaders = new G4ReactionProductVector;
176 
177  G4LorentzVector pspectators=SortResult(result,spectators,cascaders);
178 
179  // pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom);
180 
181  std::vector<G4ReactionProduct *>::iterator iter;
182 
183  // G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl;
184  // if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0)
185  // {
186  // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" <<
187  // spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl;
188  // }
189  delete result;
190  result=0;
191  G4LorentzVector momentum(pInitialState-pFinalState);
192  G4int loopcount(0);
193  //G4cout << "BLIC: momentum, pspectators : " << momentum << " / " << pspectators << G4endl;
194  while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) /* Loop checking, 31.08.2015, G.Folger */
195  // see if on loopcount
196  {
197  G4LorentzVector pCorrect(pInitialState-pspectators);
198  //G4cout << "BLIC:: BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl;
199  // Correct outgoing casacde particles.... to have momentum of (initial state - spectators)
200  G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
201  if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
202  {
203  G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl;
204  }
205  pFinalState=G4LorentzVector(0,0,0,0);
206  for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
207  {
208  pFinalState += G4LorentzVector( (*iter)->GetMomentum(), (*iter)->GetTotalEnergy() );
209  }
210  momentum=pInitialState-pFinalState;
211  if (++loopcount > 10 )
212  {
213  if ( momentum.vect().mag() - momentum.e()> 10*keV )
214  {
215  G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl;
216  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
217  } else {
218  break;
219  }
220 
221  }
222  }
223 
224  if (spectatorA > 0 )
225  {
226  // check spectator momentum
227  if ( momentum.vect().mag() - momentum.e()> 10*keV )
228  {
229 
230  for (iter=spectators->begin();iter!=spectators->end();iter++)
231  {
232  delete *iter;
233  }
234  delete spectators;
235  for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
236  {
237  delete *iter;
238  }
239  delete cascaders;
240 
241  G4cout << "G4BinaryLightIonReaction.cc: mom check: " << momentum
242  << " 3.mag "<< momentum.vect().mag() << G4endl
243  << " .. pInitialState/pFinalState/spectators " << pInitialState <<" "
244  << pFinalState << " " << pspectators << G4endl
245  << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl;
246  G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl;
247  G4cout << " Primary " << aTrack.GetDefinition()
248  << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
249  << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
250  << ", kinetic energy " << aTrack.GetKineticEnergy()
251  << G4endl;
252  G4cout << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt()
253  << "," << targetNucleus.GetZ_asInt() << ")" << G4endl;
254  G4cout << " if frequent, please submit above information as bug report"
255  << G4endl << G4endl;
256 #ifdef debug_G4BinaryLightIonReaction
258  ed << "G4BinaryLightIonreaction: Terminate for above error" << G4endl;
259  G4Exception("G4BinaryLightIonreaction::ApplyYourSelf()", "BLIC001", FatalException,
260  ed);
261 
262 #endif
263  theResult.Clear();
264  theResult.SetStatusChange(isAlive);
265  theResult.SetEnergyChange(aTrack.GetKineticEnergy());
266  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
267  return &theResult;
268  }
269 
270  // DeExciteSpectatorNucleus() also handles also case of A=1, Z=0,1
271  DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
272  } else {
273  delete spectators;
274  }
275  }
276  // Rotate to lab
277  G4LorentzRotation toZ;
278  toZ.rotateZ(-1*mom.phi());
279  toZ.rotateY(-1*mom.theta());
280  G4LorentzRotation toLab(toZ.inverse());
281 
282  // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped.
283  // theResult.Clear();
284  theResult.Clear();
285  theResult.SetStatusChange(stopAndKill);
286  G4LorentzVector ptot(0);
287  G4ReactionProductVector::iterator iter;
288  #ifdef debug_BLIR_result
289  G4LorentzVector p_raw;
290  #endif
291  //G4int i=0;
292 
293  for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
294  {
295  if((*iter)->GetNewlyAdded())
296  {
297  G4DynamicParticle * aNew =
298  new G4DynamicParticle((*iter)->GetDefinition(),
299  (*iter)->GetTotalEnergy(),
300  (*iter)->GetMomentum() );
301  G4LorentzVector tmp = aNew->Get4Momentum();
302  #ifdef debug_BLIR_result
303  p_raw+= tmp;
304  #endif
305  if(swapped)
306  {
307  tmp*=toBreit.inverse();
308  tmp.setVect(-tmp.vect());
309  }
310  tmp *= toLab;
311  aNew->Set4Momentum(tmp);
312  theResult.AddSecondary(aNew);
313  ptot += tmp;
314  //G4cout << "BLIC: Secondary " << aNew->GetDefinition()->GetParticleName()
315  // <<" "<< aNew->GetMomentum()<<" "<< aNew->GetTotalEnergy() << G4endl;
316  }
317  delete *iter;
318  }
319  delete cascaders;
320 
321 #ifdef debug_BLIR_result
322  //G4cout << "Result analysis, secondaries " << theResult.GetNumberOfSecondaries() << G4endl;
323  //G4cout << "p_tot_raw " << p_raw << " sum p final " << ptot << G4endl;
325  GetIonMass(targetNucleus.GetZ_asInt(),targetNucleus.GetA_asInt());
326  // delete? tZ=targetNucleus.GetZ_asInt();
327 
328  //G4cout << "BLIC Energy conservation initial/primary/nucleus/final/delta(init-final) "
329  // << aTrack.GetTotalEnergy() + m_nucl <<" "<< aTrack.GetTotalEnergy() <<" "<< m_nucl <<" "<<ptot.e()
330  // <<" "<< aTrack.GetTotalEnergy() + m_nucl - ptot.e() << G4endl;
331  G4cout << "BLIC momentum conservation " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl)
332  << " ptot " << ptot << " delta " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot
333  << " 3mom.mag() " << (aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot).vect().mag() << G4endl;
334 #endif
335 
336  if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### " << G4endl;
337 
338  return &theResult;
339 }
340 
341 //--------------------------------------------------------------------------------
342 
343 //****************************************************************************
344 G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector(
345  G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom)
346 //****************************************************************************
347 {
348  const int nAttemptScale = 2500;
349  const double ErrLimit = 1.E-6;
350  if (Output->empty())
351  return TRUE;
352  G4LorentzVector SumMom(0,0,0,0);
353  G4double SumMass = 0;
354  G4double TotalCollisionMass = TotalCollisionMom.m();
355  size_t i = 0;
356  // Calculate sum hadron 4-momenta and summing hadron mass
357  for(i = 0; i < Output->size(); i++)
358  {
359  SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
360  SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
361  }
362  // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass "
363  // << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl;
364  if (SumMass > TotalCollisionMass) return FALSE;
365  SumMass = SumMom.m2();
366  if (SumMass < 0) return FALSE;
367  SumMass = std::sqrt(SumMass);
368 
369  // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
370  G4ThreeVector Beta = -SumMom.boostVector();
371  //G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl;
372  //--old Output->Boost(Beta);
373  for(i = 0; i < Output->size(); i++)
374  {
375  G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
376  mom *= Beta;
377  (*Output)[i]->SetMomentum(mom.vect());
378  (*Output)[i]->SetTotalEnergy(mom.e());
379  }
380 
381  // Scale total c.m.s. hadron energy (hadron system mass).
382  // It should be equal interaction mass
383  G4double Scale = 0,OldScale=0;
384  G4double factor = 1.;
385  G4int cAttempt = 0;
386  G4double Sum = 0;
387  G4bool success = false;
388  for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
389  {
390  Sum = 0;
391  for(i = 0; i < Output->size(); i++)
392  {
393  G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
394  HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
395  G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
396  HadronMom.setE(E);
397  (*Output)[i]->SetMomentum(HadronMom.vect());
398  (*Output)[i]->SetTotalEnergy(HadronMom.e());
399  Sum += E;
400  }
401  OldScale=Scale;
402  Scale = TotalCollisionMass/Sum - 1;
403  // G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl;
404  if (std::abs(Scale) <= ErrLimit
405  || OldScale == Scale) // protect 'frozen' situation and divide by 0 in calculating new factor below
406  {
407  if (debug_G4BinaryLightIonReactionResults) G4cout << "E/p corrector: " << cAttempt << G4endl;
408  success = true;
409  break;
410  }
411  if ( cAttempt > 10 )
412  {
413  // G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl;
414  factor=std::max(1.,G4Log(std::abs(OldScale/(OldScale-Scale))));
415  // G4cout << " ? factor ? " << factor << G4endl;
416  }
417  }
418 
419  if( (!success) && debug_G4BinaryLightIonReactionResults)
420  {
421  G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl;
422  G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
423  G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
424  }
425 
426  // Compute c.m.s. interaction velocity and KTV back boost
427  Beta = TotalCollisionMom.boostVector();
428  //--old Output->Boost(Beta);
429  for(i = 0; i < Output->size(); i++)
430  {
431  G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
432  mom *= Beta;
433  (*Output)[i]->SetMomentum(mom.vect());
434  (*Output)[i]->SetTotalEnergy(mom.e());
435  }
436  return TRUE;
437 }
438 G4bool G4BinaryLightIonReaction::SetLighterAsProjectile(G4LorentzVector & mom,const G4LorentzRotation & toBreit)
439 {
440  G4bool swapped = false;
441  if(tA<pA)
442  {
443  swapped = true;
444  G4int tmp(0);
445  tmp = tA; tA=pA; pA=tmp;
446  tmp = tZ; tZ=pZ; pZ=tmp;
448  G4LorentzVector it(m1, G4ThreeVector(0,0,0));
449  mom = toBreit*it;
450  }
451  return swapped;
452 }
453 G4ReactionProductVector * G4BinaryLightIonReaction::FuseNucleiAndPrompound(const G4LorentzVector & mom)
454 {
455  // Check if kinematically nuclei can fuse.
458  G4LorentzVector pCompound(mom.e()+mTarget,mom.vect());
459  G4double m2Compound=pCompound.m2();
460  if (m2Compound < sqr(mFused) ) {
461  //G4cout << "G4BLIC: projectile p, mTarget, mFused, mCompound, delta: " <<mom << " " << mTarget << " " << mFused
462  // << " " << sqrt(m2Compound)<< " " << sqrt(m2Compound) - mFused << G4endl;
463  return 0;
464  }
465 
466  G4Fragment aPreFrag;
467  aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA);
468  aPreFrag.SetNumberOfParticles(pA);
469  aPreFrag.SetNumberOfCharged(pZ);
470  aPreFrag.SetNumberOfHoles(0);
471  //GF FIXME: whyusing plop in z direction? this will not conserve momentum?
472  //G4ThreeVector plop(0.,0., mom.vect().mag());
473  //G4LorentzVector aL(mom.t()+mTarget, plop);
474  G4LorentzVector aL(mom.t()+mTarget,mom.vect());
475  aPreFrag.SetMomentum(aL);
476 
477 
478  //G4cout << "Fragment INFO "<< pA+tA <<" "<<pZ+tZ<<" "
479  // << aL <<" "<<G4endl << aPreFrag << G4endl;
480  G4ReactionProductVector * cascaders = theProjectileFragmentation->DeExcite(aPreFrag);
481  //G4double tSum = 0;
482  for(size_t count = 0; count<cascaders->size(); count++)
483  {
484  cascaders->operator[](count)->SetNewlyAdded(true);
485  //tSum += cascaders->operator[](count)->GetKineticEnergy();
486  }
487  // G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl;
488  return cascaders;
489 }
490 G4ReactionProductVector * G4BinaryLightIonReaction::Interact(G4LorentzVector & mom, const G4LorentzRotation & toBreit)
491 {
493  G4double projectileMass(0);
494  G4LorentzVector it;
495 
496  G4int tryCount(0);
497  do
498  {
499  ++tryCount;
500  projectile3dNucleus = new G4Fancy3DNucleus;
501  projectile3dNucleus->Init(pA, pZ);
502  projectile3dNucleus->CenterNucleons();
504  projectile3dNucleus->GetCharge(),projectile3dNucleus->GetMassNumber());
505  it=toBreit * G4LorentzVector(projectileMass,G4ThreeVector(0,0,0));
506 
507  target3dNucleus = new G4Fancy3DNucleus;
508  target3dNucleus->Init(tA, tZ);
509  G4double impactMax = target3dNucleus->GetOuterRadius()+projectile3dNucleus->GetOuterRadius();
510  // G4cout << "out radius - nucleus - projectile " << target3dNucleus->GetOuterRadius()/fermi << " - " << projectile3dNucleus->GetOuterRadius()/fermi << G4endl;
511  G4double aX=(2.*G4UniformRand()-1.)*impactMax;
512  G4double aY=(2.*G4UniformRand()-1.)*impactMax;
513  G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi);
514 
515  G4KineticTrackVector * initalState = new G4KineticTrackVector;
516  projectile3dNucleus->StartLoop();
517  G4Nucleon * aNuc;
518  G4LorentzVector tmpV(0,0,0,0);
519  #ifdef debug_BLIR_finalstate
520  G4LorentzVector pinitial;
521  #endif
522  G4LorentzVector nucleonMom(1./pA*mom);
523  nucleonMom.setZ(nucleonMom.vect().mag());
524  nucleonMom.setX(0);
525  nucleonMom.setY(0);
526  theFermi.Init(pA,pZ);
527  while( (aNuc=projectile3dNucleus->GetNextNucleon()) ) /* Loop checking, 31.08.2015, G.Folger */
528  {
529  G4LorentzVector p4 = aNuc->GetMomentum();
530  tmpV+=p4;
531  G4ThreeVector nucleonPosition(aNuc->GetPosition());
532  G4double density=(projectile3dNucleus->GetNuclearDensity())->GetDensity(nucleonPosition);
533  nucleonPosition += pos;
534  G4KineticTrack * it1 = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom );
536  G4double pfermi= theFermi.GetFermiMomentum(density);
537  G4double mass = aNuc->GetDefinition()->GetPDGMass();
538  G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass;
539  it1->SetProjectilePotential(-Efermi);
540  initalState->push_back(it1);
541  #ifdef debug_BLIR_finalstate
542  pinitial += it1->Get4Momentum();
543  #endif
544  }
545 
546  result=theModel->Propagate(initalState, target3dNucleus);
547  #ifdef debug_BLIR_finalstate
548  if( result && result->size()>0)
549  {
550  G4LorentzVector presult;
551  G4ReactionProductVector::iterator iter;
553  for (iter=result->begin(); iter !=result->end(); ++iter)
554  {
555  presult += G4LorentzVector((*iter)->GetMomentum(),(*iter)->GetTotalEnergy());
556  }
557 
558  G4cout << "BLIC check result : initial " << pinitial << " mass tgt " << target3dNucleus->GetMass()
559  << " final " << presult
560  << " IF - FF " << pinitial +G4LorentzVector(target3dNucleus->GetMass()) - presult << G4endl;
561 
562  }
563  #endif
564  if( result && result->size()==0)
565  {
566  delete result;
567  result=0;
568  }
569  if ( ! result )
570  {
571  delete target3dNucleus;
572  delete projectile3dNucleus;
573  }
574 
575  // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>());
576  // delete initalState;
577 
578  } while (! result && tryCount< 150); /* Loop checking, 31.08.2015, G.Folger */
579  return result;
580 }
581 G4double G4BinaryLightIonReaction::GetProjectileExcitation()
582 {
583 
584  G4Nucleon * aNuc;
585  // the projectileNucleus excitation energy estimate...
586  G4double theStatisticalExEnergy = 0;
587  projectile3dNucleus->StartLoop();
588  while( (aNuc=projectile3dNucleus->GetNextNucleon()) ) /* Loop checking, 31.08.2015, G.Folger */
589  {
590  //G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl;
591  if(aNuc->AreYouHit()) {
592  G4ThreeVector aPosition(aNuc->GetPosition());
593  G4double localDensity = projectile3dNucleus->GetNuclearDensity()->GetDensity(aPosition);
594  G4double localPfermi = theFermi.GetFermiMomentum(localDensity);
595  G4double nucMass = aNuc->GetDefinition()->GetPDGMass();
596  G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
597  G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag());
598  theStatisticalExEnergy += deltaE;
599  }
600  }
601  return theStatisticalExEnergy;
602 }
603 
604 G4LorentzVector G4BinaryLightIonReaction::SortResult(G4ReactionProductVector * result, G4ReactionProductVector * spectators,G4ReactionProductVector * cascaders)
605 {
606  unsigned int i(0);
607  spectatorA=spectatorZ=0;
608  G4LorentzVector pspectators(0,0,0,0);
609  pFinalState=G4LorentzVector(0,0,0,0);
610  for(i=0; i<result->size(); i++)
611  {
612  if( (*result)[i]->GetNewlyAdded() )
613  {
614  pFinalState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
615  cascaders->push_back((*result)[i]);
616  }
617  else {
618  // G4cout <<" spectator ... ";
619  pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
620  spectators->push_back((*result)[i]);
621  spectatorA++;
622  spectatorZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()/eplus);
623  }
624 
625  // G4cout << (*result)[i]<< " "
626  // << (*result)[i]->GetDefinition()->GetParticleName() << " "
627  // << (*result)[i]->GetMomentum()<< " "
628  // << (*result)[i]->GetTotalEnergy() << G4endl;
629  }
630  //G4cout << "pFinalState / pspectators, (A,Z), p " << pFinalState << " / " << spectators->size()
631  // << " (" << spectatorA << ", "<< spectatorZ << "), 4-mom: " << pspectators << G4endl;
632 
633  return pspectators;
634 }
635 
636 void G4BinaryLightIonReaction::DeExciteSpectatorNucleus(G4ReactionProductVector * spectators, G4ReactionProductVector * cascaders,
637  G4double theStatisticalExEnergy, G4LorentzVector & pSpectators)
638 {
639  // call precompound model
640  G4ReactionProductVector * proFrag = 0;
641  G4LorentzVector pFragment(0.,0.,0.,0.);
642  // G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
643  G4LorentzRotation boost_fragments;
644  // G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
645  // G4LorentzRotation boost_spectator_mom(-momentum.boostVector());
646  // G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl;
647  G4LorentzVector pFragments(0,0,0,0);
648 
649  if(spectatorZ>0 && spectatorA>1)
650  {
651  // Make the fragment
652  G4Fragment aProRes;
653  aProRes.SetZandA_asInt(spectatorZ, spectatorA);
654  aProRes.SetNumberOfParticles(0);
655  aProRes.SetNumberOfCharged(0);
656  aProRes.SetNumberOfHoles(pA-spectatorA);
657  G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(spectatorZ,spectatorA);
658  pFragment=G4LorentzVector(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) );
659  aProRes.SetMomentum(pFragment);
660 
661  proFrag = theHandler->BreakItUp(aProRes);
662 
663  boost_fragments = G4LorentzRotation(pSpectators.boostVector());
664 
665  // G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE "
666  // << spectatorA <<" "<< spectatorZ <<" "<< mFragment <<" "
667  // << momentum.mag() <<" "<< momentum.mag() - mFragment
668  // << " "<<theStatisticalExEnergy
669  // << " "<< boost_fragments*pFragment<< G4endl;
670  G4ReactionProductVector::iterator ispectator;
671  for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
672  {
673  delete *ispectator;
674  }
675  }
676  else if(spectatorA!=0)
677  {
678  G4ReactionProductVector::iterator ispectator;
679  for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
680  {
681  (*ispectator)->SetNewlyAdded(true);
682  cascaders->push_back(*ispectator);
683  pFinalState+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy());
684  //G4cout << "BLIC: spectatorA>0, Z=0 from spectator "
685  // << (*ispectator)->GetDefinition()->GetParticleName() << " "
686  // << (*ispectator)->GetMomentum()<< " "
687  // << (*ispectator)->GetTotalEnergy() << G4endl;
688  }
689 
690  }
691  // / if (spectators)
692  delete spectators;
693 
694  // collect the evaporation part and boost to spectator frame
695  G4ReactionProductVector::iterator ii;
696  if(proFrag)
697  {
698  for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
699  {
700  (*ii)->SetNewlyAdded(true);
701  G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy());
702  tmp *= boost_fragments;
703  (*ii)->SetMomentum(tmp.vect());
704  (*ii)->SetTotalEnergy(tmp.e());
705  // result->push_back(*ii);
706  pFragments += tmp;
707  }
708  }
709 
710  // G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum
711  // <<" "<< pFragments-momentum << G4endl;
712 
713  // correct p/E of Cascade secondaries
714  G4LorentzVector pCas=pInitialState - pFragments;
715 
716  //G4cout <<"BLIC: Going to correct from " << pFinalState << " to " << pCas << G4endl;
717  // the creation of excited fragment did violate E/p, so correct cascaders to get overall conservation.
718  G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas);
719  if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
720  {
721  G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl;
722  }
723 
724  // Add deexcitation secondaries
725  if(proFrag)
726  {
727  for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
728  {
729  cascaders->push_back(*ii);
730  }
731  delete proFrag;
732  }
733  //G4cout << "EnergyIsCorrect? " << EnergyIsCorrect << G4endl;
734  if ( ! EnergyIsCorrect )
735  {
736  // G4cout <<" ! EnergyIsCorrect " << pFinalState << " to " << pInitialState << G4endl;
737  if (! EnergyAndMomentumCorrector(cascaders,pInitialState))
738  {
739  if(debug_G4BinaryLightIonReactionResults)
740  G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl;
741  }
742  }
743 
744 }
745 
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ping.hh:33
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
Hep3Vector boostVector() const
G4BinaryLightIonReaction(G4VPreCompoundModel *ptr=0)
virtual void ModelDescription(std::ostream &) const
const G4LorentzVector & GetMomentum() const
Definition: G4Nucleon.hh:71
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:375
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
double getT() const
void SetProjectilePotential(const G4double aPotential)
int G4int
Definition: G4Types.hh:78
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
HepLorentzRotation & rotateY(double delta)
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
tuple b
Definition: test.py:12
G4ExcitationHandler * GetExcitationHandler() const
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
CascadeState SetState(const CascadeState new_state)
const G4ParticleDefinition * GetDefinition() const
double mag() const
bool G4bool
Definition: G4Types.hh:79
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1335
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:312
G4double GetKineticEnergy() const
#define FALSE
Definition: globals.hh:52
static constexpr double eplus
Definition: G4SIunits.hh:199
void Init(G4int theA, G4int theZ)
G4double GetFermiMomentum(G4double density)
#define TRUE
Definition: globals.hh:55
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:384
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *, G4V3DNucleus *)
const G4LorentzVector & Get4Momentum() const
G4LorentzVector Get4Momentum() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
HepLorentzRotation & rotateZ(double delta)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
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
double m2() const
static G4HadronicInteractionRegistry * Instance()
Hep3Vector unit() const
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:276
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
double mag2() const
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4double GetOuterRadius()
HepLorentzRotation inverse() const
G4LorentzVector operator()(G4LorentzVector a, G4ReactionProduct *b)
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:389
T sqr(const T &x)
Definition: templates.hh:145
void setVect(const Hep3Vector &)
const G4VNuclearDensity * GetNuclearDensity() const
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
double mag() const
static constexpr double keV
Definition: G4SIunits.hh:216
G4Nucleon * GetNextNucleon()
void SetMomentumChange(const G4ThreeVector &aV)
static const G4double pos
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector