GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GiBUURESFormFactor.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <cassert>
12 #include <sstream>
13 #include <string>
14 
15 #include <TSystem.h>
16 #include <TTree.h>
17 #include <TMath.h>
18 
24 
25 using std::ostringstream;
26 using std::istream;
27 using std::ios;
28 using std::cout;
29 using std::endl;
30 
31 using namespace genie;
32 using namespace genie::utils;
33 
34 //____________________________________________________________________________
36 //____________________________________________________________________________
38 {
39  this->LoadTables();
40  fInstance = 0;
41 }
42 //____________________________________________________________________________
44 {
45  delete fFormFactors;
46  fFormFactors = 0;
47 }
48 //____________________________________________________________________________
50 {
51  if(fInstance == 0) {
52  LOG("GiBUURESFormFactor", pINFO) << "GiBUURESFormFactor late initialization";
53  static GiBUURESFormFactor::Cleaner cleaner;
55  fInstance = new GiBUURESFormFactor;
56  }
57  return fInstance;
58 }
59 //____________________________________________________________________________
61 {
62  fFormFactors = new FormFactors;
63  fFormFactors->LoadTables();
64 }
65 //____________________________________________________________________________
67 {
68  return *fFormFactors;
69 }
70 //____________________________________________________________________________
71 //
72 // FormFactors nested class
73 //
74 //____________________________________________________________________________
75 double GiBUURESFormFactor::FormFactors::fMinQ2 = 0.0; // GeV^2
76 double GiBUURESFormFactor::FormFactors::fMaxQ2 = 4.0; // GeV^2
77 //____________________________________________________________________________
79 {
80 
81 }
82 //____________________________________________________________________________
84 {
85  if(!gAbortingInErr) {
86  cout << "GiBUURESFormFactor singleton dtor: Deleting all f/f splines" << endl;
87  }
88 
89  // resonance form factor splines
90  for(int r=0; r<kNRes; r++) {
91  for(int i=0; i<kNCurr; i++) {
92  for(int j=0; j<kNHitNuc; j++) {
93  for(int k=0; k<kNFFRes; k++) {
94  if (fFFRes[r][i][j][k]) {
95  delete fFFRes[r][i][j][k];
96  fFFRes[r][i][j][k] = 0;
97  }
98  }//k
99  }//j
100  }//i
101  }//r
102 }
103 //____________________________________________________________________________
105 {
106 // Loads hadronic x-section data
107 
108  for(int r=0; r<kNRes; r++) {
109  for(int i=0; i<kNCurr; i++) {
110  for(int j=0; j<kNHitNuc; j++) {
111  for(int k=0; k<kNFFRes; k++) {
112  fFFRes[r][i][j][k] = 0;
113  }//k
114  }//j
115  }//i
116  }//r
117 
118  string data_dir = string(gSystem->Getenv("GENIE")) +
119  string("/data/evgen/gibuu");
120 
121  LOG("GiBUURESFormFactor", pNOTICE) << "Loading GiBUU data from: " << data_dir;
122 
123  //
124  // load resonance form factor data
125  //
126  for(int r=0; r<kNRes; r++) {
127  for(int i=0; i<kNCurr; i++) {
128  for(int j=0; j<kNHitNuc; j++) {
129 
130  bool skip = false;
131 
132  Resonance_t resonance = (Resonance_t)r;
133 
134  ostringstream datafile;
135  datafile << data_dir << "/form_factors/";
136 
137  switch(r){
138  case ( 0): datafile << "P33_1232"; break;
139  case ( 1): datafile << "S11_1535"; break;
140  case ( 2): datafile << "D13_1520"; break;
141  case ( 3): datafile << "S11_1650"; break;
142  case ( 5): datafile << "D15_1675"; break;
143  case ( 6): datafile << "S31_1620"; break;
144  case ( 7): datafile << "D33_1700"; break;
145  case ( 8): datafile << "P11_1440"; break;
146  case (10): datafile << "P13_1720"; break;
147  case (11): datafile << "F15_1680"; break;
148  case (12): datafile << "P31_1910"; break;
149  case (14): datafile << "F35_1905"; break;
150  case (15): datafile << "F37_1950"; break;
151  default : skip = true;
152 
153  }
154  switch(i){
155  case (0): datafile << "_CC"; break;
156  case (1): datafile << "_NC"; break;
157  case (2): datafile << "_EM"; break;
158  default : skip = true;
159  }
160  switch(j){
161  case (0): datafile << "_neutron"; break;
162  case (1): datafile << "_proton"; break;
163  default : skip = true;
164  }
165  datafile << "_FFres.dat";
166 
167  if(skip) continue;
168 
169  //-- Make sure that all data file exists
170  assert( ! gSystem->AccessPathName(datafile.str().c_str()) );
171 
172  //-- Read the data and fill a tree
173 
174  // The GiBUU files have the data organized in 9 columns:
175  // 0 1 2 3 4 5 6 7 8
176  // 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
177  // I=1/2 res: Qs F_1^V F_2^V ----- ----- F_A F_P ----- ----
178 
179  TTree data_ffres;
180  data_ffres.ReadFile(datafile.str().c_str(),
181  "Q2/D:f1/D:f2/D:f3/D:f4/D:f5/D:f6/D:f7/D:f8/D");
182 
183  LOG("GiBUURESFormFactor", pDEBUG)
184  << "Number of data rows: " << data_ffres.GetEntries();
185 
186  //
187  // I=3/2 resonances
188  //
189  if(res::IsDelta(resonance)) {
190  fFFRes[r][i][j][0] = new Spline(&data_ffres, "Q2:f1"); // F1V = f(Q2)
191  fFFRes[r][i][j][1] = new Spline(&data_ffres, "Q2:f2"); // F2V = f(Q2)
192  fFFRes[r][i][j][2] = new Spline(&data_ffres, "Q2:f5"); // FA = f(Q2)
193  fFFRes[r][i][j][3] = new Spline(&data_ffres, "Q2:f6"); // FP = f(Q2)
194  } // Delta res
195  else
196  //
197  // I=1/2 resonances
198  //
199  if(res::IsN(resonance)) {
200  fFFRes[r][i][j][4] = new Spline(&data_ffres, "Q2:f1"); // C3V = f(Q2)
201  fFFRes[r][i][j][5] = new Spline(&data_ffres, "Q2:f2"); // C4V = f(Q2)
202  fFFRes[r][i][j][6] = new Spline(&data_ffres, "Q2:f3"); // C5V = f(Q2)
203  fFFRes[r][i][j][7] = new Spline(&data_ffres, "Q2:f4"); // C6V = f(Q2)
204  fFFRes[r][i][j][8] = new Spline(&data_ffres, "Q2:f5"); // C3A = f(Q2)
205  fFFRes[r][i][j][9] = new Spline(&data_ffres, "Q2:f6"); // C4A = f(Q2)
206  fFFRes[r][i][j][10] = new Spline(&data_ffres, "Q2:f7"); // C5A = f(Q2)
207  fFFRes[r][i][j][11] = new Spline(&data_ffres, "Q2:f8"); // C6A = f(Q2)
208  } //N res
209 
210  }//j
211  }//i
212  }//r
213 
214 
215  for(int r=0; r<kNRes; r++) {
216  for(int i=0; i<kNCurr; i++) {
217  for(int j=0; j<kNHitNuc; j++) {
218  for(int k=0; k<kNFFRes; k++) {
219  if(fFFRes[r][i][j][k]) fFFRes[r][i][j][k]->YCanBeNegative(true);
220  }//k
221  }//j
222  }//i
223  }//r
224 
225  LOG("GiBUURESFormFactor", pINFO)
226  << "Done loading all resonance form factor files...";
227 }
228 //____________________________________________________________________________
230  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
231 {
232  if(!res::IsN(res)) return 0.;
233  return this->FFRes(Q2,res,hit_nucleon_pdg,it,4);
234 }
235 //____________________________________________________________________________
237  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
238 {
239  if(!res::IsN(res)) return 0.;
240  return this->FFRes(Q2,res,hit_nucleon_pdg,it,5);
241 }
242 //____________________________________________________________________________
244  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
245 {
246  if(!res::IsN(res)) return 0.;
247  return this->FFRes(Q2,res,hit_nucleon_pdg,it,6);
248 }
249 //____________________________________________________________________________
251  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
252 {
253  if(!res::IsN(res)) return 0.;
254  return this->FFRes(Q2,res,hit_nucleon_pdg,it,7);
255 }
256 //____________________________________________________________________________
258  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
259 {
260  if(!res::IsN(res)) return 0.;
261  return this->FFRes(Q2,res,hit_nucleon_pdg,it,8);
262 }
263 //____________________________________________________________________________
265  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
266 {
267  if(!res::IsN(res)) return 0.;
268  return this->FFRes(Q2,res,hit_nucleon_pdg,it,9);
269 }
270 //____________________________________________________________________________
272  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
273 {
274  if(!res::IsN(res)) return 0.;
275  return this->FFRes(Q2,res,hit_nucleon_pdg,it,10);
276 }
277 //____________________________________________________________________________
279  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
280 {
281  if(!res::IsN(res)) return 0.;
282  return this->FFRes(Q2,res,hit_nucleon_pdg,it,11);
283 }
284 //____________________________________________________________________________
286  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
287 {
288  if(!res::IsDelta(res)) return 0.;
289  return this->FFRes(Q2,res,hit_nucleon_pdg,it,0);
290 }
291 //____________________________________________________________________________
293  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
294 {
295  if(!res::IsDelta(res)) return 0.;
296  return this->FFRes(Q2,res,hit_nucleon_pdg,it,1);
297 }
298 //____________________________________________________________________________
300  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
301 {
302  if(!res::IsDelta(res)) return 0.;
303  return this->FFRes(Q2,res,hit_nucleon_pdg,it,2);
304 }
305 //____________________________________________________________________________
307  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it) const
308 {
309  if(!res::IsDelta(res)) return 0.;
310  return this->FFRes(Q2,res,hit_nucleon_pdg,it,3);
311 }
312 //____________________________________________________________________________
314  double Q2, Resonance_t res, int hit_nucleon_pdg, InteractionType_t it,
315  int ffresid) const
316 {
317  if(Q2 < fMinQ2 || Q2 > fMaxQ2) return 0.;
318 
319  int r = -1, i = -1, j = -1;
320 
321  if(ffresid<0 || ffresid >= kNFFRes) return 0.;
322 
323  r = (int)res;
324  if(r<0 || r >= kNRes) return 0.;
325 
326  if (it == kIntWeakCC) { i = 0; }
327  else if (it == kIntWeakNC) { i = 1; }
328  else if (it == kIntEM) { i = 2; }
329 
330  if (hit_nucleon_pdg == kPdgNeutron) { j = 0; }
331  else if (hit_nucleon_pdg == kPdgProton ) { j = 1; }
332 
333  const Spline * spl = fFFRes[r][i][j][ffresid];
334 
335  if(!spl) return 0;
336  else {
337  return spl->Evaluate(Q2);
338  }
339 }
340 //____________________________________________________________________________
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
Singleton to load and serve data tables provided by the GiBUU group.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
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:363
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:96
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
#define pINFO
Definition: Messenger.h:62
double F1V(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
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
const int kPdgProton
Definition: PDGCodes.h:81
#define pNOTICE
Definition: Messenger.h:61
double FP(double Q2, Resonance_t res, int nucleon_pdg, InteractionType_t it) const
bool gAbortingInErr
Definition: Messenger.cxx:34
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...
enum genie::EInteractionType InteractionType_t
void LoadTables(void)
load all form factor data tables
#define pDEBUG
Definition: Messenger.h:63