ZExpAxialFormFactorModel.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7 \author Aaron Meyer <asmeyer2012 \at uchicago.edu>
8  based off DipoleELFormFactorsModel by
9  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
10  STFC, Rutherford Appleton Laboratory
11 
12  For the class documentation see the corresponding header file.
13 
14 */
15 //____________________________________________________________________________
16 
17 #include <TMath.h>
18 #include <sstream>
19 
25 
26 using std::ostringstream;
27 
28 using namespace genie;
29 
30 //____________________________________________________________________________
32 AxialFormFactorModelI("genie::ZExpAxialFormFactorModel")
33 {
34 
35 }
36 //____________________________________________________________________________
38 AxialFormFactorModelI("genie::ZExpAxialFormFactorModel", config)
39 {
40 
41 }
42 //____________________________________________________________________________
44 {
45  delete[] fZ_An;
46 }
47 //____________________________________________________________________________
48 double ZExpAxialFormFactorModel::FA(const Interaction * interaction) const
49 {
50  // calculate and return FA
51  double q2 = interaction->KinePtr()->q2();
52  double zparam = this->CalculateZ(q2);
53  if (zparam != zparam) // checks for nan
54  {
55  LOG("ZExpAxialFormFactorModel",pWARN) << "Undefined expansion parameter";
56  return 0.;
57  }
58  double fa = 0.;
59  for (int ki=0;ki<=fKmax+(fQ4limit ? 4 : 0);ki++)
60  {
61  fa = fa + TMath::Power(zparam,ki) * fZ_An[ki];
62  }
63 
64  return fa;
65 }
66 //____________________________________________________________________________
68 {
69 
70  // calculate z expansion parameter
71  double znum = TMath::Sqrt(fTcut - q2) - TMath::Sqrt(fTcut - fT0);
72  double zden = TMath::Sqrt(fTcut - q2) + TMath::Sqrt(fTcut - fT0);
73 
74  return znum/zden;
75 }
76 //____________________________________________________________________________
78 {
79  //if (fKmax < 1 ) { fKmax = 1; }
80  //else if (fKmax > 10) { fKmax = 10; }
81 
82  if (fQ4limit) this->FixQ4Limit();
83  else this->FixA0();
84 }
85 //____________________________________________________________________________
87 {
88  // Function to fix form factor such that FA(q2=0) = gA
89  // For T0 = 0, this will set A0 = gA
90  double zparam = this->CalculateZ(0.);
91  double fa = 0.;
92  for (int ki=1;ki<=fKmax;ki++)
93  {
94  fa = fa + TMath::Power(zparam,ki) * fZ_An[ki];
95  }
96  fZ_An[0] = fFA0 - fa;
97 
98 }
99 //____________________________________________________________________________
101 {
102  // fixes parameters such that the model limits to 1/Q^4 at large t
103  // -- requires at least 5 parameters to do so --
104  // 4 parameters for Q^4 behavior, 1 for normalization to FA(q2=0)=gA
105 
106  // will use A_0 and A_Kmax through A_Kmax-3 to do the fitting
107  // calculate some useful numbers (redundancy for readability)
108  double kp4 = (double)fKmax+4;
109  double kp3 = (double)fKmax+3;
110  double kp2 = (double)fKmax+2;
111  double kp1 = (double)fKmax+1;
112  double kp0 = (double)fKmax ;
113  //double km5 = (double)fKmax-5;
114  double z0 = this->CalculateZ(0.);
115  double zkp4 = TMath::Power(z0,(int)kp4);
116  double zkp3 = TMath::Power(z0,(int)kp3);
117  double zkp2 = TMath::Power(z0,(int)kp2);
118  double zkp1 = TMath::Power(z0,(int)kp1);
119 
120  // denominator
121  double denom = \
122  6. - kp4*kp3*kp2*zkp1 + 3.*kp4*kp3*kp1*zkp2 \
123  - 3.*kp4*kp2*kp1*zkp3 + kp3*kp2*kp1*zkp4;
124 
125  // extra parameters (effectively constants)
126  // number refers to the number of derivatives
127  double b0 = 0.;
128  for (int ki = 1;ki <= fKmax;ki++)
129  {
130  b0 = b0 + fZ_An[ki];
131  }
132 
133  double b0z = -fFA0;
134  for (int ki = 1;ki <= fKmax;ki++)
135  {
136  b0z = b0z + fZ_An[ki]*TMath::Power(z0,ki);
137  }
138 
139  double b1 = 0.;
140  for (int ki = 1;ki <= fKmax;ki++)
141  {
142  b1 = b1 + ki*fZ_An[ki];
143  }
144 
145  double b2 = 0.;
146  for (int ki = 1;ki <= fKmax;ki++)
147  {
148  b2 = b2 + ki*(ki-1)*fZ_An[ki];
149  }
150 
151  double b3 = 0.;
152  for (int ki = 1;ki <= fKmax;ki++)
153  {
154  b3 = b3 + ki*(ki-1)*(ki-2)*fZ_An[ki];
155  }
156 
157  // Assign new parameters
158  fZ_An[(int)kp4] = (1./denom) * \
159  ( (b0-b0z)*kp3*kp2*kp1 \
160  + b3*( -1. + .5*kp3*kp2*zkp1 - kp3*kp1*zkp2 \
161  +.5*kp2*kp1*zkp3 ) \
162  + b2*( 3.*kp1 - kp3*kp2*kp1*zkp1 \
163  +kp3*kp1*(2*fKmax+1)*zkp2 - kp2*kp1*kp0*zkp3 ) \
164  + b1*( -3.*kp2*kp1 + .5*kp3*kp2*kp2*kp1*zkp1 \
165  -kp3*kp2*kp1*kp0*zkp2 + .5*kp2*kp1*kp1*kp0*zkp3) );
166 
167  fZ_An[(int)kp3] = (1./denom) * \
168  ( -3.*(b0-b0z)*kp4*kp2*kp1 \
169  + b3*( 3. - kp4*kp2*zkp1 + (3./2.)*kp4*kp1*zkp2 \
170  -.5*kp2*kp1*zkp4 ) \
171  + b2*( -3.*(3*fKmax+4) + kp4*kp2*(2*fKmax+3)*zkp1 \
172  -3.*kp4*kp1*kp1*zkp2 + kp2*kp1*kp0*zkp4 ) \
173  + b1*( 3.*kp1*(3*fKmax+8) - kp4*kp3*kp2*kp1*zkp1 \
174  +(3./2.)*kp4*kp3*kp1*kp0*zkp2 - .5*kp2*kp1*kp1*kp0*zkp4) );
175 
176  fZ_An[(int)kp2] = (1./denom) * \
177  ( 3.*(b0-b0z)*kp4*kp3*kp1 \
178  + b3*( -3. + .5*kp4*kp3*zkp1 - (3./2.)*kp4*kp1*zkp3 \
179  +kp3*kp1*zkp4 ) \
180  + b2*( 3.*(3*fKmax+5) - kp4*kp3*kp2*zkp1 \
181  +3.*kp4*kp1*kp1*zkp3 - kp3*kp1*(2*fKmax+1)*zkp4) \
182  + b1*( -3.*kp3*(3*fKmax+4) + .5*kp4*kp3*kp3*kp2*zkp1 \
183  -(3./2.)*kp4*kp3*kp1*kp0*zkp3 + kp3*kp2*kp1*kp0*zkp4) );
184 
185  fZ_An[(int)kp1] = (1./denom) * \
186  ( -(b0-b0z)*kp4*kp3*kp2 \
187  + b3*( 1. - .5*kp4*kp3*zkp2 + kp4*kp2*zkp3 \
188  -.5*kp3*kp2*zkp4 ) \
189  + b2*( -3.*kp2 + kp4*kp3*kp2*zkp2 \
190  -kp4*kp2*(2*fKmax+3)*zkp3 + kp3*kp2*kp1*zkp4) \
191  + b1*( 3.*kp3*kp2 - .5*kp4*kp3*kp3*kp2*zkp2 \
192  +kp4*kp3*kp2*kp1*zkp3 - .5*kp3*kp2*kp2*kp1*zkp4) );
193 
194  fZ_An[0] = (1./denom) * \
195  ( -6.*b0z \
196  + b0*( kp4*kp3*kp2*zkp1 - 3.*kp4*kp3*kp1*zkp2 \
197  +3.*kp4*kp2*kp1*zkp3 - kp3*kp2*kp1*zkp4 ) \
198  + b3*( -zkp1 + 3.*zkp2 - 3.*zkp3 + zkp4 ) \
199  + b2*( 3.*kp2*zkp1 - 3.*(3*fKmax+5)*zkp2 \
200  +3.*(3*fKmax+4)*zkp3 - 3.*kp1*zkp4 ) \
201  + b1*( -3.*kp3*kp2*zkp1 + 3.*kp3*(3*fKmax+4)*zkp2 \
202  -3.*kp1*(3*fKmax+8)*zkp3 + 3.*kp2*kp1*zkp4 ) );
203 }
204 //____________________________________________________________________________
206 {
207  Algorithm::Configure(config);
208  this->LoadConfig();
209 }
210 //____________________________________________________________________________
212 {
213  Algorithm::Configure(param_set);
214  this->LoadConfig();
215 }
216 //____________________________________________________________________________
218 {
219 // get config options from the configuration registry or set defaults
220 // from the global parameter list
221 
222  GetParam( "QEL-Q4limit", fQ4limit ) ;
223  GetParam( "QEL-Kmax", fKmax ) ;
224 
225  GetParam( "QEL-T0", fT0 ) ;
226  GetParam( "QEL-T0", fT0 ) ;
227  GetParam( "QEL-Tcut", fTcut ) ;
228 
229  GetParam( "QEL-FA0", fFA0 ) ;
230  assert(fKmax > 0);
231 
232  // z expansion coefficients
233  if (fQ4limit) fZ_An = new double [fKmax+5];
234  else fZ_An = new double [fKmax+1];
235 
236  // load the user-defined coefficient values
237  // -- A0 and An for n<fKmax are calculated from other means
238  for (int ip=1;ip<fKmax+1;ip++) {
239  ostringstream alg_key;
240  alg_key << "QEL-Z_A" << ip;
241  GetParam( alg_key.str(), fZ_An[ip] ) ;
242  }
243 
244  this->FixCoeffs();
245 }
246 //____________________________________________________________________________
247 
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Pure abstract base class. Defines the AxialFormFactorModelI interface to be implemented by LlewellynS...
double FA(const Interaction *interaction) const
Compute the axial form factor.
Definition: config.py:1
TString ip
Definition: loadincs.C:5
Double_t q2[12][num]
Definition: f2_nu.C:137
Summary information for an interaction.
Definition: Interaction.h:56
double q2(bool selected=false) const
Definition: Kinematics.cxx:151
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
void Configure(const Registry &config)
#define pWARN
Definition: Messenger.h:61
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
assert(nhit_max >=nhit_nbins)
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
TFile fa("Li7.root")