ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WilsonAblationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4WilsonAblationModel.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 // * *
21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 17191/03/NL/LvH (Aurora Programme). *
29 // * *
30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Software license. *
34 // ********************************************************************
35 //
36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 //
38 // MODULE: G4WilsonAblationModel.cc
39 //
40 // Version: 1.0
41 // Date: 08/12/2009
42 // Author: P R Truscott
43 // Organisation: QinetiQ Ltd, UK
44 // Customer: ESA/ESTEC, NOORDWIJK
45 // Contract: 17191/03/NL/LvH
46 //
47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 //
49 // CHANGE HISTORY
50 // --------------
51 //
52 // 6 October 2003, P R Truscott, QinetiQ Ltd, UK
53 // Created.
54 //
55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56 // Beta release
57 //
58 // 08 December 2009, P R Truscott, QinetiQ Ltd, UK
59 // Ver 1.0
60 // Updated as a result of changes in the G4Evaporation classes. These changes
61 // affect mostly SelectSecondariesByEvaporation, and now you have variables
62 // associated with the evaporation model which can be changed:
63 // OPTxs to select the inverse cross-section
64 // OPTxs = 0 => Dostrovski's parameterization
65 // OPTxs = 1 or 2 => Chatterjee's paramaterization
66 // OPTxs = 3 or 4 => Kalbach's parameterization
67 // useSICB => use superimposed Coulomb Barrier for inverse cross
68 // sections
69 // Other problem found with G4Fragment definition using Lorentz vector and
70 // **G4ParticleDefinition**. This does not allow A and Z to be defined for the
71 // fragment for some reason. Now the fragment is defined more explicitly:
72 // G4Fragment *fragment = new G4Fragment(A, Z, lorentzVector);
73 // to avoid this quirk. Bug found in SelectSecondariesByDefault: *type is now
74 // equated to evapType[i] whereas previously it was equated to fragType[i].
75 //
76 // 06 August 2015, A. Ribon, CERN
77 // Migrated std::exp and std::pow to the faster G4Exp and G4Pow.
78 //
79 // 09 June 2017, C. Mancini Terracciano, INFN
80 // Fixed bug on the initialization of Photon Evaporation model
81 //
82 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84 //
85 #include <iomanip>
86 #include <numeric>
87 
88 #include "G4WilsonAblationModel.hh"
89 #include "G4PhysicalConstants.hh"
90 #include "G4SystemOfUnits.hh"
91 #include "Randomize.hh"
92 #include "G4ParticleTable.hh"
93 #include "G4IonTable.hh"
94 #include "G4Alpha.hh"
95 #include "G4He3.hh"
96 #include "G4Triton.hh"
97 #include "G4Deuteron.hh"
98 #include "G4Proton.hh"
99 #include "G4Neutron.hh"
106 #include "G4PhotonEvaporation.hh"
107 #include "G4LorentzVector.hh"
108 #include "G4VEvaporationChannel.hh"
109 
110 #include "G4Exp.hh"
111 #include "G4Pow.hh"
112 
113 
115 //
117 {
118 //
119 //
120 // Send message to stdout to advise that the G4Abrasion model is being used.
121 //
123 //
124 //
125 // Set the default verbose level to 0 - no output.
126 //
127  verboseLevel = 0;
128 //
129 //
130 // Set the binding energy per nucleon .... did I mention that this is a crude
131 // model for nuclear de-excitation?
132 //
133  B = 10.0 * MeV;
134 //
135 //
136 // It is possuble to switch off secondary particle production (other than the
137 // final nuclear fragment). The default is on.
138 //
139  produceSecondaries = true;
140 //
141 //
142 // Now we need to define the decay modes. We're using the G4Evaporation model
143 // to help determine the kinematics of the decay.
144 //
145  nFragTypes = 6;
146  fragType[0] = G4Alpha::Alpha();
147  fragType[1] = G4He3::He3();
148  fragType[2] = G4Triton::Triton();
150  fragType[4] = G4Proton::Proton();
152  for(G4int i=0; i<200; ++i) { fSig[i] = 0.0; }
153 //
154 //
155 // Set verboseLevel default to no output.
156 //
157  verboseLevel = 0;
160 //
161 //
162 // Set defaults for evaporation classes. These can be overridden by user
163 // "set" methods.
164 //
165  OPTxs = 3;
166  useSICB = false;
167  fragmentVector = 0;
168 }
170 //
172 {}
173 
175 //
177  (const G4Fragment &theNucleus)
178 {
179 //
180 //
181 // Initilise the pointer to the G4FragmentVector used to return the information
182 // about the breakup.
183 //
184  fragmentVector = new G4FragmentVector;
185  fragmentVector->clear();
186 //
187 //
188 // Get the A, Z and excitation of the nucleus.
189 //
190  G4int A = theNucleus.GetA_asInt();
191  G4int Z = theNucleus.GetZ_asInt();
192  G4double ex = theNucleus.GetExcitationEnergy();
193  if (verboseLevel >= 2)
194  {
195  G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
196  <<"oooooooooooooooooooooooooooooooooooooooo"
197  <<G4endl;
198  G4cout.precision(6);
199  G4cout <<"IN G4WilsonAblationModel" <<G4endl;
200  G4cout <<"Initial prefragment A=" <<A
201  <<", Z=" <<Z
202  <<", excitation energy = " <<ex/MeV <<" MeV"
203  <<G4endl;
204  }
205 //
206 //
207 // Check that there is a nucleus to speak of. It's possible there isn't one
208 // or its just a proton or neutron. In either case, the excitation energy
209 // (from the Lorentz vector) is not used.
210 //
211  if (A == 0)
212  {
213  if (verboseLevel >= 2)
214  {
215  G4cout <<"No nucleus to decay" <<G4endl;
216  G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
217  <<"oooooooooooooooooooooooooooooooooooooooo"
218  <<G4endl;
219  }
220  return fragmentVector;
221  }
222  else if (A == 1)
223  {
224  G4LorentzVector lorentzVector = theNucleus.GetMomentum();
225  lorentzVector.setE(lorentzVector.e()-ex+10.0*eV);
226  if (Z == 0)
227  {
228  G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron());
229  fragmentVector->push_back(fragment);
230  }
231  else
232  {
233  G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton());
234  fragmentVector->push_back(fragment);
235  }
236  if (verboseLevel >= 2)
237  {
238  G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl;
239  G4cout <<(*fragmentVector)[0] <<G4endl;
240  G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
241  <<"oooooooooooooooooooooooooooooooooooooooo"
242  <<G4endl;
243  }
244  return fragmentVector;
245  }
246 //
247 //
248 // Then the number of nucleons ablated (either as nucleons or light nuclear
249 // fragments) is based on a simple argument for the binding energy per nucleon.
250 //
251  G4int DAabl = (G4int) (ex / B);
252  if (DAabl > A) DAabl = A;
253 // The following lines are no longer accurate given we now treat the final fragment
254 // if (verboseLevel >= 2)
255 // G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;
256 
257 //
258 //
259 // Determine the nuclear fragment from the ablation process by sampling the
260 // Rudstam equation.
261 //
262  G4int AF = A - DAabl;
263  G4int ZF = 0;
264 
265  if (AF > 0)
266  {
267  G4Pow* g4calc = G4Pow::GetInstance();
268  G4double AFd = (G4double) AF;
269  G4double R = 11.8 / g4calc->powZ(AF, 0.45);
270  G4int minZ = std::max(1, Z - DAabl);
271 //
272 //
273 // Here we define an integral probability distribution based on the Rudstam
274 // equation assuming a constant AF.
275 //
276  G4int zmax = std::min(199, Z);
277  G4double sum = 0.0;
278  for (ZF=minZ; ZF<=zmax; ++ZF)
279  {
280  sum += G4Exp(-R*g4calc->powA(std::abs(ZF - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
281  fSig[ZF] = sum;
282  }
283 //
284 //
285 // Now sample that distribution to determine a value for ZF.
286 //
287  sum *= G4UniformRand();
288  for (ZF=minZ; ZF<=zmax; ++ZF) {
289  if(sum <= fSig[ZF]) { break; }
290  }
291  }
292  G4int DZabl = Z - ZF;
293 //
294 //
295 // Now determine the nucleons or nuclei which have bee ablated. The preference
296 // is for the production of alphas, then other nuclei in order of decreasing
297 // binding energy. The energies assigned to the products of the decay are
298 // provisional for the moment (the 10eV is just to avoid errors with negative
299 // excitation energies due to rounding).
300 //
301  G4double totalEpost = 0.0;
302  evapType.clear();
303  for (G4int ift=0; ift<nFragTypes; ift++)
304  {
305  G4ParticleDefinition *type = fragType[ift];
306  G4double n = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10);
307  G4double n1 = 1.0E+10;
308  if (fragType[ift]->GetPDGCharge() > 0.0)
309  n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10);
310  if (n > n1) n = n1;
311  if (n > 0.0)
312  {
313  G4double mass = type->GetPDGMass();
314  for (G4int j=0; j<(G4int) n; j++)
315  {
316  totalEpost += mass;
317  evapType.push_back(type);
318  }
319  DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
320  DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
321  }
322  }
323 //
324 //
325 // Determine the properties of the final nuclear fragment. Note that if
326 // the final fragment is predicted to have a nucleon number of zero, then
327 // really it's the particle last in the vector evapType which becomes the
328 // final fragment. Therefore delete this from the vector if this is the
329 // case.
330 //
331  G4double massFinalFrag = 0.0;
332  if (AF > 0)
333  massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()->
334  GetIonMass(ZF,AF);
335  else
336  {
337  G4ParticleDefinition *type = evapType[evapType.size()-1];
338  AF = type->GetBaryonNumber();
339  ZF = (G4int) (type->GetPDGCharge() + 1.0E-10);
340  evapType.erase(evapType.end()-1);
341  }
342  totalEpost += massFinalFrag;
343 //
344 //
345 // Provide verbose output on the nuclear fragment if requested.
346 //
347  if (verboseLevel >= 2)
348  {
349  G4cout <<"Final fragment A=" <<AF
350  <<", Z=" <<ZF
351  <<G4endl;
352  for (G4int ift=0; ift<nFragTypes; ift++)
353  {
354  G4ParticleDefinition *type = fragType[ift];
355  G4int n = std::count(evapType.begin(),evapType.end(),type);
356  if (n > 0)
357  G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
358  <<", number of particles emitted = " <<n <<G4endl;
359  }
360  }
361 //
362 // Add the total energy from the fragment. Note that the fragment is assumed
363 // to be de-excited and does not undergo photo-evaporation .... I did mention
364 // this is a bit of a crude model?
365 //
366  G4double massPreFrag = theNucleus.GetGroundStateMass();
367  G4double totalEpre = massPreFrag + ex;
368  G4double excess = totalEpre - totalEpost;
369 // G4Fragment *resultNucleus(theNucleus);
370  G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum());
371  G4ThreeVector boost(0.0,0.0,0.0);
372  G4int nEvap = 0;
373  if (produceSecondaries && evapType.size()>0)
374  {
375  if (excess > 0.0)
376  {
377  SelectSecondariesByEvaporation (resultNucleus);
378  nEvap = fragmentVector->size();
379  boost = resultNucleus->GetMomentum().findBoostToCM();
380  if (evapType.size() > 0)
381  SelectSecondariesByDefault (boost);
382  }
383  else
384  SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
385  }
386 
387  if (AF > 0)
388  {
390  GetIonMass(ZF,AF);
391  G4double e = mass + 10.0*eV;
392  G4double p = std::sqrt(e*e-mass*mass);
393  G4ThreeVector direction(0.0,0.0,1.0);
394  G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
395  lorentzVector.boost(-boost);
396  G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector);
397  fragmentVector->push_back(frag);
398  }
399  delete resultNucleus;
400 //
401 //
402 // Provide verbose output on the ablation products if requested.
403 //
404  if (verboseLevel >= 2)
405  {
406  if (nEvap > 0)
407  {
408  G4cout <<"----------------------" <<G4endl;
409  G4cout <<"Evaporated particles :" <<G4endl;
410  G4cout <<"----------------------" <<G4endl;
411  }
412  G4int ie = 0;
413  G4FragmentVector::iterator iter;
414  for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++)
415  {
416  if (ie == nEvap)
417  {
418 // G4cout <<*iter <<G4endl;
419  G4cout <<"---------------------------------" <<G4endl;
420  G4cout <<"Particles from default emission :" <<G4endl;
421  G4cout <<"---------------------------------" <<G4endl;
422  }
423  G4cout <<*iter <<G4endl;
424  }
425  G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
426  <<"oooooooooooooooooooooooooooooooooooooooo"
427  <<G4endl;
428  }
429 
430  return fragmentVector;
431 }
433 //
435  (G4Fragment *intermediateNucleus)
436 {
437  G4Fragment theResidualNucleus = *intermediateNucleus;
438  G4bool evaporate = true;
439  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
440  while (evaporate && evapType.size() != 0)
441  {
442 //
443 //
444 // Here's the cheaky bit. We're hijacking the G4Evaporation model, in order to
445 // more accurately sample to kinematics, but the species of the nuclear
446 // fragments will be the ones of our choosing as above.
447 //
448  std::vector <G4VEvaporationChannel*> theChannels1;
449  theChannels1.clear();
450  std::vector <G4VEvaporationChannel*>::iterator i;
451  VectorOfFragmentTypes::iterator iter;
452  std::vector <VectorOfFragmentTypes::iterator> iters;
453  iters.clear();
454  iter = std::find(evapType.begin(), evapType.end(), G4Alpha::Alpha());
455  if (iter != evapType.end())
456  {
457  theChannels1.push_back(new G4AlphaEvaporationChannel);
458  i = theChannels1.end() - 1;
459  (*i)->SetOPTxs(OPTxs);
460  (*i)->UseSICB(useSICB);
461 // (*i)->Initialize(theResidualNucleus);
462  iters.push_back(iter);
463  }
464  iter = std::find(evapType.begin(), evapType.end(), G4He3::He3());
465  if (iter != evapType.end())
466  {
467  theChannels1.push_back(new G4He3EvaporationChannel);
468  i = theChannels1.end() - 1;
469  (*i)->SetOPTxs(OPTxs);
470  (*i)->UseSICB(useSICB);
471 // (*i)->Initialize(theResidualNucleus);
472  iters.push_back(iter);
473  }
474  iter = std::find(evapType.begin(), evapType.end(), G4Triton::Triton());
475  if (iter != evapType.end())
476  {
477  theChannels1.push_back(new G4TritonEvaporationChannel);
478  i = theChannels1.end() - 1;
479  (*i)->SetOPTxs(OPTxs);
480  (*i)->UseSICB(useSICB);
481 // (*i)->Initialize(theResidualNucleus);
482  iters.push_back(iter);
483  }
484  iter = std::find(evapType.begin(), evapType.end(), G4Deuteron::Deuteron());
485  if (iter != evapType.end())
486  {
487  theChannels1.push_back(new G4DeuteronEvaporationChannel);
488  i = theChannels1.end() - 1;
489  (*i)->SetOPTxs(OPTxs);
490  (*i)->UseSICB(useSICB);
491 // (*i)->Initialize(theResidualNucleus);
492  iters.push_back(iter);
493  }
494  iter = std::find(evapType.begin(), evapType.end(), G4Proton::Proton());
495  if (iter != evapType.end())
496  {
497  theChannels1.push_back(new G4ProtonEvaporationChannel);
498  i = theChannels1.end() - 1;
499  (*i)->SetOPTxs(OPTxs);
500  (*i)->UseSICB(useSICB);
501 // (*i)->Initialize(theResidualNucleus);
502  iters.push_back(iter);
503  }
504  iter = std::find(evapType.begin(), evapType.end(), G4Neutron::Neutron());
505  if (iter != evapType.end())
506  {
507  theChannels1.push_back(new G4NeutronEvaporationChannel);
508  i = theChannels1.end() - 1;
509  (*i)->SetOPTxs(OPTxs);
510  (*i)->UseSICB(useSICB);
511 // (*i)->Initialize(theResidualNucleus);
512  iters.push_back(iter);
513  }
514  G4int nChannels = theChannels1.size();
515 
516  G4double totalProb = 0.0;
517  G4int ich = 0;
518  G4double probEvapType[6] = {0.0};
519  std::vector<G4VEvaporationChannel*>::iterator iterEv;
520  for (iterEv=theChannels1.begin(); iterEv!=theChannels1.end(); iterEv++) {
521  totalProb += (*iterEv)->GetEmissionProbability(intermediateNucleus);
522  probEvapType[ich] = totalProb;
523  ++ich;
524  }
525  if (totalProb > 0.0) {
526 //
527 //
528 // The emission probability for at least one of the evaporation channels is
529 // positive, therefore work out which one should be selected and decay
530 // the nucleus.
531 //
532  G4double xi = totalProb*G4UniformRand();
533  G4int ii = 0;
534  for (ii=0; ii<nChannels; ii++) {
535  if (xi < probEvapType[ii]) { break; }
536  }
537  if (ii >= nChannels) { ii = nChannels - 1; }
538  G4FragmentVector *evaporationResult = theChannels1[ii]->
539  BreakUpFragment(intermediateNucleus);
540  fragmentVector->push_back((*evaporationResult)[0]);
541  intermediateNucleus = (*evaporationResult)[1];
542  delete evaporationResult;
543  }
544  else
545  {
546 //
547 //
548 // Probability for further evaporation is nil so have to escape from this
549 // routine and set the energies of the secondaries to 10eV.
550 //
551  evaporate = false;
552  }
553  }
554 
555  return;
556 }
558 //
560 {
561  for (unsigned i=0; i<evapType.size(); i++)
562  {
563  G4ParticleDefinition *type = evapType[i];
564  G4double mass = type->GetPDGMass();
565  G4double e = mass + 10.0*eV;
566  G4double p = std::sqrt(e*e-mass*mass);
567  G4double costheta = 2.0*G4UniformRand() - 1.0;
568  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
570  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
571  G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
572  lorentzVector.boost(-boost);
573 // Possibility that the following line is not correctly carrying over A and Z
574 // from particle definition. Force values. PRT 03/12/2009.
575 // G4Fragment *fragment =
576 // new G4Fragment(lorentzVector, type);
577  G4int A = type->GetBaryonNumber();
578  G4int Z = (G4int) (type->GetPDGCharge() + 1.0E-10);
579  G4Fragment *fragment =
580  new G4Fragment(A, Z, lorentzVector);
581 
582  fragmentVector->push_back(fragment);
583  }
584 }
586 //
588 {
589  G4cout <<G4endl;
590  G4cout <<" *****************************************************************"
591  <<G4endl;
592  G4cout <<" Nuclear ablation model for nuclear-nuclear interactions activated"
593  <<G4endl;
594  G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
595  <<G4endl;
596  G4cout <<" !!! WARNING: This model is not well validation and should not be used for accurate simulation !!!"
597  <<G4endl;
598  G4cout <<" *****************************************************************"
599  <<G4endl;
600  G4cout << G4endl;
601 
602  return;
603 }
605 //