MCReweight_service.cc
Go to the documentation of this file.
1 /*
2  * MCReweight_service.cc:
3  * Reweight service.
4  *
5  * Created on: Apr 20, 2016
6  * Author: J. Wolcott <jwolcott@fnal.gov>
7  */
8 
9 #include "TDatabasePDG.h"
10 #include "TParticlePDG.h"
11 
12 #include "cetlib_except/exception.h"
13 
16 
17 #include "NOvARwgt/interfaces/NuToolsInterface.h"
18 #include "NOvARwgt/rwgt/genie/GenieSystKnob.h"
19 #include "NOvARwgt/rwgt/genie/QE/RPASysts.h"
20 #include "NOvARwgt/rwgt/genie/QE/RPAWeights.h"
21 #include "NOvARwgt/rwgt/genie/QE/MAQESysts.h"
22 #include "NOvARwgt/rwgt/genie/DIS/Nonres1piWeights.h"
23 #include "NOvARwgt/rwgt/genie/DIS/HighWDISWeight.h"
24 #include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECFixups.h"
25 #include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECTuneSA.h"
26 #include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECTune2017.h"
27 #include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECSysts2017.h"
28 #include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECSysts2018.h"
29 #include "NOvARwgt/rwgt/tunes/TunesSA.h"
30 #include "NOvARwgt/rwgt/tunes/Tunes2017.h"
31 #include "NOvARwgt/rwgt/tunes/Tunes2018.h"
32 
34 
35 #include "Utilities/func/MathUtil.h"
36 
37 #include "MCReweight/MCReweight.h"
38 
39 namespace rwgt
40 {
41 
43  {
44  // clear the stored values before every event
45  registry.sPreProcessEvent.watch ( this, &MCReweight::Reset );
46  }
47 
48  //----------------------------------------------------------------------------
51  {
52  this->fWgts.clear();
53  this->fHaveOldGENIEMEC = false;
54 
55  // TODO: fix up the parameter set checking below for current files. For
56  // 2017 analysis the fHaveOldGENIEMEC should be false, so just return
57  // the function at this point.
58 
59  return;
60 
61  // For now, this is an ugly hack which can only tell the difference between
62  // an "older" and a "newer" version of GENIE,
63  // based on https://cdcvs.fnal.gov/redmine/projects/nutools/wiki/GENIE_Configuration_Files.
64  // In the future, further improvements to the generator record
65  // within the datastream should be used to fill this field properly.
67  evt.getProcessParameterSet("Genie", ps);
68 
69  MF_LOG_VERBATIM("MCReweight")
70  << ps.to_string();
71 
72  std::string GENIEGeneratorList = ps.get<fhicl::ParameterSet>("physics")
73  .get<fhicl::ParameterSet>("producers")
74  .get<fhicl::ParameterSet>("generator")
75  .get<std::string>("EventGeneratorList");
76 
77  if (GENIEGeneratorList == "DefaultPlusCCMEC")
78  this->fHaveOldGENIEMEC = false;
79  else if (GENIEGeneratorList == "Default+MEC" || GENIEGeneratorList == "Default+CCMEC ")
80  this->fHaveOldGENIEMEC = true;
81  // 'Default' doesn't identify anything. Assume it's old, I guess?
82  else if (GENIEGeneratorList == "Default")
83  this->fHaveOldGENIEMEC = true;
84  else
85  throw cet::exception("MCReweight")
86  << "Could not recognize GENIE event generator list name: "
87  << GENIEGeneratorList;
88 
89  }
90 
91  //////////
92 
93  /// Do the dirty work of matching the WeightType \a rwType to a function, then call it with \a neutrino as a parameter
95  const simb::GTruth* gtruth,
96  WeightType rwType) const
97  {
98  const novarwgt::EventRecord & rwRec = this->GetNOvARwgtEvRec(mctruth, gtruth);
99 
100  // individual weighter calls first
101  const novarwgt::IWeightGenerator * wgtr = nullptr;
102  if(rwType == kRPACCQE_SA)
104  else if(rwType == kRPACCQE_2017)
106  else if(rwType == kRPACCRES_2017)
108  else if(rwType == kMaCCQE_2018)
110  else if(rwType == kNonRes1Pi)
112  else if(rwType == kHighWDIS_2018)
114  else if(rwType == kDytmanMEC_FixItlState)
116  else if(rwType == kDytmanMEC_FixXsecEdep)
118  else if(rwType == kTufts2p2h_SA)
120  else if(rwType == kEmpiricalMEC_2017)
122  else if(rwType == kEmpiricalMEC_2018)
124 
125  if (wgtr)
126  return wgtr->GetWeight(rwRec, {});
127 
128  // the tunes have a slightly different interface
129  if(rwType == kTuftsCC_SA)
130  return novarwgt::kCVTuneSA.EventWeight(rwRec);
131  else if(rwType == kXSecCVWgt_2017)
132  return novarwgt::kCVTune2017.EventWeight(rwRec);
133  else if(rwType == kXSecCVWgt_2018)
134  return novarwgt::kCVTune2018.EventWeight(rwRec);
135 
136 
137  throw cet::exception( Form("MCReweight::DispatchToWeighter(): Unknown weight type enum value: %d", rwType) );
138 
139  } // MCReweight::DispatchToWeighter()
140 
141  //////////
142 
143  const novarwgt::EventRecord & MCReweight::GetNOvARwgtEvRec(const simb::MCTruth* mctruth,
144  const simb::GTruth* gtruth) const
145  {
146  const simb::MCNeutrino *neutrino = &(mctruth->GetNeutrino());
147  if (neutrino != this->fNOvARwgtEvtRecCache.first)
148  {
149  this->fNOvARwgtEvtRecCache.second = novarwgt::ConvertNuToolsEvent(mctruth, gtruth, {});
150  this->fNOvARwgtEvtRecCache.first = neutrino;
151  }
152 
153  return this->fNOvARwgtEvtRecCache.second;
154 
155  }
156 
157  //////////
158 
159  double MCReweight::GetWeight (const simb::MCTruth* mctruth,
160  WeightType rwType,
161  const simb::GTruth* gtruth) const
162  {
163  const simb::MCNeutrino* neutrino = &(mctruth->GetNeutrino());
164 
165  // check if this weight has already been calculated.
166  if (this->fWgts.find(neutrino) == this->fWgts.end() ||
167  this->fWgts.at(neutrino).find(rwType) == this->fWgts.at(neutrino).end())
168  {
169  // calculate if not.
170  this->fWgts[neutrino][rwType] = this->DispatchToWeighter(mctruth, gtruth, rwType);
171  }
172 
173  return this->fWgts.at(neutrino).at(rwType);
174  }
175 
176  ///////////////////////////////
177  //
178  // systematics functions
179  //
180  ///////////////////////////////
181 
183  const simb::MCTruth* mctruth,
184  const simb::GTruth* gtruth ) const
185  {
186  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
187  return novarwgt::kMAQEGenieReducedSyst2018->GetWeight(sigma, evt);
188  }
189 
191  const simb::MCTruth* mctruth,
192  const simb::GTruth* gtruth ) const
193  {
194  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
195  return novarwgt::kRPACCQESuppSyst2018->GetWeight(sigma, evt);
196  }
197 
198  double MCReweight::RPACCQEEnhSyst2018( const double sigma,
199  const simb::MCTruth* mctruth,
200  const simb::GTruth* gtruth ) const
201  {
202  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
203  return novarwgt::kRPACCQEEnhSyst2018->GetWeight(sigma, evt);
204  }
205 
206  double MCReweight::RPARESSyst2018( const double sigma,
207  const simb::MCTruth* mctruth,
208  const simb::GTruth* gtruth ) const
209  {
210  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
211  return novarwgt::kRPARESSyst2018->GetWeight(sigma, evt);
212  }
213 
215  const simb::MCTruth* mctruth,
216  const simb::GTruth* gtruth ) const
217  {
218  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
219  return novarwgt::kMECQ0Q3RespSystNu2018->GetWeight(sigma, evt);
220  }
221 
222 
224  const simb::MCTruth* mctruth,
225  const simb::GTruth* gtruth ) const
226  {
227  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
228  return novarwgt::kMECEnuShapeSyst2017_NuOnly->GetWeight(sigma, evt)
229  * novarwgt::kMECEnuShapeSyst2017_NubarOnly->GetWeight(sigma, evt);
230  }
231 
233  const simb::MCTruth* mctruth,
234  const simb::GTruth* gtruth ) const
235  {
236  const auto & evt = this->GetNOvARwgtEvRec(mctruth, gtruth);
237  return novarwgt::kMECInitStateNPFracSyst2017_NuOnly->GetWeight(sigma, evt)
239  }
240 
241  double MCReweight::GetGENIEWeight( const double sigma,
242  const simb::MCTruth* mctruth,
243  const simb::GTruth* gtruth,
244  rwgt::EReweightLabel rw_label ) const
245  {
246  if( sigma == 0 ) return 1;
247 
248  auto knob = novarwgt::GetGenieSystKnob(novarwgt::kNugenKnobTranslationTable.at(rw_label));
249  return knob->CalcWeight(sigma, novarwgt::ConvertNuToolsEvent(mctruth, gtruth, {}));
250  }
251 
252 
254 
255 } /* namespace rwgt */
const MAQEGenieReducedSyst2018 * kMAQEGenieReducedSyst2018
Definition: MAQESysts.cxx:19
double MECq0q3ShapeSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
const MECEnuShapeSyst2017 * kMECEnuShapeSyst2017_NubarOnly
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
double RPACCQEEnhSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
const novarwgt::Tune kCVTuneSA({{"Nonres1pi", novarwgt::GetWeighter< Nonres1PiWgt >(true, false)},{"MEC", novarwgt::GetWeighter< Tufts2p2hWgtSA >()},},{})
#define DEFINE_ART_SERVICE(svc)
Definition: ServiceMacros.h:88
double MECInitStateNPFracSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
double GetGENIEWeight(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth, rwgt::EReweightLabel rw_label) const
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const std::unordered_map< rwgt::ReweightLabel_t, novarwgt::ReweightKnob > kNugenKnobTranslationTable
void Reset(const art::Event &evt, art::ScheduleContext)
Clears stored weights.
bool getProcessParameterSet(std::string const &process, fhicl::ParameterSet &) const
const EmpiricalMECWgt2017 * kEmpiricalMECWgt2017
novarwgt::EventRecord ConvertNuToolsEvent(const simb::MCTruth *mctruth, const simb::GTruth *gtruth, const novarwgt::ReweightList &rwList)
const RPACCQEshapeSyst2018 * kRPACCQESuppSyst2018
Definition: RPASysts.cxx:21
WeightType
Definition: MCReweight.h:33
const Nonres1PiWgt * kNonres1PiWgt_Typo
const RPARESSyst2018 * kRPARESSyst2018
Definition: RPASysts.cxx:31
const RPAWeightQ2_2017 * kRPAWeightRES2017
Definition: RPAWeights.cxx:54
Use NuReweight to compute +/-1,2sigma shifts for all systematics.
Definition: FillTruth.h:18
const HighWDISWgt_2018 * kHighWDISWgt_2018
std::unordered_map< const simb::MCNeutrino *, std::unordered_map< WeightType, double, std::hash< int > >, NeutrinoHasher > fWgts
cache of already-calculated weights for each neutrino
Definition: MCReweight.h:122
const DytmanMECFixItlStateWgt * kDytmanMECFixItlStateWgt
const MECInitStateNPFracSyst2017 * kMECInitStateNPFracSyst2017_NuOnly
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
const novarwgt::Tune kCVTune2018({{"MA_QE", novarwgt::GetWeighter< MAQEWeight_2018 >()}, {"RPA_QE", novarwgt::GetWeighter< RPAWeightCCQE_2017 >("CV", novarwgt::kScQuasiElastic, true, false, true)},{"RPA_RES", novarwgt::GetWeighter< RPAWeightQ2_2017 >(novarwgt::kRxnCC, novarwgt::kScResonant, true)},{"Nonres1pi", novarwgt::GetWeighter< Nonres1PiWgt >(false, true)},{"HighW", novarwgt::GetWeighter< HighWDISWgt_2018 >()},{"MEC", novarwgt::GetWeighter< EmpiricalMECWgt2018 >()},}, k2018GENIEKnobs|k2018MECKnobs|k2018OtherCustomKnobs)
CV tune used for 2018 analysis.
MCReweight(const fhicl::ParameterSet &pset, art::ActivityRegistry &reg)
const Tufts2p2hWgtSA * kTufts2p2hWgtSA
double GetWeight(const simb::MCTruth *mctruth, WeightType rwType, const simb::GTruth *gtruth=nullptr) const
Retrieves a weight by name. Calculates on demand, but only calculates once per neutrino.
double MAQEGenieReducedSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
Systematics functions.
const RPAWeightCCQE_2017 * kRPAWeightCCQE2017
Definition: RPAWeights.cxx:22
GlobalSignal< detail::SignalResponseType::FIFO, void(Event const &, ScheduleContext)> sPreProcessEvent
const MECEnuShapeSyst2017 * kMECEnuShapeSyst2017_NuOnly
double sigma(TH1F *hist, double percentile)
std::pair< const simb::MCNeutrino *, novarwgt::EventRecord > fNOvARwgtEvtRecCache
Definition: MCReweight.h:124
const novarwgt::Tune kCVTune2017({{"RPA_QE", novarwgt::GetWeighter< RPAWeightCCQE_2017 >("CV", novarwgt::kScNull, true, true, false)},{"Nonres1pi", novarwgt::GetWeighter< Nonres1PiWgt >(true, false)},{"MEC", novarwgt::GetWeighter< EmpiricalMECWgt2017 >()},}, k2017GENIEKnobs|k2017CustomKnobs)
CV MC tune used in 2017 analyses.
#define MF_LOG_VERBATIM(category)
const RPAWeightCCQESA * kRPAWeightCCQESA
Definition: RPAWeights.cxx:18
const MECQ0Q3RespSyst2018 * kMECQ0Q3RespSystNu2018
double RPACCQESuppSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
const EmpiricalMECWgt2018 * kEmpiricalMECWgt2018
double MECEnuShapeSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
const DytmanMECFixXsecEdepWgt * kDytmanMECFixXsecEdepWgt
const MAQEWeight_2018 * kMAQEWeight_2018
Definition: MAQEWgts.cxx:17
Event generator information.
Definition: MCTruth.h:32
const MECInitStateNPFracSyst2017 * kMECInitStateNPFracSyst2017_NubarOnly
const RPACCQEshapeSyst2018 * kRPACCQEEnhSyst2018
Definition: RPASysts.cxx:20
Event generator information.
Definition: MCNeutrino.h:18
double RPARESSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
double DispatchToWeighter(const simb::MCTruth *mctruth, const simb::GTruth *gtruth, WeightType rwType) const
Finds the right reweight function for the given weight type.
const novarwgt::EventRecord & GetNOvARwgtEvRec(const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
std::string to_string() const
Definition: ParameterSet.h:137
enum BeamMode string