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::PPFXWeightCalc Class Reference
Inheritance diagram for evwgh::PPFXWeightCalc:
evwgh::WeightCalc

Public Member Functions

 PPFXWeightCalc ()
 
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 PPFXWeightCalc.cxx.

Constructor & Destructor Documentation

evwgh::PPFXWeightCalc::PPFXWeightCalc ( )

Definition at line 37 of file PPFXWeightCalc.cxx.

37 {}

Member Function Documentation

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

Implements evwgh::WeightCalc.

Definition at line 39 of file PPFXWeightCalc.cxx.

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

40  {
41  //get configuration for this function
42  fhicl::ParameterSet const& pset = p.get<fhicl::ParameterSet>(GetName());
43  fGenieModuleLabel = p.get<std::string>("genie_module_label");
44 
45  //ppfx setup
46  fInputLabels = pset.get<std::vector<std::string>>("input_labels");
47  fPPFXMode = pset.get<std::string>("ppfx_mode");
48  fVerbose = pset.get<int>("verbose");
49  fMode = pset.get<std::string>("mode");
50  fHorn = pset.get<std::string>("horn_curr");
51  fTarget = pset.get<std::string>("target_config");
52  fSeed = pset.get<int>("random_seed", -1);
53 
54  gSystem->Setenv("MODE", fPPFXMode.c_str());
55 
56  fPPFXrw = NeutrinoFluxReweight::MakeReweight::getInstance();
57  std::cout << "PPFX instance " << fPPFXrw << std::endl;
58 
59  std::string inputOptions; // Full path.
60  std::string mode_file = "inputs_" + fPPFXMode + ".xml"; // Just file name.
61  cet::search_path sp("FW_SEARCH_PATH");
62  sp.find_file(mode_file, inputOptions);
63 
64  if (fSeed != -1) fPPFXrw->setBaseSeed(fSeed); // Set the random seed
65  std::cout << "is PPFX setup : " << fPPFXrw->AlreadyInitialized() << std::endl;
66  std::cout << "Setting PPFX, inputs: " << inputOptions << std::endl;
67  std::cout << "Setting Horn Current Configuration to: " << fHorn << std::endl;
68  std::cout << "Setting Target Configuration to: " << fTarget << std::endl;
69  if (!(fPPFXrw->AlreadyInitialized())) { fPPFXrw->SetOptions(inputOptions); }
70  std::cout << "PPFX just set with mode: " << fPPFXMode << std::endl;
71  }
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::PPFXWeightCalc::fInputLabels
std::vector< std::string > fInputLabels
Definition: PPFXWeightCalc.cxx:25
evwgh::PPFXWeightCalc::fGenieModuleLabel
std::string fGenieModuleLabel
Definition: PPFXWeightCalc.cxx:23
evwgh::PPFXWeightCalc::fPPFXrw
NeutrinoFluxReweight::MakeReweight * fPPFXrw
Definition: PPFXWeightCalc.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(), Configure(), evwgh::PPFXCVWeightCalc::Configure(), evwgh::PPFXThinKaonWeightCalc::Configure(), evwgh::PPFXMIPPPionWeightCalc::Configure(), evwgh::PPFXOtherWeightCalc::Configure(), evwgh::PPFXTargAttenWeightCalc::Configure(), evwgh::PPFXMIPPKaonWeightCalc::Configure(), evwgh::PPFXThinPionWeightCalc::Configure(), evwgh::PPFXThinNeutronPionWeightCalc::Configure(), evwgh::PPFXThinNucAWeightCalc::Configure(), evwgh::PPFXThinNucWeightCalc::Configure(), evwgh::PPFXThinMesonWeightCalc::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::PPFXWeightCalc::GetWeight ( art::Event &  e)
overridevirtual

Implements evwgh::WeightCalc.

Definition at line 73 of file PPFXWeightCalc.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.

74  {
75  std::vector<std::vector<double>> weight;
76  evgb::MCTruthAndFriendsItr mcitr(e, fInputLabels);
77  //calculate weight(s) here
78 
79  int nmctruth = 0, nmatched = 0;
80  bool flag = true;
81  int ievt = -1;
82  std::vector<art::Handle<std::vector<bsim::Dk2Nu>>> dk2nus2 =
83  e.getMany<std::vector<bsim::Dk2Nu>>();
84  for (size_t dk2 = 0; dk2 < dk2nus2.size(); ++dk2) {
85  art::Handle<std::vector<bsim::Dk2Nu>> dk2nus = dk2nus2[dk2];
86  }
87  while ((flag = mcitr.Next())) {
88  std::string label = mcitr.GetLabel();
89  const simb::MCTruth* pmctruth = mcitr.GetMCTruth();
90  // const simb::GTruth* pgtruth = mcitr.GetGTruth();
91  //not-used//const simb::MCFlux* pmcflux = mcitr.GetMCFlux();
92  const bsim::Dk2Nu* pdk2nu = mcitr.GetDk2Nu();
93  //not-used//const bsim::NuChoice* pnuchoice = mcitr.GetNuChoice();
94  // art::Ptr<simb::MCTruth> mctruthptr = mcitr.GetMCTruthPtr();
95 
96  ++ievt;
97  ++nmctruth;
98  if (fVerbose > 0) {
99  std::cout << "FluxWeightCalculator [" << std::setw(4) << ievt << "] "
100  << " label \"" << label << "\" MCTruth@ " << pmctruth << " Dk2Nu@ " << pdk2nu
101  << std::endl;
102  }
103  if (!pdk2nu) continue;
104  ++nmatched;
105 
106  //RWH//bsim::Dk2Nu* tmp_dk2nu = fTmpDK2NUConversorAlg.construct_toy_dk2nu( &mctruth,&mcflux);
107  //RWH//bsim::DkMeta* tmp_dkmeta = fTmpDK2NUConversorAlg.construct_toy_dkmeta(&mctruth,&mcflux);
108  //RWH// those appear to have been memory leaks
109  //RWH// herein begins the replacment for the "construct_toy_dkmeta"
110 
111  static bsim::DkMeta
112  dkmeta_obj; //RWH// create this on stack (destroyed when out-of-scope) ... or static
113  dkmeta_obj.tgtcfg = fTarget;
114  dkmeta_obj.horncfg = fHorn;
115  (dkmeta_obj.vintnames).push_back("Index_Tar_In_Ancestry");
116  (dkmeta_obj.vintnames).push_back("Playlist");
117  bsim::DkMeta* tmp_dkmeta = &dkmeta_obj;
118  //RWH// herein ends this block
119 
120  // sigh ....
121  //RWH// this is the signature in NeutrinoFluxReweight::MakeReweight :
122  // //! create an interaction chain from the new dk2nu(dkmeta) format
123  // void calculateWeights(bsim::Dk2Nu* nu, bsim::DkMeta* meta);
124  //RWH// these _should_ be "const <class>*" because we don't need to change them
125  //RWH// and the pointers we get out of the ART record are going to be const.
126  bsim::Dk2Nu* tmp_dk2nu = const_cast<bsim::Dk2Nu*>(pdk2nu); // remove const-ness
127  try {
128  fPPFXrw->calculateWeights(tmp_dk2nu, tmp_dkmeta);
129  }
130  catch (...) {
131  std::cerr << "Failed to calcualte weight" << std::endl;
132  continue;
133  }
134  //Get weights:
135  if (fMode == "reweight") {
136  double ppfx_cv_wgt = fPPFXrw->GetCVWeight();
137  std::vector<double> wvec = {ppfx_cv_wgt};
138  weight.push_back(wvec);
139  }
140  else {
141  std::vector<double> vmipppion = fPPFXrw->GetWeights("MIPPNumiPionYields");
142  std::vector<double> vmippkaon = fPPFXrw->GetWeights("MIPPNumiKaonYields");
143  std::vector<double> vtgtatt = fPPFXrw->GetWeights("TargetAttenuation");
144  std::vector<double> vabsorp = fPPFXrw->GetWeights("TotalAbsorption");
145  std::vector<double> vttpcpion = fPPFXrw->GetWeights("ThinTargetpCPion");
146  std::vector<double> vttpckaon = fPPFXrw->GetWeights("ThinTargetpCKaon");
147  std::vector<double> vttpcnucleon = fPPFXrw->GetWeights("ThinTargetpCNucleon");
148  std::vector<double> vttncpion = fPPFXrw->GetWeights("ThinTargetnCPion");
149  std::vector<double> vttnucleona = fPPFXrw->GetWeights("ThinTargetnucleonA");
150  std::vector<double> vttmesoninc = fPPFXrw->GetWeights("ThinTargetMesonIncident");
151  std::vector<double> vothers = fPPFXrw->GetWeights("Other");
152 
153  std::vector<double> tmp_vhptot;
154  for (unsigned int iuniv = 0; iuniv < vtgtatt.size(); iuniv++) {
155  tmp_vhptot.push_back(float(vmipppion[iuniv] * vmippkaon[iuniv] * vtgtatt[iuniv] *
156  vabsorp[iuniv] * vttpcpion[iuniv] * vttpckaon[iuniv] *
157  vttpcnucleon[iuniv] * vttncpion[iuniv] * vttnucleona[iuniv] *
158  vttmesoninc[iuniv] * vothers[iuniv]));
159  }
160  weight.push_back(tmp_vhptot);
161  }
162  }
163  return weight;
164  }
evwgh::PPFXWeightCalc::fInputLabels
std::vector< std::string > fInputLabels
Definition: PPFXWeightCalc.cxx:25
weight
double weight
Definition: plottest35.C:25
evwgh::PPFXWeightCalc::fPPFXrw
NeutrinoFluxReweight::MakeReweight * fPPFXrw
Definition: PPFXWeightCalc.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::PPFXWeightCalc::fGenieModuleLabel
private

Definition at line 23 of file PPFXWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::PPFXWeightCalc::fHorn
private

Definition at line 28 of file PPFXWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

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

Definition at line 25 of file PPFXWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::string evwgh::PPFXWeightCalc::fMode
private

Definition at line 27 of file PPFXWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::string evwgh::PPFXWeightCalc::fPPFXMode
private

Definition at line 26 of file PPFXWeightCalc.cxx.

Referenced by Configure().

NeutrinoFluxReweight::MakeReweight* evwgh::PPFXWeightCalc::fPPFXrw
private

Definition at line 32 of file PPFXWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

int evwgh::PPFXWeightCalc::fSeed
private

Definition at line 30 of file PPFXWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::PPFXWeightCalc::fTarget
private

Definition at line 29 of file PPFXWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

int evwgh::PPFXWeightCalc::fVerbose
private

Definition at line 31 of file PPFXWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

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