GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestNaturalIsotopes.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestNaturalIsotopes
5 
6 \brief A simple test program to test the natural isotopes table
7 
8 \author Jim Dobson <j.dobson07@imperial.ac.uk>
9  Imperial College London
10 
11 \created May 30, 2008
12 
13 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
14  For the full text of the license visit http://copyright.genie-mc.org
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <iomanip>
20 #include <string>
21 
25 #include "TParticlePDG.h"
26 #include "TMath.h"
27 
28 using namespace genie;
29 using namespace std;
30 
31 int getA(int pdg);
32 int setA(int pdg, int a);
33 std::string PDGcheck(int pdg);
34 
35 int main(int argc, char** argv)
36 {
37  bool docheckpdg = false;
38  bool calcavg = false;
39 
40  for (int iarg=1; iarg < argc; ++iarg) {
41  // argv[0] is program name, skip it
42  // std::cout << "[" << iarg << "] = \""
43  // << argv[iarg] << "\"" << std::endl;
44  std::string argvstr = std::string(argv[iarg]);
45  if ( argvstr == "--check-pdg" ) docheckpdg = true;
46  else if ( argvstr == "--avg-a" ) calcavg = true;
47  else {
48  cout << "Usage: " << argv[0] << " [--check-pdg] [--avg-a]" << endl;
49  exit(1);
50  }
51  }
52 
53  LOG("test", pNOTICE)
54  << "Testing the NaturalIsotopes utility";
55 
56  /* PDGLibrary * pdglib = */ PDGLibrary::Instance(); // get msg out
58 
59 
60  if ( docheckpdg ) {
61  LOG("test", pNOTICE)
62  << "Check on PDG in PDGLibrary: [PDG,PDG-proton,PDG-neutron]";
63  }
64 
65  for(int Z = 1; Z < 104; Z++){
66  int nel = ni->NElements(Z);
67  LOG("test", pNOTICE) << "** Z = " << Z
68  << " Number of elements = " << nel;
69 
70  int pdg;
71  double abund, avg_a = 0;
72  std::string extra;
73  for(int n = 0; n < nel; n++){
74  const NaturalIsotopeElementData * el = ni->ElementData(Z,n);
75  pdg = el->PdgCode();
76  abund = el->Abundance();
77  avg_a += double(getA(pdg)) * abund;
78  extra = ( docheckpdg ) ? PDGcheck(pdg) : "";
79  LOG("test", pNOTICE)
80  << " -- Element: " << n
81  << ", PdgCode = " << pdg
82  << extra
83  << ", Abundance = " << abund;
84  } // n
85  if ( calcavg && ( nel > 1 ) ) {
86  avg_a /= 100.; // abundances were in %
87  pdg = setA(pdg,TMath::Nint(avg_a));
88  extra = ( docheckpdg ) ? PDGcheck(pdg) : "";
89  LOG("test", pNOTICE)
90  << " "
91  << " PdgCode = " << pdg
92  << extra
93  << ", average <A> " << avg_a ;
94  }
95  } // Z
96 
97  return 0;
98 }
99 
100 int getA(int pdg)
101 { return (int(pdg/10)%1000); }
102 
103 int setA(int pdg, int a)
104 {
105  return int(pdg/10000)*10000 + 10*a;
106 }
107 
108 std::string PDGcheck(int pdg)
109 {
110  const int minus_p = 10010; // Z-1, A-1
111  const int minus_n = 10; // A-1
112 
113  PDGLibrary * pdglib = PDGLibrary::Instance();
114  int pdg_minus_p = pdg - minus_p;
115  int pdg_minus_n = pdg - minus_n;
116  const TParticlePDG * part = pdglib->Find(pdg);
117  const TParticlePDG * part_minus_p = pdglib->Find(pdg_minus_p);
118  const TParticlePDG * part_minus_n = pdglib->Find(pdg_minus_n);
119  string codeck = " [";
120  codeck += ( (part) ? "ok" : "NO" );
121  if ( pdg != 1000010010 ) {
122  codeck += ",";
123  codeck += ( (part_minus_p) ? "ok" : "NO" );
124  codeck += ",";
125  codeck += ( (part_minus_n) ? "ok" : "NO" );
126  codeck += "]";
127  } else { // special case for H1
128  codeck += ",--,--]";
129  }
130 
131  return codeck;
132 }
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const double a
std::string PDGcheck(int pdg)
int NElements(int Z) const
static NaturalIsotopes * Instance(void)
int setA(int pdg, int a)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
Singleton class to load &amp; serve tables of natural occurring isotopes.
int getA(int pdg)
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
const NaturalIsotopeElementData * ElementData(int Z, int ielement) const