ECCE @ EIC Software
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EvaporationChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EvaporationChannel.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 //J.M. Quesada (August2008). Based on:
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara (Oct 1998)
30 //
31 // Modified:
32 // 03-09-2008 J.M. Quesada for external choice of inverse cross section option
33 // 06-09-2008 J.M. Quesada Also external choices have been added for superimposed
34 // Coulomb barrier (if useSICB is set true, by default is false)
35 // 17-11-2010 V.Ivanchenko in constructor replace G4VEmissionProbability by
36 // G4EvaporationProbability and do not new and delete probability
37 // object at each call; use G4Pow
38 
39 #include "G4EvaporationChannel.hh"
40 #include "G4EvaporationProbability.hh"
41 #include "G4CoulombBarrier.hh"
42 #include "G4NuclearLevelData.hh"
43 #include "G4NucleiProperties.hh"
44 #include "G4Pow.hh"
45 #include "G4Log.hh"
46 #include "G4Exp.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "Randomize.hh"
50 #include "G4RandomDirection.hh"
51 #include "G4Alpha.hh"
52 
53 G4EvaporationChannel::G4EvaporationChannel(G4int anA, G4int aZ,
54  G4EvaporationProbability* aprob):
55  G4VEvaporationChannel(),
56  theA(anA),
57  theZ(aZ),
58  theProbability(aprob),
59  theCoulombBarrier(new G4CoulombBarrier(anA, aZ))
60 {
61  resA = resZ = 0;
62  mass = resMass = 0.0;
63  evapMass = G4NucleiProperties::GetNuclearMass(theA, theZ);
64  //G4cout << "G4EvaporationChannel: Z= " << theZ << " A= " << theA
65  // << " M(GeV)= " << evapMass/GeV << G4endl;
66  evapMass2 = evapMass*evapMass;
67  theLevelData = G4NuclearLevelData::GetInstance();
68 }
69 
70 G4EvaporationChannel::~G4EvaporationChannel()
71 {
72  delete theCoulombBarrier;
73 }
74 
75 void G4EvaporationChannel::Initialise()
76 {
77  theProbability->Initialise();
78  G4VEvaporationChannel::Initialise();
79 }
80 
81 G4double G4EvaporationChannel::GetEmissionProbability(G4Fragment* fragment)
82 {
83  theProbability->ResetProbability();
84  G4int fragA = fragment->GetA_asInt();
85  G4int fragZ = fragment->GetZ_asInt();
86  resA = fragA - theA;
87  resZ = fragZ - theZ;
88 
89  // Only channels which are physically allowed are taken into account
90  if(resA < theA || resA < resZ || resZ < 0 || (resA == theA && resZ < theZ)
91  || ((resA > 1) && (resA == resZ || resZ == 0)))
92  { return 0.0; }
93 
94  G4double exEnergy = fragment->GetExcitationEnergy();
95  G4double delta0 = theLevelData->GetPairingCorrection(fragZ,fragA);
96  /*
97  G4cout << "G4EvaporationChannel::Initialize Z= "<<theZ<<" A= "<<theA
98  << " FragZ= " << fragZ << " FragA= " << fragA
99  << " exEnergy= " << exEnergy << " d0= " << delta0 << G4endl;
100  */
101  if(exEnergy < delta0) { return 0.0; }
102 
103  G4double fragMass = fragment->GetGroundStateMass();
104  mass = fragMass + exEnergy;
105 
106  resMass = G4NucleiProperties::GetNuclearMass(resA, resZ);
107  G4double bCoulomb = 0.0;
108  G4double elim = 0.0;
109  if(theZ > 0) {
110  bCoulomb = theCoulombBarrier->GetCoulombBarrier(resA,resZ,exEnergy);
111 
112  // for OPTxs >0 penetration under the barrier is taken into account
113  const G4double dCB = 3.5*CLHEP::MeV;
114  elim = (0 != OPTxs) ?
115  std::max(bCoulomb*0.5, bCoulomb - dCB*theZ) : bCoulomb;
116  }
117  /*
118  G4cout << "exEnergy= " << exEnergy << " Ec= " << bCoulomb
119  << " d0= " << delta0
120  << " Free= " << mass - resMass - evapMass
121  << G4endl;
122  */
123  if(mass <= resMass + evapMass + elim) { return 0.0; }
124 
125  G4double twoMass = mass + mass;
126  G4double ekinmax =
127  ((mass-resMass)*(mass+resMass) + evapMass2)/twoMass - evapMass;
128  G4double ekinmin = 0.0;
129  if(elim > 0.0) {
130  G4double resM = std::max(mass - evapMass - elim, resMass);
131  ekinmin =
132  std::max(((mass-resM)*(mass+resM) + evapMass2)/twoMass - evapMass,0.0);
133  }
134  /*
135  G4cout << "Emin= " <<ekinmin<<" Emax= "<<ekinmax
136  << " mass= " << mass << " resM= " << resMass
137  << " evapM= " << evapMass << G4endl;
138  */
139  if(ekinmax <= ekinmin) { return 0.0; }
140 
141  theProbability->SetDecayKinematics(resZ, resA, resMass, mass);
142  G4double prob = theProbability->TotalProbability(*fragment, ekinmin,
143  ekinmax, bCoulomb,
144  exEnergy - delta0);
145  /*
146  G4cout<<"G4EvaporationChannel: prob= "<< prob << " Z= " << theZ
147  << " A= " << theA << " E1= " << ekinmin << " E2= " << ekinmax
148  << G4endl;
149  */
150  return prob;
151 }
152 
153 G4Fragment* G4EvaporationChannel::EmittedFragment(G4Fragment* theNucleus)
154 {
155  G4double ekin;
156  // assumed, that TotalProbability(...) was already called
157  // if value iz zero no possiblity to sample final state
158  if(resA <= 4 || theProbability->GetProbability() == 0.0) {
159  ekin = 0.5*(mass*mass - resMass*resMass + evapMass2)/mass - evapMass;
160  } else {
161  ekin = theProbability->SampleEnergy();
162  }
163  ekin = std::max(ekin, 0.0);
164  G4LorentzVector lv0 = theNucleus->GetMomentum();
165  G4LorentzVector lv(std::sqrt(ekin*(ekin + 2.0*evapMass))*G4RandomDirection(),
166  ekin + evapMass);
167  lv.boost(lv0.boostVector());
168 
169  G4Fragment* evFragment = new G4Fragment(theA, theZ, lv);
170  lv0 -= lv;
171  theNucleus->SetZandA_asInt(resZ, resA);
172  theNucleus->SetMomentum(lv0);
173 
174  //G4cout << "Residual: Z= " << resZ << " A= " << resA << " Eex= "
175  // << theNucleus->GetExcitationEnergy() << G4endl;
176  return evFragment;
177 }