LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PPFXTargAttenWeightCalc.cxx
Go to the documentation of this file.
1 
4 
5 #include "CLHEP/Random/RandGaussQ.h"
6 
7 #include "MakeReweight.h"
8 #include "cetlib/search_path.h"
9 #include "dk2nu/tree/dk2nu.h"
10 #include "dk2nu/tree/dkmeta.h"
12 
13 #include "TSystem.h"
14 
15 namespace evwgh {
17  public:
19  void Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine) override;
20  std::vector<std::vector<double>> GetWeight(art::Event& e) override;
21 
22  private:
23  std::string fGenieModuleLabel;
24 
25  std::vector<std::string> fInputLabels;
26  std::string fPPFXMode;
27  std::string fMode;
28  std::string fHorn;
29  std::string fTarget;
30  int fSeed;
31  int fVerbose;
32  NeutrinoFluxReweight::MakeReweight* fPPFXrw;
33 
35  };
36 
38 
40  CLHEP::HepRandomEngine& engine)
41  {
42  //get configuration for this function
44  fGenieModuleLabel = p.get<std::string>("genie_module_label");
45 
46  //ppfx setup
47  fInputLabels = pset.get<std::vector<std::string>>("input_labels");
48  fPPFXMode = pset.get<std::string>("ppfx_mode");
49  fVerbose = pset.get<int>("verbose");
50  fMode = pset.get<std::string>("mode");
51  fHorn = pset.get<std::string>("horn_curr");
52  fTarget = pset.get<std::string>("target_config");
53  fSeed = pset.get<int>("random_seed", -1);
54 
55  gSystem->Setenv("MODE", fPPFXMode.c_str());
56 
57  fPPFXrw = NeutrinoFluxReweight::MakeReweight::getInstance();
58  std::cout << "PPFX instance " << fPPFXrw << std::endl;
59 
60  std::string inputOptions; // Full path.
61  std::string mode_file = "inputs_" + fPPFXMode + ".xml"; // Just file name.
62  cet::search_path sp("FW_SEARCH_PATH");
63  sp.find_file(mode_file, inputOptions);
64 
65  if (fSeed != -1) fPPFXrw->setBaseSeed(fSeed); // Set the random seed
66  std::cout << "is PPFX setup : " << fPPFXrw->AlreadyInitialized() << std::endl;
67  std::cout << "Setting PPFX, inputs: " << inputOptions << std::endl;
68  std::cout << "Setting Horn Current Configuration to: " << fHorn << std::endl;
69  std::cout << "Setting Target Configuration to: " << fTarget << std::endl;
70  if (!(fPPFXrw->AlreadyInitialized())) { fPPFXrw->SetOptions(inputOptions); }
71  std::cout << "PPFX just set with mode: " << fPPFXMode << std::endl;
72  }
73 
74  std::vector<std::vector<double>> PPFXTargAttenWeightCalc::GetWeight(art::Event& e)
75  {
76  std::vector<std::vector<double>> weight;
78  //calculate weight(s) here
79 
80  int nmctruth = 0, nmatched = 0;
81  bool flag = true;
82  int ievt = -1;
83  std::vector<art::Handle<std::vector<bsim::Dk2Nu>>> dk2nus2 =
84  e.getMany<std::vector<bsim::Dk2Nu>>();
85  for (size_t dk2 = 0; dk2 < dk2nus2.size(); ++dk2) {
86  art::Handle<std::vector<bsim::Dk2Nu>> dk2nus = dk2nus2[dk2];
87  }
88  while ((flag = mcitr.Next())) {
89  std::string label = mcitr.GetLabel();
90  const simb::MCTruth* pmctruth = mcitr.GetMCTruth();
91  // const simb::GTruth* pgtruth = mcitr.GetGTruth();
92  //not-used//const simb::MCFlux* pmcflux = mcitr.GetMCFlux();
93  const bsim::Dk2Nu* pdk2nu = mcitr.GetDk2Nu();
94  //not-used//const bsim::NuChoice* pnuchoice = mcitr.GetNuChoice();
95  // art::Ptr<simb::MCTruth> mctruthptr = mcitr.GetMCTruthPtr();
96 
97  ++ievt;
98  ++nmctruth;
99  if (fVerbose > 0) {
100  std::cout << "FluxWeightCalculator [" << std::setw(4) << ievt << "] "
101  << " label \"" << label << "\" MCTruth@ " << pmctruth << " Dk2Nu@ " << pdk2nu
102  << std::endl;
103  }
104  if (!pdk2nu) continue;
105  ++nmatched;
106 
107  //RWH//bsim::Dk2Nu* tmp_dk2nu = fTmpDK2NUConversorAlg.construct_toy_dk2nu( &mctruth,&mcflux);
108  //RWH//bsim::DkMeta* tmp_dkmeta = fTmpDK2NUConversorAlg.construct_toy_dkmeta(&mctruth,&mcflux);
109  //RWH// those appear to have been memory leaks
110  //RWH// herein begins the replacment for the "construct_toy_dkmeta"
111 
112  static bsim::DkMeta
113  dkmeta_obj; //RWH// create this on stack (destroyed when out-of-scope) ... or static
114  dkmeta_obj.tgtcfg = fTarget;
115  dkmeta_obj.horncfg = fHorn;
116  (dkmeta_obj.vintnames).push_back("Index_Tar_In_Ancestry");
117  (dkmeta_obj.vintnames).push_back("Playlist");
118  bsim::DkMeta* tmp_dkmeta = &dkmeta_obj;
119  //RWH// herein ends this block
120 
121  // sigh ....
122  //RWH// this is the signature in NeutrinoFluxReweight::MakeReweight :
123  // //! create an interaction chain from the new dk2nu(dkmeta) format
124  // void calculateWeights(bsim::Dk2Nu* nu, bsim::DkMeta* meta);
125  //RWH// these _should_ be "const <class>*" because we don't need to change them
126  //RWH// and the pointers we get out of the ART record are going to be const.
127  bsim::Dk2Nu* tmp_dk2nu = const_cast<bsim::Dk2Nu*>(pdk2nu); // remove const-ness
128  try {
129  fPPFXrw->calculateWeights(tmp_dk2nu, tmp_dkmeta);
130  }
131  catch (...) {
132  std::cerr << "Failed to calcualte weight" << std::endl;
133  continue;
134  }
135  //Get weights:
136  if (fMode == "reweight") {
137  double ppfx_cv_wgt = fPPFXrw->GetCVWeight();
138  std::vector<double> wvec = {ppfx_cv_wgt};
139  weight.push_back(wvec);
140  }
141  else {
142 
143  std::vector<double> vtgtatt = fPPFXrw->GetWeights("TargetAttenuation");
144 
145  std::vector<double> tmp_vhptot;
146  for (unsigned int iuniv = 0; iuniv < vtgtatt.size(); iuniv++) {
147  tmp_vhptot.push_back(float(vtgtatt[iuniv]));
148  }
149  weight.push_back(tmp_vhptot);
150  }
151  }
152  return weight;
153  }
155 }
std::vector< std::string > fInputLabels
std::vector< std::vector< double > > GetWeight(art::Event &e) override
const simb::MCTruth * GetMCTruth() const
void Configure(fhicl::ParameterSet const &p, CLHEP::HepRandomEngine &engine) override
std::string GetName()
Definition: WeightCalc.h:26
#define DECLARE_WEIGHTCALC(wghcalc)
T get(std::string const &key) const
Definition: ParameterSet.h:314
NeutrinoFluxReweight::MakeReweight * fPPFXrw
double weight
Definition: plottest35.C:25
Event generator information.
Definition: MCTruth.h:32
Float_t e
Definition: plot.C:35
std::string GetLabel() const
const bsim::Dk2Nu * GetDk2Nu() const
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
#define REGISTER_WEIGHTCALC(wghcalc)