LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
GenieWeightCalc.cxx
Go to the documentation of this file.
1 // GenieWeightCalc.cxx
2 //
3 // Handles event weights for GENIE systematics studies
4 //
5 // Updated by Marco Del Tutto on Feb 18 2017
6 // Ported from uboonecode to larsim on Feb 14 2017
7 // by Marco Del Tutto <marco.deltutto@physics.ox.ac.uk>
8 
9 #include "larsim/EventWeight/Base/WeightCalcCreator.h"
10 #include "larsim/EventWeight/Base/WeightCalc.h"
11 
12 #include "art/Framework/Services/Registry/ServiceHandle.h"
13 #include "art/Framework/Services/Optional/RandomNumberGenerator.h"
14 
15 #include "CLHEP/Random/RandGaussQ.h"
16 
17 #include "nutools/NuReweight/art/NuReweight.h" //GENIEReweight.h"
18 
19 #include "nusimdata/SimulationBase/MCFlux.h"
20 #include "nusimdata/SimulationBase/MCTruth.h"
21 #include "nusimdata/SimulationBase/GTruth.h"
22 
23 namespace evwgh {
24  class GenieWeightCalc : public WeightCalc
25  {
26  public:
27  GenieWeightCalc();
28  void Configure(fhicl::ParameterSet const& pset);
29  std::vector<std::vector<double> > GetWeight(art::Event & e);
30 
31  private:
32  // The reweighting utility class:
33  std::vector<rwgt::NuReweight *> reweightVector;
34 
35  CLHEP::RandGaussQ *fGaussRandom;
36  std::string fGenieModuleLabel;
37 
38  // What follows is the list of sereighting parameters present in LArSoft.
39  // Parameters with a (*) contains more that one reweighing parameter at the same time.
40  // They can only be modified changing the relative method in $NUTOOLS_DIR/source/NuReweight/GENIEReweight.cxx
41 
42  enum EReweight {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}; // ?
84 
85  DECLARE_WEIGHTCALC(GenieWeightCalc)
86  };
87 
88  GenieWeightCalc::GenieWeightCalc()
89  {}
90 
91  void GenieWeightCalc::Configure(fhicl::ParameterSet const& p)
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  }
353 
354  std::vector<std::vector<double> > GenieWeightCalc::GetWeight(art::Event & e)
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  }
388  REGISTER_WEIGHTCALC(GenieWeightCalc)
389 }
DECLARE_WEIGHTCALC
#define DECLARE_WEIGHTCALC(wghcalc)
Definition: WeightCalcCreator.h:25
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
evwgh::GenieWeightCalc::GetWeight
std::vector< std::vector< double > > GetWeight(art::Event &e)
Definition: GenieWeightCalc.cxx:354
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
MCFlux.h
object containing MC flux information
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
REGISTER_WEIGHTCALC
#define REGISTER_WEIGHTCALC(wghcalc)
Definition: WeightCalcCreator.h:29
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
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
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
e
Float_t e
Definition: plot.C:34
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
NuReweight.h
Wrapper for reweightings neutrino interactions within the ART framework.