GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestDISSF.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestDISSF
5 
6 \brief Program used for testing / debugging the GENIE DIS SF models
7 
8  Syntax :
9  gtestDISSF -a model -c config [-m mode] [-x x] [-q Q2]
10 
11  Options :
12  -a DIS SF model (algorithm name, eg genie::BYStructureFuncModel)
13  -c DIS SF model configuration
14  -m mode (1: make std SF ntuple, 2: vertical slice)
15  [default:1]
16  -x Specify Bjorken x to be used at the vertical slice
17  -q Specify mom. transfer Q2(>0) to be used at the vertical slice
18 
19 
20 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
21  University of Liverpool
22 
23 \created June 20, 2004
24 
25 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
26  For the full text of the license visit http://copyright.genie-mc.org
27 
28 */
29 //____________________________________________________________________________
30 
31 #include <string>
32 
33 #include <TFile.h>
34 #include <TTree.h>
35 
43 
44 using namespace genie;
45 using std::string;
46 
47 void GetCommandLineArgs(int argc, char ** argv);
48 void PrintSyntax(void);
49 
50 void BuildStdNtuple (void);
51 void VerticalSlice (void);
52 
53 int gMode = 1;
54 double gX = 0;
55 double gQ2 = 0;
56 string gDISSFAlg = "";
57 string gDISSFConfig = "";
58 
59 //__________________________________________________________________________
60 int main(int argc, char ** argv)
61 {
62  // -- parse the command line arguments (model choice,...)
63  GetCommandLineArgs(argc,argv);
64 
65  if(gMode==1) BuildStdNtuple();
66  if(gMode==2) VerticalSlice ();
67 
68  return 0;
69 }
70 //__________________________________________________________________________
71 void BuildStdNtuple(void)
72 {
73  // -- define initial states & x,Q2 values to compute the DIS SFs
74  const int kNNu = 6;
75  const int kNNuc = 2;
76  const int kNTgt = 3;
77  const int kNX = 12;
78  const int kNQ2 = 60;
79 
80  int neutrino [kNNu] = { kPdgNuE, kPdgAntiNuE, kPdgNuMu, kPdgAntiNuMu, kPdgNuTau, kPdgAntiNuTau };
81  int hit_nucleon [kNNuc] = { kPdgProton, kPdgNeutron };
82  int target [kNTgt] = { kPdgTgtFreeP, kPdgTgtFreeN, kPdgTgtFe56 };
83 
84  double xbj[kNX] = {
85  0.0150, 0.0450, 0.0800, 0.1250, 0.1750, 0.2250,
86  0.2750, 0.3500, 0.4500, 0.5500, 0.6500, 0.7500
87  };
88  double Q2[kNQ2] = {
89  0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0,
90  3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
91  13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0,
92  23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 42.0,
93  33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, 40.0, 41.0, 52.0,
94  43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0, 50.0, 51.0, 62.0
95  };
96 
97  // -- request the specified DIS SF model
99  const DISStructureFuncModelI * dissf_model =
100  dynamic_cast<const DISStructureFuncModelI *> (
102  assert(dissf_model);
103 
104  // -- create a structure functions objects and attach the specified model
105  DISStructureFunc dissf;
106  dissf.SetModel(dissf_model);
107 
108  // -- define the outout TTree
109  TTree * dissf_nt = new TTree("dissf_nt","DIS SF");
110  int brNu; // event number
111  int brNuc; // hit nucleon PDG code
112  int brTgt; // neutrino PDG code
113  double brXbj; //
114  double brQ2; //
115  double brF1; //
116  double brF2; //
117  double brF3; //
118  double brF4; //
119  double brF5; //
120  dissf_nt->Branch("nu", &brNu, "nu/I");
121  dissf_nt->Branch("nuc", &brNuc, "nuc/I");
122  dissf_nt->Branch("tgt", &brTgt, "tgt/I");
123  dissf_nt->Branch("x", &brXbj, "x/D");
124  dissf_nt->Branch("Q2", &brQ2, "Q2/D");
125  dissf_nt->Branch("F1", &brF1, "F1/D");
126  dissf_nt->Branch("F2", &brF2, "F2/D");
127  dissf_nt->Branch("F3", &brF3, "F3/D");
128  dissf_nt->Branch("F4", &brF4, "F4/D");
129  dissf_nt->Branch("F5", &brF5, "F5/D");
130 
131  // -- loop over the initial states
132  for(int inu=0; inu<kNNu; inu++) {
133  for(int inuc=0; inuc<kNNuc; inuc++) {
134  for(int itgt=0; itgt<kNTgt; itgt++) {
135 
136  // -- get a DIS interaction object & access its kinematics
137  Interaction * interaction = Interaction::DISCC(
138  target[itgt],hit_nucleon[inuc],neutrino[inu]);
139  Kinematics * kine = interaction->KinePtr();
140 
141  // -- loop over x,Q2
142  for(int ix=0; ix<kNX; ix++) {
143  for(int iq=0; iq<kNQ2; iq++) {
144 
145  // update the interaction kinematics
146  kine->Setx(xbj[ix]);
147  kine->SetQ2(Q2[iq]);
148 
149  // calculate the structure functions for the input interaction
150  dissf.Calculate(interaction);
151 
152  // update all tree branches & save
153  brNu = neutrino[inu];
154  brNuc = hit_nucleon[inuc];
155  brTgt = target[itgt];
156  brXbj = xbj[ix];
157  brQ2 = Q2[iq];
158  brF1 = dissf.F1();
159  brF2 = dissf.F2();
160  brF3 = dissf.F3();
161  brF4 = dissf.F4();
162  brF5 = dissf.F5();
163  dissf_nt->Fill();
164  }//Q2
165  }//x
166 
167  delete interaction;
168 
169  }//tgt
170  }//nuc
171  }//nu
172 
173  TFile f("./dissf.root","recreate");
174  dissf_nt->Write();
175  f.Close();
176 }
177 //__________________________________________________________________________
178 void VerticalSlice(void)
179 {
180  // manual override of Registry mesg level for more verbose output
181  Messenger * msg = Messenger::Instance();
182  msg->SetPriorityLevel("DISSF", pDEBUG);
183  msg->SetPriorityLevel("BodekYang", pDEBUG);
184  msg->SetPriorityLevel("LHAPDF5", pDEBUG);
185 
186  // -- request the specified DIS SF model
187  AlgFactory * algf = AlgFactory::Instance();
188  const DISStructureFuncModelI * dissf_model =
189  dynamic_cast<const DISStructureFuncModelI *> (
191  assert(dissf_model);
192 
193  // -- create a structure functions objects and attach the specified model
194  DISStructureFunc dissf;
195  dissf.SetModel(dissf_model);
196 
197  const int kNNu = 2;
198  const int kNNuc = 2;
199  const int kNTgt = 1;
200 
201  int neutrino [kNNu] = { kPdgNuMu, kPdgAntiNuMu };
202  int hit_nucleon [kNNuc] = { kPdgProton, kPdgNeutron };
203  int target [kNTgt] = { kPdgTgtFe56 };
204 
205  for(int inu=0; inu<kNNu; inu++) {
206  for(int inuc=0; inuc<kNNuc; inuc++) {
207  for(int itgt=0; itgt<kNTgt; itgt++) {
208 
209  // -- get a DIS interaction object & access its kinematics
210  Interaction * interaction = Interaction::DISCC(
211  target[itgt],hit_nucleon[inuc],neutrino[inu]);
212  Kinematics * kine = interaction->KinePtr();
213  kine->Setx(gX);
214  kine->SetQ2(gQ2);
215 
216  LOG("test", pNOTICE) << "*** Vertical SF slice for: ";
217  LOG("test", pNOTICE) << *interaction;
218 
219  // calculate the structure functions for the input interaction
220  dissf.Calculate(interaction);
221 
222  LOG("test", pNOTICE) << dissf;
223 
224  delete interaction;
225  }
226  }
227  }
228 
229  // print algorithms & heir configurations
230  LOG("test", pNOTICE) << *algf;
231 }
232 //__________________________________________________________________________
233 void GetCommandLineArgs(int argc, char ** argv)
234 {
235 // Parse the command line arguments
236 
237  CmdLnArgParser parser(argc,argv);
238 
239  // DIS SF alg:
240  if( parser.OptionExists('a') ) {
241  LOG("test", pINFO) << "Reading DIS SF algorithm name";
242  gDISSFAlg = parser.ArgAsString('a');
243  } else {
244  LOG("test", pINFO) << "No DIS SF algorithm was specified";
245  PrintSyntax();
246  exit(1);
247  }
248 
249  // DIS SF config:
250  if( parser.OptionExists('c') ) {
251  LOG("test", pINFO) << "Reading DIS SF algorithm config name";
252  gDISSFConfig = parser.ArgAsString('c');
253  } else {
254  LOG("test", pINFO) << "No DIS SF algorithm config was specified";
255  PrintSyntax();
256  exit(1);
257  }
258 
259  // testDISSF mode:
260  if( parser.OptionExists('m') ) {
261  LOG("test", pINFO) << "Reading testDISSF mode";
262  gMode = parser.ArgAsInt('m');
263  } else {
264  LOG("test", pINFO) << "No testDISSF was specified. Using default";
265  }
266 
267  // x,Q2 for vertical slice mode
268  if(gMode==2) {
269  // x:
270  if( parser.OptionExists('x') ) {
271  LOG("test", pINFO) << "Reading x";
272  gX = parser.ArgAsDouble('x');
273  } else {
274  LOG("test", pINFO)
275  << "No Bjorken x was specified for vertical slice";
276  PrintSyntax();
277  exit(1);
278  }
279  if( parser.OptionExists('q') ) {
280  LOG("test", pINFO) << "Reading Q2";
281  gQ2 = parser.ArgAsDouble('q');
282  } else {
283  LOG("test", pINFO)
284  << "No momentum transfer Q2 was specified for vertical slice";
285  PrintSyntax();
286  exit(1);
287  }
288  }//mode=2
289 }
290 //__________________________________________________________________________
291 void PrintSyntax(void)
292 {
293  LOG("test", pNOTICE)
294  << "\n\n" << "Syntax:" << "\n"
295  << " testDISSF -a model -c config [-m mode] [-x x] [-q Q2]\n";
296 }
297 //____________________________________________________________________________
298 
299 
void SetModel(const DISStructureFuncModelI *model)
Attach an algorithm.
Pure Abstract Base Class. Defines the DISStructureFuncModelI interface to be implemented by any algor...
string gDISSFAlg
Definition: gtestDISSF.cxx:56
const int kPdgNuE
Definition: PDGCodes.h:28
double F2(void) const
Get the computed structure function F2.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const int kPdgNuMu
Definition: PDGCodes.h:30
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
double F4(void) const
Get the computed structure function F4.
string gDISSFConfig
Definition: gtestDISSF.cxx:57
Summary information for an interaction.
Definition: Interaction.h:56
double F1(void) const
Get the computed structure function F1.
double gX
Definition: gtestDISSF.cxx:54
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static Messenger * Instance(void)
Definition: Messenger.cxx:49
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
const int kPdgTgtFreeN
Definition: PDGCodes.h:200
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
#define pINFO
Definition: Messenger.h:62
void VerticalSlice(void)
Definition: gtestDISSF.cxx:178
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
double F5(void) const
Get the computed structure function F5.
A more convenient interface to the log4cpp Message Service.
Definition: Messenger.h:258
void BuildStdNtuple(void)
Definition: gtestDISSF.cxx:71
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
A class holding Deep Inelastic Scattering (DIS) Form Factors (invariant structure funstions) ...
double F3(void) const
Get the computed structure function F3.
const int kPdgNuTau
Definition: PDGCodes.h:32
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
void SetPriorityLevel(const char *stream, log4cpp::Priority::Value p)
Definition: Messenger.cxx:88
static Interaction * DISCC(int tgt, int nuc, int probe, double E=0)
double gQ2
Definition: gtestDISSF.cxx:55
const int kPdgTgtFe56
Definition: PDGCodes.h:205
void Calculate(const Interaction *interaction)
Calculate the S/F&#39;s for the input interaction using the attached algorithm.
const int kPdgProton
Definition: PDGCodes.h:81
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
const int kPdgNeutron
Definition: PDGCodes.h:83
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
void PrintSyntax(void)
int gMode
Definition: gtestDISSF.cxx:53
#define pDEBUG
Definition: Messenger.h:63