GiBUURESFormFactor.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (C) 2003-2015, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  May 30, 2009
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ May 30, 2009 - CA
14  Was first added in v2.5.1
15  @ Sep 29, 2009 - CA
16  Add FormFactors nested class & re-organize code
17  @ May 01, 2012 - CA
18  Pick data from new location ($GENIE/data/evgen/gibuu/)
19 */
20 //____________________________________________________________________________
21 
22 #include <cassert>
23 #include <sstream>
24 #include <string>
25 
26 #include <TSystem.h>
27 #include <TTree.h>
28 #include <TMath.h>
29 
35 
36 using std::ostringstream;
37 using std::istream;
38 using std::ios;
39 using std::cout;
40 using std::endl;
41 
42 using namespace genie;
43 using namespace genie::utils;
44 
45 //____________________________________________________________________________
47 //____________________________________________________________________________
49 {
50  this->LoadTables();
51  fInstance = 0;
52 }
53 //____________________________________________________________________________
55 {
56  delete fFormFactors;
57  fFormFactors = 0;
58 }
59 //____________________________________________________________________________
61 {
62  if(fInstance == 0) {
63  LOG("GiBUURESFormFactor", pINFO) << "GiBUURESFormFactor late initialization";
64  static GiBUURESFormFactor::Cleaner cleaner;
66  fInstance = new GiBUURESFormFactor;
67  }
68  return fInstance;
69 }
70 //____________________________________________________________________________
72 {
73  fFormFactors = new FormFactors;
74  fFormFactors->LoadTables();
75 }
76 //____________________________________________________________________________
78 {
79  return *fFormFactors;
80 }
81 //____________________________________________________________________________
82 //
83 // FormFactors nested class
84 //
85 //____________________________________________________________________________
86 double GiBUURESFormFactor::FormFactors::fMinQ2 = 0.0; // GeV^2
87 double GiBUURESFormFactor::FormFactors::fMaxQ2 = 4.0; // GeV^2
88 //____________________________________________________________________________
90 {
91 
92 }
93 //____________________________________________________________________________
95 {
96  if(!gAbortingInErr) {
97  cout << "GiBUURESFormFactor singleton dtor: Deleting all f/f splines" << endl;
98  }
99 
100  // resonance form factor splines
101  for(int r=0; r<kNRes; r++) {
102  for(int i=0; i<kNCurr; i++) {
103  for(int j=0; j<kNHitNuc; j++) {
104  for(int k=0; k<kNFFRes; k++) {
105  if (fFFRes[r][i][j][k]) {
106  delete fFFRes[r][i][j][k];
107  fFFRes[r][i][j][k] = 0;
108  }
109  }//k
110  }//j
111  }//i
112  }//r
113 }
114 //____________________________________________________________________________
116 {
117 // Loads hadronic x-section data
118 
119  for(int r=0; r<kNRes; r++) {
120  for(int i=0; i<kNCurr; i++) {
121  for(int j=0; j<kNHitNuc; j++) {
122  for(int k=0; k<kNFFRes; k++) {
123  fFFRes[r][i][j][k] = 0;
124  }//k
125  }//j
126  }//i
127  }//r
128 
129  string data_dir = string(gSystem->Getenv("GENIE")) +
130  string("/data/evgen/gibuu");
131 
132  LOG("GiBUURESFormFactor", pNOTICE) << "Loading GiBUU data from: " << data_dir;
133 
134  //
135  // load resonance form factor data
136  //
137  for(int r=0; r<kNRes; r++) {
138  for(int i=0; i<kNCurr; i++) {
139  for(int j=0; j<kNHitNuc; j++) {
140 
141  bool skip = false;
142 
143  Resonance_t resonance = (Resonance_t)r;
144 
145  ostringstream datafile;
146  datafile << data_dir << "/form_factors/";
147 
148  switch(r){
149  case ( 0): datafile << "P33_1232"; break;
150  case ( 1): datafile << "S11_1535"; break;
151  case ( 2): datafile << "D13_1520"; break;
152  case ( 3): datafile << "S11_1650"; break;
153  case ( 5): datafile << "D15_1675"; break;
154  case ( 6): datafile << "S31_1620"; break;
155  case ( 7): datafile << "D33_1700"; break;
156  case ( 8): datafile << "P11_1440"; break;
157  case (10): datafile << "P13_1720"; break;
158  case (11): datafile << "F15_1680"; break;
159  case (12): datafile << "P31_1910"; break;
160  case (14): datafile << "F35_1905"; break;
161  case (15): datafile << "F37_1950"; break;
162  default : skip = true;
163 
164  }
165  switch(i){
166  case (0): datafile << "_CC"; break;
167  case (1): datafile << "_NC"; break;
168  case (2): datafile << "_EM"; break;
169  default : skip = true;
170  }
171  switch(j){
172  case (0): datafile << "_neutron"; break;
173  case (1): datafile << "_proton"; break;
174  default : skip = true;
175  }
176  datafile << "_FFres.dat";
177 
178  if(skip) continue;
179 
180  //-- Make sure that all data file exists
181  assert( ! gSystem->AccessPathName(datafile.str().c_str()) );
182 
183  //-- Read the data and fill a tree
184 
185  // The GiBUU files have the data organized in 9 columns:
186  // 0 1 2 3 4 5 6 7 8
187  // I=3/2 res: Qs C_3^V C_4^V C_5^V C_6^V C_3^A C_4^A C_5^A C_6^A
188  // I=1/2 res: Qs F_1^V F_2^V ----- ----- F_A F_P ----- ----
189 
190  TTree data_ffres;
191  data_ffres.ReadFile(datafile.str().c_str(),
192  "Q2/D:f1/D:f2/D:f3/D:f4/D:f5/D:f6/D:f7/D:f8/D");
193 
194  LOG("GiBUURESFormFactor", pDEBUG)
195  << "Number of data rows: " << data_ffres.GetEntries();
196 
197  //
198  // I=3/2 resonances
199  //
200  if(res::IsDelta(resonance)) {
201  fFFRes[r][i][j][0] = new Spline(&data_ffres, "Q2:f1"); // F1V = f(Q2)
202  fFFRes[r][i][j][1] = new Spline(&data_ffres, "Q2:f2"); // F2V = f(Q2)
203  fFFRes[r][i][j][2] = new Spline(&data_ffres, "Q2:f5"); // FA = f(Q2)
204  fFFRes[r][i][j][3] = new Spline(&data_ffres, "Q2:f6"); // FP = f(Q2)
205  } // Delta res
206  else
207  //
208  // I=1/2 resonances
209  //
210  if(res::IsN(resonance)) {
211  fFFRes[r][i][j][4] = new Spline(&data_ffres, "Q2:f1"); // C3V = f(Q2)
212  fFFRes[r][i][j][5] = new Spline(&data_ffres, "Q2:f2"); // C4V = f(Q2)
213  fFFRes[r][i][j][6] = new Spline(&data_ffres, "Q2:f3"); // C5V = f(Q2)
214  fFFRes[r][i][j][7] = new Spline(&data_ffres, "Q2:f4"); // C6V = f(Q2)
215  fFFRes[r][i][j][8] = new Spline(&data_ffres, "Q2:f5"); // C3A = f(Q2)
216  fFFRes[r][i][j][9] = new Spline(&data_ffres, "Q2:f6"); // C4A = f(Q2)
217  fFFRes[r][i][j][10] = new Spline(&data_ffres, "Q2:f7"); // C5A = f(Q2)
218  fFFRes[r][i][j][11] = new Spline(&data_ffres, "Q2:f8"); // C6A = f(Q2)
219  } //N res
220 
221  }//j
222  }//i
223  }//r
224 
225 
226  for(int r=0; r<kNRes; r++) {
227  for(int i=0; i<kNCurr; i++) {
228  for(int j=0; j<kNHitNuc; j++) {
229  for(int k=0; k<kNFFRes; k++) {
230  if(fFFRes[r][i][j][k]) fFFRes[r][i][j][k]->YCanBeNegative(true);
231  }//k
232  }//j
233  }//i
234  }//r
235 
236  LOG("GiBUURESFormFactor", pINFO)
237  << "Done loading all resonance form factor files...";
238 }
239 //____________________________________________________________________________
241  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
242 {
243  if(!res::IsN(res)) return 0.;
244  return this->FFRes(Q2,res,hit_nucleon_pdg,it,4);
245 }
246 //____________________________________________________________________________
248  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
249 {
250  if(!res::IsN(res)) return 0.;
251  return this->FFRes(Q2,res,hit_nucleon_pdg,it,5);
252 }
253 //____________________________________________________________________________
255  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
256 {
257  if(!res::IsN(res)) return 0.;
258  return this->FFRes(Q2,res,hit_nucleon_pdg,it,6);
259 }
260 //____________________________________________________________________________
262  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
263 {
264  if(!res::IsN(res)) return 0.;
265  return this->FFRes(Q2,res,hit_nucleon_pdg,it,7);
266 }
267 //____________________________________________________________________________
269  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
270 {
271  if(!res::IsN(res)) return 0.;
272  return this->FFRes(Q2,res,hit_nucleon_pdg,it,8);
273 }
274 //____________________________________________________________________________
276  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
277 {
278  if(!res::IsN(res)) return 0.;
279  return this->FFRes(Q2,res,hit_nucleon_pdg,it,9);
280 }
281 //____________________________________________________________________________
283  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
284 {
285  if(!res::IsN(res)) return 0.;
286  return this->FFRes(Q2,res,hit_nucleon_pdg,it,10);
287 }
288 //____________________________________________________________________________
290  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
291 {
292  if(!res::IsN(res)) return 0.;
293  return this->FFRes(Q2,res,hit_nucleon_pdg,it,11);
294 }
295 //____________________________________________________________________________
297  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
298 {
299  if(!res::IsDelta(res)) return 0.;
300  return this->FFRes(Q2,res,hit_nucleon_pdg,it,0);
301 }
302 //____________________________________________________________________________
304  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
305 {
306  if(!res::IsDelta(res)) return 0.;
307  return this->FFRes(Q2,res,hit_nucleon_pdg,it,1);
308 }
309 //____________________________________________________________________________
311  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
312 {
313  if(!res::IsDelta(res)) return 0.;
314  return this->FFRes(Q2,res,hit_nucleon_pdg,it,2);
315 }
316 //____________________________________________________________________________
318  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
319 {
320  if(!res::IsDelta(res)) return 0.;
321  return this->FFRes(Q2,res,hit_nucleon_pdg,it,3);
322 }
323 //____________________________________________________________________________
325  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it,
326  int ffresid) const
327 {
328  if(Q2 < fMinQ2 || Q2 > fMaxQ2) return 0.;
329 
330  int r = -1, i = -1, j = -1;
331 
332  if(ffresid<0 || ffresid >= kNFFRes) return 0.;
333 
334  r = (int)res;
335  if(r<0 || r >= kNRes) return 0.;
336 
337  if (it == kIntWeakCC) { i = 0; }
338  else if (it == kIntWeakNC) { i = 1; }
339  else if (it == kIntEM) { i = 2; }
340 
341  if (hit_nucleon_pdg == kPdgNeutron) { j = 0; }
342  else if (hit_nucleon_pdg == kPdgProton ) { j = 1; }
343 
344  const Spline * spl = fFFRes[r][i][j][ffresid];
345 
346  if(!spl) return 0;
347  else {
348  return spl->Evaluate(Q2);
349  }
350 }
351 //____________________________________________________________________________
bool IsDelta(Resonance_t res)
is it a Delta resonance?
double F2V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
double C3A(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Singleton to load and serve data tables provided by the GiBUU group.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
static double fMinQ2
min Q2 for which resonance f/f data are given
static GiBUURESFormFactor * Instance(void)
double Evaluate(double x) const
Definition: Spline.cxx:362
double C6V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
enum genie::EResonance Resonance_t
double C4V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
bool IsN(Resonance_t res)
is it an N resonance?
const FormFactors & FF(void) const
double C5V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
double FA(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
double C6A(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double FFRes(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it, int ffid) const
func to retrieve interpolated form factor values
static double fMaxQ2
max Q2 for which resonance f/f data are given
const double j
Definition: BetheBloch.cxx:29
#define pINFO
Definition: Messenger.h:63
double F1V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
OStream cout
Definition: OStream.cxx:6
double C3V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
double C5A(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
static GiBUURESFormFactor * fInstance
double C4A(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
const int kPdgProton
Definition: PDGCodes.h:65
#define pNOTICE
Definition: Messenger.h:62
double FP(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
const int kNRes
Definition: bwnorm.C:23
bool gAbortingInErr
Definition: Messenger.cxx:56
const int kPdgNeutron
Definition: PDGCodes.h:67
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
enum genie::EInteractionType InteractionType_t
void LoadTables(void)
load all form factor data tables
Root of GENIE utility namespaces.
#define pDEBUG
Definition: Messenger.h:64
enum BeamMode string