LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Static Public Member Functions | Private Attributes | List of all members
evwgh::PPFXMIPPPionWeightCalc Class Reference
Inheritance diagram for evwgh::PPFXMIPPPionWeightCalc:
evwgh::WeightCalc

Public Member Functions

 PPFXMIPPPionWeightCalc ()
 
void Configure (fhicl::ParameterSet const &p, CLHEP::HepRandomEngine &engine) override
 
std::vector< std::vector< double > > GetWeight (art::Event &e) override
 
void SetName (std::string name)
 
std::string GetName ()
 

Static Public Member Functions

static std::vector< std::vector< double > > MultiGaussianSmearing (std::vector< double > const &centralValues, std::vector< std::vector< double >> const &inputCovarianceMatrix, int n_multisims, CLHEP::RandGaussQ &GaussRandom)
 Applies Gaussian smearing to a set of data. More...
 
static std::vector< double > MultiGaussianSmearing (std::vector< double > const &centralValue, TMatrixD *const &inputCovarianceMatrix, std::vector< double > rand)
 Over load of the above function that only returns a single varied parameter set. More...
 
static std::vector< double > MultiGaussianSmearing (std::vector< double > const &centralValue, TMatrixD *const &LowerTriangleCovarianceMatrix, bool isDecomposed, std::vector< double > rand)
 

Private Attributes

std::string fGenieModuleLabel
 
std::vector< std::string > fInputLabels
 
std::string fPPFXMode
 
std::string fMode
 
std::string fHorn
 
std::string fTarget
 
int fSeed
 
int fVerbose
 
NeutrinoFluxReweight::MakeReweight * fPPFXrw
 

Detailed Description

Definition at line 16 of file PPFXMIPPPionWeightCalc.cxx.

Constructor & Destructor Documentation

evwgh::PPFXMIPPPionWeightCalc::PPFXMIPPPionWeightCalc ( )

Definition at line 37 of file PPFXMIPPPionWeightCalc.cxx.

37 {}

Member Function Documentation

void evwgh::PPFXMIPPPionWeightCalc::Configure ( fhicl::ParameterSet const &  p,
CLHEP::HepRandomEngine &  engine 
)
overridevirtual

Implements evwgh::WeightCalc.

Definition at line 39 of file PPFXMIPPPionWeightCalc.cxx.

References fGenieModuleLabel, fHorn, fInputLabels, fMode, fPPFXMode, fPPFXrw, fSeed, fTarget, fVerbose, fhicl::ParameterSet::get(), and evwgh::WeightCalc::GetName().

41  {
42  //get configuration for this function
43  fhicl::ParameterSet const& pset = p.get<fhicl::ParameterSet>(GetName());
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  }
evwgh::PPFXMIPPPionWeightCalc::fInputLabels
std::vector< std::string > fInputLabels
Definition: PPFXMIPPPionWeightCalc.cxx:25
evwgh::WeightCalc::GetName
std::string GetName()
Definition: WeightCalc.h:26
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:314
evwgh::PPFXMIPPPionWeightCalc::fPPFXrw
NeutrinoFluxReweight::MakeReweight * fPPFXrw
Definition: PPFXMIPPPionWeightCalc.cxx:32
std::string evwgh::WeightCalc::GetName ( )
inlineinherited

Definition at line 26 of file WeightCalc.h.

References lar::dump::vector().

Referenced by evwgh::GenieWeightCalc::CheckForIncompatibleSystematics(), evwgh::PPFXWeightCalc::Configure(), evwgh::PPFXCVWeightCalc::Configure(), evwgh::PPFXMIPPKaonWeightCalc::Configure(), Configure(), evwgh::PPFXOtherWeightCalc::Configure(), evwgh::PPFXTargAttenWeightCalc::Configure(), evwgh::PPFXThinKaonWeightCalc::Configure(), evwgh::PPFXThinMesonWeightCalc::Configure(), evwgh::PPFXThinNeutronPionWeightCalc::Configure(), evwgh::PPFXThinNucAWeightCalc::Configure(), evwgh::PPFXThinNucWeightCalc::Configure(), evwgh::PPFXThinPionWeightCalc::Configure(), evwgh::PPFXTotAbsorpWeightCalc::Configure(), and evwgh::GenieWeightCalc::Configure().

26 { return fName; }
evwgh::WeightCalc::fName
std::string fName
Definition: WeightCalc.h:54
std::vector< std::vector< double > > evwgh::PPFXMIPPPionWeightCalc::GetWeight ( art::Event &  e)
overridevirtual

Implements evwgh::WeightCalc.

Definition at line 74 of file PPFXMIPPPionWeightCalc.cxx.

References fHorn, fInputLabels, fMode, fPPFXrw, fTarget, fVerbose, evgb::MCTruthAndFriendsItr::GetDk2Nu(), evgb::MCTruthAndFriendsItr::GetLabel(), art::ProductRetriever::getMany(), evgb::MCTruthAndFriendsItr::GetMCTruth(), evgb::MCTruthAndFriendsItr::Next(), REGISTER_WEIGHTCALC, and weight.

75  {
76  std::vector<std::vector<double>> weight;
77  evgb::MCTruthAndFriendsItr mcitr(e, fInputLabels);
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  std::vector<double> vmipppion = fPPFXrw->GetWeights("MIPPNumiPionYields");
143 
144  std::vector<double> tmp_vhptot;
145  for (unsigned int iuniv = 0; iuniv < vmipppion.size(); iuniv++) {
146  tmp_vhptot.push_back(float(vmipppion[iuniv]));
147  }
148  weight.push_back(tmp_vhptot);
149  }
150  }
151  return weight;
152  }
evwgh::PPFXMIPPPionWeightCalc::fInputLabels
std::vector< std::string > fInputLabels
Definition: PPFXMIPPPionWeightCalc.cxx:25
weight
double weight
Definition: plottest35.C:25
evwgh::PPFXMIPPPionWeightCalc::fPPFXrw
NeutrinoFluxReweight::MakeReweight * fPPFXrw
Definition: PPFXMIPPPionWeightCalc.cxx:32
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
art::ProductRetriever::getMany
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: ProductRetriever.h:303
std::vector< std::vector< double > > evwgh::WeightCalc::MultiGaussianSmearing ( std::vector< double > const &  centralValues,
std::vector< std::vector< double >> const &  inputCovarianceMatrix,
int  n_multisims,
CLHEP::RandGaussQ &  GaussRandom 
)
staticinherited

Applies Gaussian smearing to a set of data.

Parameters
centralValuesthe values to be smeared
inputCovarianceMatrixcovariance matrix for smearing
n_multisimsnumber of sets of smeared values to be produced
Returns
a set of n_multisims value sets smeared from the central value

If centralValues is of dimension N, inputCovarianceMatrix needs to be NxN, and each of the returned data sets will be also of dimension N.

Definition at line 14 of file WeightCalc.cxx.

References col, art::errors::Configuration, n, and art::errors::StdException.

19  {
20 
21  std::vector<std::vector<double>> setOfSmearedCentralValues;
22 
23  //Check that covarianceMatrix is of compatible dimension with the central values
24  unsigned int covarianceMatrix_dim = centralValue.size();
25 
26  if (inputCovarianceMatrix.size() != covarianceMatrix_dim) {
27  throw art::Exception(art::errors::Configuration)
28  << "inputCovarianceMatrix rows " << inputCovarianceMatrix.size()
29  << " not equal to entries of centralValue[]: " << covarianceMatrix_dim;
30  }
31 
32  for (auto const& row : inputCovarianceMatrix) { // Range-for (C++11).
33  if (row.size() != covarianceMatrix_dim) {
34  throw art::Exception(art::errors::Configuration)
35  << "inputCovarianceMatrix columns " << row.size()
36  << " not equal to entries of centralValue[]: " << covarianceMatrix_dim;
37  }
38  }
39 
40  //change covariance matrix into a TMartixD object
41  int dim = int(inputCovarianceMatrix.size());
42  TMatrixD covarianceMatrix(dim, dim);
43  int i = 0;
44  for (auto const& row : inputCovarianceMatrix) {
45  int j = 0;
46  for (auto const& element : row) {
47  covarianceMatrix[i][j] = element;
48  ++j;
49  }
50  ++i;
51  }
52 
53  //perform Choleskey Decomposition
54  TDecompChol dc = TDecompChol(covarianceMatrix);
55  if (!dc.Decompose()) {
56  throw art::Exception(art::errors::StdException)
57  << "Cannot decompose covariance matrix to begin smearing.";
58  return setOfSmearedCentralValues;
59  }
60 
61  //Get upper triangular matrix. This maintains the relations in the
62  //covariance matrix, but simplifies the structure.
63  TMatrixD U = dc.GetU();
64 
65  //for every multisim
66  for (int n = 0; n < n_multisims; ++n) {
67 
68  //get a gaussian random number for every central value element
69  int dim = centralValue.size();
70 
71  std::vector<double> rands(dim);
72  GaussRandom.fireArray(dim, rands.data());
73 
74  //compute the smeared central values
75  std::vector<double> smearedCentralValues;
76  for (int col = 0; col < dim; ++col) {
77  //find the weight of each col of the upper triangular cov. matrix
78  double weightFromU = 0.;
79  for (int row = 0; row < col + 1; ++row) {
80  weightFromU += U(row, col) * rands[row];
81  }
82  //multiply this weight by each col of the central values to obtain
83  //the gaussian smeared and constrained central values
84  // smearedCentralValues.push_back(weightFromU * centralValue[col]);
85  smearedCentralValues.push_back(weightFromU + centralValue[col]);
86  }
87 
88  //collect the smeared central values into a set
89  setOfSmearedCentralValues.push_back(smearedCentralValues);
90  } //loop over multisims
91  return setOfSmearedCentralValues;
92  } // WeightCalc::MultiGaussianSmearing() - gives a vector of sets of smeared parameters
col
Int_t col[ntarg]
Definition: Style.C:29
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
n
Char_t n[5]
Definition: comparison_ascii.C:55
std::vector< double > evwgh::WeightCalc::MultiGaussianSmearing ( std::vector< double > const &  centralValue,
TMatrixD *const &  inputCovarianceMatrix,
std::vector< double >  rand 
)
staticinherited

Over load of the above function that only returns a single varied parameter set.

Definition at line 97 of file WeightCalc.cxx.

References col, and art::errors::StdException.

100  {
101 
102  //compute the smeared central values
103  std::vector<double> smearedCentralValues;
104 
105  //
106  //perform Choleskey Decomposition
107  //
108  // Best description I have found
109  // is in the PDG (Monte Carlo techniques, Algorithms, Gaussian distribution)
110  //
111  // http://pdg.lbl.gov/2016/reviews/rpp2016-rev-monte-carlo-techniques.pdf (Page 5)
112  //
113  //
114  TDecompChol dc = TDecompChol(*(inputCovarianceMatrix));
115  if (!dc.Decompose()) {
116  throw art::Exception(art::errors::StdException)
117  << "Cannot decompose covariance matrix to begin smearing.";
118  return smearedCentralValues;
119  }
120 
121  //Get upper triangular matrix. This maintains the relations in the
122  //covariance matrix, but simplifies the structure.
123  TMatrixD U = dc.GetU();
124 
125  for (unsigned int col = 0; col < centralValue.size(); ++col) {
126  //find the weight of each col of the upper triangular cov. matrix
127  double weightFromU = 0.;
128 
129  for (unsigned int row = 0; row < col + 1; ++row) {
130  weightFromU += U(row, col) * rand[row];
131  }
132 
133  //multiply this weight by each col of the central values to obtain
134  //the gaussian smeared and constrained central values
135  // smearedCentralValues.push_back(weightFromU * centralValue[col]);
136  smearedCentralValues.push_back(weightFromU + centralValue[col]);
137  }
138  return smearedCentralValues;
139  } // WeightCalc::MultiGaussianSmearing() - Gives a single set of smeared parameters
col
Int_t col[ntarg]
Definition: Style.C:29
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::vector< double > evwgh::WeightCalc::MultiGaussianSmearing ( std::vector< double > const &  centralValue,
TMatrixD *const &  LowerTriangleCovarianceMatrix,
bool  isDecomposed,
std::vector< double >  rand 
)
staticinherited

Definition at line 151 of file WeightCalc.cxx.

References col, and art::errors::StdException.

156  {
157 
158  //compute the smeared central values
159  std::vector<double> smearedCentralValues;
160 
161  // This guards against accidentally
162  if (!isDecomposed) {
163  throw art::Exception(art::errors::StdException)
164  << "Must supply the decomposed lower triangular covariance matrix.";
165  return smearedCentralValues;
166  }
167 
168  for (unsigned int col = 0; col < centralValue.size(); ++col) {
169  //find the weight of each col of the upper triangular cov. matrix
170  double weightFromU = 0.;
171 
172  for (unsigned int row = 0; row < col + 1; ++row) {
173  weightFromU += LowerTriangleCovarianceMatrix[0][row][col] * rand[row];
174  }
175 
176  //multiply this weight by each col of the central values to obtain
177  //the gaussian smeared and constrained central values
178  // smearedCentralValues.push_back(weightFromU * centralValue[col]);
179  smearedCentralValues.push_back(weightFromU + centralValue[col]);
180  }
181  return smearedCentralValues;
182  } // WeightCalc::MultiGaussianSmearing() - Gives a single set of smeared parameters
col
Int_t col[ntarg]
Definition: Style.C:29
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
void evwgh::WeightCalc::SetName ( std::string  name)
inlineinherited

Definition at line 25 of file WeightCalc.h.

25 { fName = name; }
evwgh::WeightCalc::fName
std::string fName
Definition: WeightCalc.h:54

Member Data Documentation

std::string evwgh::PPFXMIPPPionWeightCalc::fGenieModuleLabel
private

Definition at line 23 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::PPFXMIPPPionWeightCalc::fHorn
private

Definition at line 28 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::vector<std::string> evwgh::PPFXMIPPPionWeightCalc::fInputLabels
private

Definition at line 25 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::string evwgh::PPFXMIPPPionWeightCalc::fMode
private

Definition at line 27 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::string evwgh::PPFXMIPPPionWeightCalc::fPPFXMode
private

Definition at line 26 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure().

NeutrinoFluxReweight::MakeReweight* evwgh::PPFXMIPPPionWeightCalc::fPPFXrw
private

Definition at line 32 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

int evwgh::PPFXMIPPPionWeightCalc::fSeed
private

Definition at line 30 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::PPFXMIPPPionWeightCalc::fTarget
private

Definition at line 29 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

int evwgh::PPFXMIPPPionWeightCalc::fVerbose
private

Definition at line 31 of file PPFXMIPPPionWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

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