ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StatMF.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4StatMF.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 //
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara
30 
31 #include "G4StatMF.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Pow.hh"
35 
36 // Default constructor
37 G4StatMF::G4StatMF() : _theEnsemble(0) {}
38 
39 
40 // Destructor
41 G4StatMF::~G4StatMF() {} //{if (_theEnsemble != 0) delete _theEnsemble;}
42 
43 
45 {
46  // G4FragmentVector * theResult = new G4FragmentVector;
47 
48  if (theFragment.GetExcitationEnergy() <= 0.0) {
49  //G4FragmentVector * theResult = new G4FragmentVector;
50  //theResult->push_back(new G4Fragment(theFragment));
51  return 0;
52  }
53 
54 
55  // Maximun average multiplicity: M_0 = 2.6 for A ~ 200
56  // and M_0 = 3.3 for A <= 110
57  G4double MaxAverageMultiplicity =
59 
60 
61  // We'll use two kinds of ensembles
62  G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
63  G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
64 
65  //-------------------------------------------------------
66  // Direct simulation part (Microcanonical ensemble)
67  //-------------------------------------------------------
68 
69  // Microcanonical ensemble initialization
70  theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
71 
72  G4int Iterations = 0;
73  G4int IterationsLimit = 100000;
74  G4double Temperature = 0.0;
75 
76  G4bool FirstTime = true;
77  G4StatMFChannel * theChannel = 0;
78 
79  G4bool ChannelOk;
80  do { // Try to de-excite as much as IterationLimit permits
81  do {
82 
83  G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
84  if (theMeanMult <= MaxAverageMultiplicity) {
85  // G4cout << "MICROCANONICAL" << G4endl;
86  // Choose fragments atomic numbers and charges from direct simulation
87  theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
88  _theEnsemble = theMicrocanonicalEnsemble;
89  } else {
90  //-----------------------------------------------------
91  // Non direct simulation part (Macrocanonical Ensemble)
92  //-----------------------------------------------------
93  if (FirstTime) {
94  // Macrocanonical ensemble initialization
95  theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
96  _theEnsemble = theMacrocanonicalEnsemble;
97  FirstTime = false;
98  }
99  // G4cout << "MACROCANONICAL" << G4endl;
100  // Select calculated fragment total multiplicity,
101  // fragment atomic numbers and fragment charges.
102  theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
103  }
104 
105  ChannelOk = theChannel->CheckFragments();
106  if (!ChannelOk) delete theChannel;
107 
108  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
109  } while (!ChannelOk);
110 
111 
112  if (theChannel->GetMultiplicity() <= 1) {
113  G4FragmentVector * theResult = new G4FragmentVector;
114  theResult->push_back(new G4Fragment(theFragment));
115  delete theMicrocanonicalEnsemble;
116  if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
117  delete theChannel;
118  return theResult;
119  }
120 
121  //--------------------------------------
122  // Second part of simulation procedure.
123  //--------------------------------------
124 
125  // Find temperature of breaking channel.
126  Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature
127 
128  if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
129 
130  // Do not forget to delete this unusable channel, for which we failed to find the temperature,
131  // otherwise for very proton-reach nuclei it would lead to memory leak due to large
132  // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ()
133 
134  // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl;
135 
136  delete theChannel;
137 
138  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
139  } while (Iterations++ < IterationsLimit );
140 
141 
142 
143  // If Iterations >= IterationsLimit means that we couldn't solve for temperature
144  if (Iterations >= IterationsLimit)
145  throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
146 
147 
148  G4FragmentVector * theResult = theChannel->
149  GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature);
150 
151 
152 
153  // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!!
154  // Original nucleus 4-momentum in CM system
155  G4LorentzVector InitialMomentum(theFragment.GetMomentum());
156  InitialMomentum.boost(-InitialMomentum.boostVector());
157  G4double ScaleFactor = 0.0;
158  G4double SavedScaleFactor = 0.0;
159  do {
160  G4double FragmentsEnergy = 0.0;
161  G4FragmentVector::iterator j;
162  for (j = theResult->begin(); j != theResult->end(); j++)
163  FragmentsEnergy += (*j)->GetMomentum().e();
164  SavedScaleFactor = ScaleFactor;
165  ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
166  G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
167  for (j = theResult->begin(); j != theResult->end(); j++) {
168  ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
169  G4double Mass = (*j)->GetMomentum().m();
170  G4LorentzVector NewMomentum;
171  NewMomentum.setVect(ScaledMomentum);
172  NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
173  (*j)->SetMomentum(NewMomentum);
174  }
175  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
176  } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
177  // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178 
179  // Perform Lorentz boost
180  G4FragmentVector::iterator i;
181  for (i = theResult->begin(); i != theResult->end(); i++) {
182  G4LorentzVector FourMom = (*i)->GetMomentum();
183  FourMom.boost(theFragment.GetMomentum().boostVector());
184  (*i)->SetMomentum(FourMom);
185  }
186 
187  // garbage collection
188  delete theMicrocanonicalEnsemble;
189  if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
190  delete theChannel;
191 
192  return theResult;
193 }
194 
195 
197  const G4StatMFChannel * aChannel,
198  G4double & Temperature)
199  // This finds temperature of breaking channel.
200 {
201  G4int A = theFragment.GetA_asInt();
202  G4int Z = theFragment.GetZ_asInt();
203  G4double U = theFragment.GetExcitationEnergy();
204 
205  G4double T = std::max(Temperature,0.0012*MeV);
206  G4double Ta = T;
207  G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T);
208 
209  G4double Da = (U - TotalEnergy)/U;
210  G4double Db = 0.0;
211 
212  // bracketing the solution
213  if (Da == 0.0) {
214  Temperature = T;
215  return true;
216  } else if (Da < 0.0) {
217  do {
218  T *= 0.5;
219  if (T < 0.001*MeV) return false;
220 
221  TotalEnergy = CalcEnergy(A,Z,aChannel,T);
222 
223  Db = (U - TotalEnergy)/U;
224  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
225  } while (Db < 0.0);
226 
227  } else {
228  do {
229  T *= 1.5;
230 
231  TotalEnergy = CalcEnergy(A,Z,aChannel,T);
232 
233  Db = (U - TotalEnergy)/U;
234  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
235  } while (Db > 0.0);
236  }
237 
238  G4double eps = 1.0e-14 * std::abs(T-Ta);
239  //G4double eps = 1.0e-3 ;
240 
241  // Start the bisection method
242  for (G4int j = 0; j < 1000; j++) {
243  G4double Tc = (Ta+T)*0.5;
244  if (std::abs(Ta-Tc) <= eps) {
245  Temperature = Tc;
246  return true;
247  }
248 
249  T = Tc;
250 
251  TotalEnergy = CalcEnergy(A,Z,aChannel,T);
252 
253  G4double Dc = (U - TotalEnergy)/U;
254 
255  if (Dc == 0.0) {
256  Temperature = Tc;
257  return true;
258  }
259 
260  if (Da*Dc < 0.0) {
261  T = Tc;
262  Db = Dc;
263  } else {
264  Ta = Tc;
265  Da = Dc;
266  }
267  }
268 
269  Temperature = (Ta+T)*0.5;
270  return false;
271 }
272 
274  G4double T)
275 {
276  G4double MassExcess0 = G4NucleiProperties::GetMassExcess(A,Z);
277  G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T);
278  return -MassExcess0 + G4StatMFParameters::GetCoulomb() + ChannelEnergy;
279 }
280 
281 
282