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