| Geant4 Cross Reference |
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 // Class Description 26 // Class Description
27 // Cross-section data set for a high precision 27 // Cross-section data set for a high precision (based on JENDL_HE evaluated data
28 // libraries) description of elastic scatterin 28 // libraries) description of elastic scattering 20 MeV ~ 3 GeV;
29 // Class Description - End 29 // Class Description - End
30 30
31 // 15-Nov-06 First Implementation is done by T 31 // 15-Nov-06 First Implementation is done by T. Koi (SLAC/SCCS)
32 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 // 33 //
34 #include "G4ParticleHPJENDLHEData.hh" 34 #include "G4ParticleHPJENDLHEData.hh"
35 << 35 #include "G4SystemOfUnits.hh"
>> 36 #include "G4LPhysicsFreeVector.hh"
36 #include "G4ElementTable.hh" 37 #include "G4ElementTable.hh"
37 #include "G4ParticleHPData.hh" 38 #include "G4ParticleHPData.hh"
38 #include "G4PhysicsFreeVector.hh" <<
39 #include "G4Pow.hh" <<
40 #include "G4SystemOfUnits.hh" <<
41 39
42 G4bool G4ParticleHPJENDLHEData::IsApplicable(c << 40 G4bool G4ParticleHPJENDLHEData::IsApplicable(const G4DynamicParticle*aP, const G4Element* anE)
43 { 41 {
44 G4bool result = true; <<
45 G4double eKin = aP->GetKineticEnergy(); <<
46 // if(eKin>20*MeV||aP->GetDefinition()!=G4Ne <<
47 if (eKin < 20 * MeV || 3 * GeV < eKin || aP- <<
48 result = false; <<
49 } <<
50 // Element Check <<
51 else if (!(vElement[anE->GetIndex()])) <<
52 result = false; <<
53 42
54 return result; << 43 G4bool result = true;
>> 44 G4double eKin = aP->GetKineticEnergy();
>> 45 //if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
>> 46 if ( eKin < 20*MeV || 3*GeV < eKin || aP->GetDefinition()!=G4Neutron::Neutron() )
>> 47 {
>> 48 result = false;
>> 49 }
>> 50 // Element Check
>> 51 else if ( !(vElement[ anE->GetIndex() ]) ) result = false;
>> 52
>> 53 return result;
>> 54
55 } 55 }
56 56
>> 57
>> 58
57 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa 59 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEData()
58 { 60 {
59 for (auto& itZ : mIsotope) { << 61 for ( std::map< G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itZ = mIsotope.begin();
60 std::map<G4int, G4PhysicsVector*>* pointer << 62 itZ != mIsotope.end(); ++itZ ) {
61 if (pointer_map != nullptr) { << 63 std::map< G4int , G4PhysicsVector* >* pointer_map = itZ->second;
62 for (auto& itA : *pointer_map) { << 64 if ( pointer_map ) {
63 G4PhysicsVector* pointerPhysicsVector << 65 for ( std::map< G4int , G4PhysicsVector* >::iterator itA = pointer_map->begin();
64 if (pointerPhysicsVector != nullptr) { << 66 itA != pointer_map->end() ; ++itA ) {
>> 67 G4PhysicsVector* pointerPhysicsVector = itA->second;
>> 68 if ( pointerPhysicsVector ) {
65 delete pointerPhysicsVector; 69 delete pointerPhysicsVector;
66 itA.second = NULL; << 70 itA->second = NULL;
67 } 71 }
68 } 72 }
69 delete pointer_map; 73 delete pointer_map;
70 itZ.second = NULL; << 74 itZ->second = NULL;
71 } 75 }
72 } 76 }
73 mIsotope.clear(); 77 mIsotope.clear();
74 } 78 }
>> 79
75 80
76 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa << 81
77 : G4VCrossSectionDataSet("JENDLHE" + reactio << 82 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEData( G4String reaction , G4ParticleDefinition* pd )
>> 83 :G4VCrossSectionDataSet( "JENDLHE"+reaction+"CrossSection" )
78 { 84 {
79 reactionName = reaction; << 85 reactionName = reaction;
80 BuildPhysicsTable(*pd); << 86 BuildPhysicsTable( *pd );
81 } 87 }
82 88
83 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHED <<
84 89
85 void G4ParticleHPJENDLHEData::BuildPhysicsTabl << 90
86 { << 91 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHEData()
87 particleName = aP.GetParticleName(); << 92 {
88 << 93 ;
89 const G4String& baseName = G4FindDataDir("G4 << 94 //delete theCrossSections;
90 const G4String& dirName = baseName + "/JENDL <<
91 const G4String& aFSType = "/CrossSection/"; <<
92 G4ParticleHPNames theNames; <<
93 <<
94 G4String filename; <<
95 <<
96 // Create JENDL_HE data <<
97 // Create map element or isotope <<
98 <<
99 std::size_t numberOfElements = G4Element::Ge <<
100 <<
101 // make a PhysicsVector for each element <<
102 <<
103 auto theElementTable = G4Element::GetElement <<
104 vElement.clear(); <<
105 vElement.resize(numberOfElements); <<
106 for (std::size_t i = 0; i < numberOfElements <<
107 G4Element* theElement = (*theElementTable) <<
108 vElement[i] = false; <<
109 <<
110 // isotope <<
111 auto nIso = (G4int)(*theElementTable)[i]-> <<
112 auto Z = (G4int)(*theElementTable)[i]->Get <<
113 for (G4int i1 = 0; i1 < nIso; ++i1) { <<
114 G4int A = theElement->GetIsotope(i1)->Ge <<
115 <<
116 if (isThisNewIsotope(Z, A)) { <<
117 std::stringstream ss; <<
118 ss << dirName << aFSType << Z << "_" << A << <<
119 filename = ss.str(); <<
120 std::fstream file; <<
121 file.open(filename, std::fstream::in); <<
122 G4int dummy; <<
123 file >> dummy; <<
124 if (file.good()) { <<
125 vElement[i] = true; <<
126 <<
127 // read the file <<
128 G4PhysicsVector* aPhysVec = readAFile(&fil <<
129 registAPhysicsVector(Z, A, aPhysVec); <<
130 } <<
131 file.close(); <<
132 } <<
133 } <<
134 } <<
135 } 95 }
>> 96
136 97
137 void G4ParticleHPJENDLHEData::DumpPhysicsTable <<
138 {} <<
139 98
140 G4double G4ParticleHPJENDLHEData::GetCrossSect << 99 void G4ParticleHPJENDLHEData::BuildPhysicsTable( const G4ParticleDefinition& aP )
141 <<
142 { 100 {
143 // Primary energy >20MeV << 101
144 // Thus not taking into account of Doppler b << 102 // if ( &aP != G4Neutron::Neutron() )
145 // also not taking into account of Target th << 103 // throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
146 << 104 particleName = aP.GetParticleName();
147 G4double result = 0; << 105
148 << 106 G4String baseName = getenv( "G4NEUTRONHPDATA" );
149 G4double ek = aP->GetKineticEnergy(); << 107 G4String dirName = baseName+"/JENDL_HE/"+particleName+"/"+reactionName ;
150 << 108 G4String aFSType = "/CrossSection/";
151 auto nIso = (G4int)anE->GetNumberOfIsotopes( << 109 G4ParticleHPNames theNames;
152 auto Z = (G4int)anE->GetZ(); << 110
153 for (G4int i1 = 0; i1 < nIso; ++i1) { << 111 G4String filename;
154 G4int A = anE->GetIsotope(i1)->GetN(); << 112
155 G4double frac = anE->GetRelativeAbundanceV << 113 // Create JENDL_HE data
156 // This case does NOT request "*perCent". << 114 // Create map element or isotope
157 result += frac * getXSfromThisIsotope(Z, A << 115
158 } << 116 size_t numberOfElements = G4Element::GetNumberOfElements();
159 return result; << 117 //theCrossSections = new G4PhysicsTable( numberOfElements );
>> 118
>> 119 // make a PhysicsVector for each element
>> 120
>> 121 static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
>> 122 vElement.clear();
>> 123 vElement.resize( numberOfElements );
>> 124 for ( size_t i = 0; i < numberOfElements; ++i )
>> 125 {
>> 126
>> 127 G4Element* theElement = (*theElementTable)[i];
>> 128 vElement[i] = false;
>> 129
>> 130 // isotope
>> 131 G4int nIso = (*theElementTable)[i]->GetNumberOfIsotopes();
>> 132 G4int Z = static_cast<G4int> ((*theElementTable)[i]->GetZ());
>> 133 if ( nIso!=0 )
>> 134 {
>> 135 G4bool found_at_least_one = false;
>> 136 for ( G4int i1 = 0; i1 < nIso; i1++ )
>> 137 {
>> 138 G4int A = theElement->GetIsotope(i1)->GetN();
>> 139
>> 140 if ( isThisNewIsotope( Z , A ) )
>> 141 {
>> 142
>> 143 std::stringstream ss;
>> 144 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
>> 145 filename = ss.str();
>> 146 std::fstream file;
>> 147 file.open ( filename , std::fstream::in );
>> 148 G4int dummy;
>> 149 file >> dummy;
>> 150 if ( file.good() )
>> 151 {
>> 152
>> 153 //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
>> 154 found_at_least_one = true;
>> 155
>> 156 // read the file
>> 157 G4PhysicsVector* aPhysVec = readAFile ( &file );
>> 158
>> 159 //Regist
>> 160
>> 161 registAPhysicsVector( Z , A , aPhysVec );
>> 162
>> 163 }
>> 164 else
>> 165 {
>> 166 //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
>> 167 }
>> 168
>> 169 file.close();
>> 170
>> 171 }
>> 172 else
>> 173 {
>> 174 found_at_least_one = TRUE;
>> 175 }
>> 176 }
>> 177
>> 178 if ( found_at_least_one ) vElement[i] = true;
>> 179
>> 180 }
>> 181 else
>> 182 {
>> 183 G4StableIsotopes theStableOnes;
>> 184 G4int first = theStableOnes.GetFirstIsotope( Z );
>> 185 G4bool found_at_least_one = FALSE;
>> 186 for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(theElement->GetZ() ) ); i1++)
>> 187 {
>> 188 G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
>> 189
>> 190 if ( isThisNewIsotope( Z , A ) )
>> 191 {
>> 192
>> 193 std::stringstream ss;
>> 194 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
>> 195 filename = ss.str();
>> 196
>> 197 std::fstream file;
>> 198 file.open ( filename , std::fstream::in );
>> 199 G4int dummy;
>> 200 file >> dummy;
>> 201 if ( file.good() )
>> 202 {
>> 203 //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
>> 204 found_at_least_one = TRUE;
>> 205 //Read the file
>> 206
>> 207 G4PhysicsVector* aPhysVec = readAFile ( &file );
>> 208
>> 209 //Regist the PhysicsVector
>> 210 registAPhysicsVector( Z , A , aPhysVec );
>> 211
>> 212 }
>> 213 else
>> 214 {
>> 215 //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
>> 216 }
>> 217
>> 218 file.close();
>> 219 }
>> 220 else
>> 221 {
>> 222 found_at_least_one = TRUE;
>> 223 }
>> 224 }
>> 225
>> 226 if ( found_at_least_one ) vElement[i] = true;
>> 227
>> 228 }
>> 229
>> 230 }
>> 231
160 } 232 }
161 233
162 G4PhysicsVector* G4ParticleHPJENDLHEData::read << 234
>> 235
>> 236 void G4ParticleHPJENDLHEData::DumpPhysicsTable(const G4ParticleDefinition& aP)
163 { 237 {
164 G4int dummy; << 238 if(&aP!=G4Neutron::Neutron())
165 G4int len; << 239 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
166 *file >> dummy; << 240 // G4cout << "G4ParticleHPJENDLHEData::DumpPhysicsTable still to be implemented"<<G4endl;
167 *file >> len; << 241 }
168 242
169 std::vector<G4double> v_e; <<
170 std::vector<G4double> v_xs; <<
171 243
172 for (G4int i = 0; i < len; ++i) { <<
173 G4double e; <<
174 G4double xs; <<
175 244
176 *file >> e; << 245 G4double G4ParticleHPJENDLHEData::
177 *file >> xs; << 246 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double )
178 // data are written in eV and barn. << 247 // aTemp
179 v_e.push_back(e * eV); << 248 {
180 v_xs.push_back(xs * barn); <<
181 } <<
182 249
183 auto aPhysVec = new G4PhysicsFreeVector(stat << 250 // Primary energy >20MeV
>> 251 // Thus
>> 252 // Not take account of Doppler broadening
>> 253 // also
>> 254 // Not take account of Target thermal motions
>> 255
>> 256 G4double result = 0;
>> 257
>> 258 G4double ek = aP->GetKineticEnergy();
>> 259
>> 260 G4int nIso = anE->GetNumberOfIsotopes();
>> 261 G4int Z = static_cast<G4int> ( anE->GetZ() );
>> 262 if ( nIso!=0 )
>> 263 {
>> 264 for ( G4int i1 = 0; i1 < nIso; i1++ )
>> 265 {
184 266
185 for (G4int i = 0; i < len; ++i) { << 267 G4int A = anE->GetIsotope(i1)->GetN();
186 aPhysVec->PutValues(static_cast<std::size_ << 268 G4double frac = anE->GetRelativeAbundanceVector()[ i1 ]; // This case do NOT request "*perCent".
187 } << 269
>> 270 result += frac * getXSfromThisIsotope( Z , A , ek );
>> 271 //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl;
>> 272
>> 273 }
>> 274 }
>> 275 else
>> 276 {
>> 277
>> 278 G4StableIsotopes theStableOnes;
>> 279 G4int first = theStableOnes.GetFirstIsotope( Z );
>> 280 for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(anE->GetZ() ) ); i1++)
>> 281 {
>> 282
>> 283 G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
>> 284 G4double frac = theStableOnes.GetAbundance( first+i1 )*perCent; // This case request "*perCent".
>> 285
>> 286 result += frac * getXSfromThisIsotope( Z , A , ek );
>> 287 //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl;
>> 288
>> 289 }
>> 290 }
>> 291 return result;
188 292
189 return aPhysVec; <<
190 } 293 }
191 294
192 G4bool G4ParticleHPJENDLHEData::isThisInMap(G4 << 295
>> 296
>> 297 G4PhysicsVector* G4ParticleHPJENDLHEData::readAFile ( std::fstream* file )
193 { 298 {
194 if (mIsotope.find(z) == mIsotope.end()) retu << 299
195 if (mIsotope.find(z)->second->find(a) == mIs << 300 G4int dummy;
196 return true; << 301 G4int len;
>> 302 *file >> dummy;
>> 303 *file >> len;
>> 304
>> 305 std::vector< G4double > v_e;
>> 306 std::vector< G4double > v_xs;
>> 307
>> 308 for ( G4int i = 0 ; i < len ; i++ )
>> 309 {
>> 310 G4double e;
>> 311 G4double xs;
>> 312
>> 313 *file >> e;
>> 314 *file >> xs;
>> 315 // data are written in eV and barn.
>> 316 v_e.push_back( e*eV );
>> 317 v_xs.push_back( xs*barn );
>> 318 }
>> 319
>> 320 G4LPhysicsFreeVector* aPhysVec = new G4LPhysicsFreeVector( static_cast< size_t >( len ) , v_e.front() , v_e.back() );
>> 321
>> 322 for ( G4int i = 0 ; i < len ; i++ )
>> 323 {
>> 324 aPhysVec->PutValues( static_cast< size_t >( i ) , v_e[ i ] , v_xs[ i ] );
>> 325 }
>> 326
>> 327 return aPhysVec;
197 } 328 }
198 329
199 void G4ParticleHPJENDLHEData::registAPhysicsVe << 330
>> 331
>> 332 G4bool G4ParticleHPJENDLHEData::isThisInMap( G4int z , G4int a )
200 { 333 {
201 std::pair<G4int, G4PhysicsVector*> aPair = s << 334 if ( mIsotope.find ( z ) == mIsotope.end() ) return false;
202 auto itm = mIsotope.find(Z); << 335 if ( mIsotope.find ( z ) -> second->find ( a ) == mIsotope.find ( z ) -> second->end() ) return false;
203 if (itm != mIsotope.cend()) { << 336 return true;
204 itm->second->insert(aPair); <<
205 } <<
206 else { <<
207 auto aMap = new std::map<G4int, G4PhysicsV <<
208 aMap->insert(aPair); <<
209 mIsotope.insert(std::pair<G4int, std::map< <<
210 } <<
211 } 337 }
212 338
213 G4double G4ParticleHPJENDLHEData::getXSfromThi << 339
>> 340
>> 341 void G4ParticleHPJENDLHEData::registAPhysicsVector( G4int Z , G4int A , G4PhysicsVector* aPhysVec )
214 { 342 {
215 G4double aXSection = 0.0; <<
216 343
217 G4PhysicsVector* aPhysVec; << 344 std::pair< G4int , G4PhysicsVector* > aPair = std::pair < G4int , G4PhysicsVector* > ( A , aPhysVec );
218 if (mIsotope.find(Z)->second->find(A) != mIs << 345
219 aPhysVec = mIsotope.find(Z)->second->find( << 346 std::map < G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itm;
220 aXSection = aPhysVec->Value(ek); << 347 itm = mIsotope.find ( Z );
221 } << 348 if ( itm != mIsotope.end() )
222 else { << 349 {
223 // Select closest one in the same Z << 350 itm->second->insert ( aPair );
224 G4int delta0 = 99; // no mean for 99 << 351 }
225 for (auto it = mIsotope.find(Z)->second->c << 352 else
226 { 353 {
227 G4int delta = std::abs(A - it->first); << 354 std::map< G4int , G4PhysicsVector* >* aMap = new std::map< G4int , G4PhysicsVector* >;
228 if (delta < delta0) delta0 = delta; << 355 aMap->insert ( aPair );
>> 356 mIsotope.insert( std::pair< G4int , std::map< G4int , G4PhysicsVector* >* > ( Z , aMap ) );
229 } 357 }
230 358
231 // Randomize of selection larger or smalle << 359 }
232 if (G4UniformRand() < 0.5) delta0 *= -1; <<
233 G4int A1 = A + delta0; <<
234 if (mIsotope.find(Z)->second->find(A1) != <<
235 aPhysVec = mIsotope.find(Z)->second->fin <<
236 } <<
237 else { <<
238 A1 = A - delta0; <<
239 aPhysVec = mIsotope.find(Z)->second->fin <<
240 } <<
241 360
242 aXSection = aPhysVec->Value(ek); <<
243 // X^(2/3) factor <<
244 aXSection *= G4Pow::GetInstance()->A23(1.0 <<
245 } <<
246 361
247 return aXSection; << 362
>> 363 G4double G4ParticleHPJENDLHEData::getXSfromThisIsotope( G4int Z , G4int A , G4double ek )
>> 364 {
>> 365
>> 366 G4double aXSection = 0.0;
>> 367 G4bool outOfRange;
>> 368
>> 369 G4PhysicsVector* aPhysVec;
>> 370 if ( mIsotope.find ( Z )->second->find ( A ) != mIsotope.find ( Z )->second->end() )
>> 371 {
>> 372
>> 373 aPhysVec = mIsotope.find ( Z )->second->find ( A )->second;
>> 374 aXSection = aPhysVec->GetValue( ek , outOfRange );
>> 375
>> 376 }
>> 377 else
>> 378 {
>> 379
>> 380 //Select closest one in the same Z
>> 381 std::map < G4int , G4PhysicsVector* >::iterator it;
>> 382 G4int delta0 = 99; // no mean for 99
>> 383 for ( it = mIsotope.find ( Z )->second->begin() ; it != mIsotope.find ( Z )->second->end() ; it++ )
>> 384 {
>> 385 G4int delta = std::abs( A - it->first );
>> 386 if ( delta < delta0 ) delta0 = delta;
>> 387 }
>> 388
>> 389 // Randomize of selection larger or smaller than A
>> 390 if ( G4UniformRand() < 0.5 ) delta0 *= -1;
>> 391 G4int A1 = A + delta0;
>> 392 if ( mIsotope.find ( Z )->second->find ( A1 ) != mIsotope.find ( Z )->second->end() )
>> 393 {
>> 394 aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
>> 395 }
>> 396 else
>> 397 {
>> 398 A1 = A - delta0;
>> 399 aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
>> 400 }
>> 401
>> 402 aXSection = aPhysVec->GetValue( ek , outOfRange );
>> 403 // X^(2/3) factor
>> 404 aXSection *= std::pow ( 1.0*A/ A1 , 2.0 / 3.0 );
>> 405
>> 406 }
>> 407
>> 408 return aXSection;
248 } 409 }
249 410