GammaProductionAna.h

Go to the documentation of this file.

#ifndef _GAMMA_PRODUCTION_ANALYSIS_H_
#define _GAMMA_PRODUCTION_ANALYSIS_H_
//ROOT includes
#include "TH1.h"
#include "SpectrumAnalysis.h"
#include "TH2.h"
#include "TMatrixD.h"
#include "TFitResult.h"
#include "TRandom3.h"
#include "TALYSUtils.h"
#include "FluxMatrixAna.h"
#include "NSDPhysicsEvents.h"
#include "TChain.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 DataValuePair
{
public:
    DataValuePair();
    double m_xVal;
    double m_xUnc;
    double m_xUncUp=0.;
    double m_yVal;
    double m_yUnc;
  friend std::ostream& operator<<(std::ostream& os, const DataValuePair& a_dataPair);
  friend std::istream& operator>>(std::istream& a_istream, DataValuePair& a_dataPair);
};
 
class GammaProductionAna
{
public:
    GammaProductionAna();
    GammaProductionAna(int a_ZID, int a_AID, string a_talysPath);
// ------------------
    void setupTALYSUtils(int a_ZID, int a_AID, string a_talysPath);
    void setTalysPath(std::string a_path);
    TH2F* buildFluxMatrix(std::string a_stofFluxFileName,
                                                double a_minEnergy,
                                                double a_maxEnergy,
                                                double a_flightPath,    // cm
                                                double a_timeResolution, // ns [sigma]
                                                double a_RFPeriod,
                                                double a_timeBinWidth=-1);
    void addFluxMatrix(FluxAna a_fluxAna);
    void readGammaEffParams(std::string a_effParamFile, 
                                                    std::string a_fitType="Debertin",
                                                    bool a_readEffCovMat = false);
    void readEGammaLevels(std::string a_lvlFile);
    void loadExpTOFvsEGamma(vector<std::string> a_expFiles);
    void loadExpTOFvsEGamma(std::string a_expFile);
  // loads tree without making histograms allowing you to load several
  // files as a seperate step
  void addTree(std::string a_expTree);
  void addTreeBKG(std::string a_expTree);
  // builds hists from the loaded trees
  void buildHists(vector<int> a_cloverList,
                              bool a_saveFile,
                              string a_outFileName);
  void buildHistsBKG(vector<int> a_cloverList,
                              bool a_saveFile,
                              string a_outFileName);
    void loadExpBKGData(vector<std::string> a_expFile);
    void loadExpBKGData(std::string a_expFile);
    void setBackgroundScalar(double a_scalar);
    void addEXFORCrossSection(std::string a_fileName,
                                                        double a_gammaEnergy,
                                                        bool a_draw=false,
                                                        std::string a_dataSetName = "");
// ------------------
    void setVerbosity(int a_verbosity);
    void setTotalCharge(double a_charge);
    void setTargetRhoR(double a_rhoR);
    void setTargetSolidAngle(double a_solidAngle);
  void setDetSolidAngle(double a_solidAngle);
    void setRFPeriod(double a_RFPeriod);
    void setTimeResolution(double a_timeResolution);
    void setFlightPath(double a_flightPath);
// ------------------
    vector<double> calcGammaEff(int a_detID, int a_elementID, double a_energy);
 
    vector<Eigen::MatrixXf> getFluxMatrix(double a_minEnergy);
    vector<Eigen::MatrixXf> getFluxMatrix(double a_gammaEnergy,
                                                                                double a_thresholdAdder);
    TALYSUtils* getTALYSUtils();
    TH2F* getFluxHist();
// ------------------   
    PeakInfo getPeakArea(double a_tof, double a_tofWidth,
                                                            int a_detID, int a_elementID,
                                                            double a_energy,
                                                            bool a_draw=0,
                                                            bool a_fitBkg=1,
                                                            bool a_drawBkg=0,
                                                            bool a_saveFit=1);
    void setNewPeakFit(double a_gammaEnergy,
                                        int a_detID,
                                        int a_elementID, 
                                        int a_tof, 
                                        PeakInfo a_peakInfo);
    void writePeakFits(double a_gammaEnergy,
                                            int a_detID,
                                            int a_elementID);
    void readPeakFits(double a_gammaEnergy,
                                            int a_detID,
                                            int a_elementID);
    map<double, std::pair<double,double>> getTotExpYield(double a_gammaEnergy,
                                                                                                    int a_detID, int a_elementID,
                                                                                                    double a_tofWidth=-1,
                                                                                                    bool a_draw=true,
                                                                                                    int a_plotType=0);
    
    Eigen::VectorXf calculateYieldUnc(Eigen::MatrixXf a_fluxUncMat,
                                                                        Eigen::VectorXf a_crossSection);
    Eigen::VectorXf calculateYieldUnc(vector<Eigen::MatrixXf> a_fluxMats,
                                                                        Eigen::VectorXf a_crossSectionUnc,
                                                                        Eigen::VectorXf a_crossSection);
    Eigen::VectorXf convolveXSec(double a_gammaEnergy, 
                                                int a_detID, 
                                                int a_elementID,
                                                bool a_calcUnc=false,
                                                bool a_useEXFOR=false,
                                                double a_thresholdAdder=0.);
// ------------------
    double plotComp(double a_gammaEnergy,
                                int a_detID,
                                int a_elementID,
                                bool a_useEXFOR=false,
                                double a_normalize=1,
                                bool a_drawEnAxis=true,
                                bool a_draw=true);
    double plotRatio(double a_energy1, 
                                    double a_energy2, 
                                    int a_detID, 
                                    int a_elementID,
                                    bool a_useEXFOR=false,
                                    bool a_drawEnAxis=true,
                                    bool a_draw=true);
    void compDetectorElements(double a_gammaEnergy,
                                                        int a_detID1,
                                                        int a_elementID1,
                                                        int a_detID2,
                                                        int a_elementID2);
    vector<DataValuePair> calcUnwrappedXSec(double a_gammaEnergy,
                                                int a_detID,
                                                int a_elementID,
                                                double a_thresholdAdder=0.,
                                                double a_normalize=0,
                                                bool a_draw=true);
    vector<DataValuePair> calcUnwrappedXSec(map<double, std::pair<double,double>> a_expYields,
                                                double a_gammaEnergy,
                                                double a_threshold,
                                                bool a_normalize=0,
                                                bool a_draw=true);
    //overloaded method to input gamma data from radware
  vector<DataValuePair> calcUnwrappedXSec(string a_expYields,
                        double a_gammaEnergy,
                        double a_threshold,
                        bool a_normalize=0,
                        bool a_draw=true);
 
    void testUnwrapping(double a_gammaEnergy,
                                            int a_detID,
                                            int a_elementID);
    void noExpUnwrapTest(double a_gammaEnergy,
                                                int a_detID=0,
                                                int a_elementID=2,
                                                double a_alpha=1.);
    void plotSVD();
    double testForwardModel(double a_gammaEnergy,
                                                int a_detID,
                                                int a_elementID,
                                                double a_smoothness=1.e-2,
                                                bool a_draw=true,
                                                bool a_useEXFOR=false);
    double forwardModelFit(const double* a_params);
    void setupTALYSParams(std::string a_paramFile);
    void talysForwardModel(vector<double> a_gammaEnergy,
                                                int a_detID,
                                                int a_elementID,
                                                bool a_scale=false,
                                                bool a_runECIS=true);
    double talysModelFit(const double* a_params);
    vector<std::pair<double,double>> scanTALYSParameter(int a_parameterNum,
                                                    vector<double> a_gammaEnergy,
                                                    bool a_runECIS=true,
                                                    int a_detID=1,
                                                    int a_elementID=2);
    vector<std::pair<double,double>> scanTALYSParameter(std::string a_paramName,
                                                    vector<double> a_gammaEnergy,
                                                    bool a_runECIS=true,
                                                    int a_detID=1,
                                                    int a_elementID=2); 
    void setSaveTALYSFits(bool a_saveFits);
    vector<double> getEnergySpecLowerBinEdges();
    Eigen::VectorXf getEXFORCrossSection(double a_gammaEnergy,
                                                                            bool a_calcUnc=false,
                                                                            double a_thresholdAdder=0);
    SpectrumAnalysis m_specAna;
    TH1F* m_lastHist;
 
    std::map<string, int> m_talysParams;
    map<int,map<int,TH2F*>> m_tofEgamma;
 
    void writeNonWrappedXSec(double a_gammaE, int a_detID, int a_elemID);
    map<int, PeakInfo> getPeakFits(double a_gammaE, int a_detID, int a_elemID);
private:
    TH2F* m_fluxHist;
    Eigen::MatrixXf m_fullFluxMat;
    Eigen::MatrixXf m_fullFluxMatErr;
    double m_bkgScalar;
    vector<double> m_ener0BinEdges;
    vector<double> m_enerBinWidths;
    vector<double> m_enerBinCenters;
    // map<int,map<int,TH2F*>> m_tofEgamma;
    map<int,map<int,TH2F*>> m_tofEgamma0;
    map<int,map<int,TH2F*>> m_tofEgammaBKG;
    map<int,map<int,vector<double>>> m_effParams;
    map<double,vector<std::pair<std::string,std::string>>> m_eGammaLevels;
    map<double, vector<double>> m_eGammaThreshold;
    std::string m_talysPath;
    double m_rhoR = 5.283e21;
    double m_solidAngle = 3.67e-5;
    double m_solidAngleDet = 1;
    double m_totalCharge;
    double m_RFPeriod;
    double m_timeResolution;
    double m_flightPath;
    map<double, map<int, map<int, map<int, PeakInfo>>>> m_peakFits;
    map<int, map<int, Eigen::MatrixXf>> m_gammaEffCovMats;
    map<double, std::pair<vector<double>, vector<double>>> m_exforData;
    map<double, std::string> m_exforDataSetNames;
    map<string, map<double, vector<DataValuePair>>> m_exforDataAll;
    bool m_hasGammaEffUnc;
    bool m_hasFluxUncMat;
    bool m_useTALYSscaler;
    bool m_runECIS;
 
    TALYSUtils m_talysUtil; 
    // vector<double> getEnergyPoints(double a_gammaEnergy);
 
    map<int, map<int,map<double, vector<DataValuePair>>>> m_unwrappedXSecs;
    map<double, vector<DataValuePair>> m_currUnwrappedData;
    map<double, vector<int>> m_currUnwrappedDataIndex;
 
    TRandom3* m_rand_gen;
 
    FluxAna m_fluxMatAna;
public:
    vector<double> calcNonWrappedEners(vector<double> a_cumulativeCounts,                                                                       vector<double> a_eners);
private:    
int m_verbosity;
public:
vector<double> getEnergyPoints(double a_gammaEnergy, int a_coinGammaID=0);
    bool m_saveFits;
    bool m_ratio;
    int setupTALYSMinimizer(vector<double> a_gammaEnergy,
                                                    int a_detID,
                                                    int a_elementID);
    vector<Eigen::MatrixXf> m_currFluxMat;
    vector<Eigen::MatrixXf> m_currFluxUncMat;
    vector<Eigen::VectorXf> m_currData;
    vector<Eigen::VectorXf> m_currDataUnc;
    std::pair<int ,int> m_currDetElemID;
    // std::map<string, std::pair<double,std::pair<double,double>>> m_talysParams;
    vector<double> m_tempGammaEnergy;
    ROOT::Math::Minimizer* m_minimum;
 
    double m_smoothness;
    // Double_t fTOFWrapper(Double_t* x, Double_t* par);
  //adding calibrated tree support
  //storing loaded trees
  TChain* m_trees;
  std::vector<string> m_loadedTrees;
  bool m_hasFile;
  NSDPhysicsEvents::CloverEvent* m_currentEvent;
 
  TChain* m_treesBKG;
  std::vector<string> m_loadedTreesBKG;
  bool m_hasFileBKG;
  NSDPhysicsEvents::CloverEvent* m_currentEventBKG;
};
 
#endif