Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Private Attributes | List of all members
rwgt::MakeGENIEReweightTable Class Reference

Use NuReweight to compute +/-1,2sigma shifts for all systematics. More...

Inheritance diagram for rwgt::MakeGENIEReweightTable:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 MakeGENIEReweightTable (const fhicl::ParameterSet &pset)
 
virtual ~MakeGENIEReweightTable ()
 
void produce (art::Event &evt) override
 
void beginJob () override
 
void endJob () override
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Attributes

const art::ProductToken< std::vector< simb::MCTruth > > fGeneratorToken
 
std::string fGeneratorLabel
 
bool fConsiderNDRock
 Include rock interaction in the ND? More...
 
bool fConsiderFDRock
 Include rock interaction in the FD? More...
 
bool fTimeKnobs
 Keep tabs on the reweight knobs & measure the time they take? More...
 
bool fSkipCCQEWgtsForHyperons
 
std::vector< std::stringfKnobNamesToSkip
 
std::unordered_map< std::string, std::unordered_set< rwgt::ReweightLabel_t, std::hash< int > > > fKnobsToSkip
 
std::vector< rwgt::NuReweight * > fRWsMinus2
 
std::vector< rwgt::NuReweight * > fRWsMinus1
 
std::vector< rwgt::NuReweight * > fRWsPlus1
 
std::vector< rwgt::NuReweight * > fRWsPlus2
 
std::map< rwgt::ReweightLabel_t, TH1D * > fKnobTimes
 
art::ServiceHandle< art::TFileServicefFileService
 

Detailed Description

Use NuReweight to compute +/-1,2sigma shifts for all systematics.

Definition at line 50 of file MakeGENIEReweightTable_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

rwgt::MakeGENIEReweightTable::MakeGENIEReweightTable ( const fhicl::ParameterSet pset)
explicit

Definition at line 85 of file MakeGENIEReweightTable_module.cc.

References febshutoff_auto::end, allTimeWatchdog::endl, fKnobsToSkip, fhicl::ParameterSet::get(), and rwgt::REWEIGHT_KNOB_INDICES.

86  : fGeneratorToken(consumes<std::vector<simb::MCTruth>>(pset.get<std::string>("GeneratorLabel"))),
87  fGeneratorLabel(pset.get<std::string>("GeneratorLabel")),
88  fConsiderNDRock(pset.get<bool>("ConsiderNDRock")),
89  fConsiderFDRock(pset.get<bool>("ConsiderFDRock")),
90  fTimeKnobs(pset.get<bool>("TimeKnobs")),
91  fSkipCCQEWgtsForHyperons(pset.get<bool>("SkipCCQEWgtsForHyperons"))
92 // fKnobNamesToSkip(pset.get<std::unordered_map<std::string, std::vector<std::string>>>("KnobsToSkip"))
93  {
94  produces<std::vector<rwgt::GENIEReweightTable>>();
95 
96  produces<art::Assns<rwgt::GENIEReweightTable, simb::MCTruth>>();
97 
98  mf::LogInfo("MakeGENIEReweightTable") << "Skipping the following reweight knobs:";
99  auto knobsToSkipLists = pset.get<fhicl::ParameterSet>("KnobsToSkip");
100  if (!knobsToSkipLists.has_key("default"))
101  {
102  mf::LogWarning("MakeGENIEReweightTable") << "No 'default' key in KnobsToSkip table. Assuming empty list...";
103  fKnobsToSkip["default"]; // default-initializes...
104  }
105  for (const auto & knobListKey : knobsToSkipLists.get_names())
106  {
107  mf::LogInfo("MakeGENIEReweightTable") << " for tune '" << knobListKey << "'" << std::endl;
108  for (const auto &knobName : knobsToSkipLists.get<std::vector<std::string>>(knobListKey))
109  {
111  throw cet::exception("MakeGENIEReweightTable",
112  "Unknown GENIE reweight knob name: " + knobName);
113 
114  this->fKnobsToSkip[knobListKey].insert(rwgt::REWEIGHT_KNOB_INDICES.at(knobName));
115  mf::LogInfo("MakeGENIEReweightTable")
116  << " " << knobName
117  << " (index: " << rwgt::REWEIGHT_KNOB_INDICES.at(knobName) << ")";
118  }
119  }
120  }
bool fTimeKnobs
Keep tabs on the reweight knobs & measure the time they take?
bool fConsiderNDRock
Include rock interaction in the ND?
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const std::unordered_map< std::string, ReweightLabel_t > REWEIGHT_KNOB_INDICES
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const art::ProductToken< std::vector< simb::MCTruth > > fGeneratorToken
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::unordered_map< std::string, std::unordered_set< rwgt::ReweightLabel_t, std::hash< int > > > fKnobsToSkip
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
ProductToken< T > consumes(InputTag const &)
bool fConsiderFDRock
Include rock interaction in the FD?
enum BeamMode string
rwgt::MakeGENIEReweightTable::~MakeGENIEReweightTable ( )
virtual

Definition at line 123 of file MakeGENIEReweightTable_module.cc.

References fRWsMinus1, fRWsMinus2, fRWsPlus1, and fRWsPlus2.

124  {
125  for(rwgt::NuReweight* rw: fRWsMinus2) delete rw;
126  for(rwgt::NuReweight* rw: fRWsMinus1) delete rw;
127  for(rwgt::NuReweight* rw: fRWsPlus1) delete rw;
128  for(rwgt::NuReweight* rw: fRWsPlus2) delete rw;
129  }
std::vector< rwgt::NuReweight * > fRWsPlus1
std::vector< rwgt::NuReweight * > fRWsMinus2
std::vector< rwgt::NuReweight * > fRWsPlus2
std::vector< rwgt::NuReweight * > fRWsMinus1

Member Function Documentation

void rwgt::MakeGENIEReweightTable::beginJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 132 of file MakeGENIEReweightTable_module.cc.

References ERROR, fRWsMinus1, fRWsMinus2, fRWsPlus1, fRWsPlus2, MECModelEnuComparisons::i, genie::Messenger::Instance(), rwgt::kMaxReweightIdx, stan::math::resize(), genie::Messenger::SetPriorityLevel(), and sigma().

133  {
134  mf::LogInfo("MakeGENIEReweightTable") << "Initializing GENIE reweight calculators...";
135 
136  // for GENIE 3, need to set the tune before you do anything.
137  // (this function is a no-op for GENIE 2 so it's safe anyway.)
138  // note that this sets the tune to the one that the genie_xsec UPS product
139  // is configured to load. 99.9999% of the time that's the one you want,
140  // but if you're trying to create reweights for MC that was generated
141  // in a different release than the one you're using,
142  // you may need to add some FCL handles to pass options to this function.
143  evgb::SetEventGeneratorListAndTune();
144 
145  // Don't need to see every tiny thing these modules are thinking...
147  log4cpp::Priority::WARN);
149  log4cpp::Priority::WARN);
151  log4cpp::Priority::WARN);
152 
153  // The way NuReweight constructs the EventRecord it inevitably trips this
156 
157  genie::Messenger::Instance()->SetPriorityLevel("GENIEReweight",
158  log4cpp::Priority::WARN);
159 
160  fRWsMinus2.resize(kMaxReweightIdx+1, 0);
161  fRWsMinus1.resize(kMaxReweightIdx+1, 0);
164 
165  // Start at 1. Zero is kReweightNull. Not worth dealing with shifting all
166  // the indices by one everywhere, just leave it empty.
167  for(int i = 1; i <= kMaxReweightIdx; ++i){
168  for(int sigma: {-2, -1, +1, +2}){
169  rwgt::NuReweight* rw = new rwgt::NuReweight;
170 
171  switch(sigma){
172  case -2: fRWsMinus2[i] = rw; break;
173  case -1: fRWsMinus1[i] = rw; break;
174  case +1: fRWsPlus1 [i] = rw; break;
175  case +2: fRWsPlus2 [i] = rw; break;
176  }
177  double shift = sigma;
178 
179  //Catch the cases where shape parameters are called instead of rate parameters.
180  switch(i) {
181  case rwgt::fReweightNormCCQE:
182  case rwgt::fReweightNormCCQEenu:
183  case rwgt::fReweightMaCCQEshape:
184  rw->MaQEshape();
185  break;
186 
187  case rwgt::fReweightNormCCRES:
188  case rwgt::fReweightMaCCRESshape:
189  case rwgt::fReweightMvCCRESshape:
190  rw->CCRESshape();
191  break;
192 
193  case rwgt::fReweightNormNCRES:
194  case rwgt::fReweightMaNCRESshape:
195  case rwgt::fReweightMvNCRESshape:
196  rw->NCRESshape();
197  break;
198 
199  case rwgt::fReweightAhtBYshape:
200  case rwgt::fReweightBhtBYshape:
201  case rwgt::fReweightCV1uBYshape:
202  case rwgt::fReweightCV2uBYshape:
203  rw->DIS_BYshape();
204  break;
205 
206  case rwgt::fReweightCCQEMomDistroFGtoSF:
207  case rwgt::fReweightTheta_Delta2Npi:
208  // These are continuous knobs.
209  // Encode -2sigma = .25, -1sigma = .5, +1sigma = .75, +2sigma = 1
210  shift = (shift < 0) ? (shift+3.0)/4.0 : (shift + 2)/4.0;
211  break;
212 
213  case rwgt::fReweightRnubarnuCC:
214 // case rwgt::fReweightDISNuclMod: // disabled in GENIE v3. may return later?
215  // These are binary switches. Make the negative sigmas be "off" and
216  // the positive be "on".
217  shift = (shift < 0) ? 0 : 1;
218  break;
219  }
220 
221  rw->AddReweightValue(rwgt::ReweightLabel_t(i), shift);
222  rw->Configure();
223  }
224  }
225 
226  mf::LogInfo("MakeGENIEReweightTable") << "Done initializing GENIE reweight calculators...";
227 
228  }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< rwgt::NuReweight * > fRWsPlus1
std::vector< rwgt::NuReweight * > fRWsMinus2
std::vector< rwgt::NuReweight * > fRWsPlus2
std::vector< rwgt::NuReweight * > fRWsMinus1
void resize(T &x, std::vector< int > dims)
Definition: resize.hpp:41
static Messenger * Instance(void)
Definition: Messenger.cxx:71
#define ERROR
double sigma(TH1F *hist, double percentile)
void SetPriorityLevel(const char *stream, log4cpp::Priority::Value p)
Definition: Messenger.cxx:110
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
void rwgt::MakeGENIEReweightTable::endJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 230 of file MakeGENIEReweightTable_module.cc.

References bin, om::cout, allTimeWatchdog::endl, fKnobTimes, getGoodRuns4SAM::n, cet::pow(), sum, submit_syst::x, and submit_syst::y.

231  {
232  std::cout << "Statistics regarding time spent in each reweight knob:" << std::endl;
233  std::cout << "knob # | total time (s) | mean time (s) | median time (s)" << std::endl;
234  std::cout << "---------------------------------------------------------" << std::endl;
235 
236  for (auto & knobTimePair : this->fKnobTimes)
237  {
238  std::cout << std::setw(7) << knobTimePair.first << "|";
239  TH1D * knobHist = knobTimePair.second;
240 
241  double sum = 0;
242  for (int bin = 1; bin < knobHist->GetNbinsX()+1; bin++)
243  sum += pow(10, knobHist->GetBinCenter(bin)) * knobHist->GetBinContent(bin);
244  std::cout << std::setw(16) << sum << "|";
245 
246  std::cout << std::setw(15);
247  if (knobHist->GetEntries() > 0)
248  std::cout << pow(10, knobHist->GetMean());
249  else
250  std::cout << "-";
251  std::cout << "|";
252 
253  std::cout << std::setw(16);
254  int n = knobHist->GetXaxis()->GetNbins();
255  std::vector<double> x(n);
256  knobHist->GetXaxis()->GetCenter( &x[0] );
257  const double * y = knobHist->GetArray();
258  std::cout << pow(10, TMath::Median(n, &x[0], &y[1]));
259 
260  std::cout << std::endl;
261  }
262  }
constexpr T pow(T x)
Definition: pow.h:75
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
std::map< rwgt::ReweightLabel_t, TH1D * > fKnobTimes
Double_t sum
Definition: plot.C:31
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void rwgt::MakeGENIEReweightTable::produce ( art::Event evt)
overridevirtual

Implements art::EDProducer.

Definition at line 265 of file MakeGENIEReweightTable_module.cc.

References ana::assert(), genie::utils::res::AsString(), util::CreateAssn(), DEFINE_ART_MODULE(), geo::GeometryBase::DetId(), art::Handle< T >::failedToGet(), fConsiderFDRock, fConsiderNDRock, fFileService, fGeneratorLabel, fGeneratorToken, simb::GTruth::fgQ2, simb::GTruth::fIsCharm, simb::GTruth::fIsStrange, fKnobsToSkip, fKnobTimes, fRWsMinus1, fRWsMinus2, fRWsPlus1, fRWsPlus2, fSkipCCQEWgtsForHyperons, fTimeKnobs, geom(), art::DataViewImpl::getByToken(), geo::GeometryBase::isInsideDetectorBigBox(), simb::kCC, novadaq::cnv::kFARDET, simb::kGENIE, rwgt::kMaxReweightIdx, novadaq::cnv::kNEARDET, simb::kQE, genie::kScQuasiElastic, std::log10(), art::TFileDirectory::make(), art::Event::put(), rwgt::GENIEReweightTable::SetWeights(), febshutoff_auto::start, string, and getGoodRuns4SAM::table.

266  {
268 
269  std::unique_ptr<std::vector<rwgt::GENIEReweightTable>> tablecol(new std::vector<rwgt::GENIEReweightTable>);
270 
271  std::unique_ptr<art::Assns<rwgt::GENIEReweightTable, simb::MCTruth>> assns(new art::Assns<rwgt::GENIEReweightTable, simb::MCTruth>);
272 
273 
275  evt.getByToken(fGeneratorToken, mctruths);
276  assert(!mctruths.failedToGet());
277 
278  art::FindManyP<simb::GTruth> fmp(mctruths, evt, fGeneratorLabel);
279  assert(fmp.isValid());
280 
281  for(unsigned int truthIdx = 0; truthIdx < mctruths->size(); ++truthIdx){
282  const std::vector<art::Ptr<simb::GTruth>> gtruths = fmp.at(truthIdx);
283  assert(gtruths.size() == 1);
284 
285  const simb::MCTruth& mctruth = (*mctruths)[truthIdx];
286 
287  // obviously we won't have a reweight table if this isn't a GENIE event
288  if (mctruth.GeneratorInfo().generator != simb::Generator_t::kGENIE)
289  continue;
290 
291  const simb::GTruth& gtruth = *gtruths[0];
292 
293  mf::LogDebug("MakeGENIEReweightTable") << "Reweighting event:\n"
294  << "Interaction type: " << (mctruth.GetNeutrino().CCNC() == simb::kCC ? "CC" : "NC") << " " << ", mode = " << mctruth.GetNeutrino().Mode();
295 
296  // If we want to skip rock events in this detector
297  if((geom->DetId() == novadaq::cnv::kNEARDET && !fConsiderNDRock) ||
298  (geom->DetId() == novadaq::cnv::kFARDET && !fConsiderFDRock)){
299  // Where's the true vertex?
300  const TVector3 vtx = mctruth.GetNeutrino().Nu().Position().Vect();
301 
302  // If the vertex is outside the detector, skip this event
303  if(!geom->isInsideDetectorBigBox(vtx.X(), vtx.Y(), vtx.Z())) continue;
304  }
305 
306  // Signal of a nu/e interaction, which can cause trouble
307  if(gtruth.fgQ2 < 0) continue;
308 
309  GENIEReweightTable table;
310 
311  std::string tuneName = "default";
312  const auto * knobsToSkip = &this->fKnobsToSkip.at("default");
313  if (!mctruth.GeneratorInfo().generatorVersion.empty()
314  && (mctruth.GeneratorInfo().generatorVersion != "unknown")
315  && std::atoi(&mctruth.GeneratorInfo().generatorVersion[0]) >= 3
316  && mctruth.GeneratorInfo().generatorConfig.find("tune") != mctruth.GeneratorInfo().generatorConfig.end())
317  {
318  tuneName = mctruth.GeneratorInfo().generatorConfig.at("tune");
319  if (this->fKnobsToSkip.find(tuneName) == this->fKnobsToSkip.end())
320  {
321  mf::LogWarning("MakeGENIEReweightTable") << "No tune-specific knob list for tune '" << tuneName
322  << "'. Using default list.";
323  tuneName = "default";
324  }
325  else
326  knobsToSkip = &this->fKnobsToSkip.at(tuneName);
327  }
328  mf::LogDebug("MakeGENIEReweightTable") << "Using knob list for tune: '" << tuneName << "'";
329  bool hyperonCCQE = ((gtruth.fIsStrange || gtruth.fIsCharm)
330  && mctruth.GetNeutrino().Mode() == simb::kQE
331  && mctruth.GetNeutrino().CCNC() == simb::kCC);
332  for(int systIdx = 1; systIdx <= kMaxReweightIdx; ++systIdx){
333  if (knobsToSkip->find(ReweightLabel_t(systIdx)) != knobsToSkip->end())
334  continue;
335 
336  // CCQE knobs sometimes barf on CCQE-strange antinu events. assume it probably does the same on CCQE-charm nu events too
337  if (hyperonCCQE && fSkipCCQEWgtsForHyperons)
338  {
339  genie::rew::GReWeight * wgtCalc = fRWsPlus1[systIdx]->WeightCalculator();
340  bool skipThisKnob = false;
341  for (const auto & wgtrName : wgtCalc->WghtCalcNames())
342  {
343  const genie::rew::GReWeightI * wgtr = wgtCalc->WghtCalc(wgtrName);
344  if (wgtr->AppliesTo(genie::kScQuasiElastic, true))
345  {
346  skipThisKnob = true;
347  break;
348  }
349  }
350  if (skipThisKnob)
351  {
352  mf::LogWarning("MakeGENIEReweightTable") << "Skipping knob "
353  << genie::rew::GSyst::AsString(static_cast<genie::rew::GSyst_t>(systIdx))
354  << " for QES event:";
355  mf::LogWarning("MakeGENIEReweightTable") << gtruth;
356  continue;
357  }
358  }
359 
360  std::clock_t start = std::clock();
361 
362  const double wm2 = fRWsMinus2[systIdx]->CalcWeight(mctruth, gtruth);
363  const double wm1 = fRWsMinus1[systIdx]->CalcWeight(mctruth, gtruth);
364  const double wp1 = fRWsPlus1 [systIdx]->CalcWeight(mctruth, gtruth);
365  const double wp2 = fRWsPlus2 [systIdx]->CalcWeight(mctruth, gtruth);
366 
367  double duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
368 
369  rwgt::EReweightLabel systLabel = rwgt::EReweightLabel(systIdx);
370  if (this->fTimeKnobs)
371  {
372  if (this->fKnobTimes.find(systLabel) == this->fKnobTimes.end())
373  this->fKnobTimes[systLabel] = this->fFileService->make<TH1D>( Form("knob_%d", systIdx), Form(";Log_{10}(Time reweight knob %d used in event) (s)", systIdx), 100, -7, 3);
374 
375  if (duration <= 0)
376  this->fKnobTimes[systLabel]->Fill(-7);
377  else
378  this->fKnobTimes[systLabel]->Fill(log10(duration));
379  }
380 
381  table.SetWeights(systIdx, wm2, wm1, wp1, wp2);
382  } // end for systIdx
383 
384  tablecol->push_back(table);
385 
386  util::CreateAssn(*this, evt, *tablecol,
387  art::Ptr<simb::MCTruth>(mctruths, truthIdx),
388  *assns);
389  } // end for truthIdx
390 
391 
392  evt.put(std::move(tablecol));
393  evt.put(std::move(assns));
394  }
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
bool fTimeKnobs
Keep tabs on the reweight knobs & measure the time they take?
bool fConsiderNDRock
Include rock interaction in the ND?
std::vector< rwgt::NuReweight * > fRWsPlus1
std::vector< rwgt::NuReweight * > fRWsMinus2
bool fIsStrange
strange production // added version 13
Definition: GTruth.h:73
std::vector< rwgt::NuReweight * > fRWsPlus2
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const art::ProductToken< std::vector< simb::MCTruth > > fGeneratorToken
::xsd::cxx::tree::duration< char, simple_type > duration
Definition: Database.h:188
Far Detector at Ash River, MN.
std::vector< rwgt::NuReweight * > fRWsMinus1
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
bool fIsCharm
did the interaction produce a charmed hadron?
Definition: GTruth.h:71
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
T * make(ARGS...args) const
std::map< rwgt::ReweightLabel_t, TH1D * > fKnobTimes
T log10(T number)
Definition: d0nt_math.hpp:120
std::unordered_map< std::string, std::unordered_set< rwgt::ReweightLabel_t, std::hash< int > > > fKnobsToSkip
double fgQ2
Definition: GTruth.h:61
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
const char * AsString(Resonance_t res)
resonance id -> string
Event generator information.
Definition: MCTruth.h:32
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
bool fConsiderFDRock
Include rock interaction in the FD?
bool failedToGet() const
Definition: Handle.h:196
art::ServiceHandle< art::TFileService > fFileService
enum BeamMode string
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Member Data Documentation

bool rwgt::MakeGENIEReweightTable::fConsiderFDRock
private

Include rock interaction in the FD?

Definition at line 65 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().

bool rwgt::MakeGENIEReweightTable::fConsiderNDRock
private

Include rock interaction in the ND?

Definition at line 64 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().

art::ServiceHandle<art::TFileService> rwgt::MakeGENIEReweightTable::fFileService
private

Definition at line 81 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().

std::string rwgt::MakeGENIEReweightTable::fGeneratorLabel
private

Definition at line 62 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().

const art::ProductToken<std::vector<simb::MCTruth> > rwgt::MakeGENIEReweightTable::fGeneratorToken
private

Definition at line 61 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().

std::vector<std::string> rwgt::MakeGENIEReweightTable::fKnobNamesToSkip
private

Definition at line 71 of file MakeGENIEReweightTable_module.cc.

std::unordered_map<std::string, std::unordered_set<rwgt::ReweightLabel_t, std::hash<int> > > rwgt::MakeGENIEReweightTable::fKnobsToSkip
private

Definition at line 73 of file MakeGENIEReweightTable_module.cc.

Referenced by MakeGENIEReweightTable(), and produce().

std::map<rwgt::ReweightLabel_t, TH1D*> rwgt::MakeGENIEReweightTable::fKnobTimes
private

Definition at line 80 of file MakeGENIEReweightTable_module.cc.

Referenced by endJob(), and produce().

std::vector<rwgt::NuReweight*> rwgt::MakeGENIEReweightTable::fRWsMinus1
private

Definition at line 76 of file MakeGENIEReweightTable_module.cc.

Referenced by beginJob(), produce(), and ~MakeGENIEReweightTable().

std::vector<rwgt::NuReweight*> rwgt::MakeGENIEReweightTable::fRWsMinus2
private

Definition at line 75 of file MakeGENIEReweightTable_module.cc.

Referenced by beginJob(), produce(), and ~MakeGENIEReweightTable().

std::vector<rwgt::NuReweight*> rwgt::MakeGENIEReweightTable::fRWsPlus1
private

Definition at line 77 of file MakeGENIEReweightTable_module.cc.

Referenced by beginJob(), produce(), and ~MakeGENIEReweightTable().

std::vector<rwgt::NuReweight*> rwgt::MakeGENIEReweightTable::fRWsPlus2
private

Definition at line 78 of file MakeGENIEReweightTable_module.cc.

Referenced by beginJob(), produce(), and ~MakeGENIEReweightTable().

bool rwgt::MakeGENIEReweightTable::fSkipCCQEWgtsForHyperons
private

Definition at line 69 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().

bool rwgt::MakeGENIEReweightTable::fTimeKnobs
private

Keep tabs on the reweight knobs & measure the time they take?

Definition at line 67 of file MakeGENIEReweightTable_module.cc.

Referenced by produce().


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