gtestMuELoss.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestMuELoss
5 
6 \brief Program for the MuELoss utility package. The program saves the
7  computed data in an output ROOT ntuple.
8 
9  Syntax :
10  gtestMuELoss -m materials
11 
12  Options :
13  -m specifies a comma seperated list of materials
14  (the material ids correspond to the enumeration that can be seen
15  in $GENIE/src/MuELoss/MuELMaterial.h)
16 
17 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
18  University of Liverpool & STFC Rutherford Appleton Lab
19 
20 \created March 10, 2006
21 
22 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
23  For the full text of the license visit http://copyright.genie-mc.org
24  or see $GENIE/LICENSE
25 */
26 //____________________________________________________________________________
27 
28 #include <cassert>
29 #include <cstdlib>
30 #include <string>
31 #include <vector>
32 
33 #include <TFile.h>
34 #include <TNtuple.h>
35 
36 #include "Algorithm/AlgFactory.h"
37 #include "Conventions/Units.h"
38 #include "MuELoss/MuELossI.h"
39 #include "MuELoss/MuELMaterial.h"
40 #include "MuELoss/MuELProcess.h"
41 #include "Messenger/Messenger.h"
42 #include "Utils/StringUtils.h"
43 #include "Utils/CmdLnArgParser.h"
44 
45 using std::string;
46 using std::vector;
47 using namespace genie;
48 using namespace genie::utils;
49 using namespace genie::mueloss;
50 
51 void GetCommandLineArgs (int argc, char ** argv);
52 
54 
55 //____________________________________________________________________________
56 int main(int argc, char ** argv)
57 {
58  GetCommandLineArgs(argc, argv);
59 
60  const int N = 14;
61  double E[N] = {1,5,10,15,20,30,50,100,200,500, 1000,2000, 5000, 9000}; //GeV
62 
63  // split the comma separated list of materials
64  vector<string> mtv = utils::str::Split(gOptMaterials, ",");
65 
66  // get muon energy loss algorithms
68 
69  const MuELossI * betheBloch =
70  dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
71  "genie::mueloss::BetheBlochModel","Default"));
72 
73  const MuELossI * petrukhinShestakov =
74  dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
75  "genie::mueloss::PetrukhinShestakovModel","Default"));
76 
77  const MuELossI * kokoulinPetroukhin =
78  dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
79  "genie::mueloss::KokoulinPetrukhinModel","Default"));
80 
81  const MuELossI * bezroukovBugaev =
82  dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
83  "genie::mueloss::BezrukovBugaevModel","Default"));
84 
85  assert ( betheBloch );
86  assert ( petrukhinShestakov );
87  assert ( kokoulinPetroukhin );
88  assert ( bezroukovBugaev );
89 
90  double myunits_conversion = units::GeV/(units::g/units::cm2);
91  string myunits_name = " GeV/(gr/cm^2)";
92 
93  // open a ROOT file and define the output ntuple.
94  TFile froot("./genie-mueloss.root", "RECREATE");
95  TNtuple muntp("muntp","muon dE/dx", "material:E:ion:brem:pair:pnucl");
96 
97  //loop over materials
98  vector<string>::iterator iter;
99  for(iter = mtv.begin(); iter != mtv.end(); ++iter) {
100 
101  MuELMaterial_t mt = (MuELMaterial_t) atoi(iter->c_str());
102 
103  LOG("test", pINFO)
104  << "---------- Computing/Printing muon energy losses in "
105  << MuELMaterial::AsString(mt) << " ----------";
106 
107  // loop over energies
108  for(int i=0; i<N; i++) {
109 
110  // ionization
111  double ion = betheBloch->dE_dx(E[i],mt) / myunits_conversion;
112 
113  LOG("test", pINFO)
114  << "Process: " << MuELProcess::AsString(betheBloch->Process())
115  << ", Model: " << betheBloch->Id().Key()
116  << " : \n -dE/dx(E=" << E[i] << ") = " << ion << myunits_name;
117 
118  // bremsstrahlung
119  double brem = petrukhinShestakov->dE_dx(E[i],mt) / myunits_conversion;
120 
121  LOG("test", pINFO)
122  << "Process: " << MuELProcess::AsString(petrukhinShestakov->Process())
123  << ", Model: " << petrukhinShestakov->Id().Key()
124  << " : \n -dE/dx(E=" << E[i] << ") = " << brem << myunits_name;
125 
126  // e-e+ pair production
127  double pair = kokoulinPetroukhin->dE_dx(E[i],mt) / myunits_conversion;
128 
129  LOG("test", pINFO)
130  << "Process: " << MuELProcess::AsString(kokoulinPetroukhin->Process())
131  << ", Model: " << kokoulinPetroukhin->Id().Key()
132  << " : \n -dE/dx(E=" << E[i] << ") = " << pair << myunits_name;
133 
134  // photonuclear interactions
135  double pnucl = bezroukovBugaev->dE_dx(E[i],mt) / myunits_conversion;
136 
137  LOG("test", pINFO)
138  << "Process: " << MuELProcess::AsString(bezroukovBugaev->Process())
139  << ", Model: " << bezroukovBugaev->Id().Key()
140  << " : \n -dE/dx(E=" << E[i] << ") = " << pnucl << myunits_name
141  << "\n\n";
142 
143  muntp.Fill( (int)mt,E[i],ion,brem,pair,pnucl);
144  }//e
145  }//m
146 
147  muntp.Write();
148  froot.Close();
149 
150  return 0;
151 }
152 //____________________________________________________________________________
153 void GetCommandLineArgs(int argc, char ** argv)
154 {
155  LOG("test", pNOTICE) << "Parsing command line arguments";
156 
157  CmdLnArgParser parser(argc,argv);
158 
159  if ( parser.OptionExists('m') ) {
160  LOG("test", pINFO) << "Reading material ids";
161  gOptMaterials = parser.ArgAsString('m');
162  } else {
163  LOG("test", pINFO) << "Unspecified material ids - Exiting";
164  exit(1);
165  }
166 }
167 //____________________________________________________________________________
168 
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
virtual MuELProcess_t Process(void) const =0
Double_t ion
Definition: plot.C:18
The MuELoss utility package that computes muon energy losses in the energy range from 1 GeV to 10 TeV...
static const double g
Definition: Units.h:145
static const double cm2
Definition: Units.h:77
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
#define pINFO
Definition: Messenger.h:63
void GetCommandLineArgs(int argc, char **argv)
int main(int argc, char **argv)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
string gOptMaterials
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
exit(0)
assert(nhit_max >=nhit_nbins)
const char * AsString(Resonance_t res)
resonance id -> string
virtual double dE_dx(double E, MuELMaterial_t m) const =0
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
static const double GeV
Definition: Units.h:29
bool OptionExists(char opt)
was option set?
string Key(void) const
Definition: AlgId.h:47
Root of GENIE utility namespaces.
enum genie::mueloss::EMuELMaterial MuELMaterial_t