Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DPMJET2_5CrossSection.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 // * *
21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 19770/06/NL/JD (Technology Research Programme). *
29 // * *
30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Software license. *
34 // ********************************************************************
35 //
38 //
39 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40 //
41 // MODULE: G4DPMJET2_5CrossSection.cc
42 //
43 // Version: 0.A
44 // Date: 02/04/08
45 // Author: P R Truscott
46 // Organisation: QinetiQ Ltd, UK
47 // Customer: ESA/ESTEC, NOORDWIJK
48 // Contract: 19770/06/NL/JD
49 //
50 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 //
53 #ifdef G4_USE_DPMJET
54 
55 
56 #include "G4DPMJET2_5CrossSection.hh"
57 #include "G4ParticleTable.hh"
58 #include "G4DynamicParticle.hh"
59 #include "G4IonTable.hh"
60 
61 #include "G4HadronicException.hh"
62 #include "G4StableIsotopes.hh"
63 #include "G4HadTmpUtil.hh"
64 
65 #include "globals.hh"
66 
67 #include <iomanip>
68 #include <fstream>
69 #include <sstream>
70 
71 #include "G4DynamicParticle.hh"
72 
73 using namespace std;
74 
76 //
77 G4DPMJET2_5CrossSection::G4DPMJET2_5CrossSection ():
78  upperLimit ( 1000.0 * TeV ), lowerLimit ( 5.0 * GeV ), maxA(240)
79 {
80  theCrossSectionIndex.clear();
81  Initialise();
82 //
83 //
84 // vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
85 // This next bit is provisional, stating that this cross-section estimator
86 // is applicable to hydrogen targets. However, the cross-section will be
87 // set to zero.
88 //
89  ATmin = 1;
90 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
91 //
92 }
94 //
95 G4DPMJET2_5CrossSection::~G4DPMJET2_5CrossSection ()
96 {
97 //
98 // Go through the list of cross-section fit parameters and delete the arrays.
99 //
100  G4cout << "G4DPMJET2_5CrossSection::~G4DPMJET2_5CrossSection" << G4endl;
101  G4cout << "Size: " << theCrossSectionIndex.size() << G4endl;
102  /*
103  if(theCrossSectionIndex.size() > 0) {
104 
105  G4DPMJET2_5CrossSectionIndex::iterator it;
106  for (it=theCrossSectionIndex.begin(); it!=theCrossSectionIndex.end(); ++it)
107  {
108  G4DPMJET2_5CrossSectionParamSet *ptr = it->second;
109  for (G4DPMJET2_5CrossSectionParamSet *ptr1=ptr; ptr1<ptr+maxA; ptr1++)
110  { delete ptr1; }
111  }
112  }
113  */
114 }
116 //
117 G4bool G4DPMJET2_5CrossSection::IsApplicable
118  (const G4DynamicParticle* theProjectile, const G4Element* theTarget)
119 {
120  static G4StableIsotopes theDefaultIsotopes; // natural abundances and
121  // stable isotopes
122 //
123 //
124 // Determine first if the user has defined her own isotopic composition for the
125 // element.
126 //
127  G4int nIso = theTarget->GetNumberOfIsotopes();
128  G4bool result = true;
129 
130  if (nIso) {
131 //
132 //
133 // Determine whether we have the necessary data loaded for the user-defined
134 // cross-section.
135 //
136  G4IsotopeVector* isoVector = theTarget->GetIsotopeVector();
137  G4double ZZ = 0.0;
138  G4double AA = 0.0;
139  G4int i = 0;
140 
141  do {
142  ZZ = G4double( (*isoVector)[i]->GetZ() );
143  AA = G4double( (*isoVector)[i]->GetN() );
144  result = IsZAApplicable(theProjectile, ZZ, AA);
145  } while (result && ++i < nIso);
146  } else {
147 //
148 //
149 // Determine whether we have the necessary data loaded for the natural
150 // abundance composition of the element.
151 //
152  G4int ZZ = G4lrint(theTarget->GetZ());
153  nIso = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
154  G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
155  G4double AA = 0.0;
156  G4int i = 0;
157  do {
158  AA = G4double( theDefaultIsotopes.GetIsotopeNucleonCount(index+i) );
159  result = IsZAApplicable(theProjectile, G4double(ZZ), AA);
160  } while (result && ++i < nIso);
161  }
162  G4cout << "G4DPMJET2_5CrossSection::IsApplicable E(GeV)= "
163  << theProjectile->GetKineticEnergy()/GeV << " off "
164  << theTarget->GetName() << " - " << result << G4endl;
165  return result;
166 }
168 //
169 G4bool G4DPMJET2_5CrossSection::IsZAApplicable
170  (const G4DynamicParticle* theProjectile, G4double , G4double AA)
171 {
172  const G4int AT = G4lrint(AA);
173  const G4int AP = G4lrint(theProjectile->GetDefinition()->GetBaryonNumber());
174  G4double EPN = theProjectile->GetKineticEnergy()/
175  theProjectile->GetDefinition()->GetBaryonNumber();
176  G4bool result = EPN >= lowerLimit && EPN <= upperLimit &&
177  AT >= ATmin && AT <= ATmax &&
178  AP >= APmin && AP <= APmax;
179  return result;
180 }
182 //
183 G4double G4DPMJET2_5CrossSection::GetIsoZACrossSection
184  (const G4DynamicParticle* theProjectile, G4double ZZ, G4double AA,
185  G4double /*theTemperature*/)
186 {
187 //
188 // Initialise the result.
189  G4double result = 0.0;
190 //
191 //
192 // Get details of the projectile and target (nucleon number, atomic number,
193 // kinetic enery and energy/nucleon.
194 //
195  const G4int AT = G4lrint(AA);
196  G4int AP = G4lrint(theProjectile->GetDefinition()->GetBaryonNumber());
197  const G4double TP = theProjectile->GetKineticEnergy();
198  G4double EPN = TP / AP;
199 
200  if (AT < ATmin || AT > ATmax || AP < APmin || AP > APmax ||
201  EPN < lowerLimit || EPN > upperLimit)
202  {
203  G4cout <<G4endl;
204  G4cout <<"ERROR IN G4DPMJET2_5CrossSection::GetIsoZACrossSection" <<G4endl;
205  G4cout <<"ATTEMPT TO USE CROSS-SECTION OUTSIDE OF RANGE" <<G4endl;
206  G4cout <<"NUCLEON NUMBER OF PROJECTILE = " <<AP <<G4endl;
207  G4cout <<"NUCLEON NUMBER OF TARGET = " <<AT <<G4endl;
208  G4cout <<"ENERGY PER NUCLEON = " <<EPN*MeV <<G4endl;
209  G4cout <<"ACCEPTABLE RANGE FOR AP = " <<APmin
210  <<" TO " <<APmax <<G4endl;
211  G4cout <<"ACCEPTABLE RANGE FOR AT = " <<ATmin
212  <<" TO " <<ATmax <<G4endl;
213  G4cout <<"ACCEPTABLE RANGE FOR ENERGY = " <<lowerLimit
214  <<" MeV/n TO " <<upperLimit
215  <<" MeV/n" <<G4endl;
216  G4cout <<G4endl;
217  return result;
218  }
219 //
220 //
221 // vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
222 // This next bit is provisional, stating that this cross-section hydrogen
223 // targets is zero.
224 //
225  if ( AT == 1 ) return 0.0;
226 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
227 //
228 //
229 //
230 // Results are parameterised as a function of the natural logarithm of the
231 // centre of mass energy of the projectile and target system.
232 //
233  G4double sigma = 0.0;
234  G4double mT = G4ParticleTable::GetParticleTable()
235  ->GetIonTable()
236  ->GetIonMass(static_cast<G4int>(ZZ), static_cast<G4int>(AT));
237  G4double EP = theProjectile->GetTotalEnergy();
238  G4double mP = EP - TP;
239  G4double lnECM = std::log(std::sqrt(mP*mP + mT*mT + 2.0*mT*EP));
240  G4DPMJET2_5CrossSectionIndex::iterator it = theCrossSectionIndex.find(AT);
241  if (it != theCrossSectionIndex.end())
242  {
243  G4DPMJET2_5CrossSectionParamSet *ptr = (it->second) + AP;
244  G4double cc0 = (*ptr)[0];
245  G4double cc1 = (*ptr)[1];
246  G4double cc2 = (*ptr)[2];
247  sigma = cc0 + cc1*lnECM + cc2*lnECM*lnECM;
248  sigma = sigma * millibarn;
249  if (verboseLevel >= 2) {
250  G4cout <<"***************************************************************"
251  <<G4endl;
252  G4cout <<"G4DPMJET2_5CrossSection::GetIsoZACrossSection" <<G4endl;
253  G4cout <<"PROJECTILE = "
254  <<theProjectile->GetDefinition()->GetParticleName() <<G4endl;
255  G4cout <<"TARGET (A,Z) = (" <<AA <<"," <<ZZ <<")" <<G4endl;
256  G4cout <<"K. ENERGY/NUC = " <<EPN/MeV <<" MeV/n" <<G4endl;
257  G4cout <<"CROSS SECTION = " <<sigma/millibarn <<" MILLIBARNS" <<G4endl;
258  G4cout <<"***************************************************************"
259  <<G4endl;
260  }
261  }
262  else
263  {
264  G4cout <<G4endl;
265  G4cout <<"ERROR IN G4DPMJET2_5CrossSection::GetIsoZACrossSection" <<G4endl;
266  G4cout <<"NO CROSS-SECTION FIT DATA LOADED FOR AT = " <<AT <<G4endl;
267  G4cout <<G4endl;
268  }
269 
270  return sigma;
271 }
273 //
274 G4double G4DPMJET2_5CrossSection::GetCrossSection
275  (const G4DynamicParticle* theProjectile, const G4Element* theTarget,
276  G4double theTemperature)
277 {
278 // DumpPhysicsTable(*(theProjectile->GetDefinition()));
279  static G4StableIsotopes theDefaultIsotopes; // natural abundances and
280  // stable isotopes
281 //
282 //
283 // Determine first if the user has defined her own isotopic composition for the
284 // element.
285 //
286  G4int nIso = theTarget->GetNumberOfIsotopes();
287  G4double xsection = 0;
288 
289  if (nIso) {
290 //
291 //
292 // User-defined cross-section.
293 //
294  G4double sig = 0.0;
295  G4IsotopeVector* isoVector = theTarget->GetIsotopeVector();
296  G4double* abundVector = theTarget->GetRelativeAbundanceVector();
297  G4double ZZ = 0.0;
298  G4double AA = 0.0;
299 
300  for (G4int i = 0; i < nIso; i++) {
301  ZZ = G4double( (*isoVector)[i]->GetZ() );
302  AA = G4double( (*isoVector)[i]->GetN() );
303  sig = GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
304  xsection += sig*abundVector[i];
305  }
306 
307  } else {
308 //
309 //
310 // Calculation of cross-section for natural abundance.
311 //
312  G4double sig = 0.0;
313  G4int ZZ = G4lrint(theTarget->GetZ());
314  nIso = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
315  G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
316  G4double AA = 0.0;
317  G4double ab = 0.0;
318  for (G4int i = 0; i < nIso; i++) {
319  AA = G4double( theDefaultIsotopes.GetIsotopeNucleonCount(index+i) );
320  ab = theDefaultIsotopes.GetAbundance(index+i)/100.0;
321  sig = GetIsoZACrossSection(theProjectile, G4double(ZZ), AA,
322  theTemperature);
323  xsection += sig*ab;
324  }
325  }
326 
327  if (verboseLevel >= -2) {
328  G4int AP = G4lrint(theProjectile->GetDefinition()->GetBaryonNumber());
329  const G4double TP = theProjectile->GetKineticEnergy();
330  G4double EPN = TP / AP;
331  G4cout <<"***************************************************************"
332  <<G4endl;
333  G4cout <<"G4DPMJET2_5CrossSection::GetCrossSection" <<G4endl;
334  G4cout <<"PROJECTILE = "
335  <<theProjectile->GetDefinition()->GetParticleName() <<G4endl;
336  G4cout <<"TARGET = " <<theTarget->GetName() <<G4endl;
337  G4cout <<"K. ENERGY/NUC = " <<EPN/MeV <<" MeV/n" <<G4endl;
338  G4cout <<"CROSS SECTION = " <<xsection/millibarn <<" MILLIBARNS" <<G4endl;
339  G4cout <<"***************************************************************"
340  <<G4endl;
341  }
342 
343  return xsection;
344 }
346 //
347 void G4DPMJET2_5CrossSection::Initialise ()
348 {
349  static G4StableIsotopes theDefaultIsotopes; // natural abundances and
350  // stable isotopes
351  //verboseLevel = 2;
352 //
353 //
354 // Determine first if the environment variable G4DPMJET2_5DATA is set. If not
355 // then ask for it to be set and call exception.
356 //
357  if ( !getenv("G4DPMJET2_5DATA") )
358  {
359  G4cout <<"ENVIRONMENT VARIABLE G4DPMJET2_5DATA NOT SET " <<G4endl;
360  throw G4HadronicException(__FILE__, __LINE__,
361  "Please setenv G4DPMJET2_5DATA to point to the dpmjet2.5 data files.");
362  }
363 
364  G4String filename = G4String(getenv("G4DPMJET2_5DATA")) + "/" +
365  "GlauberCrossSections.dat";
366 
367  std::ifstream glauberXSFile(filename);
368  if (glauberXSFile) {
369 //
370 //
371 // Glaubercross-section file does exist, so read in maximum and minimum A
372 // for target and projectile.
373 //
374  glauberXSFile >>APmin >>APmax >>ATmin >>ATmax;
375 //
376 //
377 // Determine the list of targets based on the G4ElementList. The list of
378 // target nucleon numbers is stored as a ket to the map theCrossSectionIndex.
379 // G4double[240][3] array objects are created to allow storage of the
380 // cross-section fit parameters.
381 //
382  const G4ElementTable *theElementTable = G4Element::GetElementTable();
383  G4ElementTable::const_iterator it;
384  for (it=theElementTable->begin(); it!=theElementTable->end(); it++)
385  {
386  G4int nIso = (*it)->GetNumberOfIsotopes();
387  if (nIso) {
388 //
389 //
390 // The user has defined her own isotopes associated with this element. Read
391 // the nucleon numbers.
392 //
393  G4IsotopeVector* isoVector = (*it)->GetIsotopeVector();
394  for (G4int i = 0; i < nIso; i++)
395  {
396  G4int AA = (*isoVector)[i]->GetN();
397  if (theCrossSectionIndex.count(AA) == 0 && AA >= ATmin && AA <= ATmax)
398  {
399 //
400 //
401 // Whilst the use of std::map should eliminate duplication of keys, we need to
402 // know whether isotope's with the same nucleon number have been declared before
403 // creating the large arrays, hence the use of the "count" member function.
404 //
405  G4DPMJET2_5CrossSectionParamSet *a =
406  new G4DPMJET2_5CrossSectionParamSet[maxA];
407  theCrossSectionIndex.insert(
408  G4DPMJET2_5CrossSectionIndex::value_type(AA,a));
409  }
410  }
411  } else {
412 //
413 //
414 // Determine the natural isotopic abundances for this element.
415 //
416  G4int ZZ = G4lrint((*it)->GetZ());
417  nIso = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
418  G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
419  G4int AA = 0;
420  for (G4int i = 0; i < nIso; i++) {
421  AA = theDefaultIsotopes.GetIsotopeNucleonCount(index+i);
422  if (theCrossSectionIndex.count(AA) == 0 && AA >= ATmin && AA <= ATmax)
423  {
424  G4DPMJET2_5CrossSectionParamSet *a =
425  new G4DPMJET2_5CrossSectionParamSet[maxA];
426  theCrossSectionIndex.insert(
427  G4DPMJET2_5CrossSectionIndex::value_type(AA,a));
428  }
429  }
430  }
431  }
432 //
433 //
434 // Now proceed to read in the remainder of the GlauberCrossSection.dat file,
435 // loading into theCrossSectionIndex any relevant fitting parameters to the
436 // target nuclei.
437 //
438  char inputChars[80]={' '};
439  G4String inputLine;
440  while (-glauberXSFile.getline(inputChars, 80).eof() != EOF)
441  {
442  inputLine = inputChars;
443  if (inputLine.length() != 0)
444  {
445  std::istringstream tmpStream(inputLine);
446  G4int AP, AT;
447  G4double cc0, cc1, cc2;
448  tmpStream >>AP >>AT >>cc0 >>cc1 >>cc2;
449  G4DPMJET2_5CrossSectionIndex::iterator it = theCrossSectionIndex.find(AT);
450  if (it != theCrossSectionIndex.end())
451  {
452  G4DPMJET2_5CrossSectionParamSet *ptr = (it->second) + AP;
453  *ptr = G4DPMJET2_5CrossSectionParamSet(cc0,cc1,cc2);
454  }
455  }
456  }
457 
458  glauberXSFile.close();
459  G4cout << "G4DPMJET2_5CrossSection::Initialise () done!" << G4endl;
460  }
461  else {
462  G4cout <<"GlauberCrossSections.dat DOES NOT EXIST" <<G4endl;
463  throw G4HadronicException(__FILE__, __LINE__,
464  "GlauberCrossSections.dat should be located in $G4DPMJET2_5DATA directory.");
465  }
466 }
468 //
469 void G4DPMJET2_5CrossSection::BuildPhysicsTable(const G4ParticleDefinition&)
470 {;}
472 //
473 void G4DPMJET2_5CrossSection::DumpPhysicsTable(const G4ParticleDefinition
474  &theProjectile)
475 {
476  const G4int AP = G4lrint(theProjectile.GetBaryonNumber());
477  G4cout <<G4endl;
478  G4cout <<"G4DPMJET2_5CrossSection::DumpPhysicsTable" <<G4endl;
479  G4cout <<"DUMPING CROSS-SECTION FITTING COEFFICIENTS FOR AP = "
480  <<AP <<G4endl;
481  G4cout <<G4endl;
482  G4cout <<" AT"
483  <<" c0"
484  <<" c1"
485  <<" c2"
486  <<G4endl;
487  G4DPMJET2_5CrossSectionIndex::iterator it;
488  for (it=theCrossSectionIndex.begin(); it!=theCrossSectionIndex.end(); it++)
489  {
490  G4cout.unsetf(std::ios::scientific);
491  G4cout.setf(std::ios::fixed|std::ios::right|std::ios::adjustfield);
492  G4cout.precision(0);
493  G4cout <<std::setw(5) <<it->first;
494 
495  G4cout.unsetf(std::ios::fixed);
496  G4cout.setf(std::ios::scientific|std::ios::right|std::ios::adjustfield);
497  G4cout.precision(7);
498  G4DPMJET2_5CrossSectionParamSet *ptr = (it->second) + AP;
499  G4cout <<std::setw(15) <<(*ptr)[0]
500  <<std::setw(15) <<(*ptr)[1]
501  <<std::setw(15) <<(*ptr)[2]
502  <<G4endl;
503  }
504  G4cout.setf(std::ios::fixed);
505 }
506 #endif
Generated on Sat May 25 2013 14:32:56 for Geant4 by   doxygen 1.8.4