Inheritance diagram for evwgh::PPFXThinNeutronPionWeightCalc:

Public Member Functions
	PPFXThinNeutronPionWeightCalc ()

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 PPFXThinNeutronPionWeightCalc.cxx.

Constructor & Destructor Documentation

evwgh::PPFXThinNeutronPionWeightCalc::PPFXThinNeutronPionWeightCalc ( )

Definition at line 37 of file PPFXThinNeutronPionWeightCalc.cxx.

37 {}

Member Function Documentation

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

overridevirtual

Implements evwgh::WeightCalc.

Definition at line 39 of file PPFXThinNeutronPionWeightCalc.cxx.

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

   {
     //get configuration for this function
     fhicl::ParameterSet const& pset = p.get<fhicl::ParameterSet>(GetName());
     fGenieModuleLabel = p.get<std::string>("genie_module_label");
 
     //ppfx setup
     fInputLabels = pset.get<std::vector<std::string>>("input_labels");
     fPPFXMode = pset.get<std::string>("ppfx_mode");
     fVerbose = pset.get<int>("verbose");
     fMode = pset.get<std::string>("mode");
     fHorn = pset.get<std::string>("horn_curr");
     fTarget = pset.get<std::string>("target_config");
     fSeed = pset.get<int>("random_seed", -1);
 
     gSystem->Setenv("MODE", fPPFXMode.c_str());
 
     fPPFXrw = NeutrinoFluxReweight::MakeReweight::getInstance();
     std::cout << "PPFX instance " << fPPFXrw << std::endl;
 
     std::string inputOptions;                               // Full path.
     std::string mode_file = "inputs_" + fPPFXMode + ".xml"; // Just file name.
     cet::search_path sp("FW_SEARCH_PATH");
     sp.find_file(mode_file, inputOptions);
 
     if (fSeed != -1) fPPFXrw->setBaseSeed(fSeed); // Set the random seed
     std::cout << "is PPFX setup : " << fPPFXrw->AlreadyInitialized() << std::endl;
     std::cout << "Setting PPFX, inputs: " << inputOptions << std::endl;
     std::cout << "Setting Horn Current Configuration to: " << fHorn << std::endl;
     std::cout << "Setting Target Configuration to: " << fTarget << std::endl;
     if (!(fPPFXrw->AlreadyInitialized())) { fPPFXrw->SetOptions(inputOptions); }
     std::cout << "PPFX just set with mode: " << fPPFXMode << std::endl;
   }

std::string evwgh::WeightCalc::GetName ( )

inlineinherited

std::vector< std::vector< double > > evwgh::PPFXThinNeutronPionWeightCalc::GetWeight ( art::Event & e )

overridevirtual

Implements evwgh::WeightCalc.

Definition at line 74 of file PPFXThinNeutronPionWeightCalc.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.

   {
     std::vector<std::vector<double>> weight;
     evgb::MCTruthAndFriendsItr mcitr(e, fInputLabels);
     //calculate weight(s) here
 
     int nmctruth = 0, nmatched = 0;
     bool flag = true;
     int ievt = -1;
     std::vector<art::Handle<std::vector<bsim::Dk2Nu>>> dk2nus2 =
       e.getMany<std::vector<bsim::Dk2Nu>>();
     for (size_t dk2 = 0; dk2 < dk2nus2.size(); ++dk2) {
       art::Handle<std::vector<bsim::Dk2Nu>> dk2nus = dk2nus2[dk2];
     }
     while ((flag = mcitr.Next())) {
       std::string label = mcitr.GetLabel();
       const simb::MCTruth* pmctruth = mcitr.GetMCTruth();
       // const simb::GTruth*  pgtruth  = mcitr.GetGTruth();
       //not-used//const simb::MCFlux*      pmcflux   = mcitr.GetMCFlux();
       const bsim::Dk2Nu* pdk2nu = mcitr.GetDk2Nu();
       //not-used//const bsim::NuChoice*    pnuchoice = mcitr.GetNuChoice();
       // art::Ptr<simb::MCTruth>  mctruthptr = mcitr.GetMCTruthPtr();
 
       ++ievt;
       ++nmctruth;
       if (fVerbose > 0) {
         std::cout << "FluxWeightCalculator [" << std::setw(4) << ievt << "] "
                   << " label \"" << label << "\" MCTruth@ " << pmctruth << " Dk2Nu@ " << pdk2nu
                   << std::endl;
       }
       if (!pdk2nu) continue;
       ++nmatched;
 
       //RWH//bsim::Dk2Nu*  tmp_dk2nu  = fTmpDK2NUConversorAlg.construct_toy_dk2nu( &mctruth,&mcflux);
       //RWH//bsim::DkMeta* tmp_dkmeta = fTmpDK2NUConversorAlg.construct_toy_dkmeta(&mctruth,&mcflux);
       //RWH// those appear to have been memory leaks
       //RWH// herein begins the replacment for the "construct_toy_dkmeta"
 
       static bsim::DkMeta
         dkmeta_obj; //RWH// create this on stack (destroyed when out-of-scope)  ... or static
       dkmeta_obj.tgtcfg = fTarget;
       dkmeta_obj.horncfg = fHorn;
       (dkmeta_obj.vintnames).push_back("Index_Tar_In_Ancestry");
       (dkmeta_obj.vintnames).push_back("Playlist");
       bsim::DkMeta* tmp_dkmeta = &dkmeta_obj;
       //RWH// herein ends this block
 
       // sigh ....
       //RWH// this is the signature in NeutrinoFluxReweight::MakeReweight :
       //      //! create an interaction chain from the new dk2nu(dkmeta) format
       //      void calculateWeights(bsim::Dk2Nu* nu, bsim::DkMeta* meta);
       //RWH// these _should_ be "const <class>*" because we don't need to change them
       //RWH// and the pointers we get out of the ART record are going to be const.
       bsim::Dk2Nu* tmp_dk2nu = const_cast<bsim::Dk2Nu*>(pdk2nu); // remove const-ness
       try {
         fPPFXrw->calculateWeights(tmp_dk2nu, tmp_dkmeta);
       }
       catch (...) {
         std::cerr << "Failed to calcualte weight" << std::endl;
         continue;
       }
       //Get weights:
       if (fMode == "reweight") {
         double ppfx_cv_wgt = fPPFXrw->GetCVWeight();
         std::vector<double> wvec = {ppfx_cv_wgt};
         weight.push_back(wvec);
       }
       else {
         std::vector<double> vttncpion = fPPFXrw->GetWeights("ThinTargetnCPion");
 
         std::vector<double> tmp_vhptot;
         for (unsigned int iuniv = 0; iuniv < vttncpion.size(); iuniv++) {
           tmp_vhptot.push_back(float(vttncpion[iuniv]));
         }
         weight.push_back(tmp_vhptot);
       }
     }
     return weight;
   }

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

centralValues	the values to be smeared
inputCovarianceMatrix	covariance matrix for smearing
n_multisims	number 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.

   {
 
     std::vector<std::vector<double>> setOfSmearedCentralValues;
 
     //Check that covarianceMatrix is of compatible dimension with the central values
     unsigned int covarianceMatrix_dim = centralValue.size();
 
     if (inputCovarianceMatrix.size() != covarianceMatrix_dim) {
       throw art::Exception(art::errors::Configuration)
         << "inputCovarianceMatrix rows " << inputCovarianceMatrix.size()
         << " not equal to entries of centralValue[]: " << covarianceMatrix_dim;
     }
 
     for (auto const& row : inputCovarianceMatrix) { // Range-for (C++11).
       if (row.size() != covarianceMatrix_dim) {
         throw art::Exception(art::errors::Configuration)
           << "inputCovarianceMatrix columns " << row.size()
           << " not equal to entries of centralValue[]: " << covarianceMatrix_dim;
       }
     }
 
     //change covariance matrix into a TMartixD object
     int dim = int(inputCovarianceMatrix.size());
     TMatrixD covarianceMatrix(dim, dim);
     int i = 0;
     for (auto const& row : inputCovarianceMatrix) {
       int j = 0;
       for (auto const& element : row) {
         covarianceMatrix[i][j] = element;
         ++j;
       }
       ++i;
     }
 
     //perform Choleskey Decomposition
     TDecompChol dc = TDecompChol(covarianceMatrix);
     if (!dc.Decompose()) {
       throw art::Exception(art::errors::StdException)
         << "Cannot decompose covariance matrix to begin smearing.";
       return setOfSmearedCentralValues;
     }
 
     //Get upper triangular matrix. This maintains the relations in the
     //covariance matrix, but simplifies the structure.
     TMatrixD U = dc.GetU();
 
     //for every multisim
     for (int n = 0; n < n_multisims; ++n) {
 
       //get a gaussian random number for every central value element
       int dim = centralValue.size();
 
       std::vector<double> rands(dim);
       GaussRandom.fireArray(dim, rands.data());
 
       //compute the smeared central values
       std::vector<double> smearedCentralValues;
       for (int col = 0; col < dim; ++col) {
         //find the weight of each col of the upper triangular cov. matrix
         double weightFromU = 0.;
         for (int row = 0; row < col + 1; ++row) {
           weightFromU += U(row, col) * rands[row];
         }
         //multiply this weight by each col of the central values to obtain
         //the gaussian smeared and constrained central values
         // smearedCentralValues.push_back(weightFromU * centralValue[col]);
         smearedCentralValues.push_back(weightFromU + centralValue[col]);
       }
 
       //collect the smeared central values into a set
       setOfSmearedCentralValues.push_back(smearedCentralValues);
     } //loop over multisims
     return setOfSmearedCentralValues;
   } // WeightCalc::MultiGaussianSmearing() - gives a vector of sets of smeared parameters

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.

   {
 
     //compute the smeared central values
     std::vector<double> smearedCentralValues;
 
     //
     //perform Choleskey Decomposition
     //
     // Best description I have found
     //     is in the PDG (Monte Carlo techniques, Algorithms, Gaussian distribution)
     //
     //  http://pdg.lbl.gov/2016/reviews/rpp2016-rev-monte-carlo-techniques.pdf (Page 5)
     //
     //
     TDecompChol dc = TDecompChol(*(inputCovarianceMatrix));
     if (!dc.Decompose()) {
       throw art::Exception(art::errors::StdException)
         << "Cannot decompose covariance matrix to begin smearing.";
       return smearedCentralValues;
     }
 
     //Get upper triangular matrix. This maintains the relations in the
     //covariance matrix, but simplifies the structure.
     TMatrixD U = dc.GetU();
 
     for (unsigned int col = 0; col < centralValue.size(); ++col) {
       //find the weight of each col of the upper triangular cov. matrix
       double weightFromU = 0.;
 
       for (unsigned int row = 0; row < col + 1; ++row) {
         weightFromU += U(row, col) * rand[row];
       }
 
       //multiply this weight by each col of the central values to obtain
       //the gaussian smeared and constrained central values
       // smearedCentralValues.push_back(weightFromU * centralValue[col]);
       smearedCentralValues.push_back(weightFromU + centralValue[col]);
     }
     return smearedCentralValues;
   } // WeightCalc::MultiGaussianSmearing() - Gives a single set of smeared parameters

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.

   {
 
     //compute the smeared central values
     std::vector<double> smearedCentralValues;
 
     // This guards against accidentally
     if (!isDecomposed) {
       throw art::Exception(art::errors::StdException)
         << "Must supply the decomposed lower triangular covariance matrix.";
       return smearedCentralValues;
     }
 
     for (unsigned int col = 0; col < centralValue.size(); ++col) {
       //find the weight of each col of the upper triangular cov. matrix
       double weightFromU = 0.;
 
       for (unsigned int row = 0; row < col + 1; ++row) {
         weightFromU += LowerTriangleCovarianceMatrix[0][row][col] * rand[row];
       }
 
       //multiply this weight by each col of the central values to obtain
       //the gaussian smeared and constrained central values
       // smearedCentralValues.push_back(weightFromU * centralValue[col]);
       smearedCentralValues.push_back(weightFromU + centralValue[col]);
     }
     return smearedCentralValues;
   } // WeightCalc::MultiGaussianSmearing() - Gives a single set of smeared parameters