Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4StatMF Class Reference

#include <G4StatMF.hh>

Inheritance diagram for G4StatMF:
Collaboration diagram for G4StatMF:

Public Member Functions

 G4StatMF ()
 
 ~G4StatMF ()
 
G4FragmentVectorBreakItUp (const G4Fragment &theNucleus)
 
- Public Member Functions inherited from G4VMultiFragmentation
 G4VMultiFragmentation ()
 
virtual ~G4VMultiFragmentation ()
 

Detailed Description

Definition at line 47 of file G4StatMF.hh.

Constructor & Destructor Documentation

G4StatMF::G4StatMF ( )

Definition at line 38 of file G4StatMF.cc.

38 : _theEnsemble(0) {}
G4StatMF::~G4StatMF ( )

Definition at line 42 of file G4StatMF.cc.

42 {} //{if (_theEnsemble != 0) delete _theEnsemble;}

Member Function Documentation

G4FragmentVector * G4StatMF::BreakItUp ( const G4Fragment theNucleus)
virtual

Implements G4VMultiFragmentation.

Definition at line 45 of file G4StatMF.cc.

46 {
47  // G4FragmentVector * theResult = new G4FragmentVector;
48 
49  if (theFragment.GetExcitationEnergy() <= 0.0) {
50  //G4FragmentVector * theResult = new G4FragmentVector;
51  //theResult->push_back(new G4Fragment(theFragment));
52  return 0;
53  }
54 
55 
56  // Maximun average multiplicity: M_0 = 2.6 for A ~ 200
57  // and M_0 = 3.3 for A <= 110
58  G4double MaxAverageMultiplicity =
59  G4StatMFParameters::GetMaxAverageMultiplicity(theFragment.GetA_asInt());
60 
61 
62  // We'll use two kinds of ensembles
63  G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
64  G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
65 
66  //-------------------------------------------------------
67  // Direct simulation part (Microcanonical ensemble)
68  //-------------------------------------------------------
69 
70  // Microcanonical ensemble initialization
71  theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
72 
73  G4int Iterations = 0;
74  G4int IterationsLimit = 100000;
75  G4double Temperature = 0.0;
76 
77  G4bool FirstTime = true;
78  G4StatMFChannel * theChannel = 0;
79 
80  G4bool ChannelOk;
81  do { // Try to de-excite as much as IterationLimit permits
82  do {
83 
84  G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
85  if (theMeanMult <= MaxAverageMultiplicity) {
86  // G4cout << "MICROCANONICAL" << G4endl;
87  // Choose fragments atomic numbers and charges from direct simulation
88  theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
89  _theEnsemble = theMicrocanonicalEnsemble;
90  } else {
91  //-----------------------------------------------------
92  // Non direct simulation part (Macrocanonical Ensemble)
93  //-----------------------------------------------------
94  if (FirstTime) {
95  // Macrocanonical ensemble initialization
96  theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
97  _theEnsemble = theMacrocanonicalEnsemble;
98  FirstTime = false;
99  }
100  // G4cout << "MACROCANONICAL" << G4endl;
101  // Select calculated fragment total multiplicity,
102  // fragment atomic numbers and fragment charges.
103  theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
104  }
105 
106  ChannelOk = theChannel->CheckFragments();
107  if (!ChannelOk) delete theChannel;
108 
109  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
110  } while (!ChannelOk);
111 
112 
113  if (theChannel->GetMultiplicity() <= 1) {
114  G4FragmentVector * theResult = new G4FragmentVector;
115  theResult->push_back(new G4Fragment(theFragment));
116  delete theMicrocanonicalEnsemble;
117  if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
118  delete theChannel;
119  return theResult;
120  }
121 
122  //--------------------------------------
123  // Second part of simulation procedure.
124  //--------------------------------------
125 
126  // Find temperature of breaking channel.
127  Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature
128 
129  if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
130 
131  // Do not forget to delete this unusable channel, for which we failed to find the temperature,
132  // otherwise for very proton-reach nuclei it would lead to memory leak due to large
133  // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ()
134 
135  // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl;
136 
137  delete theChannel;
138 
139  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
140  } while (Iterations++ < IterationsLimit );
141 
142 
143 
144  // If Iterations >= IterationsLimit means that we couldn't solve for temperature
145  if (Iterations >= IterationsLimit)
146  throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
147 
148 
149  G4FragmentVector * theResult = theChannel->
150  GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature);
151 
152 
153 
154  // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!!
155  // Original nucleus 4-momentum in CM system
156  G4LorentzVector InitialMomentum(theFragment.GetMomentum());
157  InitialMomentum.boost(-InitialMomentum.boostVector());
158  G4double ScaleFactor = 0.0;
159  G4double SavedScaleFactor = 0.0;
160  do {
161  G4double FragmentsEnergy = 0.0;
162  G4FragmentVector::iterator j;
163  for (j = theResult->begin(); j != theResult->end(); j++)
164  FragmentsEnergy += (*j)->GetMomentum().e();
165  SavedScaleFactor = ScaleFactor;
166  ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
167  G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
168  for (j = theResult->begin(); j != theResult->end(); j++) {
169  ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
170  G4double Mass = (*j)->GetMomentum().m();
171  G4LorentzVector NewMomentum;
172  NewMomentum.setVect(ScaledMomentum);
173  NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
174  (*j)->SetMomentum(NewMomentum);
175  }
176  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
177  } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
178  // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
179 
180  // Perform Lorentz boost
181  G4FragmentVector::iterator i;
182  for (i = theResult->begin(); i != theResult->end(); i++) {
183  G4LorentzVector FourMom = (*i)->GetMomentum();
184  FourMom.boost(theFragment.GetMomentum().boostVector());
185  (*i)->SetMomentum(FourMom);
186  }
187 
188  // garbage collection
189  delete theMicrocanonicalEnsemble;
190  if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
191  delete theChannel;
192 
193  return theResult;
194 }
G4double GetMeanTemperature(void) const
size_t GetMultiplicity(void)
int G4int
Definition: G4Types.hh:78
G4StatMFChannel * ChooseAandZ(const G4Fragment &theFragment)
static G4double GetMaxAverageMultiplicity(G4int A)
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
void setVect(const Hep3Vector &)
double G4double
Definition: G4Types.hh:76
G4bool CheckFragments(void)

Here is the call graph for this function:


The documentation for this class was generated from the following files: