Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFChannel.cc (Version 11.3.0) and /processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFChannel.cc (Version 8.2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
                                                   >>  27 // $Id: G4StatMFChannel.cc,v 1.6 2006/06/29 20:24:45 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-08-02 $
 27 //                                                 29 //
 28 // Hadronic Process: Nuclear De-excitations        30 // Hadronic Process: Nuclear De-excitations
 29 // by V. Lara                                      31 // by V. Lara
 30 //                                             << 
 31 // Modified:                                   << 
 32 // 25.07.08 I.Pshenichnov (in collaboration wi << 
 33 //          Mishustin (FIAS, Frankfurt, INR, M << 
 34 //          Moscow, pshenich@fias.uni-frankfur << 
 35                                                << 
 36 #include <numeric>                             << 
 37                                                    32 
 38 #include "G4StatMFChannel.hh"                      33 #include "G4StatMFChannel.hh"
 39 #include "G4PhysicalConstants.hh"              << 
 40 #include "G4HadronicException.hh"                  34 #include "G4HadronicException.hh"
 41 #include "Randomize.hh"                        <<  35 #include <numeric>
 42 #include "G4Pow.hh"                            <<  36 
 43 #include "G4Exp.hh"                            <<  37 class SumCoulombEnergy : public std::binary_function<G4double,G4double,G4double>
 44 #include "G4RandomDirection.hh"                <<  38 {
 45                                                <<  39 public:
 46 G4StatMFChannel::G4StatMFChannel() :           <<  40   SumCoulombEnergy() : total(0.0) {}
 47   _NumOfNeutralFragments(0),                   <<  41   G4double operator() (G4double& , G4StatMFFragment*& frag)
 48   _NumOfChargedFragments(0)                    <<  42   { 
 49 {                                              <<  43       total += frag->GetCoulombEnergy();
 50   Pos.resize(8);                               <<  44       return total;
 51   Vel.resize(8);                               <<  45     }
 52   Accel.resize(8);                             <<  46     
 53 }                                              <<  47   G4double GetTotal() { return total; }
 54                                                <<  48 public:
 55 G4StatMFChannel::~G4StatMFChannel()            <<  49   G4double total;  
 56 {                                              <<  50 };
 57   if (!_theFragments.empty()) {                <<  51 
 58     std::for_each(_theFragments.begin(),_theFr <<  52 
 59       DeleteFragment());                       <<  53 
 60   }                                            <<  54 
                                                   >>  55 
                                                   >>  56 // Copy constructor
                                                   >>  57 G4StatMFChannel::G4StatMFChannel(const G4StatMFChannel & )
                                                   >>  58 {
                                                   >>  59     throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::copy_constructor meant to not be accessable");
 61 }                                                  60 }
 62                                                    61 
                                                   >>  62 // Operators
                                                   >>  63 
                                                   >>  64 G4StatMFChannel & G4StatMFChannel::
                                                   >>  65 operator=(const G4StatMFChannel & )
                                                   >>  66 {
                                                   >>  67     throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator= meant to not be accessable");
                                                   >>  68     return *this;
                                                   >>  69 }
                                                   >>  70 
                                                   >>  71 
                                                   >>  72 G4bool G4StatMFChannel::operator==(const G4StatMFChannel & ) const
                                                   >>  73 {
                                                   >>  74     //  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator== meant to not be accessable");
                                                   >>  75     return false;
                                                   >>  76 }
                                                   >>  77  
                                                   >>  78 
                                                   >>  79 G4bool G4StatMFChannel::operator!=(const G4StatMFChannel & ) const 
                                                   >>  80 {
                                                   >>  81     //  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator!= meant to not be accessable");
                                                   >>  82     return true;
                                                   >>  83 }
                                                   >>  84 
                                                   >>  85 
 63 G4bool G4StatMFChannel::CheckFragments(void)       86 G4bool G4StatMFChannel::CheckFragments(void)
 64 {                                                  87 {
 65   std::deque<G4StatMFFragment*>::iterator i;   <<  88     std::deque<G4StatMFFragment*>::iterator i;
 66   for (i = _theFragments.begin();              <<  89     for (i = _theFragments.begin(); 
 67        i != _theFragments.end(); ++i)          <<  90    i != _theFragments.end(); ++i) 
 68     {                                          <<  91       {
 69       G4int A = (*i)->GetA();                  <<  92   G4int A = static_cast<G4int>((*i)->GetA());
 70       G4int Z = (*i)->GetZ();                  <<  93   G4int Z = static_cast<G4int>((*i)->GetZ());
 71       if ( (A > 1 && (Z > A || Z <= 0)) || (A= <<  94   if (A > 1 && (Z >= A || Z <= 0) || (A==1 && Z > A) || A <= 0) return false;
 72     }                                              95     }
                                                   >>  96     
 73     return true;                                   97     return true;
 74 }                                                  98 }
 75                                                    99 
 76 void G4StatMFChannel::CreateFragment(G4int A,  << 100 
 77 // Create a new fragment.                      << 101 
 78 // Fragments are automatically sorted: first c << 102 
 79 // then neutral ones.                          << 103 void G4StatMFChannel::CreateFragment(const G4double A, const G4double Z)
 80 {                                              << 104     // Create a new fragment.
 81   if (Z <= 0.5) {                              << 105     // Fragments are automatically sorted: first charged fragments, 
 82     _theFragments.push_back(new G4StatMFFragme << 106     // then neutral ones.
 83     _NumOfNeutralFragments++;                  << 107 {
 84   } else {                                     << 108     if (Z <= 0.5) {
 85     _theFragments.push_front(new G4StatMFFragm << 109   _theFragments.push_back(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z)));
 86     _NumOfChargedFragments++;                  << 110   _NumOfNeutralFragments++;
 87   }                                            << 111     } else {
                                                   >> 112   _theFragments.push_front(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z)));
                                                   >> 113   _NumOfChargedFragments++;
                                                   >> 114     }
 88                                                   115   
 89   return;                                      << 116     return;
 90 }                                                 117 }
 91                                                   118 
                                                   >> 119 
 92 G4double G4StatMFChannel::GetFragmentsCoulombE    120 G4double G4StatMFChannel::GetFragmentsCoulombEnergy(void)
 93 {                                                 121 {
 94   G4double Coulomb =                           << 122     G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
 95     std::accumulate(_theFragments.begin(),_the << 123            0.0,SumCoulombEnergy());
 96                     0.0,                       << 124 //      G4double Coulomb = 0.0;
 97                     [](const G4double& running << 125 //      for (unsigned int i = 0;i < _theFragments.size(); i++)
 98                        G4StatMFFragment*& frag << 126 //    Coulomb += _theFragments[i]->GetCoulombEnergy();
 99                     {                          << 127     return Coulomb;
100                       return running_total + f << 
101                     } );                       << 
102   //      G4double Coulomb = 0.0;              << 
103   //      for (unsigned int i = 0;i < _theFrag << 
104   //    Coulomb += _theFragments[i]->GetCoulom << 
105   return Coulomb;                              << 
106 }                                                 128 }
107                                                   129 
108 G4double G4StatMFChannel::GetFragmentsEnergy(G << 130 
                                                   >> 131 
                                                   >> 132 G4double G4StatMFChannel::GetFragmentsEnergy(const G4double T) const
109 {                                                 133 {
110   G4double Energy = 0.0;                       << 134     G4double Energy = 0.0;
111                                                   135   
112   G4double TranslationalEnergy = 1.5*T*_theFra << 136     G4double TranslationalEnergy = (3./2.)*T*static_cast<G4double>(_theFragments.size());
113                                                   137 
114   std::deque<G4StatMFFragment*>::const_iterato << 138     std::deque<G4StatMFFragment*>::const_iterator i;
115   for (i = _theFragments.begin(); i != _theFra << 139     for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
116     {                                          << 140       {
117       Energy += (*i)->GetEnergy(T);            << 141   Energy += (*i)->GetEnergy(T);
118     }                                          << 142       }
119   return Energy + TranslationalEnergy;         << 143     return Energy + TranslationalEnergy;  
120 }                                                 144 }
121                                                   145 
122 G4FragmentVector * G4StatMFChannel::GetFragmen << 146 G4FragmentVector * G4StatMFChannel::GetFragments(const G4double anA, 
123              G4int anZ,                        << 147              const G4double anZ,
124              G4double T)                       << 148              const G4double T)
                                                   >> 149     //
125 {                                                 150 {
126   // calculate momenta of charged fragments    << 151     // calculate momenta of charged fragments  
127   CoulombImpulse(anA,anZ,T);                   << 152     CoulombImpulse(anA,anZ,T);
128                                                   153   
129   // calculate momenta of neutral fragments    << 154     // calculate momenta of neutral fragments
130   FragmentsMomenta(_NumOfNeutralFragments, _Nu << 155     FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
                                                   >> 156 
131                                                   157 
132   G4FragmentVector * theResult = new G4Fragmen << 158     G4FragmentVector * theResult = new G4FragmentVector;
133   std::deque<G4StatMFFragment*>::iterator i;   << 159     std::deque<G4StatMFFragment*>::iterator i;
134   for (i = _theFragments.begin(); i != _theFra << 160     for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
135     theResult->push_back((*i)->GetFragment(T)) << 161   theResult->push_back((*i)->GetFragment(T));
                                                   >> 162 
                                                   >> 163     return theResult;
136                                                   164 
137   return theResult;                            << 
138 }                                                 165 }
139                                                   166 
140 void G4StatMFChannel::CoulombImpulse(G4int anA << 167 
141 // Aafter breakup, fragments fly away under Co << 168 
142 // This method calculates asymptotic fragments << 169 void G4StatMFChannel::CoulombImpulse(const G4double anA, const G4double anZ, const G4double T)
                                                   >> 170     // Aafter breakup, fragments fly away under Coulomb field.
                                                   >> 171     // This method calculates asymptotic fragments momenta.
143 {                                                 172 {
144   // First, we have to place the fragments ins << 173     // First, we have to place the fragments inside of the original nucleus volume
145   PlaceFragments(anA);                         << 174     PlaceFragments(anA);
146                                                   175   
147   // Second, we sample initial charged fragmen << 176   // Second, we sample initial charged fragments momenta. There are
148   // _NumOfChargedFragments charged fragments  << 177   // _NumOfChargedFragments charged fragments and they start at the begining
149   // of the vector _theFragments (i.e. 0)      << 178   // of the vector _theFragments (i.e. 0) 
150   FragmentsMomenta(_NumOfChargedFragments, 0,  << 179     FragmentsMomenta(_NumOfChargedFragments, 0, T);
151                                                   180 
152   // Third, we have to figure out the asymptot << 181     // Third, we have to figure out the asymptotic momenta of charged fragments 
153   // For taht we have to solve equations of mo << 182     // For taht we have to solve equations of motion for fragments
154   SolveEqOfMotion(anA,anZ,T);                  << 183     SolveEqOfMotion(anA,anZ,T);
155                                                   184 
156   return;                                      << 185     return;
157 }                                                 186 }
158                                                   187 
159 void G4StatMFChannel::PlaceFragments(G4int anA << 
160 // This gives the position of fragments at the << 
161 // Fragments positions are sampled inside prol << 
162 {                                              << 
163   G4Pow* g4calc = G4Pow::GetInstance();        << 
164   const G4double R0 = G4StatMFParameters::Getr << 
165   G4double Rsys = 2.0*R0*g4calc->Z13(anA);     << 
166                                                   188 
167   G4bool TooMuchIterations;                    << 189 
168   do                                           << 190 void G4StatMFChannel::PlaceFragments(const G4double anA)
169     {                                          << 191     // This gives the position of fragments at the breakup instant. 
170       TooMuchIterations = false;               << 192     // Fragments positions are sampled inside prolongated ellipsoid.
                                                   >> 193 {
                                                   >> 194     const G4double R0 = G4StatMFParameters::Getr0();
                                                   >> 195     const G4double Rsys = 2.0*R0*std::pow(anA,1./3.);
                                                   >> 196 
                                                   >> 197     G4bool TooMuchIterations;
                                                   >> 198     do 
                                                   >> 199       {
                                                   >> 200   TooMuchIterations = false;
171                                                   201   
172       // Sample the position of the first frag << 202   // Sample the position of the first fragment
173       G4double R = (Rsys - R0*g4calc->Z13(_the << 203   G4double R = (Rsys - R0*std::pow(_theFragments[0]->GetA(),1./3.))*
174   g4calc->A13(G4UniformRand());                << 204     std::pow(G4UniformRand(),1./3.);
175       _theFragments[0]->SetPosition(R*G4Random << 205   _theFragments[0]->SetPosition(IsotropicVector(R));
176                                                   206   
177                                                   207   
178       // Sample the position of the remaining  << 208   // Sample the position of the remaining fragments
179       G4bool ThereAreOverlaps = false;         << 209   G4bool ThereAreOverlaps = false;
180       std::deque<G4StatMFFragment*>::iterator  << 210   std::deque<G4StatMFFragment*>::iterator i;
181       for (i = _theFragments.begin()+1; i != _ << 211   for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i) 
182   {                                            << 212     {
183     G4int counter = 0;                         << 213       G4int counter = 0;
184     do                                         << 214       do 
185       {                                        << 215         {
186         R = (Rsys - R0*g4calc->Z13((*i)->GetA( << 216     R = (Rsys - R0*std::pow((*i)->GetA(),1./3.))*std::pow(G4UniformRand(),1./3.);
187         (*i)->SetPosition(R*G4RandomDirection( << 217     (*i)->SetPosition(IsotropicVector(R));
188                                                   218     
189         // Check that there are not overlappin << 219     // Check that there are not overlapping fragments
190         std::deque<G4StatMFFragment*>::iterato << 220     std::deque<G4StatMFFragment*>::iterator j;
191         for (j = _theFragments.begin(); j != i << 221     for (j = _theFragments.begin(); j != i; ++j) 
192     {                                          << 222       {
193       G4ThreeVector FragToFragVector =         << 223         G4ThreeVector FragToFragVector = (*i)->GetPosition() - (*j)->GetPosition();
194         (*i)->GetPosition() - (*j)->GetPositio << 224         G4double Rmin = R0*(std::pow((*i)->GetA(),1./3.) +
195       G4double Rmin = R0*(g4calc->Z13((*i)->Ge << 225           std::pow((*j)->GetA(),1./3));
196               g4calc->Z13((*j)->GetA()));      << 226         if (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)) break;
197       if ( (ThereAreOverlaps = (FragToFragVect << 227       }
198         { break; }                             << 228     counter++;
199     }                                          << 229         } while (ThereAreOverlaps && counter < 1000);
200         counter++;                             << 
201         // Loop checking, 05-Aug-2015, Vladimi << 
202       } while (ThereAreOverlaps && counter < 1 << 
203                                                   230       
204     if (counter >= 1000)                       << 231       if (counter >= 1000) 
205       {                                        << 232         {
206         TooMuchIterations = true;              << 233     TooMuchIterations = true;
207         break;                                 << 234     break;
208       }                                        << 235         } 
209   }                                            << 236     }
210       // Loop checking, 05-Aug-2015, Vladimir  << 
211     } while (TooMuchIterations);                  237     } while (TooMuchIterations);
                                                   >> 238     
212     return;                                       239     return;
213 }                                                 240 }
214                                                   241 
215 void G4StatMFChannel::FragmentsMomenta(G4int N << 242 
216                G4double T)                     << 243 void G4StatMFChannel::FragmentsMomenta(const G4int NF, const G4int idx,
217 // Calculate fragments momenta at the breakup  << 244                const G4double T)
218 // Fragment kinetic energies are calculated ac << 245     // Calculate fragments momenta at the breakup instant
219 // Boltzmann distribution at given temperature << 246     // Fragment kinetic energies are calculated according to the 
220 // NF is number of fragments                   << 247     // Boltzmann distribution at given temperature.
221 // idx is index of first fragment              << 248     // NF is number of fragments
                                                   >> 249     // idx is index of first fragment
222 {                                                 250 {
223   G4double KinE = 1.5*T*NF;                    << 251     G4double KinE = (3./2.)*T*static_cast<G4double>(NF);
224   G4ThreeVector p(0.,0.,0.);                   << 
225                                                   252   
226   if (NF <= 0) return;                         << 253     G4ThreeVector p;
227   else if (NF == 1)                            << 254   
228     {                                          << 255     if (NF <= 0) return;
229       // We have only one fragment to deal wit << 256     else if (NF == 1) 
230       p = std::sqrt(2.0*_theFragments[idx]->Ge << 257       {
231   *G4RandomDirection();                        << 258   // We have only one fragment to deal with
232       _theFragments[idx]->SetMomentum(p);      << 259   p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE));
233     }                                          << 260   _theFragments[idx]->SetMomentum(p);
234   else if (NF == 2)                            << 261       } 
235     {                                          << 262     else if (NF == 2) 
236       // We have only two fragment to deal wit << 263       {
237       G4double M1 = _theFragments[idx]->GetNuc << 264   // We have only two fragment to deal with
238       G4double M2 = _theFragments[idx+1]->GetN << 265   G4double M1 = _theFragments[idx]->GetNuclearMass();
239       p = std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))* << 266   G4double M2 = _theFragments[idx+1]->GetNuclearMass();
240       _theFragments[idx]->SetMomentum(p);      << 267   p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2)));   
241       _theFragments[idx+1]->SetMomentum(-p);   << 268   _theFragments[idx]->SetMomentum(p);
242     }                                          << 269   _theFragments[idx+1]->SetMomentum(-p);
243   else                                         << 270       } 
244     {                                          << 271     else 
245       // We have more than two fragments       << 272       {
246       G4double AvailableE;                     << 273   // We have more than two fragments
247       G4int i1,i2;                             << 274   G4double AvailableE;
248       G4double SummedE;                        << 275   G4int i1,i2;
249       G4ThreeVector SummedP(0.,0.,0.);         << 276   G4double SummedE;
250       do                                       << 277   G4ThreeVector SummedP;
251   {                                            << 278   do 
252     // Fisrt sample momenta of NF-2 fragments  << 279     {
253     // according to Boltzmann distribution     << 280       // Fisrt sample momenta of NF-2 fragments 
254     AvailableE = 0.0;                          << 281       // according to Boltzmann distribution
255     SummedE = 0.0;                             << 282       AvailableE = 0.0;
256     SummedP.setX(0.0);SummedP.setY(0.0);Summed << 283       SummedE = 0.0;
257     for (G4int i = idx; i < idx+NF-2; ++i)     << 284       SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
258       {                                        << 285       for (G4int i = idx; i < idx+NF-2; i++) 
259         G4double E;                            << 286         {
260         G4double RandE;                        << 287     G4double E;
261         do                                     << 288     G4double RandE;
262     {                                          << 289     G4double Boltzmann;
263       E = 9.0*G4UniformRand();                 << 290     do 
264       RandE = std::sqrt(0.5/E)*G4Exp(E-0.5)*G4 << 291       {
265     }                                          << 292         E = 9.0*T*G4UniformRand();
266         // Loop checking, 05-Aug-2015, Vladimi << 293         Boltzmann = std::sqrt(E)*std::exp(-E/T);
267         while (RandE > 1.0);                   << 294         RandE = std::sqrt(T/2.)*std::exp(-0.5)*G4UniformRand();
268         E *= T;                                << 295       } 
269         p = std::sqrt(2.0*E*_theFragments[i]-> << 296     while (RandE > Boltzmann);
270     *G4RandomDirection();                      << 297     p = IsotropicVector(std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass()));
271         _theFragments[i]->SetMomentum(p);      << 298     _theFragments[i]->SetMomentum(p);
272         SummedE += E;                          << 299     SummedE += E;
273         SummedP += p;                          << 300     SummedP += p;
274       }                                        << 301         } 
275     // Calculate momenta of last two fragments << 302       // Calculate momenta of last two fragments in such a way
276     // that constraints are satisfied          << 303       // that constraints are satisfied
277     i1 = idx+NF-2;  // before last fragment in << 304       i1 = idx+NF-2;  // before last fragment index
278     i2 = idx+NF-1;  // last fragment index     << 305       i2 = idx+NF-1;  // last fragment index
279     p = -SummedP;                              << 306       p = -SummedP;
280     AvailableE = KinE - SummedE;               << 307       AvailableE = KinE - SummedE;
281     // Available Kinetic Energy should be shar << 308       // Available Kinetic Energy should be shared between two last fragments
282   }                                            << 309     } 
283       // Loop checking, 05-Aug-2015, Vladimir  << 310   while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
284       while (AvailableE <= p.mag2()/(2.0*(_the << 311               _theFragments[i2]->GetNuclearMass())));   
285             _theFragments[i2]->GetNuclearMass( << 312     
286       G4double H = 1.0 + _theFragments[i2]->Ge << 313   G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()/_theFragments[i1]->GetNuclearMass();
287   /_theFragments[i1]->GetNuclearMass();        << 314   G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()*AvailableE/p.mag2());
288       G4double CTM12 = H*(1.0 - 2.0*_theFragme << 315   G4double CosTheta1;
289         *AvailableE/p.mag2());                 << 316   G4double Sign;
290       G4double CosTheta1;                      << 
291       G4double Sign;                           << 
292                                                << 
293       if (CTM12 > 1.) {CosTheta1 = 1.;}        << 
294       else {                                   << 
295   do                                              317   do 
296     {                                             318     {
297       do                                          319       do 
298         {                                         320         {
299     CosTheta1 = 1.0 - 2.0*G4UniformRand();        321     CosTheta1 = 1.0 - 2.0*G4UniformRand();
300         }                                         322         } 
301       // Loop checking, 05-Aug-2015, Vladimir  << 
302       while (CosTheta1*CosTheta1 < CTM12);        323       while (CosTheta1*CosTheta1 < CTM12);
303     }                                             324     }
304   // Loop checking, 05-Aug-2015, Vladimir Ivan << 
305   while (CTM12 >= 0.0 && CosTheta1 < 0.0);        325   while (CTM12 >= 0.0 && CosTheta1 < 0.0);
306       }                                        << 
307                                                   326 
308       if (CTM12 < 0.0) Sign = 1.0;             << 327   if (CTM12 < 0.0) Sign = 1.0;
309       else if (G4UniformRand() <= 0.5) Sign =  << 328   else if (G4UniformRand() <= 0.5) Sign = -1.0;
310       else Sign = 1.0;                         << 329   else Sign = 1.0;
311                                                   330     
312       G4double P1 = (p.mag()*CosTheta1+Sign*st << 
313                   *(CosTheta1*CosTheta1-CTM12) << 
314       G4double P2 = std::sqrt(P1*P1+p.mag2() - << 
315       G4double Phi = twopi*G4UniformRand();    << 
316       G4double SinTheta1 = std::sqrt(1.0 - Cos << 
317       G4double CosPhi1 = std::cos(Phi);        << 
318       G4double SinPhi1 = std::sin(Phi);        << 
319       G4double CosPhi2 = -CosPhi1;             << 
320       G4double SinPhi2 = -SinPhi1;             << 
321       G4double CosTheta2 = (p.mag2() + P2*P2 - << 
322       G4double SinTheta2 = 0.0;                << 
323       if (CosTheta2 > -1.0 && CosTheta2 < 1.0) << 
324   SinTheta2 = std::sqrt(1.0 - CosTheta2*CosThe << 
325       }                                        << 
326       G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1 << 
327       G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2 << 
328       G4ThreeVector b(1.0,0.0,0.0);            << 
329                                                << 
330       p1 = RotateMomentum(p,b,p1);             << 
331       p2 = RotateMomentum(p,b,p2);             << 
332                                                << 
333       SummedP += p1 + p2;                      << 
334       SummedE += p1.mag2()/(2.0*_theFragments[ << 
335   p2.mag2()/(2.0*_theFragments[i2]->GetNuclear << 
336                                                   331     
337       _theFragments[i1]->SetMomentum(p1);      << 332   G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()*(CosTheta1*CosTheta1-CTM12)))/H;
338       _theFragments[i2]->SetMomentum(p2);      << 333   G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
                                                   >> 334   G4double Phi = twopi*G4UniformRand();
                                                   >> 335   G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
                                                   >> 336   G4double CosPhi1 = std::cos(Phi);
                                                   >> 337   G4double SinPhi1 = std::sin(Phi);
                                                   >> 338   G4double CosPhi2 = -CosPhi1;
                                                   >> 339   G4double SinPhi2 = -SinPhi1;
                                                   >> 340   G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
                                                   >> 341   G4double SinTheta2 = 0.0;
                                                   >> 342   if (CosTheta2 > -1.0 && CosTheta2 < 1.0) SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
                                                   >> 343   
                                                   >> 344   G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
                                                   >> 345   G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
                                                   >> 346   G4ThreeVector b(1.0,0.0,0.0);
                                                   >> 347   
                                                   >> 348   p1 = RotateMomentum(p,b,p1);
                                                   >> 349   p2 = RotateMomentum(p,b,p2);
                                                   >> 350   
                                                   >> 351   SummedP += p1 + p2;
                                                   >> 352   SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) + 
                                                   >> 353       p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());    
339                                                   354     
340     }                                          << 355   _theFragments[i1]->SetMomentum(p1);
341   return;                                      << 356   _theFragments[i2]->SetMomentum(p2);
                                                   >> 357     
                                                   >> 358       }
                                                   >> 359 
                                                   >> 360     return;
342 }                                                 361 }
343                                                   362 
344 void G4StatMFChannel::SolveEqOfMotion(G4int an << 363 
345 // This method will find a solution of Newton' << 364 void G4StatMFChannel::SolveEqOfMotion(const G4double anA, const G4double anZ, const G4double T)
346 // for fragments in the self-consistent time-d << 365     // This method will find a solution of Newton's equation of motion
347 {                                              << 366     // for fragments in the self-consistent time-dependent Coulomb field
348   G4Pow* g4calc = G4Pow::GetInstance();        << 367 {
349   G4double CoulombEnergy = 0.6*CLHEP::elm_coup << 368   G4double CoulombEnergy = (3./5.)*(elm_coupling*anZ*anZ)*
350     g4calc->A13(1.0+G4StatMFParameters::GetKap << 369     std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
351     (G4StatMFParameters::Getr0()*g4calc->Z13(a << 370     (G4StatMFParameters::Getr0()*std::pow(anA,1./3.))
                                                   >> 371     - GetFragmentsCoulombEnergy();
352   if (CoulombEnergy <= 0.0) return;               372   if (CoulombEnergy <= 0.0) return;
353                                                   373   
354   G4int Iterations = 0;                           374   G4int Iterations = 0;
355   G4double TimeN = 0.0;                           375   G4double TimeN = 0.0;
356   G4double TimeS = 0.0;                           376   G4double TimeS = 0.0;
357   G4double DeltaTime = 10.0;                      377   G4double DeltaTime = 10.0;
358   if (_NumOfChargedFragments > (G4int)Pos.size << 378   
359     Pos.resize(_NumOfChargedFragments);        << 379   G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments];
360     Vel.resize(_NumOfChargedFragments);        << 380   G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments];
361     Accel.resize(_NumOfChargedFragments);      << 381   G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments];
362   }                                            << 
363                                                   382   
364   G4int i;                                        383   G4int i;
365   for (i = 0; i < _NumOfChargedFragments; ++i) << 384   for (i = 0; i < _NumOfChargedFragments; i++) 
366     {                                             385     {
367       Vel[i] = (1.0/(_theFragments[i]->GetNucl    386       Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
368   _theFragments[i]->GetMomentum();                387   _theFragments[i]->GetMomentum();
369       Pos[i] = _theFragments[i]->GetPosition()    388       Pos[i] = _theFragments[i]->GetPosition();
370     }                                             389     }
                                                   >> 390   
                                                   >> 391   do 
                                                   >> 392     {
371                                                   393 
372   G4ThreeVector distance(0.,0.,0.);            << 394       G4ThreeVector distance;
373   G4ThreeVector force(0.,0.,0.);               << 395       G4ThreeVector force;
374   G4ThreeVector SavedVel(0.,0.,0.);            << 396 
375   do {                                         << 397       for (i = 0; i < _NumOfChargedFragments; i++) 
376     for (i = 0; i < _NumOfChargedFragments; ++ << 398   {
377       {                                        << 399     force.setX(0.0); force.setY(0.0); force.setZ(0.0);
378   force.set(0.,0.,0.);                         << 400     for (G4int j = 0; j < _NumOfChargedFragments; j++) 
379   for (G4int j = 0; j < _NumOfChargedFragments << 401       {
380     {                                          << 402         if (i != j) 
381       if (i != j)                              << 403     {
382         {                                      << 404       distance = Pos[i] - Pos[j];
383     distance = Pos[i] - Pos[j];                << 405       force += (elm_coupling*(_theFragments[i]->GetZ()*_theFragments[j]->GetZ())/
384     force += (_theFragments[i]->GetZ()*_theFra << 406           (distance.mag2()*distance.mag()))*distance;
385         (distance.mag2()*distance.mag()))*dist << 407     }
386         }                                      << 408       }
387     }                                          << 409     Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
388   Accel[i] = CLHEP::elm_coupling*CLHEP::fermi* << 410   }
389       }                                        << 
390                                                   411 
391     TimeN = TimeS + DeltaTime;                 << 412       TimeN = TimeS + DeltaTime;
392                                                   413   
393     for ( i = 0; i < _NumOfChargedFragments; + << 414       G4ThreeVector SavedVel;
394       {                                        << 415       for ( i = 0; i < _NumOfChargedFragments; i++) 
395   SavedVel = Vel[i];                           << 416   {
396   Vel[i] += Accel[i]*(TimeN-TimeS);            << 417     SavedVel = Vel[i];
397   Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0. << 418     Vel[i] += Accel[i]*(TimeN-TimeS);
398       }                                        << 419     Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
399     TimeS = TimeN;                             << 420   }
400                                                << 421       
401     // Loop checking, 05-Aug-2015, Vladimir Iv << 422       //    if (Iterations >= 50 && Iterations < 75) DeltaTime = 4.;
402   } while (Iterations++ < 100);                << 423       //    else if (Iterations >= 75) DeltaTime = 10.;
                                                   >> 424 
                                                   >> 425       TimeS = TimeN;
                                                   >> 426 
                                                   >> 427     } 
                                                   >> 428   while (Iterations++ < 100);
403                                                   429   
404   // Summed fragment kinetic energy               430   // Summed fragment kinetic energy
405   G4double TotalKineticEnergy = 0.0;              431   G4double TotalKineticEnergy = 0.0;
406   for (i = 0; i < _NumOfChargedFragments; ++i) << 432   for (i = 0; i < _NumOfChargedFragments; i++) 
407     {                                             433     {
408       TotalKineticEnergy += _theFragments[i]->    434       TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
409   0.5*Vel[i].mag2();                              435   0.5*Vel[i].mag2();
410     }                                             436     }
411   // Scaling of fragment velocities               437   // Scaling of fragment velocities
412   G4double KineticEnergy = 1.5*_NumOfChargedFr << 438   G4double KineticEnergy = (3./2.)*static_cast<G4double>(_theFragments.size())*T;
413   G4double Eta = std::pow(( CoulombEnergy + Ki << 439   G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
414                                                << 440   for (i = 0; i < _NumOfChargedFragments; i++) 
                                                   >> 441     {
                                                   >> 442       Vel[i] *= Eta;
                                                   >> 443     }
                                                   >> 444   
415   // Finally calculate fragments momenta          445   // Finally calculate fragments momenta
416   for (i = 0; i < _NumOfChargedFragments; ++i) << 446   for (i = 0; i < _NumOfChargedFragments; i++) 
417     {                                             447     {
418       _theFragments[i]->SetMomentum((_theFragm << 448       _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
419     }                                             449     }
                                                   >> 450   
                                                   >> 451   // garbage collection
                                                   >> 452   delete [] Pos;
                                                   >> 453   delete [] Vel;
                                                   >> 454   delete [] Accel;
                                                   >> 455   
420   return;                                         456   return;
421 }                                                 457 }
422                                                   458 
                                                   >> 459 
                                                   >> 460 
423 G4ThreeVector G4StatMFChannel::RotateMomentum(    461 G4ThreeVector G4StatMFChannel::RotateMomentum(G4ThreeVector Pa,
424                 G4ThreeVector V, G4ThreeVector    462                 G4ThreeVector V, G4ThreeVector P)
425     // Rotates a 3-vector P to close momentum     463     // Rotates a 3-vector P to close momentum triangle Pa + V + P = 0
426 {                                                 464 {
427   G4ThreeVector U = Pa.unit();                    465   G4ThreeVector U = Pa.unit();
428                                                   466   
429   G4double Alpha1 = U * V;                        467   G4double Alpha1 = U * V;
430                                                   468   
431   G4double Alpha2 = std::sqrt(V.mag2() - Alpha    469   G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
432                                                   470 
433   G4ThreeVector N = (1./Alpha2)*U.cross(V);       471   G4ThreeVector N = (1./Alpha2)*U.cross(V);
434                                                   472   
435   G4ThreeVector RotatedMomentum(                  473   G4ThreeVector RotatedMomentum(
436         ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.    474         ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
437         ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.    475         ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
438         ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.    476         ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
439         );                                        477         );
440   return RotatedMomentum;                         478   return RotatedMomentum;
441 }                                                 479 }
442                                                   480 
                                                   >> 481 
                                                   >> 482 
                                                   >> 483 
                                                   >> 484 
                                                   >> 485 G4ThreeVector G4StatMFChannel::IsotropicVector(const G4double Magnitude)
                                                   >> 486     // Samples a isotropic random vector with a magnitud given by Magnitude.
                                                   >> 487     // By default Magnitude = 1
                                                   >> 488 {
                                                   >> 489     G4double CosTheta = 1.0 - 2.0*G4UniformRand();
                                                   >> 490     G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
                                                   >> 491     G4double Phi = twopi*G4UniformRand();
                                                   >> 492     G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
                                                   >> 493        Magnitude*std::cos(Phi)*CosTheta,
                                                   >> 494        Magnitude*std::sin(Phi));
                                                   >> 495     return Vector;
                                                   >> 496 }
443                                                   497