MakeG4WeightTable_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MakeG4WeightTable
3 // Plugin Type: producer (art v2_13_00)
4 // File: MakeG4WeightTable_module.cc
5 //
6 // Generated at Mon Nov 25 17:59:23 2019 by Cathal Sweeney using cetskelgen
7 // from cetlib version v3_06_01.
8 ////////////////////////////////////////////////////////////////////////
9 
10 
11 // Framework includes
21 #include "fhiclcpp/ParameterSet.h"
23 
24 #include <memory>
25 
26 // NOvASoft includes
27 #include "MCCheater/BackTracker.h"
29 #include "Geometry/Geometry.h"
30 #include "TGeoMaterial.h"
32 #include "Utilities/AssociationUtil.h"
33 #include "G4Rwgt/G4WeightTable.h"
34 #include "helperFunctions.h"
35 
36 // geant4reweight includes
37 
38 //#include "ReweightBase/G4PiPlusReweighter.hh"
39 #include "ReweightBase/G4MultiReweighter.hh"
40 #include "ReweightBase/G4ReweightTraj.hh"
41 #include "ReweightBase/G4ReweightStep.hh"
42 #include "PropBase/G4ReweightParameterMaker.hh"
43 
44 
45 // ROOT includes
46 #include "TTree.h"
47 
48 //C++ includes
49 #include <algorithm>
50 #include <cmath>
51 
52 
53 //--------------------------------------------------------------------
54 
55 namespace g4rwgt
56 {
57 
59  public:
60  explicit MakeG4WeightTable(const fhicl::ParameterSet& pset);
61  virtual ~MakeG4WeightTable() = default;
62 
63 
64  // Required functions.
65  void produce(art::Event & evt) override;
66 
67  // Selected optional functions.
68  void beginJob() override;
69  void endJob() override;
70  void reconfigure(const fhicl::ParameterSet& pset);
71 
72 
73  private:
74 
75  // Declare member data here.
77 
78  //FCL pars
81 
82  int fSeed;
83 
84  bool fConsiderNDRock; ///< Include rock interaction in the ND?
85  bool fConsiderFDRock; ///< Include rock interaction in the FD?
86 
89  //TFile fPiplusFitResults;
90 
93  //TFile fPiminusFitResults;
94 
97 
98  int fNThrows;
99  bool fDebug;
100 
102 
103  std::vector< fhicl::ParameterSet > fPionFitParSets;
104  std::vector< fhicl::ParameterSet > fProtonFitParSets;
105 
106  // These MultiReweighters will be used to get central value weights
107  // Note: in a world where central knob values are all 1.0, central
108  // value weights will all be equal to 1.0, and these nom
109  // reweighters are useless. Leaving them in because we may
110  // in the future we may use central values != 1.0
111  //
112  // N.B. I think there's some confusing naming:
113  // here "nom" means "set knobs to central values"
114  // not "set knobs to 1.0"
115  G4MultiReweighter* piplusMRW_nom;
116  G4MultiReweighter* piminusMRW_nom;
117  G4MultiReweighter* protonMRW_nom;
118 
119  // A vector of MultiReweighters where each MultiReweighter is set to
120  // a different universe. Adds a bit of memory bulk, but cuts down
121  // run time massively (MultiReweighters take *ages* to reconfigure
122  // themselves for a different universe, so we do it just once, all up front)
123  std::vector<G4MultiReweighter> piplusMRW_vec;
124  std::vector<G4MultiReweighter> piminusMRW_vec;
125  std::vector<G4MultiReweighter> protonMRW_vec;
126 
127  //Values to use when setting knob values by hand
128  std::vector<std::vector<double>> knob_values;
129  std::vector<double> throwVals;
130  };
131 
132  //--------------------------------------------------------------------
133 
135  : fGeneratorToken(consumes<std::vector<simb::MCTruth>>(pset.get<std::string>("GeneratorLabel"))),
136  fGeneratorLabel(pset.get<std::string>("GeneratorLabel")),
137  //fHitLabel(pset.get<std::string>("HitLabel")),
138  fSeed( pset.get<int>( "Seed" ) ),
139  fConsiderNDRock(pset.get<bool>("ConsiderNDRock")),
140  fConsiderFDRock(pset.get<bool>("ConsiderFDRock")),
141  fPiplusFracsFile( (pset.get< std::string >( "PiplusFracs" )).c_str(), "OPEN" ),
142  fPiplusXSecFile( (pset.get< std::string >( "PiplusXSec" )).c_str(), "OPEN" ),
143  //fPiplusFitResults( (pset.get< std::string >( "PiplusFit" )).c_str(), "OPEN" ),
144  fPiminusFracsFile( (pset.get< std::string >( "PiminusFracs" )).c_str(), "OPEN" ),
145  fPiminusXSecFile( (pset.get< std::string >( "PiminusXSec" )).c_str(), "OPEN" ),
146  //fPiminusFitResults( (pset.get< std::string >( "PiminusFit" )).c_str(), "OPEN" ),
147  fProtonFracsFile( (pset.get< std::string >( "ProtonFracs" )).c_str(), "OPEN" ),
148  fProtonXSecFile( (pset.get< std::string >( "ProtonXSec" )).c_str(), "OPEN" ),
149  fNThrows(pset.get<int>("NThrows")),
150  fDebug(pset.get<bool>("Debug")),
151  fReweightProtons(pset.get<bool>("ReweightProtons")),
152  fPionFitParSets( pset.get< std::vector<fhicl::ParameterSet> >( "pionParameterSet" )),
153  fProtonFitParSets( pset.get< std::vector<fhicl::ParameterSet> >( "protonParameterSet" ))
154  {
155 
156  //When not using Fit file order of knobs for pions is same as in MakeG4WeightTable.fcl :
157  // reac,
158  // cex,
159  // dcex,
160  // inel,
161  // abs,
162  // prod
163  //
164  //For protons there is only a single knob - reac
165  throwVals = {0.8, 0.9, 0.95, 1.05, 1.1, 1.2};
166  int nThrow_vals = throwVals.size();
167 
168  int nKnobs = fPionFitParSets.size();
169  int myNThrows = nKnobs * nThrow_vals; // +/- 5, 10, 20 % for each knob
170 
171  //When running in the hypercross configuration
172  //the user-specified fNThrows is redundant. I'm leaving it
173  //in for the now, especially as things may change in the future
174  if(myNThrows != fNThrows){
175  std::cout << "Number of throws specified in fcl not equal to (nThrow_vals x nKnobs) \n";
176  abort();
177  }
178 
179  std::vector<double> template_knobVals(nKnobs,1.0);
180  for(int i=0; i< nKnobs; i++){
181  for(int j=0; j<nThrow_vals; j++){
182  template_knobVals[i] = throwVals[j];
183  knob_values.push_back(template_knobVals);
184  template_knobVals[i] = 1.0; //Set it back to default
185  }
186  }
187 
188 
189  // dummy trajectory which is only used to set MultiReweighters below
190  G4ReweightTraj dummy_traj(0, 0, 0, 0, std::make_pair(0,0));
191 
193  = new G4MultiReweighter( 211,
197  fNThrows,
198  fSeed);
199  // we don't actually care about the weight we get here, this is just a way of
200  // setting this MultiReweighter to nominal (central value) universe
201  piplusMRW_nom->GetWeightFromNominal(dummy_traj);
202 
204  = new G4MultiReweighter( -211,
208  fNThrows,
209  fSeed);
210  piminusMRW_nom->GetWeightFromNominal(dummy_traj);
211 
212  if(fReweightProtons){
214  = new G4MultiReweighter( 2212,
218  fNThrows,
219  fSeed);
220  protonMRW_nom->GetWeightFromNominal(dummy_traj);
221  }
222 
223  for(int i=0; i<fNThrows; i++){
224 
225  G4MultiReweighter piplusMRW
226  = G4MultiReweighter( 211,
230  fNThrows,
231  fSeed);
232 
233  piplusMRW.SetAllParameterValues(knob_values[i]);
234  piplusMRW_vec.push_back(piplusMRW);
235 
236  G4MultiReweighter piminusMRW
237  = G4MultiReweighter( -211,
241  fNThrows,
242  fSeed);
243 
244  piminusMRW.SetAllParameterValues(knob_values[i]);
245  piminusMRW_vec.push_back(piminusMRW);
246 
247  // We only have one knob for protons, so need far fewer universes
248  if(fReweightProtons && (i < nThrow_vals) ){
249  G4MultiReweighter protonMRW
250  = G4MultiReweighter( 2212,
254  fNThrows,
255  fSeed);
256  std::vector proton_knobVals{throwVals[i]};
257  protonMRW.SetAllParameterValues(proton_knobVals);
258  protonMRW_vec.push_back(protonMRW);
259  }
260 
261  }
262 
263  // Call appropriate produces<>() functions here.
264  produces<std::vector< G4WeightTable >>();
265 
266  produces<art::Assns< G4WeightTable, simb::MCTruth>>();
267  }
268 
269  //--------------------------------------------------------------------
270 
272  {
273 
274  fG4ModuleLabel = pset.get< std::string >("G4ModuleLabel", "geantgen");
275 
276  return;
277  }
278 
279  //--------------------------------------------------------------------
280 
282  {
283  //Utilities
286 
287 
288  std::unique_ptr< std::vector<G4WeightTable> > g4Tablecol(new std::vector<G4WeightTable>);
289 
290  std::unique_ptr<art::Assns< G4WeightTable, simb::MCTruth>> assns(new art::Assns< G4WeightTable, simb::MCTruth>);
291 
292  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
293 
295  std::vector< art::Ptr< simb::MCTruth > > mcTruths;
296  if( evt.getByLabel("generator", mcTruthHandle) ){
297  art::fill_ptr_vector(mcTruths, mcTruthHandle);
298  }
299 
300  // Loop over nu interactions
301  for(unsigned int truthIdx = 0; truthIdx < mcTruths.size(); ++truthIdx){
302 
303  art::Ptr<simb::MCTruth> mctruth = mcTruths.at(truthIdx);
304 
305  std::vector<const sim::Particle*> parts =
306  bt->MCTruthToParticles(mctruth);
307 
308  // If we want to skip "rock" events in this detector
309  if((geom->DetId() == novadaq::cnv::kNEARDET && !fConsiderNDRock) ||
310  (geom->DetId() == novadaq::cnv::kFARDET && !fConsiderFDRock)){
311 
312  // If interaction vertex outside det, skip
313  if( !RockFilter_C(mctruth, geom) ) continue;
314 
315  }
316 
317  std::vector<double> piplus_weights_vec(fNThrows, 1.0);
318  std::vector<double> piminus_weights_vec(fNThrows, 1.0);
319  std::vector<double> proton_weights_vec(fNThrows, 1.0);
320 
321  double piplus_cv = 1.0;
322  double piminus_cv = 1.0;
323  double proton_cv = 1.0;
324 
325  for(uint i=0; i < parts.size(); i++){
326 
327  const sim::Particle *part = parts.at(i);
328 
329  //if( !isPrimary(part) ) continue;
330 
331 
332  int pdg = part->PdgCode();
333  std::string process = part->Process();
334 
335  if( abs(pdg) != 211 && (pdg != 2212 || !fReweightProtons) ) continue;
336 
337 
338  G4ReweightTraj theTraj = makeTraj(part, evt, bt);
339 
340  if(theTraj.GetNSteps() < 1) continue;
341 
342  if(pdg == 211){
343 
344  // get CV weight from MultiReweighters with nominal universe settings
345  double temp_piplus_cv = piplusMRW_nom->GetWeightFromSetParameters(theTraj);
346  piplus_cv *= temp_piplus_cv;
347 
348  for(int j=0; j < fNThrows; ++j){
349  double temp_piplus_weight = (piplusMRW_vec.at(j)).GetWeightFromSetParameters(theTraj);
350  temp_piplus_weight /= temp_piplus_cv; // Divide by central value
351  piplus_weights_vec.at(j) *= temp_piplus_weight ; // multiply weight with weights from other piplus in event
352  }
353  }
354 
355  if(pdg == -211){
356 
357  double temp_piminus_cv = piminusMRW_nom->GetWeightFromSetParameters(theTraj);
358  piminus_cv *= temp_piminus_cv;
359 
360  for(int j=0; j < fNThrows; ++j){
361  double temp_piminus_weight =(piminusMRW_vec.at(j)).GetWeightFromSetParameters(theTraj);
362  temp_piminus_weight /= temp_piminus_cv; // Divide by central value
363  piminus_weights_vec.at(j) *= temp_piminus_weight ; // multiply weight with weights from other piminus in event
364  }
365  }
366 
367  if(pdg == 2212){
368 
369  double temp_proton_cv = (protonMRW_nom)->GetWeightFromSetParameters(theTraj);
370  proton_cv *= temp_proton_cv;
371 
372  for(int j=0; j < (int)throwVals.size(); ++j){
373  double temp_proton_weight =(protonMRW_vec.at(j)).GetWeightFromSetParameters(theTraj);
374  temp_proton_weight /= temp_proton_cv; // Divide by central value
375  proton_weights_vec.at(j) *= temp_proton_weight ; // multiply weight with weights from other proton in event
376  }
377  }
378 
379 
380 
381 
382  } // End loop over particles
383 
384 
385  G4WeightTable g4Table;
386  g4Table.SetNUniv((uint)fNThrows);
387  g4Table.SetCV(211, piplus_cv);
388  g4Table.SetCV(-211, piminus_cv);
389  g4Table.SetCV(2212, proton_cv);
390  for(int i = 0; i < fNThrows; ++i){
391 
392  g4Table.SetWeight(i, 211, piplus_weights_vec.at(i));
393  g4Table.SetWeight(i, -211, piminus_weights_vec.at(i));
394  g4Table.SetWeight(i, 2212, proton_weights_vec.at(i));
395 
396  }
397  g4Tablecol->push_back(g4Table);
398 
399 
400  util::CreateAssn(*this, evt, *g4Tablecol,
401  art::Ptr<simb::MCTruth>(mcTruthHandle, truthIdx),
402  *assns);
403 
404 
405 
406  } // End loop over mcTruths (neutrino interactions)
407  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
408 
409  evt.put(std::move(g4Tablecol));
410  evt.put(std::move(assns));
411  }
412 
413 
414 
415 
417  {
418  // Implementation of optional member function here.
419  }
420 
422  {
423  // Implementation of optional member function here.
424  }
425 
426 
427 
429 
430 } // namespace
back track the reconstruction to the simulation
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.
int PdgCode() const
Definition: MCParticle.h:211
G4ReweightTraj makeTraj(const sim::Particle *part, art::Event &evt, art::ServiceHandle< cheat::BackTracker > bt)
std::vector< G4MultiReweighter > protonMRW_vec
void abs(TH1 *hist)
bool RockFilter_C(art::Ptr< simb::MCTruth > mctruth, art::ServiceHandle< geo::Geometry > geom)
Store +/- 1 sigma shifts for all geant4reweight channels.
Definition: FillTruth.h:19
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
void SetCV(int pdg, float w)
Definition: G4WeightTable.h:75
bool fConsiderFDRock
Include rock interaction in the FD?
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
Definition: MCParticle.h:214
void reconfigure(const fhicl::ParameterSet &pset)
std::vector< G4MultiReweighter > piminusMRW_vec
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
std::vector< fhicl::ParameterSet > fProtonFitParSets
TString part[npart]
Definition: Style.C:32
Far Detector at Ash River, MN.
std::vector< std::vector< double > > knob_values
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
Near Detector in the NuMI cavern.
virtual ~MakeG4WeightTable()=default
const double j
Definition: BetheBloch.cxx:29
This class describes a particle created in the detector Monte Carlo simulation.
OStream cout
Definition: OStream.cxx:6
std::vector< G4MultiReweighter > piplusMRW_vec
const art::ProductToken< std::vector< simb::MCTruth > > fGeneratorToken
void produce(art::Event &evt) override
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
bool fConsiderNDRock
Include rock interaction in the ND?
void SetNUniv(uint n)
Definition: G4WeightTable.h:95
MakeG4WeightTable(const fhicl::ParameterSet &pset)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
void SetWeight(uint i, int pdg, float w)
Definition: G4WeightTable.h:45
std::vector< fhicl::ParameterSet > fPionFitParSets
ProductToken< T > consumes(InputTag const &)
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
enum BeamMode string
unsigned int uint