EffectiveSF.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: Brian Coopersmith, University of Rochester
8 
9  For the class documentation see the corresponding header file.
10 
11 */
12 //____________________________________________________________________________
13 
14 #include <sstream>
15 #include <cstdlib>
16 #include <TMath.h>
17 
30 
31 using std::ostringstream;
32 using namespace genie;
33 using namespace genie::constants;
34 using namespace genie::utils;
35 using namespace genie::utils::config;
36 
37 //____________________________________________________________________________
39 NuclearModelI("genie::EffectiveSF")
40 {
41 
42 }
43 //____________________________________________________________________________
45 NuclearModelI("genie::EffectiveSF", config)
46 {
47 
48 }
49 //____________________________________________________________________________
50 // Cleans up memory from probability distribution maps
51 //____________________________________________________________________________
53 {
54  map<string, TH1D*>::iterator iter = fProbDistroMap.begin();
55  for( ; iter != fProbDistroMap.begin(); ++iter) {
56  TH1D * hst = iter->second;
57  if(hst) {
58  delete hst;
59  hst=0;
60  }
61  }
62  fProbDistroMap.clear();
63 }
64 //____________________________________________________________________________
65 // Set the removal energy, 3 momentum, and FermiMover interaction type
66 //____________________________________________________________________________
68 {
69  assert(target.HitNucIsSet());
71  fCurrMomentum.SetXYZ(0,0,0);
72 
73  //-- set fermi momentum vector
74  //
75 
76  if ( target.A() > 1 ) {
77  TH1D * prob = this->ProbDistro(target);
78  if(!prob) {
79  LOG("EffectiveSF", pNOTICE)
80  << "Null nucleon momentum probability distribution";
81  exit(1);
82  }
83 
84  double p = prob->GetRandom();
85 
86 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
87  LOG("EffectiveSF", pDEBUG) << "|p,nucleon| = " << p;
88 #endif
89 
91 
92  double costheta = -1. + 2. * rnd->RndGen().Rndm();
93  double sintheta = TMath::Sqrt(1.-costheta*costheta);
94  double fi = 2 * kPi * rnd->RndGen().Rndm();
95  double cosfi = TMath::Cos(fi);
96  double sinfi = TMath::Sin(fi);
97 
98  double px = p*sintheta*cosfi;
99  double py = p*sintheta*sinfi;
100  double pz = p*costheta;
101 
102  fCurrMomentum.SetXYZ(px, py, pz);
103 
104  }
105 
106  //-- set removal energy
107  //
108 
109  fCurrRemovalEnergy = this->ReturnBindingEnergy(target);
110  double f1p1h = this->Returnf1p1h(target);
111  // Since TE increases the QE peak via a 2p2h process, we decrease f1p1h
112  // in order to increase the 2p2h interaction to account for this enhancement.
113  f1p1h /= this->GetTransEnh1p1hMod(target);
114  if ( RandomGen::Instance() -> RndGen().Rndm() < f1p1h) {
116  } else if (fEjectSecondNucleon2p2h) {
118  } else {
120  }
121 
122  return true;
123 }
124 //____________________________________________________________________________
125 // Returns the probability of the bin with given momentum. I don't know what w
126 // is supposed to be, but I copied its implementation from Bodek-Ritchie.
127 // Implements the interface.
128 //____________________________________________________________________________
129 double EffectiveSF::Prob(double mom, double w, const Target & target) const
130 {
131  if(w < 0) {
132  TH1D * prob_distr = this->ProbDistro(target);
133  int bin = prob_distr->FindBin(mom);
134  double y = prob_distr->GetBinContent(bin);
135  double dx = prob_distr->GetBinWidth(bin);
136  double prob = y * dx;
137  return prob;
138  }
139  return 1;
140 }
141 //____________________________________________________________________________
142 // Check the map of nucleons to see if we have a probability distribution to
143 // compute with. If not, make one.
144 //____________________________________________________________________________
145 TH1D * EffectiveSF::ProbDistro(const Target & target) const
146 {
147  //-- return stored /if already computed/
148  map<string, TH1D*>::iterator it = fProbDistroMap.find(target.AsString());
149  if(it != fProbDistroMap.end()) return it->second;
150 
151  LOG("EffectiveSF", pNOTICE)
152  << "Computing P = f(p_nucleon) for: " << target.AsString();
153  LOG("EffectiveSF", pNOTICE)
154  << "P(cut-off) = " << fPCutOff << ", P(max) = " << fPMax;
155 
156  //-- get information for the nuclear target
157  int nucleon_pdgc = target.HitNucPdg();
158  assert( pdg::IsProton(nucleon_pdgc) || pdg::IsNeutron(nucleon_pdgc) );
159  return this->MakeEffectiveSF(target);
160 
161 }
162 //____________________________________________________________________________
163 // If transverse enhancement form factor modification is enabled, we must
164 // increase the 2p2h contribution to account for the QE peak enhancement.
165 // This gets that factor based on the target.
166 //____________________________________________________________________________
168  double transEnhMod;
170  fRangeTransEnh1p1hMods, &transEnhMod) &&
171  transEnhMod > 0) {
172  return transEnhMod;
173  }
174  // If none specified, assume no enhancement to quasielastic peak
175  return 1.0;
176 }
177 //____________________________________________________________________________
178 // Makes a momentum distribtuion for the given target using parameters
179 // from the config file.
180 //____________________________________________________________________________
182 {
183  // First check for individually specified nuclei
184  int pdgc = pdg::IonPdgCode(target.A(), target.Z());
185  map<int,vector<double> >::const_iterator it = fProbDistParams.find(pdgc);
186  if(it != fProbDistParams.end()) {
187  vector<double> v = it->second;
188  return this->MakeEffectiveSF(v[0], v[1], v[2], v[3],
189  v[4], v[5], v[6], target);
190  }
191 
192  // Then check in the ranges of A
193  map<pair<int, int>, vector<double> >::const_iterator range_it = fRangeProbDistParams.begin();
194  for(; range_it != fRangeProbDistParams.end(); ++range_it) {
195  if (target.A() >= range_it->first.first && target.A() <= range_it->first.second) {
196  vector<double> v = range_it->second;
197  return this->MakeEffectiveSF(v[0], v[1], v[2], v[3],
198  v[4], v[5], v[6], target);
199  }
200  }
201 
202  return NULL;
203 }
204 //____________________________________________________________________________
205 // Makes a momentum distribution using the factors below (see reference) and
206 // inserts it into the nucleus/momentum distribution map.
207 //____________________________________________________________________________
208 TH1D * EffectiveSF::MakeEffectiveSF(double bs, double bp, double alpha,
209  double beta, double c1, double c2,
210  double c3, const Target & target) const
211 {
212  //-- create the probability distribution
213  int npbins = (int) (1000 * fPMax);
214 
215  TH1D * prob = new TH1D("", "", npbins, 0, fPMax);
216  prob->SetDirectory(0);
217 
218  double dp = fPMax / (npbins-1);
219  for(int i = 0; i < npbins; i++) {
220  double p = i * dp;
221  double y = p / 0.197;
222  double as = c1 * exp(-pow(bs*y,2));
223  double ap = c2 * pow(bp * y, 2) * exp(-pow(bp * y, 2));
224  double at = c3 * pow(y, beta) * exp(-alpha * (y - 2));
225  double rr = (3.14159265 / 4) * (as + ap + at) * pow(y, 2) / 0.197;
226  double dP_dp = rr / 1.01691371;
227  if(p>fPCutOff)
228  dP_dp = 0;
229  assert(dP_dp >= 0);
230  // calculate probability density : dProbability/dp
231 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
232  LOG("EffectiveSF", pDEBUG) << "p = " << p << ", dP/dp = " << dP_dp;
233 #endif
234  prob->Fill(p, dP_dp);
235  }
236 
237  //-- normalize the probability distribution
238  prob->Scale( 1.0 / prob->Integral("width") );
239 
240  //-- store
241  fProbDistroMap.insert(
242  map<string, TH1D*>::value_type(target.AsString(),prob));
243  return prob;
244 }
245 //____________________________________________________________________________
246 // Returns the binding energy for a given nucleus.
247 //____________________________________________________________________________
249 {
250  double binding_en;
251  if (GetValueFromNuclearMaps(target, fNucRmvE, fRangeNucRmvE, &binding_en) &&
252  binding_en > 0) {
253  return binding_en;
254  }
255  return 0;
256 }
257 //____________________________________________________________________________
258 // Returns the fraction of 1p1h events for a given nucleus. All other events
259 // are 2p2h.
260 //____________________________________________________________________________
262 {
263  double f1p1h;
264  if (GetValueFromNuclearMaps(target, f1p1hMap, fRange1p1hMap, &f1p1h) &&
265  f1p1h >= 0 && f1p1h <= 1) {
266  return f1p1h;
267  }
268  return 1;
269 }
270 //____________________________________________________________________________
272 {
273  Algorithm::Configure(config);
274  this->LoadConfig();
275 }
276 //____________________________________________________________________________
277 void EffectiveSF::Configure(string param_set)
278 {
279  Algorithm::Configure(param_set);
280  this->LoadConfig();
281 }
282 //____________________________________________________________________________
283 // Every parameter for this comes from the config files.
284 //____________________________________________________________________________
286 {
287  this->GetParamDef("EjectSecondNucleon2p2h", fEjectSecondNucleon2p2h, false);
288 
289  this->GetParamDef("MomentumMax", fPMax, 1.0);
290  this->GetParamDef("MomentumCutOff", fPCutOff, 0.65);
291  assert(fPMax > 0 && fPCutOff > 0 && fPCutOff <= fPMax);
292 
293  // Find out if Transverse enhancement is enabled to figure out whether to load
294  // the 2p2h enhancement parameters.
295  this->GetParam("UseElFFTransverseEnhancement", fUseElFFTransEnh );
296  if (!fUseElFFTransEnh) {
297  LOG("EffectiveSF", pINFO)
298  << "Transverse enhancement not used; "
299  << "Do not increase the 2p2h cross section.";
300  }
301  else {
302  LoadAllIsotopesForKey("TransEnhf1p1hMod", "EffectiveSF",
304  LoadAllNucARangesForKey("TransEnhf1p1hMod", "EffectiveSF",
306  }
307 
308  LoadAllIsotopesForKey("BindingEnergy", "EffectiveSF", GetOwnedConfig(), &fNucRmvE);
309  LoadAllNucARangesForKey("BindingEnergy", "EffectiveSF",
311  LoadAllIsotopesForKey("f1p1h", "EffectiveSF", GetOwnedConfig(), &f1p1hMap);
312  LoadAllNucARangesForKey("f1p1h", "EffectiveSF", GetOwnedConfig(), &fRange1p1hMap);
313 
314  for (int Z = 1; Z < 140; Z++) {
315  for (int A = Z; A < 3*Z; A++) {
316  const int pdgc = pdg::IonPdgCode(A, Z);
317  double bs, bp, alpha, beta, c1, c2, c3;
318  if (GetDoubleKeyPDG("bs", pdgc, GetOwnedConfig(), &bs) &&
319  GetDoubleKeyPDG("bp", pdgc, GetOwnedConfig(), &bp) &&
320  GetDoubleKeyPDG("alpha", pdgc, GetOwnedConfig(), &alpha) &&
321  GetDoubleKeyPDG("beta", pdgc, GetOwnedConfig(), &beta) &&
322  GetDoubleKeyPDG("c1", pdgc, GetOwnedConfig(), &c1) &&
323  GetDoubleKeyPDG("c2", pdgc, GetOwnedConfig(), &c2) &&
324  GetDoubleKeyPDG("c3", pdgc, GetOwnedConfig(), &c3)) {
325  vector<double> pars = vector<double>();
326  pars.push_back(bs);
327  pars.push_back(bp);
328  pars.push_back(alpha);
329  pars.push_back(beta);
330  pars.push_back(c1);
331  pars.push_back(c2);
332  pars.push_back(c3);
333  LOG("EffectiveSF", pINFO)
334  << "Nucleus: " << pdgc << " -> using bs = " << bs << "; bp = "<< bp
335  << "; alpha = " << alpha << "; beta = "<<beta<<"; c1 = "<<c1
336  <<"; c2 = "<<c2<< "; c3 = " << c3;
337  fProbDistParams[pdgc] = pars;
338  }
339  }
340  }
341  for(int lowA = 1; lowA < 3 * 140; lowA++) {
342  for(int highA = lowA; highA < 3 * 140; highA++) {
343  double bs, bp, alpha, beta, c1, c2, c3;
344  if (GetDoubleKeyRangeNucA("bs", lowA, highA, GetOwnedConfig(), &bs) &&
345  GetDoubleKeyRangeNucA("bp", lowA, highA, GetOwnedConfig(), &bp) &&
346  GetDoubleKeyRangeNucA("alpha",lowA, highA, GetOwnedConfig(), &alpha) &&
347  GetDoubleKeyRangeNucA("beta", lowA, highA, GetOwnedConfig(), &beta) &&
348  GetDoubleKeyRangeNucA("c1", lowA, highA, GetOwnedConfig(), &c1) &&
349  GetDoubleKeyRangeNucA("c2", lowA, highA, GetOwnedConfig(), &c2) &&
350  GetDoubleKeyRangeNucA("c3", lowA, highA, GetOwnedConfig(), &c3)) {
351  vector<double> pars = vector<double>();
352  pars.push_back(bs);
353  pars.push_back(bp);
354  pars.push_back(alpha);
355  pars.push_back(beta);
356  pars.push_back(c1);
357  pars.push_back(c2);
358  pars.push_back(c3);
359  LOG("EffectiveSF", pINFO) << "For "<< lowA - 1 <<" < A < " << highA + 1
360  <<" -> using bs = " << bs << "; bp = "<< bp
361  << "; alpha = " << alpha << "; beta = "<<beta<<"; c1 = "<<c1
362  <<"; c2 = "<<c2<< "; c3 = " << c3;
363  fRangeProbDistParams[pair<int, int>(lowA, highA)] = pars;
364  }
365  }
366  }
367 }
368 //____________________________________________________________________________
const double kPi
Registry * GetOwnedConfig(void)
Definition: Algorithm.cxx:287
void Configure(const Registry &config)
Basic constants.
string AsString(void) const
Definition: Target.cxx:400
const XML_Char * target
Definition: expat.h:268
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
int HitNucPdg(void) const
Definition: Target.cxx:321
map< string, TH1D * > fProbDistroMap
Definition: EffectiveSF.h:72
int A(void) const
Definition: Target.h:71
const char * p
Definition: xmltok.h:285
Simple functions for loading and reading nucleus dependent keys from config files.
constexpr T pow(T x)
Definition: pow.h:75
double Prob(double mom, double w, const Target &t) const
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
Definition: config.py:1
Double_t beta
map< int, std::vector< double > > fProbDistParams
Definition: EffectiveSF.h:81
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
TH1D * MakeEffectiveSF(const Target &target) const
map< pair< int, int >, std::vector< double > > fRangeProbDistParams
Definition: EffectiveSF.h:88
map< int, double > fTransEnh1p1hMods
Definition: EffectiveSF.h:82
virtual ~EffectiveSF()
Definition: EffectiveSF.cxx:52
Float_t Z
Definition: plot.C:38
c2
Definition: demo5.py:33
void LoadAllIsotopesForKey(const char *key_name, const char *log_tool_name, Registry *config, map< int, double > *nuc_to_val)
map< int, double > f1p1hMap
Definition: EffectiveSF.h:80
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
double dx[NP][NC]
map< pair< int, int >, double > fRangeNucRmvE
Definition: EffectiveSF.h:86
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
FermiMoverInteractionType_t fFermiMoverInteractionType
map< pair< int, int >, double > fRange1p1hMap
Definition: EffectiveSF.h:87
double GetTransEnh1p1hMod(const Target &target) const
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int Z(void) const
Definition: Target.h:69
#define pINFO
Definition: Messenger.h:63
double ReturnBindingEnergy(const Target &target) const
float bin[41]
Definition: plottest35.C:14
double Returnf1p1h(const Target &target) const
TH1D * ProbDistro(const Target &t) const
static const double A
Definition: Units.h:82
void LoadConfig(void)
Definition: bp.py:1
bool GetValueFromNuclearMaps(const Target &target, const map< int, double > &nuc_to_val, const map< pair< int, int >, double > &nucA_range_to_val, double *val)
bool HitNucIsSet(void) const
Definition: Target.cxx:300
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
bool GetDoubleKeyRangeNucA(const char *valName, const int lowA, const int highA, Registry *config, double *val)
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
assert(nhit_max >=nhit_nbins)
c1
Definition: demo5.py:24
bool GetDoubleKeyPDG(const char *valName, const int pdgc, Registry *config, double *val)
bool fEjectSecondNucleon2p2h
Definition: EffectiveSF.h:75
map< int, double > fNucRmvE
Definition: EffectiveSF.h:79
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
bool GenerateNucleon(const Target &t) const
Definition: EffectiveSF.cxx:67
Float_t w
Definition: plot.C:20
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
void LoadAllNucARangesForKey(const char *key_name, const char *log_tool_name, Registry *config, map< pair< int, int >, double > *nuc_rangeA_to_val)
Root of GENIE utility namespaces.
map< pair< int, int >, double > fRangeTransEnh1p1hMods
Definition: EffectiveSF.h:89
#define pDEBUG
Definition: Messenger.h:64