LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Static Public Member Functions | Private Types | Private Attributes | List of all members
evwgh::GenieWeightCalc Class Reference
Inheritance diagram for evwgh::GenieWeightCalc:
evwgh::WeightCalc

Public Member Functions

 GenieWeightCalc ()
 
void Configure (fhicl::ParameterSet const &pset)
 
std::vector< std::vector< double > > GetWeight (art::Event &e)
 
void SetName (std::string name)
 
std::string GetName ()
 

Static Public Member Functions

static std::vector< std::vector< double > > MultiGaussianSmearing (std::vector< double > const &centralValues, std::vector< std::vector< double >> const &inputCovarianceMatrix, int n_multisims, CLHEP::RandGaussQ &GaussRandom)
 Applies Gaussian smearing to a set of data. More...
 
static std::vector< double > MultiGaussianSmearing (std::vector< double > const &centralValue, TMatrixD *const &inputCovarianceMatrix, std::vector< double > rand)
 Over load of the above function that only returns a single varied parameter set. More...
 
static std::vector< double > MultiGaussianSmearing (std::vector< double > const &centralValue, TMatrixD *const &LowerTriangleCovarianceMatrix, bool isDecomposed, std::vector< double > rand)
 

Private Types

enum  EReweight {
  kNCELaxial, kNCELeta, kQEMA, kQEVec,
  kCCResAxial, kCCResVector, kResGanged, kNCResAxial,
  kNCResVector, kCohMA, kCohR0, kNonResRvp1pi,
  kNonResRvbarp1pi, kNonResRvp2pi, kNonResRvbarp2pi, kResDecayGamma,
  kResDecayEta, kResDecayTheta, kNC, kDISAth,
  kDISBth, kDISCv1u, kDISCv2u, kDISnucl,
  kAGKYxF, kAGKYpT, kFormZone, kFermiGasModelKf,
  kFermiGasModelSf, kIntraNukeNmfp, kIntraNukeNcex, kIntraNukeNel,
  kIntraNukeNinel, kIntraNukeNabs, kIntraNukeNpi, kIntraNukePImfp,
  kIntraNukePIcex, kIntraNukePIel, kIntraNukePIinel, kIntraNukePIabs,
  kIntraNukePIpi, kNReWeights
}
 

Private Attributes

std::vector< rwgt::NuReweight * > reweightVector
 
CLHEP::RandGaussQ * fGaussRandom
 
std::string fGenieModuleLabel
 

Detailed Description

Definition at line 24 of file GenieWeightCalc.cxx.

Member Enumeration Documentation

enum evwgh::GenieWeightCalc::EReweight
private
Enumerator
kNCELaxial 
kNCELeta 
kQEMA 
kQEVec 
kCCResAxial 
kCCResVector 
kResGanged 
kNCResAxial 
kNCResVector 
kCohMA 
kCohR0 
kNonResRvp1pi 
kNonResRvbarp1pi 
kNonResRvp2pi 
kNonResRvbarp2pi 
kResDecayGamma 
kResDecayEta 
kResDecayTheta 
kNC 
kDISAth 
kDISBth 
kDISCv1u 
kDISCv2u 
kDISnucl 
kAGKYxF 
kAGKYpT 
kFormZone 
kFermiGasModelKf 
kFermiGasModelSf 
kIntraNukeNmfp 
kIntraNukeNcex 
kIntraNukeNel 
kIntraNukeNinel 
kIntraNukeNabs 
kIntraNukeNpi 
kIntraNukePImfp 
kIntraNukePIcex 
kIntraNukePIel 
kIntraNukePIinel 
kIntraNukePIabs 
kIntraNukePIpi 
kNReWeights 

Definition at line 42 of file GenieWeightCalc.cxx.

42  {kNCELaxial, // Axial mass for NC elastic
43  kNCELeta, // Strange axial form factor for NC elastic
44  kQEMA, // Axial mass for CC quasi-elastic
45  kQEVec, // Choice of CCQE vector form factor (sigma = 0 => BBA05; sigma = 1 => Dipole)
46  kCCResAxial, // Axial mass for CC resonance neutrino production
47  kCCResVector, // Vector mass for CC resonance neutrino production
48  kResGanged, // CC Res && NC Res (NOT ACTIVE)
49  kNCResAxial, // Axial mass for NC resonance neutrino production
50  kNCResVector, // Vector mass for NC resonance neutrino production
51  kCohMA, // Axial mass for CC and NC coherent pion production
52  kCohR0, // Nuclear size param. controlling pi absorption in Rein-Sehgal model
53  kNonResRvp1pi, // v+p and vbar + n (1 pi) type interactions (*)
54  kNonResRvbarp1pi, // v+n and vbar + p (1 pi) type interactions (*)
55  kNonResRvp2pi, // v+p and vbar + n (2 pi) type interactions (*)
56  kNonResRvbarp2pi, // v+n and vbar + p (2 pi) type interactions (*)
57  kResDecayGamma, // BR for radiative resonance decay
58  kResDecayEta, // BR for single-eta resonance decay
59  kResDecayTheta, // Pion angular distibution in Delta -> pi N (sigma = 0 => isotropic; sigma = 1 => RS)
60  kNC,
61  kDISAth, // Ath higher twist param in BY model scaling variable xi_w
62  kDISBth, // Bth higher twist param in BY model scaling variable xi_w
63  kDISCv1u, // Cv1u u valence GRV98 PDF correction param in BY model
64  kDISCv2u, // Cv2u u valence GRV98 PDF correction param in BY model
65  kDISnucl, // NOT IMPLEMENTED IN GENIE
66  kAGKYxF, // Pion Feynman x for Npi states in AGKY
67  kAGKYpT, // Pion transverse momentum for Npi states in AGKY
68  kFormZone, // Hadron Formation Zone
69  kFermiGasModelKf, // CCQE Pauli Suppression via changes in Fermi level kF
70  kFermiGasModelSf, // Choice of model (sigma = 0 => FermiGas; sigma = 1 => SF (spectral function))
71  kIntraNukeNmfp, // Nucleon mean free path (total rescattering probability)
72  kIntraNukeNcex, // Nucleon charge exchange probability
73  kIntraNukeNel, // Nucleon elastic reaction probability
74  kIntraNukeNinel, // Nucleon inelastic reaction probability
75  kIntraNukeNabs, // Nucleon absorption probability
76  kIntraNukeNpi, // Nucleon pi-production probability
77  kIntraNukePImfp, // Pi mean free path (total rescattering probability)
78  kIntraNukePIcex, // Pi charge exchange probability
79  kIntraNukePIel, // Pi elastic reaction probability
80  kIntraNukePIinel, // Pi inelastic reaction probability
81  kIntraNukePIabs, // Pi absorption probability
82  kIntraNukePIpi, // Pi pi-production probability
83  kNReWeights}; // ?

Constructor & Destructor Documentation

evwgh::GenieWeightCalc::GenieWeightCalc ( )

Definition at line 88 of file GenieWeightCalc.cxx.

89  {}

Member Function Documentation

void evwgh::GenieWeightCalc::Configure ( fhicl::ParameterSet const &  pset)
virtual

Implements evwgh::WeightCalc.

Definition at line 91 of file GenieWeightCalc.cxx.

References fGaussRandom, fGenieModuleLabel, rwgt::fReweightFrAbs_N, rwgt::fReweightFrAbs_pi, rwgt::fReweightFrCEx_N, rwgt::fReweightFrCEx_pi, rwgt::fReweightFrElas_N, rwgt::fReweightFrElas_pi, rwgt::fReweightFrInel_N, rwgt::fReweightFrInel_pi, rwgt::fReweightFrPiProd_N, rwgt::fReweightFrPiProd_pi, rwgt::fReweightMFP_N, rwgt::fReweightMFP_pi, fhicl::ParameterSet::get(), art::RandomNumberGenerator::getEngine(), evwgh::WeightCalc::GetName(), kAGKYpT, kAGKYxF, kCCResAxial, kCCResVector, kCohMA, kCohR0, kDISAth, kDISBth, kDISCv1u, kDISCv2u, kDISnucl, kFermiGasModelKf, kFermiGasModelSf, kFormZone, kIntraNukeNabs, kIntraNukeNcex, kIntraNukeNel, kIntraNukeNinel, kIntraNukeNmfp, kIntraNukeNpi, kIntraNukePIabs, kIntraNukePIcex, kIntraNukePIel, kIntraNukePIinel, kIntraNukePImfp, kIntraNukePIpi, kNC, kNCELaxial, kNCELeta, kNCResAxial, kNCResVector, kNonResRvbarp1pi, kNonResRvbarp2pi, kNonResRvp1pi, kNonResRvp2pi, kNReWeights, kQEMA, kQEVec, kResDecayEta, kResDecayGamma, kResDecayTheta, kResGanged, reweightVector, and s.

92  {
93  // Global Config
94  fGenieModuleLabel = p.get<std::string> ("genie_module_label");
95  fhicl::ParameterSet const &pset = p.get<fhicl::ParameterSet> (GetName());
96 
97  // Calc Config
98  std::vector<std::string> pars = pset.get<std::vector<std::string>> ("parameter_list");
99  std::vector<float> parsigmas = pset.get<std::vector<float>> ("parameter_sigma");
100  std::string mode = pset.get<std::string>("mode");
101 
102  if (pars.size() != parsigmas.size())
103  throw cet::exception(__PRETTY_FUNCTION__) << GetName()
104  << "::Bad fcl configuration. parameter_list and parameter_sigma need to have same number of parameters." << std::endl;
105 
106  int number_of_multisims = pset.get<int> ("number_of_multisims");
107 
108  std::vector<EReweight> erwgh;
109  for (auto & s : pars) {
110  if (s == "NCELaxial") erwgh.push_back(kNCELaxial);
111  else if (s == "NCELeta") erwgh.push_back(kNCELeta);
112  else if (s == "QEMA") erwgh.push_back(kQEMA);
113  else if (s == "QEVec") erwgh.push_back(kQEVec);
114  else if (s == "CCResAxial") erwgh.push_back(kCCResAxial);
115  else if (s == "CCResVector") erwgh.push_back(kCCResVector);
116  else if (s == "ResGanged") erwgh.push_back(kResGanged);
117  else if (s == "NCResAxial") erwgh.push_back(kNCResAxial);
118  else if (s == "NCResVector") erwgh.push_back(kNCResVector);
119  else if (s == "CohMA") erwgh.push_back(kCohMA);
120  else if (s == "CohR0") erwgh.push_back(kCohR0);
121  else if (s == "NonResRvp1pi") erwgh.push_back(kNonResRvp1pi);
122  else if (s == "NonResRvbarp1pi") erwgh.push_back(kNonResRvbarp1pi);
123  else if (s == "NonResRvp2pi") erwgh.push_back(kNonResRvp2pi);
124  else if (s == "NonResRvbarp2pi") erwgh.push_back(kNonResRvbarp2pi);
125  else if (s == "ResDecayGamma") erwgh.push_back(kResDecayGamma);
126  else if (s == "ResDecayEta") erwgh.push_back(kResDecayEta);
127  else if (s == "ResDecayTheta") erwgh.push_back(kResDecayTheta);
128  else if (s == "NC") erwgh.push_back(kNC);
129  else if (s == "DISAth") erwgh.push_back(kDISAth);
130  else if (s == "DISBth") erwgh.push_back(kDISBth);
131  else if (s == "DISCv1u") erwgh.push_back(kDISCv1u);
132  else if (s == "DISCv2u") erwgh.push_back(kDISCv2u);
133  else if (s == "DISnucl") erwgh.push_back(kDISnucl);
134  else if (s == "AGKYxF") erwgh.push_back(kAGKYxF);
135  else if (s == "AGKYpT") erwgh.push_back(kAGKYpT);
136  else if (s == "FormZone") erwgh.push_back(kFormZone);
137  else if (s == "FermiGasModelKf") erwgh.push_back(kFermiGasModelKf);
138  else if (s == "FermiGasModelSf") erwgh.push_back(kFermiGasModelSf);
139  else if (s == "IntraNukeNmfp") erwgh.push_back(kIntraNukeNmfp);
140  else if (s == "IntraNukeNcex") erwgh.push_back(kIntraNukeNcex);
141  else if (s == "IntraNukeNel") erwgh.push_back(kIntraNukeNel);
142  else if (s == "IntraNukeNinel") erwgh.push_back(kIntraNukeNinel);
143  else if (s == "IntraNukeNabs") erwgh.push_back(kIntraNukeNabs);
144  else if (s == "IntraNukeNpi") erwgh.push_back(kIntraNukeNpi);
145  else if (s == "IntraNukePImfp") erwgh.push_back(kIntraNukePImfp);
146  else if (s == "IntraNukePIcex") erwgh.push_back(kIntraNukePIcex);
147  else if (s == "IntraNukePIel") erwgh.push_back(kIntraNukePIel);
148  else if (s == "IntraNukePIinel") erwgh.push_back(kIntraNukePIinel);
149  else if (s == "IntraNukePIabs") erwgh.push_back(kIntraNukePIabs);
150  else if (s == "IntraNukePIpi") erwgh.push_back(kIntraNukePIpi);
151  else {
152  throw cet::exception(__PRETTY_FUNCTION__) << GetName()
153  << "::Physical process " << s << " you requested is not available to reweight." << std::endl;
154  }
155  }
156 
157  //Prepare sigmas
158  art::ServiceHandle<art::RandomNumberGenerator> rng;
159  fGaussRandom = new CLHEP::RandGaussQ(rng->getEngine(GetName()));
160 
161  std::vector<std::vector<float>> reweightingSigmas(erwgh.size());
162 
163  if (mode.find("pm1sigma") != std::string::npos ) {
164  number_of_multisims = 2; // only +-1 sigma if pm1sigma is specified
165  }
166  for (unsigned int i = 0; i < reweightingSigmas.size(); ++i) {
167  reweightingSigmas[i].resize(number_of_multisims);
168  for (int j = 0; j < number_of_multisims; j ++) {
169  if (mode.find("multisim") != std::string::npos )
170  reweightingSigmas[i][j] = parsigmas[i]*fGaussRandom->shoot(&rng->getEngine(GetName()),0.,1.);
171  else if (mode.find("pm1sigma") != std::string::npos )
172  reweightingSigmas[i][j] = (j == 0 ? 1.: -1.); // j==0 => 1; j==1 => -1 if pm1sigma is specified
173  else
174  reweightingSigmas[i][j] = parsigmas[i];
175  }
176  }
177 
178  reweightVector.resize(number_of_multisims);
179 
180  for (int weight_point = 0;
181  weight_point < number_of_multisims;
182  weight_point++){
183 
184  reweightVector[weight_point] = new rwgt::NuReweight;
185 
186  for (unsigned int i_reweightingKnob=0;i_reweightingKnob<erwgh.size();i_reweightingKnob++) {
187  std::cout<<GetName() << "::Setting up rwgh " <<weight_point << "\t" << i_reweightingKnob << "\t" << erwgh[i_reweightingKnob] << std::endl;
188 
189  switch (erwgh[i_reweightingKnob]){
190 
191  case kNCELaxial:
192  reweightVector[weight_point] -> ReweightNCEL(reweightingSigmas[i_reweightingKnob][weight_point], 0.);
193  break;
194  case kNCELeta:
195  reweightVector[weight_point] -> ReweightNCEL(0., reweightingSigmas[i_reweightingKnob][weight_point]);
196  break;
197 
198  case kQEMA:
199  reweightVector[weight_point]
200  -> ReweightQEMA(reweightingSigmas[i_reweightingKnob][weight_point]);
201  break;
202 
203  case kQEVec:
204  reweightVector[weight_point]
205  -> ReweightQEVec(reweightingSigmas[i_reweightingKnob][weight_point]);
206  break;
207 
208  case kResGanged:
209  //reweightVector[weight_point]
210  // -> ReweightResGanged(reweightingSigmas[i_reweightingKnob][weight_point]);
211  break;
212 
213  case kCCResAxial:
214  reweightVector[weight_point] -> ReweightCCRes(reweightingSigmas[i_reweightingKnob][weight_point], 0.);
215  break;
216  case kCCResVector:
217  reweightVector[weight_point] -> ReweightCCRes(0., reweightingSigmas[i_reweightingKnob][weight_point]);
218  break;
219 
220  case kNCResAxial:
221  reweightVector[weight_point] -> ReweightNCRes(reweightingSigmas[i_reweightingKnob][weight_point], 0.);
222  break;
223  case kNCResVector:
224  reweightVector[weight_point] -> ReweightNCRes(0., reweightingSigmas[i_reweightingKnob][weight_point]);
225  break;
226 
227  case kCohMA:
228  reweightVector[weight_point] -> ReweightCoh(reweightingSigmas[i_reweightingKnob][weight_point], 0.);
229  break;
230  case kCohR0:
231  reweightVector[weight_point] -> ReweightCoh(0., reweightingSigmas[i_reweightingKnob][weight_point]);
232  break;
233 
234 
235  case kNonResRvp1pi:
236  reweightVector[weight_point]
237  -> ReweightNonResRvp1pi(reweightingSigmas[i_reweightingKnob][weight_point]);
238  break;
239  case kNonResRvbarp1pi:
240  reweightVector[weight_point]
241  -> ReweightNonResRvbarp1pi(reweightingSigmas[i_reweightingKnob][weight_point]);
242  break;
243  case kNonResRvp2pi:
244  reweightVector[weight_point]
245  -> ReweightNonResRvp2pi(reweightingSigmas[i_reweightingKnob][weight_point]);
246  break;
247  case kNonResRvbarp2pi:
248  reweightVector[weight_point]
249  -> ReweightNonResRvbarp2pi(reweightingSigmas[i_reweightingKnob][weight_point]);
250  break;
251 
252  case kResDecayGamma:
253  reweightVector[weight_point] -> ReweightResDecay(reweightingSigmas[i_reweightingKnob][weight_point], 0., 0.);
254  break;
255  case kResDecayEta:
256  reweightVector[weight_point] -> ReweightResDecay(0., reweightingSigmas[i_reweightingKnob][weight_point], 0.);
257  break;
258  case kResDecayTheta:
259  reweightVector[weight_point] -> ReweightResDecay(0., 0., reweightingSigmas[i_reweightingKnob][weight_point]);
260  break;
261 
262  case kNC:
263  reweightVector[weight_point]
264  -> ReweightNC(reweightingSigmas[i_reweightingKnob][weight_point]);
265  break;
266 
267  case kDISAth:
268  reweightVector[weight_point] -> ReweightDIS(reweightingSigmas[i_reweightingKnob][weight_point], 0., 0., 0.);
269  break;
270  case kDISBth:
271  reweightVector[weight_point] -> ReweightDIS(0., reweightingSigmas[i_reweightingKnob][weight_point], 0., 0.);
272  break;
273  case kDISCv1u:
274  reweightVector[weight_point] -> ReweightDIS(0., 0., reweightingSigmas[i_reweightingKnob][weight_point], 0.);
275  break;
276  case kDISCv2u:
277  reweightVector[weight_point] -> ReweightDIS(0., 0., 0., reweightingSigmas[i_reweightingKnob][weight_point]);
278  break;
279 
280  case kDISnucl:
281  reweightVector[weight_point]
282  -> ReweightDISnucl(reweightingSigmas[i_reweightingKnob][weight_point]);
283  break;
284 
285  case kAGKYxF:
286  reweightVector[weight_point] -> ReweightAGKY(reweightingSigmas[i_reweightingKnob][weight_point], 0.);
287  break;
288  case kAGKYpT:
289  reweightVector[weight_point] -> ReweightAGKY(0., reweightingSigmas[i_reweightingKnob][weight_point]);
290  break;
291 
292  case kFormZone:
293  reweightVector[weight_point] -> ReweightFormZone(reweightingSigmas[i_reweightingKnob][weight_point]);
294  break;
295 
296  case kFermiGasModelKf:
297  reweightVector[weight_point] -> ReweightFGM(reweightingSigmas[i_reweightingKnob][weight_point], 0.);
298  break;
299  case kFermiGasModelSf:
300  reweightVector[weight_point] -> ReweightFGM(0., reweightingSigmas[i_reweightingKnob][weight_point]);
301  break;
302 
303  case kIntraNukeNmfp:
304  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightMFP_N, reweightingSigmas[i_reweightingKnob][weight_point]);
305  break;
306  case kIntraNukeNcex:
307  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrCEx_N, reweightingSigmas[i_reweightingKnob][weight_point]);
308  break;
309  case kIntraNukeNel:
310  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrElas_N, reweightingSigmas[i_reweightingKnob][weight_point]);
311  break;
312  case kIntraNukeNinel:
313  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrInel_N, reweightingSigmas[i_reweightingKnob][weight_point]);
314  break;
315  case kIntraNukeNabs:
316  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrAbs_N, reweightingSigmas[i_reweightingKnob][weight_point]);
317  break;
318  case kIntraNukeNpi:
319  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrPiProd_N, reweightingSigmas[i_reweightingKnob][weight_point]);
320  break;
321  case kIntraNukePImfp:
322  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightMFP_pi, reweightingSigmas[i_reweightingKnob][weight_point]);
323  break;
324  case kIntraNukePIcex:
325  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrCEx_pi, reweightingSigmas[i_reweightingKnob][weight_point]);
326  break;
327  case kIntraNukePIel:
328  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrElas_pi, reweightingSigmas[i_reweightingKnob][weight_point]);
329  break;
330  case kIntraNukePIinel:
331  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrInel_pi, reweightingSigmas[i_reweightingKnob][weight_point]);
332  break;
333  case kIntraNukePIabs:
334  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrAbs_pi, reweightingSigmas[i_reweightingKnob][weight_point]);
335  break;
336  case kIntraNukePIpi:
337  reweightVector[weight_point] -> ReweightIntraNuke(rwgt::fReweightFrPiProd_pi, reweightingSigmas[i_reweightingKnob][weight_point]);
338  break;
339 
340  case kNReWeights:
341  break;
342  }
343  }
344  } // loop over nWeights
345 
346  // Tell all of the reweight drivers to configure themselves:
347  std::cout<< GetName()<<"::Setting up "<<reweightVector.size()<<" reweightcalcs"<<std::endl;
348  for(auto & driver : reweightVector){
349  driver -> Configure();
350  }
351 
352  }
rwgt::fReweightFrInel_pi
tweak inelastic probability for pions, for given total rescattering probability
Definition: ReweightLabels.h:104
s
Float_t s
Definition: plot.C:23
rwgt::fReweightFrAbs_N
tweak absorption probability for nucleons, for given total rescattering probability ...
Definition: ReweightLabels.h:110
rwgt::fReweightFrInel_N
tweak inelastic probability for nucleons, for given total rescattering probability ...
Definition: ReweightLabels.h:109
rwgt::fReweightFrPiProd_N
tweak pion production probability for nucleons, for given total rescattering probability ...
Definition: ReweightLabels.h:111
evwgh::WeightCalc::GetName
std::string GetName()
Definition: WeightCalc.h:24
art::RandomNumberGenerator::getEngine
base_engine_t & getEngine() const
rwgt::fReweightFrElas_N
tweak elastic probability for nucleons, for given total rescattering probability
Definition: ReweightLabels.h:108
evwgh::GenieWeightCalc::fGaussRandom
CLHEP::RandGaussQ * fGaussRandom
Definition: GenieWeightCalc.cxx:35
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
rwgt::fReweightFrAbs_pi
tweak absorption probability for pions, for given total rescattering probability
Definition: ReweightLabels.h:105
rwgt::fReweightMFP_pi
tweak mean free path for pions
Definition: ReweightLabels.h:100
evwgh::GenieWeightCalc::Configure
void Configure(fhicl::ParameterSet const &pset)
Definition: GenieWeightCalc.cxx:91
rwgt::fReweightFrCEx_N
tweak charge exchange probability for nucleons, for given total rescattering probability ...
Definition: ReweightLabels.h:107
rwgt::fReweightFrPiProd_pi
tweak pion production probability for pions, for given total rescattering probability ...
Definition: ReweightLabels.h:106
rwgt::fReweightFrCEx_pi
tweak charge exchange probability for pions, for given total rescattering probability ...
Definition: ReweightLabels.h:102
rwgt::fReweightMFP_N
tweak mean free path for nucleons
Definition: ReweightLabels.h:101
rwgt::fReweightFrElas_pi
tweak elastic probability for pions, for given total rescattering probability
Definition: ReweightLabels.h:103
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
evwgh::GenieWeightCalc::reweightVector
std::vector< rwgt::NuReweight * > reweightVector
Definition: GenieWeightCalc.cxx:33
std::string evwgh::WeightCalc::GetName ( )
inlineinherited

Definition at line 24 of file WeightCalc.h.

References evwgh::WeightCalc::fName, evwgh::WeightCalc::MultiGaussianSmearing(), and lar::dump::vector().

Referenced by Configure().

24 {return fName;}
evwgh::WeightCalc::fName
std::string fName
Definition: WeightCalc.h:54
std::vector< std::vector< double > > evwgh::GenieWeightCalc::GetWeight ( art::Event &  e)
virtual

Implements evwgh::WeightCalc.

Definition at line 354 of file GenieWeightCalc.cxx.

References fGenieModuleLabel, art::fill_ptr_vector(), art::DataViewImpl::getByLabel(), REGISTER_WEIGHTCALC, reweightVector, and weight.

355  {
356  // Returns a vector of weights for each neutrino interaction in the event
357 
358  // Get the MC generator information out of the event
359  // These are all handles to mc information.
360  art::Handle< std::vector<simb::MCTruth> > mcTruthHandle;
361  art::Handle< std::vector<simb::MCFlux> > mcFluxHandle;
362  art::Handle< std::vector<simb::GTruth> > gTruthHandle;
363 
364  // Actually go and get the stuff
365  e.getByLabel(fGenieModuleLabel,mcTruthHandle);
366  e.getByLabel(fGenieModuleLabel,mcFluxHandle);
367  e.getByLabel(fGenieModuleLabel,gTruthHandle);
368 
369  std::vector<art::Ptr<simb::MCTruth> > mclist;
370  art::fill_ptr_vector(mclist, mcTruthHandle);
371 
372  std::vector<art::Ptr<simb::GTruth > > glist;
373  art::fill_ptr_vector(glist, gTruthHandle);
374 
375  // Calculate weight(s) here
376  std::vector<std::vector<double> >weight(mclist.size());
377  for ( unsigned int inu=0; inu<mclist.size();inu++) {
378  weight[inu].resize(reweightVector.size());
379  for (unsigned int i_weight = 0;
380  i_weight < reweightVector.size();
381  i_weight ++){
382  weight[inu][i_weight]= reweightVector[i_weight]-> CalcWeight(*mclist[inu],*glist[inu]);
383  }
384  }
385  return weight;
386 
387  }
weight
double weight
Definition: plottest35.C:25
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
evwgh::GenieWeightCalc::reweightVector
std::vector< rwgt::NuReweight * > reweightVector
Definition: GenieWeightCalc.cxx:33
static std::vector<std::vector<double> > evwgh::WeightCalc::MultiGaussianSmearing ( std::vector< double > const &  centralValues,
std::vector< std::vector< double >> const &  inputCovarianceMatrix,
int  n_multisims,
CLHEP::RandGaussQ &  GaussRandom 
)
staticinherited

Applies Gaussian smearing to a set of data.

Parameters
centralValuesthe values to be smeared
inputCovarianceMatrixcovariance matrix for smearing
n_multisimsnumber of sets of smeared values to be produced
Returns
a set of n_multisims value sets smeared from the central value

If centralValues is of dimension N, inputCovarianceMatrix needs to be NxN, and each of the returned data sets will be also of dimension N.

Referenced by evwgh::WeightCalc::GetName().

std::vector< double > evwgh::WeightCalc::MultiGaussianSmearing ( std::vector< double > const &  centralValue,
TMatrixD *const &  inputCovarianceMatrix,
std::vector< double >  rand 
)
staticinherited

Over load of the above function that only returns a single varied parameter set.

Definition at line 106 of file WeightCalc.cxx.

References col, and art::errors::StdException.

107  {
108 
109  //compute the smeared central values
110  std::vector<double> smearedCentralValues;
111 
112  //
113  //perform Choleskey Decomposition
114  //
115  // Best description I have found
116  // is in the PDG (Monte Carlo techniques, Algorithms, Gaussian distribution)
117  //
118  // http://pdg.lbl.gov/2016/reviews/rpp2016-rev-monte-carlo-techniques.pdf (Page 5)
119  //
120  //
121  TDecompChol dc = TDecompChol(*(inputCovarianceMatrix));
122  if(!dc.Decompose())
123  {
124  throw art::Exception(art::errors::StdException)
125  << "Cannot decompose covariance matrix to begin smearing.";
126  return smearedCentralValues;
127  }
128 
129  //Get upper triangular matrix. This maintains the relations in the
130  //covariance matrix, but simplifies the structure.
131  TMatrixD U = dc.GetU();
132 
133 
134  for(unsigned int col = 0; col < centralValue.size(); ++col)
135  {
136  //find the weight of each col of the upper triangular cov. matrix
137  double weightFromU = 0.;
138 
139  for(unsigned int row = 0; row < col+1; ++row)
140  {
141  weightFromU += U(row,col)*rand[row];
142  }
143 
144  //multiply this weight by each col of the central values to obtain
145  //the gaussian smeared and constrained central values
146  // smearedCentralValues.push_back(weightFromU * centralValue[col]);
147  smearedCentralValues.push_back(weightFromU + centralValue[col]);
148  }
149  return smearedCentralValues;
150  } // WeightCalc::MultiGaussianSmearing() - Gives a single set of smeared parameters
col
Int_t col[ntarg]
Definition: Style.C:29
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::vector< double > evwgh::WeightCalc::MultiGaussianSmearing ( std::vector< double > const &  centralValue,
TMatrixD *const &  LowerTriangleCovarianceMatrix,
bool  isDecomposed,
std::vector< double >  rand 
)
staticinherited

Definition at line 163 of file WeightCalc.cxx.

References col, and art::errors::StdException.

164  {
165 
166  //compute the smeared central values
167  std::vector<double> smearedCentralValues;
168 
169  // This guards against accidentally
170  if(!isDecomposed){
171  throw art::Exception(art::errors::StdException)
172  << "Must supply the decomposed lower triangular covariance matrix.";
173  return smearedCentralValues;
174  }
175 
176 
177  for(unsigned int col = 0; col < centralValue.size(); ++col)
178  {
179  //find the weight of each col of the upper triangular cov. matrix
180  double weightFromU = 0.;
181 
182  for(unsigned int row = 0; row < col+1; ++row)
183  {
184  weightFromU += LowerTriangleCovarianceMatrix[0][row][col]*rand[row];
185  }
186 
187  //multiply this weight by each col of the central values to obtain
188  //the gaussian smeared and constrained central values
189  // smearedCentralValues.push_back(weightFromU * centralValue[col]);
190  smearedCentralValues.push_back(weightFromU + centralValue[col]);
191  }
192  return smearedCentralValues;
193  } // WeightCalc::MultiGaussianSmearing() - Gives a single set of smeared parameters
col
Int_t col[ntarg]
Definition: Style.C:29
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
void evwgh::WeightCalc::SetName ( std::string  name)
inlineinherited

Definition at line 23 of file WeightCalc.h.

References evwgh::WeightCalc::fName.

23 {fName=name;}
evwgh::WeightCalc::fName
std::string fName
Definition: WeightCalc.h:54

Member Data Documentation

CLHEP::RandGaussQ* evwgh::GenieWeightCalc::fGaussRandom
private

Definition at line 35 of file GenieWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::GenieWeightCalc::fGenieModuleLabel
private

Definition at line 36 of file GenieWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::vector<rwgt::NuReweight *> evwgh::GenieWeightCalc::reweightVector
private

Definition at line 33 of file GenieWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

The documentation for this class was generated from the following file: