GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestAxialFormFactor.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestAxialFormFactor
5 
6 \brief Program used for testing / debugging the axial form factor
7 
8 \author Aaron Meyer <asmeyer2012 \at uchicago.edu>
9  based off testElFormFactors by
10  Costas Andreopoulos <c.andreopoulos \at cern.ch>
11  University of Liverpool
12 
13 \created August 20, 2013
14 
15 \syntax gtestAxialFormFactor [-l] [-o output_filename_prefix]
16  [-c min,max,inc[,min,max,inc...]]
17 
18  [] denotes an optionan argument
19  -l switch which evaluates over logarithmic Q^2
20  -o output filename prefix for root file output
21  -c allows scanning over z-expansion coefficients
22  -- must give a multiple of 3 numbers as arguments: min, max, increment
23  -- will scan over at most MAX_COEF coefficients and fill ntuple tree with all entries
24  -- can change MAX_COEF and recompile as necessary
25  -- to skip scanning over a coefficient, put max < min for that coefficient
26  note that Kmax in UserPhysicsOptions should match the number of scanned coeffs
27 
28 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
29  For the full text of the license visit http://copyright.genie-mc.org
30 
31 */
32 //____________________________________________________________________________
33 
34 #include <string>
35 #include <sstream>
36 
37 #include <TFile.h>
38 #include <TTree.h>
39 //#include <TNtupleD.h>
40 
55 
56 // number of coefficient values stored in tree
57 #define MAX_COEF 3
58 // Q^2 ranges (GeV^2) and number of entries for tree
59 #define Q2_LIN 4
60 #define Q2_LOG_MIN 0.1
61 #define Q2_LOG_MAX 100
62 #define N_TREE 100
63 
64 using namespace genie;
65 using std::string;
66 using std::ostringstream;
67 
68 struct tdata // addresses of data which will be added to tree
69 {
70  int* pKmax;
71  int* pmod;
72  double* pQ2;
73  double* pz;
74  double* pFA;
75  double* pcAn; // passed as just the pointer to array
76 };
77 
78 void GetCommandLineArgs (int argc, char ** argv);
80  TTree * affnt, tdata params, \
81  AlgFactory * algf, Registry * r, \
82  const Registry * gc, \
83  Interaction* interaction, double t0);//, int Kmax);
84 
85 bool IncrementCoefficients(double* coefmin, double* coefmax, double* coefinc, \
86  int kmaxinc, tdata params, AlgFactory * algf,
87  Registry* r, const Registry * gc);
88 
89 string kDefOptEvFilePrefix = "test.axialff";
92 double gOptCoeffMin[MAX_COEF] = {0.};
93 double gOptCoeffMax[MAX_COEF] = {0.};
94 double gOptCoeffInc[MAX_COEF] = {0.};
95 int gOptKmaxInc = 0;
96 
97 //__________________________________________________________________________
98 int main(int argc, char ** argv)
99 {
100 
101  GetCommandLineArgs(argc,argv);
102 
103  if (gOptKmaxInc > 0)
104  {
105  LOG("testAxialFormFactor", pINFO) << "Found coefficient ranges:";
106  for (int ik=0;ik<gOptKmaxInc;ik++)
107  {
108  LOG("testAxialFormFactor", pINFO) << "( min:" << gOptCoeffMin[ik] \
109  << " , max:" << gOptCoeffMax[ik] \
110  << " , inc:" << gOptCoeffInc[ik] << " )";
111  }
112  }
113 
114  TTree * affnt = new TTree("affnt","Axial Form Factor Test Tree");
115  AxialFormFactor axff;
116  AlgFactory * algf = AlgFactory::Instance();
118  const Registry * gc = confp->GlobalParameterList();
119 
120  AlgId id("genie::ZExpAxialFormFactorModel","Default");
121  const Algorithm * alg = algf->GetAlgorithm(id);
122  Registry * r = confp->FindRegistry(alg);
123 
124  int Kmax;
125  int mod;
126  double Q2;
127  double z;
128  double FA;
129  double cAn[MAX_COEF];
130  affnt = new TTree("axff","Axial Form Factor Test");
131  affnt->Branch("Kmax", &Kmax, "Kmax/I"); // number of coefficients
132  affnt->Branch("mod", &mod, "mod/I" ); // model
133  affnt->Branch("Q2", &Q2, "Q2/D" ); // Q^2
134  affnt->Branch("z" , &z , "z/D" ); // z
135  affnt->Branch("FA", &FA, "FA/D" ); // F_A axial form factor
136  affnt->Branch("cAn", cAn, "cAn[Kmax]/D"); // z-expansion coefficients
137 
138  // load all parameters into a struct to make passing them as arguments easier
139  tdata params; //defined tdata struct above
140  params.pKmax = &Kmax;
141  params.pmod = &mod;
142  params.pQ2 = &Q2;
143  params.pz = &z;
144  params.pFA = &FA;
145  params.pcAn = cAn;
146 
147  // flag for single output or loop
148  bool do_once = (gOptKmaxInc < 1);
149  bool do_Q4limit = r->GetBoolDef("QEL-Q4limit", gc->GetBool("QEL-Q4limit"));
150 
151  if (gOptKmaxInc < 1) Kmax = r->GetIntDef("QEL-Kmax", gc->GetInt("QEL-Kmax"));
152  else Kmax = gOptKmaxInc;
153  //if (do_Q4limit) Kmax = Kmax+4;
154  if (do_Q4limit) LOG("testAxialFormFactor",pWARN) \
155  << "Q4limit specified, note that coefficients will be altered";
156 
157  // initialize to zero
158  for (int ip=0;ip<MAX_COEF;ip++)
159  {
160  params.pcAn[ip] = 0.;
161  }
162 
163  // Get T0 from registry value
164  double t0 = r->GetDoubleDef("QEL-T0", gc->GetDouble("QEL-T0"));
165 
166  // make sure coefficients are getting incremented correctly
167  if (MAX_COEF < gOptKmaxInc ||
168  (gOptKmaxInc == 0 && MAX_COEF < Kmax) )
169  {
170  LOG("testAxialFormFactor",pWARN) \
171  << "Too many coefficient ranges, reducing to " << MAX_COEF;
172  Kmax = MAX_COEF;
174  }
175  LOG("testAxialFormFactor",pWARN) << "pKmax = " << *params.pKmax;
176 
177  // set up interaction
178  Interaction * interaction =
180 
181  if (do_once)
182  { // calculate once and be done with it
183 
184  // pre-load parameters
185  ostringstream alg_key;
186  for (int ip=1;ip<Kmax+1;ip++)
187  {
188  alg_key.str("");
189  alg_key << "QEL-Z_A" << ip;
190  params.pcAn[ip-1] =
191  r->GetDoubleDef(alg_key.str(), gc->GetDouble(alg_key.str()));
192  LOG("testAxialFormFactor",pINFO) << "Loading " << alg_key.str() \
193  << " : " << params.pcAn[ip-1];
194  }
195 
196  CalculateFormFactor(axff,affnt,params,algf,r,gc,interaction,t0);//,t0,Kmax);
197  } else {
198  // override and control coefficient values
199  r->UnLock();
200  r->Set("QEL-Kmax",*params.pKmax);
201 
202  ostringstream alg_key;
203  for (int ip=1;ip<gOptKmaxInc+1;ip++)
204  {
205  alg_key.str("");
206  alg_key << "QEL-Z_A" << ip;
207  r->Set(alg_key.str(),gOptCoeffMin[ip-1]);
208  params.pcAn[ip-1] = gOptCoeffMin[ip-1];
209  LOG("testAxialFormFactor",pWARN) << "Set parameter: " << params.pcAn[ip-1];
210  }
211  algf->ForceReconfiguration();
212 
213  do
214  { CalculateFormFactor(axff,affnt,params,algf,r,gc,interaction,t0); }//,t0,Kmax); }
216  params,algf,r,gc));
217 
218  }
219 
220  TFile f((gOptEvFilePrefix + ".root").c_str(),"recreate");
221  affnt->Write();
222  f.Close();
223 
224  return 0;
225 }
226 //__________________________________________________________________________
228  TTree * affnt, tdata params, \
229  AlgFactory * algf, Registry * r, \
230  const Registry * gc, \
231  Interaction* interaction, double t0) //, int Kmax)
232 {
233  const AxialFormFactorModelI * dipole =
234  dynamic_cast<const AxialFormFactorModelI *> (
235  algf->GetAlgorithm("genie::DipoleAxialFormFactorModel", "Default"));
236  const AxialFormFactorModelI * zexp =
237  dynamic_cast<const AxialFormFactorModelI *> (
238  algf->GetAlgorithm("genie::ZExpAxialFormFactorModel", "Default"));
239 
240  interaction->KinePtr()->SetQ2(0.);
241  double t = interaction->KinePtr()->q2();
242  double tcut = 9.0 * TMath::Power(constants::kPi0Mass, 2);
243  double Q2 = 0.;
244 
245  for(int iq=0; iq<N_TREE; iq++) {
246 
248  {
249  Q2 = TMath::Exp( TMath::Log(double(Q2_LOG_MIN))
250  + (double(iq+1)/double(N_TREE))
251  * TMath::Log(double(Q2_LOG_MAX)/double(Q2_LOG_MIN)));
252  //Q2 = TMath::Exp(((iq+1)*double(Q2_RANGE_LOG/N_TREE))*TMath::Log(10.)); // logarithmic
253  } else { Q2 = (iq+1)*double(Q2_LIN)/double(N_TREE) ; } // linear
254 
255  interaction->KinePtr()->SetQ2(Q2);
256  *params.pQ2 = Q2;
257 
258  // calculate z parameter used in expansion
259  t = interaction->KinePtr()->q2();
260  double znum = TMath::Sqrt(tcut - t) - TMath::Sqrt(tcut - t0);
261  double zden = TMath::Sqrt(tcut - t) + TMath::Sqrt(tcut - t0);
262  double zpar = znum/zden;
263 
264  if (zpar != zpar) LOG("testAxialFormFactor", pERROR) << "Undefined expansion parameter";
265  *params.pz = zpar;
266 
267  axff.SetModel(dipole);
268  axff.Calculate(interaction);
269  *params.pmod = 0;
270  *params.pFA = axff.FA();
271  affnt->Fill();
272 
273  axff.SetModel(zexp);
274  axff.Calculate(interaction);
275  *params.pmod = 1;
276  *params.pFA = axff.FA();
277  affnt->Fill();
278 
279  }
280 }
281 //__________________________________________________________________________
282 bool IncrementCoefficients(double* coefmin, double* coefmax, double* coefinc, \
283  int kmaxinc, tdata params, \
284  AlgFactory * algf, Registry* r, const Registry * gc)
285 
286 {
287 
288  if (kmaxinc < 1)
289  {
290  LOG("testAxialFormFactor",pERROR) << "No coefficients to increment";
291  return false;
292  } else {
293 
294  ostringstream alg_key;
295  int ip = -1;
296  bool stopflag;
297  do
298  {
299  if (ip > -1)
300  { // a previous iteration went over max
301  params.pcAn[ip] = coefmin[ip];
302  r->Set(alg_key.str(),params.pcAn[ip]);
303  }
304  stopflag = true;
305  alg_key.str(""); // reset sstream
306  ip++; // increment index
307  if (ip == kmaxinc) return false; // stop if gone too far
308  alg_key << "QEL-Z_A" << ip+1; // get new name
309  params.pcAn[ip] += coefinc[ip]; // increment parameter
310  r->Set(alg_key.str(),params.pcAn[ip]);
311  if (params.pcAn[ip] > coefmax[ip]) stopflag=false;
312 
313  } while (! stopflag); // loop
314 
315 
316  } // if kmaxinc < 1
317 
318  algf->ForceReconfiguration();
319  return true;
320 
321 }
322 //__________________________________________________________________________
323 void GetCommandLineArgs(int argc, char ** argv)
324 {
325  LOG("testAxialFormFactor", pINFO) << "Parsing command line arguments";
326 
327  CmdLnArgParser parser(argc,argv);
328 
329  // event file prefix
330  if( parser.OptionExists('o') ) {
331  LOG("testAxialFormFactor", pINFO) << "Reading the event filename prefix";
332  gOptEvFilePrefix = parser.ArgAsString('o');
333  } else {
334  LOG("testAxialFormFactor", pDEBUG)
335  << "Will set the default event filename prefix";
337  } //-o
338 
339  // logarithmic vs linear in q2
340  if( parser.OptionExists('l') ) {
341  LOG("testAxialFormFactor", pINFO) << "Looping over logarithmic Q2";
342  gOptDoLogarithmicQ2 = true;
343  } else {
344  LOG("testAxialFormFactor", pINFO) << "Looping over linear Q2";
345  gOptDoLogarithmicQ2 = false;
346  } //-l
347 
348  if( parser.OptionExists('c') ) {
349  LOG("testAxialFormFactor", pINFO) << "Reading Coefficient Ranges";
350  string coef = parser.ArgAsString('c');
351 
352  // split into sections of min,max,inc(rement)
353  vector<string> coefrange = utils::str::Split(coef, ",");
354  assert(coefrange.size() % 3 == 0);
355  gOptKmaxInc = coefrange.size() / 3;
356  LOG("testAxialFormFactor", pINFO) << "Number of ranges to implement : " << gOptKmaxInc;
357  for (int ik = 0;ik<gOptKmaxInc;ik++)
358  {
359  gOptCoeffMin[ik] = atof(coefrange[ik*3 ].c_str());
360  gOptCoeffMax[ik] = atof(coefrange[ik*3+1].c_str());
361  gOptCoeffInc[ik] = atof(coefrange[ik*3+2].c_str());
362  }
363 
364  }
365 
366 }
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
Definition: Registry.cxx:535
int gOptKmaxInc
double gOptCoeffMax[MAX_COEF]
string gOptEvFilePrefix
Definition: gAtmoEvGen.cxx:310
#define pERROR
Definition: Messenger.h:59
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
Pure abstract base class. Defines the AxialFormFactorModelI interface to be implemented by LlewellynS...
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:40
const int kPdgNuMu
Definition: PDGCodes.h:30
bool gOptDoLogarithmicQ2
Algorithm abstract base class.
Definition: Algorithm.h:54
#define Q2_LIN
void Calculate(const Interaction *interaction)
Calculate the form factors for the input interaction using the attached algorithm.
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:474
RgInt GetInt(RgKey key) const
Definition: Registry.cxx:467
static Interaction * QELCC(int tgt, int nuc, int probe, double E=0)
void SetModel(const AxialFormFactorModelI *model)
Attach an algorithm.
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
string kDefOptEvFilePrefix
Definition: gAtmoEvGen.cxx:320
A class holding the Axial Form Factor.
Summary information for an interaction.
Definition: Interaction.h:56
double q2(bool selected=false) const
Definition: Kinematics.cxx:141
#define Q2_LOG_MIN
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void ForceReconfiguration(bool ignore_alg_opt_out=false)
Definition: AlgFactory.cxx:131
#define N_TREE
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
double FA(void) const
Get the computed axial form factor.
#define pINFO
Definition: Messenger.h:62
RgInt GetIntDef(RgKey key, RgInt def_opt, bool set_def=true)
Definition: Registry.cxx:530
#define MAX_COEF
Program used for testing / debugging the axial form factor.
Registry * GlobalParameterList(void) const
#define pWARN
Definition: Messenger.h:60
void CalculateFormFactor(AxialFormFactor axff, TTree *affnt, tdata params, AlgFactory *algf, Registry *r, const Registry *gc, Interaction *interaction, double t0)
void UnLock(void)
unlocks the registry (doesn&#39;t unlock items)
Definition: Registry.cxx:153
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:34
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
double gOptCoeffInc[MAX_COEF]
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double gOptCoeffMin[MAX_COEF]
RgBool GetBool(RgKey key) const
Definition: Registry.cxx:460
const int kPdgTgtFe56
Definition: PDGCodes.h:205
Registry * FindRegistry(string key) const
const int kPdgProton
Definition: PDGCodes.h:81
Command line argument parser.
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
bool IncrementCoefficients(double *coefmin, double *coefmax, double *coefinc, int kmaxinc, tdata params, AlgFactory *algf, Registry *r, const Registry *gc)
void Set(RgIMapPair entry)
Definition: Registry.cxx:267
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
RgBool GetBoolDef(RgKey key, RgBool def_opt, bool set_def=true)
Definition: Registry.cxx:525
static AlgConfigPool * Instance()
#define pDEBUG
Definition: Messenger.h:63
#define Q2_LOG_MAX