GenesisNeutronSinglesAna.h

Go to the documentation of this file.

#ifndef _GEN_NEUTRON_SINGLES_ANA_H_
#define _GEN_NEUTRON_SINGLES_ANA_H_
//ROOT includes
#include "TH1.h"
#include "GammaProductionAna.h"
#include "FluxMatrixAna.h"
#include "TALYSUtils.h"
#include "TH2.h"
#include "TMatrixD.h"
#include "TFitResult.h"
#include "TRandom3.h"
#include "TALYSUtils.h"
#include "TTree.h"
//c++includes
#include <cstdint>
#include <vector>
#include <utility>
#include <deque>
using std::vector;
#include <map>
using std::map;
using std::pair;
 
// #include "../../BasicSupport/external/eigen/Eigen/Dense"
#include <Eigen/Dense>
 
class GenesisNeutronSinglesAna
{
public:
  GenesisNeutronSinglesAna();
  GenesisNeutronSinglesAna(FluxAna a_fluxMat,
                            TALYSUtils a_talysUtils);
  void addFluxMatrix(FluxAna a_fluxMat);
  void addTALYSUtils(TALYSUtils a_talysUtils);
  TALYSUtils* getTALYSUtils();
  void readTotalNeutronTOFHists(string a_fileName,
                                string a_fileNameBKG,
                                double a_backgroundScaler);
  void readTotalNeutronTOFTree(string a_fileName);
  void readTotalNeutronTOFTreeBKG(string a_fileNameBKG);
  void subtractBackground(double a_bkgScalar);
  void readNeutronResponseFunctions(string a_fileNameTarget,
                                    double a_numSim = 2.e9);
  void readAngularNeutronResponseFunctions(string a_fileNameTarget,
                                          double a_numSim = 2.e9);
  void useAngularResponseFunctions(bool a_useAngularRespFuncs);
  void setRespFuncNumEnerBins(int a_numEnerBins);
  void setRespFuncNumAngleBins(int a_numAngleBins);
  void setRespFuncNumLYBins(int a_numLYBins);
  void setRespFuncMaxEner(double a_maxEnergy);
  void setTotCharge(double a_totCharge);
  void setTargetRhoR(double a_rhoR);
  void setTargetSolidAngle(double a_solidAngle);
  void readDetDistAng(string a_configFileName);
  double getDetectorDistance(int a_detID);
  double getDetectorAngle(int a_detID);
 
  void readENDFEventTree(string a_fileName);
 
  void readYAHFCEventTree(string a_fileName);
  // should return vector (with STOF errors)
  vector<Eigen::MatrixXf> plotTALYSTotalTOFvsEnOut(int a_detID, bool a_draw=1);
  vector<Eigen::MatrixXf> plotENDFTotalTOFvsEnOut(int a_detID, bool a_draw=1);
  // vector<Eigen::MatrixXf> plotYAHFCTotalTOFvsEnOut(int a_detID, bool a_draw=1);
  vector<Eigen::MatrixXf> convolveTotalTOFResponse(int a_detID);
  vector<Eigen::MatrixXf> convolveTotalTOFResponse(double a_angle,
                                                    double a_width);
  TH2F* getModelResponseHist(int a_detID);
  TH2F* getModelResponseHist(double a_angle, double a_width);
  double calcSliceChi2(int a_detID, double a_tof);
  double calcSliceChi2(double a_angle, double a_width, double a_tof);
  double getFullChi2(int a_detID);
  double getFullChi2(double a_angle, double a_width);
  void setForwardModelDetectors(vector<int> a_dets);
  void setForwardModelAngles(vector<std::pair<double,double>> a_angles);
  void forwardModel();
  double talysModelFit(const double* a_params);
 
  int readTALYSOutput();
  void plotTotalComp(int a_detID);
  void plotTotalComp(double a_angle, double a_width);
  TH2F* getExpHists(int a_detID);
  TH2F* getExpHists(double a_angle, double a_width);
  void plot1DComp(int a_detID, double a_minLight);
  void plotSliceComp(int a_detID, double a_tof);
  void plotSliceComp(double a_angle, 
                      double a_width,
                      double a_tof);
  TH2F* getResponseHist(int a_detID);
  TH2F* getResponseHist(double a_angle,
                        double a_width);
  vector<double> getCovarianceMatrix(std::ostream& stream);
private:
  map<int, TH2F*> m_hists_bkg;
  map<int, TH2F*> m_hists_orig;
  map<int, TH2F*> m_hists_blank;
  map<int, TH2F*> m_respFuncs;
  map<int, Eigen::MatrixXf> m_respFuncsMat;
  map<int, Eigen::MatrixXf> m_respFuncsMat2;
  map<int, Eigen::MatrixXf> m_respFuncsMatUnc;
  map<int, vector<double>> m_detDistAng;
  double m_bkgScalar;
  double m_totalCharge;
  double m_targetRhoR;
  double m_targetSolidAngle;
  // vector<Eigen::MatrixXf> m_currModelMats;
 
  FluxAna m_fluxMat;
  TALYSUtils m_talysUtils;
  
  vector<double> m_flux1D;
  vector<double> m_flux1D_unc;
  vector<double> m_flux1D_binEdges;
 
  TTree* m_endfEvents;
  TTree* m_yahfcEvents;
  bool m_useAngularRespFuncs;
 
  TH2F* m_totalExpHists;
 
  int m_numLYBins;
  int m_numTTOFBins;
  int m_numEnerBins;
  int m_numAngleBins;
  double m_maxEnergy;
 
  vector<int> m_forwardModelDetIDs;
  vector<std::pair<double,double>> m_forwardModelDetAngles;
public:
  ROOT::Math::Minimizer* m_minimum;
  vector<Eigen::MatrixXf> m_currModelMats;
  std::map<string, int> m_talysParams;
 
  int m_tofWidth = 0;
};
 
#endif