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

Public Member Functions

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

Constructor & Destructor Documentation

evwgh::PPFXCVWeightCalc::PPFXCVWeightCalc ( )

Definition at line 37 of file PPFXCVWeightCalc.cxx.

37 {}

Member Function Documentation

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

Implements evwgh::WeightCalc.

Definition at line 39 of file PPFXCVWeightCalc.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 
58  std::string inputOptions; // Full path.
59  std::string mode_file = "inputs_" + fPPFXMode + ".xml"; // Just file name.
60  cet::search_path sp("FW_SEARCH_PATH");
61  sp.find_file(mode_file, inputOptions);
62 
63  if (fSeed != -1) fPPFXrw->setBaseSeed(fSeed); // Set the random seed
64  std::cout << "is PPFX setup : " << fPPFXrw->AlreadyInitialized() << std::endl;
65  std::cout << "Setting PPFX, inputs: " << inputOptions << std::endl;
66  std::cout << "Setting Horn Current Configuration to: " << fHorn << std::endl;
67  std::cout << "Setting Target Configuration to: " << fTarget << std::endl;
68  if (!(fPPFXrw->AlreadyInitialized())) { fPPFXrw->SetOptions(inputOptions); }
69  std::cout << "PPFX just set with mode: " << fPPFXMode << std::endl;
70  }
evwgh::PPFXCVWeightCalc::fPPFXrw
NeutrinoFluxReweight::MakeReweight * fPPFXrw
Definition: PPFXCVWeightCalc.cxx:32
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::PPFXCVWeightCalc::fInputLabels
std::vector< std::string > fInputLabels
Definition: PPFXCVWeightCalc.cxx:25
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(), 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::PPFXCVWeightCalc::GetWeight ( art::Event &  e)
overridevirtual

Implements evwgh::WeightCalc.

Definition at line 72 of file PPFXCVWeightCalc.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.

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

Definition at line 23 of file PPFXCVWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::PPFXCVWeightCalc::fHorn
private

Definition at line 28 of file PPFXCVWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

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

Definition at line 25 of file PPFXCVWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::string evwgh::PPFXCVWeightCalc::fMode
private

Definition at line 27 of file PPFXCVWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

std::string evwgh::PPFXCVWeightCalc::fPPFXMode
private

Definition at line 26 of file PPFXCVWeightCalc.cxx.

Referenced by Configure().

NeutrinoFluxReweight::MakeReweight* evwgh::PPFXCVWeightCalc::fPPFXrw
private

Definition at line 32 of file PPFXCVWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

int evwgh::PPFXCVWeightCalc::fSeed
private

Definition at line 30 of file PPFXCVWeightCalc.cxx.

Referenced by Configure().

std::string evwgh::PPFXCVWeightCalc::fTarget
private

Definition at line 29 of file PPFXCVWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

int evwgh::PPFXCVWeightCalc::fVerbose
private

Definition at line 31 of file PPFXCVWeightCalc.cxx.

Referenced by Configure(), and GetWeight().

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