Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CompetitiveFission.cc
Go to the documentation of this file.
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 // $Id$
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara (Oct 1998)
31 //
32 // J. M. Quesada (March 2009). Bugs fixed:
33 // - Full relativistic calculation (Lorentz boosts)
34 // - Fission pairing energy is included in fragment excitation energies
35 // Now Energy and momentum are conserved in fission
36 
37 #include "G4CompetitiveFission.hh"
38 #include "G4PairingCorrection.hh"
39 #include "G4ParticleMomentum.hh"
40 #include "G4Pow.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 
45 {
46  theFissionBarrierPtr = new G4FissionBarrier;
47  MyOwnFissionBarrier = true;
48 
49  theFissionProbabilityPtr = new G4FissionProbability;
50  MyOwnFissionProbability = true;
51 
52  theLevelDensityPtr = new G4FissionLevelDensityParameter;
53  MyOwnLevelDensity = true;
54 
55  MaximalKineticEnergy = -1000.0*MeV;
56  FissionBarrier = 0.0;
57  FissionProbability = 0.0;
58  LevelDensityParameter = 0.0;
59 }
60 
62 {
63  if (MyOwnFissionBarrier) delete theFissionBarrierPtr;
64 
65  if (MyOwnFissionProbability) delete theFissionProbabilityPtr;
66 
67  if (MyOwnLevelDensity) delete theLevelDensityPtr;
68 }
69 
71 {
72  G4int anA = fragment->GetA_asInt();
73  G4int aZ = fragment->GetZ_asInt();
74  G4double ExEnergy = fragment->GetExcitationEnergy() -
76 
77 
78  // Saddle point excitation energy ---> A = 65
79  // Fission is excluded for A < 65
80  if (anA >= 65 && ExEnergy > 0.0) {
81  FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy);
82  MaximalKineticEnergy = ExEnergy - FissionBarrier;
83  LevelDensityParameter =
84  theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
85  FissionProbability =
86  theFissionProbabilityPtr->EmissionProbability(*fragment,MaximalKineticEnergy);
87  }
88  else {
89  MaximalKineticEnergy = -1000.0*MeV;
90  LevelDensityParameter = 0.0;
91  FissionProbability = 0.0;
92  }
93  return FissionProbability;
94 }
95 
97 {
98  // Nucleus data
99  // Atomic number of nucleus
100  G4int A = theNucleus.GetA_asInt();
101  // Charge of nucleus
102  G4int Z = theNucleus.GetZ_asInt();
103  // Excitation energy (in MeV)
104  G4double U = theNucleus.GetExcitationEnergy() -
106  // Check that U > 0
107  if (U <= 0.0) {
108  G4FragmentVector * theResult = new G4FragmentVector;
109  theResult->push_back(new G4Fragment(theNucleus));
110  return theResult;
111  }
112 
113  // Atomic Mass of Nucleus (in MeV)
114  G4double M = theNucleus.GetGroundStateMass();
115 
116  // Nucleus Momentum
117  G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
118 
119  // Calculate fission parameters
120  G4FissionParameters theParameters(A,Z,U,FissionBarrier);
121 
122  // First fragment
123  G4int A1 = 0;
124  G4int Z1 = 0;
125  G4double M1 = 0.0;
126 
127  // Second fragment
128  G4int A2 = 0;
129  G4int Z2 = 0;
130  G4double M2 = 0.0;
131 
132  G4double FragmentsExcitationEnergy = 0.0;
133  G4double FragmentsKineticEnergy = 0.0;
134 
135  //JMQ 04/03/09 It will be used latter to fix the bug in energy conservation
136  G4double FissionPairingEnergy=
138 
139  G4int Trials = 0;
140  do {
141 
142  // First fragment
143  A1 = FissionAtomicNumber(A,theParameters);
144  Z1 = FissionCharge(A,Z,A1);
146 
147  // Second Fragment
148  A2 = A - A1;
149  Z2 = Z - Z1;
150  if (A2 < 1 || Z2 < 0) {
151  throw G4HadronicException(__FILE__, __LINE__,
152  "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
153  }
155 
156  // Check that fragment masses are less or equal than total energy
157  if (M1 + M2 > theNucleusMomentum.e()) {
158  throw G4HadronicException(__FILE__, __LINE__,
159  "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
160  }
161  // Maximal Kinetic Energy (available energy for fragments)
162  G4double Tmax = M + U - M1 - M2;
163 
164  FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
165  A1, Z1,
166  A2, Z2,
167  U , Tmax,
168  theParameters);
169 
170  // Excitation Energy
171  // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
172  // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
173  // fragments carry the fission pairing energy in form of
174  //excitation energy
175 
176  FragmentsExcitationEnergy =
177  Tmax - FragmentsKineticEnergy+FissionPairingEnergy;
178 
179  } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
180 
181  if (FragmentsExcitationEnergy <= 0.0) {
182  throw G4HadronicException(__FILE__, __LINE__,
183  "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
184  }
185 
186  // while (FragmentsExcitationEnergy < 0 && Trials < 100);
187 
188  // Fragment 1
189  G4double U1 = FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
190  // Fragment 2
191  G4double U2 = FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
192 
193  //JMQ 04/03/09 Full relativistic calculation is performed
194  //
195  G4double Fragment1KineticEnergy=
196  (FragmentsKineticEnergy*(FragmentsKineticEnergy+2*(M2+U2)))
197  /(2*(M1+U1+M2+U2+FragmentsKineticEnergy));
198  G4ParticleMomentum Momentum1(IsotropicVector(std::sqrt(Fragment1KineticEnergy*(Fragment1KineticEnergy+2*(M1+U1)))));
199  G4ParticleMomentum Momentum2(-Momentum1);
200  G4LorentzVector FourMomentum1(Momentum1,std::sqrt(Momentum1.mag2()+(M1+U1)*(M1+U1)));
201  G4LorentzVector FourMomentum2(Momentum2,std::sqrt(Momentum2.mag2()+(M2+U2)*(M2+U2)));
202 
203  //JMQ 04/03/09 now we do Lorentz boosts (instead of Galileo boosts)
204  FourMomentum1.boost(theNucleusMomentum.boostVector());
205  FourMomentum2.boost(theNucleusMomentum.boostVector());
206 
208  // There was vioation of energy momentum conservation
209 
210  // G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
211  // ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
212 
213  //G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
214  // G4ParticleMomentum momentum2( -momentum1 );
215 
216  // Perform a Galileo boost for fragments
217  // momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
218  // momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
219 
220 
221  // Create 4-momentum for first fragment
222  // Warning!! Energy conservation is broken
223  //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
224  // G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
225 
226  // Create 4-momentum for second fragment
227  // Warning!! Energy conservation is broken
228  //JMQ 04/03/09 ...NOT ANY MORE!! BUGS FIXED: Energy and momentum are NOW conserved
229  // G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
230 
232 
233  // Create Fragments
234  G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
235  G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2);
236 
237  // Create Fragment Vector
238  G4FragmentVector * theResult = new G4FragmentVector;
239 
240  theResult->push_back(Fragment1);
241  theResult->push_back(Fragment2);
242 
243 #ifdef debug
244  CheckConservation(theNucleus,theResult);
245 #endif
246 
247  return theResult;
248 }
249 
250 G4int
251 G4CompetitiveFission::FissionAtomicNumber(G4int A,
252  const G4FissionParameters & theParam)
253  // Calculates the atomic number of a fission product
254 {
255 
256  // For Simplicity reading code
257  const G4double A1 = theParam.GetA1();
258  const G4double A2 = theParam.GetA2();
259  const G4double As = theParam.GetAs();
260  // const G4double Sigma1 = theParam.GetSigma1();
261  const G4double Sigma2 = theParam.GetSigma2();
262  const G4double SigmaS = theParam.GetSigmaS();
263  const G4double w = theParam.GetW();
264 
265  // G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
266  // std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
267 
268  // G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
269 
270  G4double C2A = A2 + 3.72*Sigma2;
271  G4double C2S = As + 3.72*SigmaS;
272 
273  G4double C2 = 0.0;
274  if (w > 1000.0 ) C2 = C2S;
275  else if (w < 0.001) C2 = C2A;
276  else C2 = std::max(C2A,C2S);
277 
278  G4double C1 = A-C2;
279  if (C1 < 30.0) {
280  C2 = A-30.0;
281  C1 = 30.0;
282  }
283 
284  G4double Am1 = (As + A1)/2.0;
285  G4double Am2 = (A1 + A2)/2.0;
286 
287  // Get Mass distributions as sum of symmetric and asymmetric Gasussians
288  G4double Mass1 = MassDistribution(As,A,theParam);
289  G4double Mass2 = MassDistribution(Am1,A,theParam);
290  G4double Mass3 = MassDistribution(A1,A,theParam);
291  G4double Mass4 = MassDistribution(Am2,A,theParam);
292  G4double Mass5 = MassDistribution(A2,A,theParam);
293  // get maximal value among Mass1,...,Mass5
294  G4double MassMax = Mass1;
295  if (Mass2 > MassMax) MassMax = Mass2;
296  if (Mass3 > MassMax) MassMax = Mass3;
297  if (Mass4 > MassMax) MassMax = Mass4;
298  if (Mass5 > MassMax) MassMax = Mass5;
299 
300  // Sample a fragment mass number, which lies between C1 and C2
301  G4double xm;
302  G4double Pm;
303  do {
304  xm = C1+G4UniformRand()*(C2-C1);
305  Pm = MassDistribution(xm,A,theParam);
306  } while (MassMax*G4UniformRand() > Pm);
307  G4int ires = G4lrint(xm);
308 
309  return ires;
310 }
311 
312 G4double
313 G4CompetitiveFission::MassDistribution(G4double x, G4double A,
314  const G4FissionParameters & theParam)
315  // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
316  // which consist of symmetric and asymmetric sum of gaussians components.
317 {
318  G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
319  (theParam.GetSigmaS()*theParam.GetSigmaS()));
320 
321  G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
322  (theParam.GetSigma2()*theParam.GetSigma2())) +
323  std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
324  (theParam.GetSigma2()*theParam.GetSigma2())) +
325  0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
326  (theParam.GetSigma1()*theParam.GetSigma1())) +
327  0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
328  (theParam.GetSigma1()*theParam.GetSigma1()));
329 
330  if (theParam.GetW() > 1000) return Xsym;
331  else if (theParam.GetW() < 0.001) return Xasym;
332  else return theParam.GetW()*Xsym+Xasym;
333 }
334 
335 G4int G4CompetitiveFission::FissionCharge(G4double A, G4double Z,
336  G4double Af)
337  // Calculates the charge of a fission product for a given atomic number Af
338 {
339  const G4double sigma = 0.6;
340  G4double DeltaZ = 0.0;
341  if (Af >= 134.0) DeltaZ = -0.45; // 134 <= Af
342  else if (Af <= (A-134.0)) DeltaZ = 0.45; // Af <= (A-134)
343  else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0)); // (A-134) < Af < 134
344 
345  G4double Zmean = (Af/A)*Z + DeltaZ;
346 
347  G4double theZ;
348  do {
349  theZ = G4RandGauss::shoot(Zmean,sigma);
350  } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af);
351  // return static_cast<G4int>(theZ+0.5);
352  return static_cast<G4int>(theZ+0.5);
353 }
354 
355 G4double
356 G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z,
357  G4double Af1, G4double /*Zf1*/,
358  G4double Af2, G4double /*Zf2*/,
359  G4double /*U*/, G4double Tmax,
360  const G4FissionParameters & theParam)
361  // Gives the kinetic energy of fission products
362 {
363  // Find maximal value of A for fragments
364  G4double AfMax = std::max(Af1,Af2);
365  if (AfMax < (A/2.0)) AfMax = A - AfMax;
366 
367  // Weights for symmetric and asymmetric components
368  G4double Pas;
369  if (theParam.GetW() > 1000) Pas = 0.0;
370  else {
371  G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
372  (theParam.GetSigma1()*theParam.GetSigma1()));
373 
374  G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
375  (theParam.GetSigma2()*theParam.GetSigma2()));
376 
377  Pas = P1+P2;
378  }
379 
380  G4double Ps;
381  if (theParam.GetW() < 0.001) Ps = 0.0;
382  else {
383  Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
384  (theParam.GetSigmaS()*theParam.GetSigmaS()));
385  }
386  G4double Psy = Ps/(Pas+Ps);
387 
388  // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
389  G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
390  G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
391  G4double Xas = PPas / (PPas+PPsy);
392  G4double Xsy = PPsy / (PPas+PPsy);
393 
394  // Average kinetic energy for symmetric and asymmetric components
395  G4double Eaverage = 0.1071*MeV*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2*MeV;
396 
397 
398  // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
399  G4double TaverageAfMax;
400  G4double ESigma;
401  // Select randomly fission mode (symmetric or asymmetric)
402  if (G4UniformRand() > Psy) { // Asymmetric Mode
403  G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
404  G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
405  G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
406  G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
407  // scale factor
408  G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
409  theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
410  // Compute average kinetic energy for fragment with AfMax
411  TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax);
412  ESigma = 10.0*MeV; // MeV
413 
414  } else { // Symmetric Mode
415  G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
416  // scale factor
417  G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
418  // Compute average kinetic energy for fragment with AfMax
419  TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax);
420  ESigma = 8.0*MeV;
421  }
422 
423 
424  // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
425  G4double KineticEnergy;
426  G4int i = 0;
427  do {
428  KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
429  if (i++ > 100) return Eaverage;
430  } while (KineticEnergy < Eaverage-3.72*ESigma ||
431  KineticEnergy > Eaverage+3.72*ESigma ||
432  KineticEnergy > Tmax);
433 
434  return KineticEnergy;
435 }
436 
437 G4double G4CompetitiveFission::AsymmetricRatio(G4int A, G4double A11)
438 {
439  const G4double B1 = 23.5;
440  const G4double A00 = 134.0;
441  return Ratio(G4double(A),A11,B1,A00);
442 }
443 
444 G4double G4CompetitiveFission::SymmetricRatio(G4int A, G4double A11)
445 {
446  const G4double B1 = 5.32;
447  const G4double A00 = A/2.0;
448  return Ratio(G4double(A),A11,B1,A00);
449 }
450 
451 G4double G4CompetitiveFission::Ratio(G4double A, G4double A11,
452  G4double B1, G4double A00)
453 {
454  if (A == 0.0) {
455  throw G4HadronicException(__FILE__, __LINE__,
456  "G4CompetitiveFission::Ratio: A == 0!");
457  }
458  if (A11 >= A/2.0 && A11 <= (A00+10.0)) {
459  return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
460  } else {
461  return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
462  }
463 }
464 
465 G4ThreeVector G4CompetitiveFission::IsotropicVector(const G4double Magnitude)
466  // Samples a isotropic random vectorwith a magnitud given by Magnitude.
467  // By default Magnitude = 1.0
468 {
469  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
470  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
471  G4double Phi = twopi*G4UniformRand();
472  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
473  Magnitude*std::sin(Phi)*SinTheta,
474  Magnitude*CosTheta);
475  return Vector;
476 }
477 
478 #ifdef debug
479 void G4CompetitiveFission::CheckConservation(const G4Fragment & theInitialState,
480  G4FragmentVector * Result) const
481 {
482  G4double ProductsEnergy =0;
483  G4ThreeVector ProductsMomentum;
484  G4int ProductsA = 0;
485  G4int ProductsZ = 0;
486  G4FragmentVector::iterator h;
487  for (h = Result->begin(); h != Result->end(); h++) {
488  G4LorentzVector tmp = (*h)->GetMomentum();
489  ProductsEnergy += tmp.e();
490  ProductsMomentum += tmp.vect();
491  ProductsA += (*h)->GetA_asInt();
492  ProductsZ += (*h)->GetZ_asInt();
493  }
494 
495  if (ProductsA != theInitialState.GetA_asInt()) {
496  G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
497  G4cout << "G4CompetitiveFission.cc: Barionic Number Conservation test for fission fragments"
498  << G4endl;
499  G4cout << "Initial A = " << theInitialState.GetA_asInt()
500  << " Fragments A = " << ProductsA << " Diference --> "
501  << theInitialState.GetA_asInt() - ProductsA << G4endl;
502  }
503  if (ProductsZ != theInitialState.GetZ_asInt()) {
504  G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
505  G4cout << "G4CompetitiveFission.cc: Charge Conservation test for fission fragments"
506  << G4endl;
507  G4cout << "Initial Z = " << theInitialState.GetZ_asInt()
508  << " Fragments Z = " << ProductsZ << " Diference --> "
509  << theInitialState.GetZ() - ProductsZ << G4endl;
510  }
511  if (std::fabs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
512  G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
513  G4cout << "G4CompetitiveFission.cc: Energy Conservation test for fission fragments"
514  << G4endl;
515  G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
516  << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
517  << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
518  }
519  if (std::fabs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV ||
520  std::fabs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
521  std::fabs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
522  G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
523  G4cout << "G4CompetitiveFission.cc: Momentum Conservation test for fission fragments"
524  << G4endl;
525  G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
526  << " Fragments P = " << ProductsMomentum << " MeV Diference --> "
527  << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
528  }
529  return;
530 }
531 #endif
532 
533 
534 
535