FluxWeightCalculator_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file TestFluxWeightCalculator_module.cc
3 // \brief Handle PPFX Weights and store FluxWeights.
4 // \author Leo Aliaga - laliaga@fnal.gov
5 ////////////////////////////////////////////////////////////////////////
6 
7 // NOvASoft includes
8 #include "Utilities/AssociationUtil.h"
9 
11 
12 // hack-a-lishioous from ROOT
13 #include "TSystem.h"
14 
15 // Framework includes
22 #include "fhiclcpp/ParameterSet.h"
23 
24 //RWH//#include "FluxWeightCalculator/TmpDK2NUConversorAlg.h"
25 //RWH// all these were slightly modified
26 //RWH// (or should be when proper paths are included, e.g. MakeReweight
27 //RWH// is part of ppfx ... but that isn't clear)
28 //GSD// Can at least, for now point to same place for SRT and CMAKE
29 #include "MakeReweight.h"
30 #include "dk2nu/tree/dk2nu.h"
31 #include "dk2nu/tree/dkmeta.h"
32 #include "MCReweight/FluxWeights.h"
33 
34 /* RWH not necessary with above headers
35 namespace NeutrinoFluxReweight{
36  class MakeReweight;
37 }
38 namespace bsim {
39  class Dk2Nu;
40  class DkMeta;
41  class NuRay;
42  class Decay;
43  class Ancestor;
44  class TgtExit;
45  class Traj;
46  class Location;
47 }
48 */
52 
53 namespace fxwgt
54 {
56  {
57  public:
58  explicit FluxWeightCalculator(const fhicl::ParameterSet& pset);
60 
61  virtual void reconfigure(const fhicl::ParameterSet& pset);
62  virtual void produce(art::Event& evt);
63  private:
64 
65  //RWH// allow a set list of input labels to process, or none (==all)
66  std::vector<std::string> fInputLabels; /// label(s) for generator information
67  std::string fPPFXMode; /// PPFX mode.
68  int fVerbose; ///Set verbose level.
69  NeutrinoFluxReweight::MakeReweight* makerew;
70  protected:
71  //RWH//TmpDK2NUConversorAlg fTmpDK2NUConversorAlg;
72  };
73 
74  //......................................................................
76  {
77  reconfigure(pset);
78  produces<std::vector<fxwgt::FluxWeights>>();
79  //RWH ?? should this be other direction??
80  produces<art::Assns<fxwgt::FluxWeights, simb::MCTruth>>();
81  //RWH// fails because CreateAssn is stupid .. you don't just give it to ptrs
82  //produces<art::Assns<simb::MCTruth,fxwgt::FluxWeights>>();
83  }
84 
85  //......................................................................
87  {
88  }
89 
90  //......................................................................
92  {
93  std::cerr << "FluxWeightCalculator::reconfigure(const fhicl::ParameterSet& pset)" << std::endl;
94  fInputLabels = pset.get<std::vector<std::string>>("InputLabels");
95  fPPFXMode = pset.get<std::string>("PPFXMode");
96  fVerbose = pset.get<int>("Verbose");
97  const char* ppfx_Dir = std::getenv("PPFX_DIR");
98  const char* ppfx = std::getenv("PPFX");
99 
100  // internally PPFX core uses ${PPFX_DIR} to be top everything
101  // with subdirs: src include scripts uncertainties data
102 
103  // $ echo $PPFX
104  // /grid/fermiapp/products/nova/externals/ppfx/v01_09/Linux64bit+2.6-2.12-e9-r5-debug
105  // $ echo $PPFX_DIR
106  // /grid/fermiapp/products/nova/externals/ppfx/v01_09
107  //
108  // but those subdirs are in equivalent of ${PPFX}/ppfx
109 
110 
111  // write code that does the right thing _if_ things are correct as well as current state
112 
113  std::cerr << " ppfx_Dir is \"" << ppfx_Dir << "\"" << std::endl;
114  //std::cerr << " ppfx is \"" << ppfx << "\"" << std::endl; //this variable is actually used.
115 
116  std::string right_PPFX_DIR, testParamPath;
117  bool not_found;
118 
119  right_PPFX_DIR = std::string(ppfx_Dir);
120  // virtual Bool_t AccessPathName (const char *path, EAccessMode mode=kFileExists)
121  // Returns FALSE if one can access a file using the specified access mode. More...
122  testParamPath = right_PPFX_DIR + "/uncertainties/Parameters_default.xml";
123  not_found = gSystem->AccessPathName(testParamPath.c_str()); // true if NOT found
124  if ( not_found ) {
125  //RWH// hack because env variables aren't right
126  // try this path
127  std::cerr << " huh, first try didn't work using \"" << right_PPFX_DIR << "\"" << std::endl;
128  right_PPFX_DIR = std::string(ppfx) + std::string("/ppfx");
129  testParamPath = right_PPFX_DIR + "/uncertainties/Parameters_default.xml";
130  not_found = gSystem->AccessPathName(testParamPath.c_str()); // true if NOT found
131  if ( not_found ) {
132  std::cerr << " tried PPFX_DIR as \"" << right_PPFX_DIR << "\", still not right"
133  << std::endl
134  << " everything is likely to end in tears..."
135  << std::endl;
136  } else {
137  gSystem->Setenv("PPFX_DIR", right_PPFX_DIR.c_str());
138  std::cerr << " setenv PPFX_DIR to \"" << right_PPFX_DIR << "\"" << std::endl;
139  }
140  }
141  std::cerr << " using right_PPFX_DIR is \"" << right_PPFX_DIR << "\"" << std::endl;
142 
143  std::cerr << " about to getInstance()" << std::endl;
144  makerew = NeutrinoFluxReweight::MakeReweight::getInstance();
145 
146  std::string inputOptions = right_PPFX_DIR +"/scripts/inputs_"+fPPFXMode+".xml";
147  not_found = gSystem->AccessPathName(inputOptions.c_str());
148  if ( not_found ) {
149  std::cerr << " huh, try to find the location of input_defaul.xml\"" << "\"" << std::endl;
150  inputOptions = right_PPFX_DIR +"/xml/inputs_"+fPPFXMode+".xml";
151  not_found = gSystem->AccessPathName(inputOptions.c_str());
152  if ( not_found ) {
153  std::cerr << " tried to locate input_default.xml but it is not at "<<inputOptions<<std::endl;
154  }
155  else {
156  std::cout << "input_default.xml found at "<<inputOptions<<std::endl;
157  }
158  }
159 
160  std::cout << "Setting PPFX, inputs: " << inputOptions << std::endl;
161  makerew->SetOptions(inputOptions);
162  std::cout << "PPFX just set with mode: " << fPPFXMode << std::endl;
163  }
164 
165  //......................................................................
167 
168  std::unique_ptr<std::vector<fxwgt::FluxWeights>>
169  rescol(new std::vector<fxwgt::FluxWeights>);
170 
171  //std::unique_ptr<art::Assns<simb::MCTruth, fxwgt::FluxWeights>>
172  // assns(new art::Assns<simb::MCTruth, fxwgt::FluxWeights>);
173  std::unique_ptr<art::Assns<fxwgt::FluxWeights, simb::MCTruth>>
175 
176 
177  int verbose=fVerbose;
178  int nmctruth=0, nmatched=0;
179 
180  // get a fancy iterator that hides a bunch of association stuff
181  evg::MCTruthToDk2NuHackItr mcitr(evt,fInputLabels,verbose);
182 
183  bool flag = true;
184  int ievt = -1;
185  while ( ( flag = mcitr.Next() ) ) {
186  std::string label = mcitr.GetLabel();
187  const simb::MCTruth* pmctruth = mcitr.GetMCTruth();
188  // const simb::GTruth* pgtruth = mcitr.GetGTruth();
189  //not-used//const simb::MCFlux* pmcflux = mcitr.GetMCFlux();
190  const bsim::Dk2Nu* pdk2nu = mcitr.GetDk2Nu();
191  //not-used//const bsim::NuChoice* pnuchoice = mcitr.GetNuChoice();
192 
193  art::Ptr<simb::MCTruth> mctruthptr = mcitr.GetMCTruthPtr();
194 
195  ++ievt;
196  ++nmctruth;
197  if ( verbose > 0 ) {
198  std::cout << "FluxWeightCalculator [" << std::setw(4) << ievt << "] "
199  << " label \"" << label << "\" MCTruth@ " << pmctruth
200  << " Dk2Nu@ " << pdk2nu << std::endl;
201  }
202  if ( ! pdk2nu ) continue;
203  ++nmatched;
204 
205  //RWH//bsim::Dk2Nu* tmp_dk2nu = fTmpDK2NUConversorAlg.construct_toy_dk2nu( &mctruth,&mcflux);
206  //RWH//bsim::DkMeta* tmp_dkmeta = fTmpDK2NUConversorAlg.construct_toy_dkmeta(&mctruth,&mcflux);
207  //RWH// those appear to have been memory leaks
208  //RWH// herein begins the replacment for the "construct_toy_dkmeta"
209  static bsim::DkMeta dkmeta_obj; //RWH// create this on stack (destroyed when out-of-scope) ... or static
210  dkmeta_obj.tgtcfg = "me000z"; //RWH// this should come from FCL parameter
211  dkmeta_obj.horncfg = "200i"; //RWH// as should this
212  (dkmeta_obj.vintnames).push_back("Index_Tar_In_Ancestry");
213  (dkmeta_obj.vintnames).push_back("Playlist");
214  bsim::DkMeta* tmp_dkmeta = &dkmeta_obj;
215  //RWH// herein ends this block
216 
217  // sigh ....
218  //RWH// this is the signature in NeutrinoFluxReweight::MakeReweight :
219  // //! create an interaction chain from the new dk2nu(dkmeta) format
220  // void calculateWeights(bsim::Dk2Nu* nu, bsim::DkMeta* meta);
221  //RWH// these _should_ be "const <class>*" because we don't need to change them
222  //RWH// and the pointers we get out of the ART record are going to be const.
223  bsim::Dk2Nu* tmp_dk2nu = const_cast<bsim::Dk2Nu*>(pdk2nu); // remove const-ness
224 
225  makerew->calculateWeights(tmp_dk2nu,tmp_dkmeta);
226 
227  //Get weights:
228  double ppfx_cv_wgt = makerew->GetCVWeight();
229  std::vector<double> vmipppion = makerew->GetWeights("MIPPNumiPionYields");
230  std::vector<double> vmippkaon = makerew->GetWeights("MIPPNumiKaonYields");
231  std::vector<double> vtgtatt = makerew->GetWeights("TargetAttenuation");
232  std::vector<double> vabsorp = makerew->GetWeights("TotalAbsorption");
233  std::vector<double> vttpcpion = makerew->GetWeights("ThinTargetpCPion");
234  std::vector<double> vttpckaon = makerew->GetWeights("ThinTargetpCKaon");
235  std::vector<double> vttpcnucleon = makerew->GetWeights("ThinTargetpCNucleon");
236  std::vector<double> vttncpion = makerew->GetWeights("ThinTargetnCPion");
237  std::vector<double> vttnucleona = makerew->GetWeights("ThinTargetnucleonA");
238  std::vector<double> vttmesoninc = makerew->GetWeights("ThinTargetMesonIncident");
239  std::vector<double> vothers = makerew->GetWeights("Other");
240 
241  std::vector<float> tmp_vhptot;
242  std::vector<float> tmp_vbftot;
243  for(unsigned int iuniv=0;iuniv<vtgtatt.size();iuniv++){
244  tmp_vhptot.push_back(float(vmipppion[iuniv]*vmippkaon[iuniv]*vtgtatt[iuniv]*vabsorp[iuniv]*vttpcpion[iuniv]*
245  vttpckaon[iuniv]*vttpcnucleon[iuniv]*vttncpion[iuniv]*vttnucleona[iuniv]*
246  vttmesoninc[iuniv]*vothers[iuniv]));
247  tmp_vbftot.push_back(1.0);
248  }
249 
250  fxwgt::FluxWeights fluxweights;
251  fluxweights.SetHadronProductionCentralValue(float(ppfx_cv_wgt));
252  fluxweights.SetBeamFocusingCentralValue(1.0);
253  fluxweights.SetHadronProductionMultiUniverses(tmp_vhptot);
254  fluxweights.SetBeamFocusingMultiUniverses(tmp_vbftot);
255 
256  //Test:
257  /* std::cout<<"=> Test of new flux weigth product:"<<std::endl;
258  std::cout<<"=> cv, Nuniv: "<<ppfx_cv_wgt<<", "<<tmp_vhptot.size()<<std::endl;
259  std::cout<<"";for(unsigned int j=0;j<tmp_vhptot.size();j++)std::cout<<tmp_vhptot[j]<<" ";std::cout<<""<<std::endl;
260  std::cout<<" "<<std::endl;*/
261 
262  rescol->push_back(fluxweights);
263  // need pmctruth to be Ptr<> ??
264  util::CreateAssn(*this, evt, *rescol, mctruthptr, *assns);
265 
266  } // loop over MCTruthToDk2NuHackItr items
267 
268  if ( verbose > 0 ) {
269  std::cout << "FluxWeightCalculator saw " << nmctruth
270  << " MCTruth obj, matched " << nmatched
271  << " with Dk2Nu and generated "
272  << rescol->size() << " fxwgt::FluxWeight obj "
273  << std::endl;
274  }
275 
276  evt.put(std::move(rescol));
277  evt.put(std::move(assns));
278 
279  } // end of produce() method
280 
282 
283 } // namespace fxwgt
284 ////////////////////////////////////////////////////////////////////////
void SetBeamFocusingMultiUniverses(std::vector< float > vuniv)
Set a vector of beam focusing weights. Each entry represents the correction for one universe...
Definition: FluxWeights.cxx:30
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.
std::vector< std::string > vintnames
names of elements for user defined vector of integers
Definition: dkmeta.h:126
void SetBeamFocusingCentralValue(float cv)
Set beam focusing central value correction from PPFX.
Definition: FluxWeights.cxx:24
virtual void reconfigure(const fhicl::ParameterSet &pset)
OStream cerr
Definition: OStream.cxx:7
NeutrinoFluxReweight::MakeReweight * makerew
Set verbose level.
DEFINE_ART_MODULE(TestTMapFile)
object containing MC flux information
const char * label
const simb::MCTruth * GetMCTruth() const
std::string tgtcfg
target config e.g. "minos/epoch3/-10cm"
Definition: dkmeta.h:98
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Store flux weigths for neutrino correction.
Definition: FluxWeights.h:15
art::Ptr< simb::MCTruth > GetMCTruthPtr() const
std::string getenv(std::string const &name)
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
base_types push_back(int_type())
const bsim::Dk2Nu * GetDk2Nu() const
FluxWeightCalculator(const fhicl::ParameterSet &pset)
void SetHadronProductionCentralValue(float cv)
Set HP central value correction from PPFX.
Definition: FluxWeights.cxx:21
OStream cout
Definition: OStream.cxx:6
std::string horncfg
horn config e.g. "FHC/185A/LE/h1xoff=1mm"
Definition: dkmeta.h:99
std::vector< std::string > fInputLabels
virtual void produce(art::Event &evt)
Event generator information.
Definition: MCTruth.h:32
std::string fPPFXMode
label(s) for generator information
void SetHadronProductionMultiUniverses(std::vector< float > vuniv)
Set a vector of HP weights. Each entry represents the correction for one universe.
Definition: FluxWeights.cxx:27
enum BeamMode string