SpectrumAnalysis.h

Go to the documentation of this file.

#ifndef _SPECTRUM_ANALYSIS_H_
#define _SPECTRUM_ANALYSIS_H_
#include "ConfigClasses.h"
#include "TH1.h"
#include "TFile.h"
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h" 
#include <vector>
#include <map>
#include <string>
#include <iostream>
using std::map;
using std::vector;
using std::string;
 
class PeakInfo
{
public:
  PeakInfo();
  PeakInfo(TFitResultPtr& a_fitResult, double a_binWidth);
  vector<double> m_amp;
  vector<double> m_sig;
  vector<double> m_bkg;
  vector<double> m_mean;
  double m_binWidth;
  double m_chi2;
  void printPeak();
  vector<double> getArea();
  friend std::ostream& operator<<(std::ostream& os, const PeakInfo& a_peakInfo);
 
  friend std::istream& operator>>(std::istream& a_istream, PeakInfo& a_peakInfo);
 
};
 
class SpectrumAnalysis
{
public:
  using EffPeakSet = std::map<double, std::pair<double, double>>;
  
  SpectrumAnalysis();
 
  void loadHists(string a_fname, int a_modID, double a_countTime=0);
  // void loadSimHists(string a_fname);
  void loadBackgroundHists(string a_fname, int a_modID, double a_countTime);
  void addHists(string a_fname, int a_modID, double a_countTime=0);
  void addBackgroundHists(string a_fname, int a_modID, double a_countTime);
  void subtractBackgroundSpectra();
  void setResolutionParameters(vector<double> a_params);  
  void fixResolution(bool a_fixRes);
  void fixMean(bool a_fixMean);
  void setVerbosity(int a_verbosity);
  // void rebinSpectrum(int a_channel, 
  //                     int a_rebinFactor=2,
  //                     bool a_fullClover = false);
  void rebinSpectrum(DetType a_detType,
                      int a_detID,
                      int a_elemID, 
                      int a_rebinFactor=2);
  // vector<PeakInfo> FitPeaksMan(int a_channel, 
  //                         double a_lower, 
  //                         double a_higher,// int a_numPeaks,
  //                         vector<double> a_peakMeans,
  //                         bool a_fullClover=false, 
  //                         bool a_addBack=false, 
  //                         bool a_checkPeakRes=true);
  vector<PeakInfo> FitPeaksMan(DetType a_detType,
                          int a_detID,
                          int a_elemID, 
                          double a_lower, 
                          double a_higher,
                          vector<double> a_peakMeans,
                          bool a_addBack=false,
                          bool a_bgoVeto = false, 
                          bool a_checkPeakRes=true);
  vector<PeakInfo> FitPeaksMan(TH1F* a_hist,
                            double a_lower,
                            double a_higher,
                            vector<double> a_peakMeans,
                            int a_polyDegree=0);
  // vector<PeakInfo> FitPeaksManPol1(int a_channel, 
  //                             int a_lower, 
  //                             int a_higher, 
  //                             // int a_numPeaks, 
  //                             vector<double> a_peakMeans,
  //                             bool a_checkPeakRes=true,
  //                             bool a_fullClover=false, 
  //                             bool a_addBack=false);
  vector<PeakInfo> FitPeaksManPol1(DetType a_detType,
                          int a_detID,
                          int a_elemID, 
                          double a_lower, 
                          double a_higher,
                          vector<double> a_peakMeans,
                          bool a_addBack=false,
                          bool a_bgoVeto = false, 
                          bool a_checkPeakRes=true);
  void printPeakFits();
  // vector<PeakInfo> FitPeaksErf(int a_channel, 
  //                         double a_lower, 
  //                         double a_higher, 
  //                         double a_peakMean,
  //                         bool a_fullClover=false, 
  //                         bool a_addBack=false);
  vector<PeakInfo> FitPeaksErf(DetType a_detType,
                          int a_detID,
                          int a_elemID, 
                          double a_lower, 
                          double a_higher, 
                          double a_peakMean,
                          bool a_addBack=false,
                          bool a_bgoVeto = false);
  vector<PeakInfo> FitPeaksErf(TH1F* a_hist, 
                          double a_lower, 
                          double a_higher, 
                          double a_peakMean);
  // Specific function (but didn't want to muddy FitPeaksMan more) to fit the
  // 6111 keV peak for Cl-35 efficiency for clover leafs that are double-peaking
  vector<PeakInfo> FitPeaks6111Double(TH1F* a_hist);
  // EffPeakSet Eu152Efficiency(int a_channel, 
  //                            double a_activity, 
  //                            double a_countTime,
  //                            double a_simBrTot=.2101,
  //                            bool a_res=false,
  //                            bool a_fullClover=false,
  //                            bool a_addBack=false,
  //                            bool a_draw =false);
  EffPeakSet Eu152Efficiency(DetType a_detType,
                              int a_detID,
                              int a_elemID, 
                             double a_activity, 
                             double a_countTime,
                             double a_simBrTot=.2101,
                             bool a_res=false,
                             bool a_bgoVeto=false,
                             bool a_addBack=false,
                             bool a_draw =false);
  EffPeakSet Eu152Efficiency(TH1F* a_hist, 
                             double a_activity, 
                             double a_countTime,
                             double a_simBrTot=.2101,
                             bool a_res=false,
                             bool a_draw =false);
 
 
  EffPeakSet Co56Efficiency(TH1F* a_hist, 
                             double a_activity, 
                             double a_countTime,
                             double a_simBrTot=.9993999,
                             bool a_res=false,
                             bool a_draw =false);
 
  // multi_peaking int is for using a double- or triple-Gaussian fit (FitPeaksMan with two
  // or three peaks in the vector argument) for when a clover is double-peaking (or triple...)
    EffPeakSet Ga66Efficiency(TH1F* a_hist, 
                             double a_activity, 
                             double a_countTime,
                             int multi_peaking=0,
                             double a_simBrTot=.37,
                             bool a_res=false,
                             bool a_draw =false);
 
  // double peak Boolean is for using a double-Gaussian fit (FitPeaksMan with two
  // peaks in the vector argument) for when a clover is double-peaking
    EffPeakSet Cl35Efficiency(TH1F* a_hist, 
                             double a_activity, 
                             double a_countTime,
                             bool double_peak=false,
                             double a_simBrTot=1.0,
                             bool a_res=false,
                             bool a_draw =false);
 
  // Creates an EffPeakSet object from a text file that can then be
  // passed to other functions.
  // The purpose is so that you can combine multiple isotope's efficiencies
  // into one EffPeakSet object for doing things like calcEffCurve. 
  // The text file must be in a three-column, space-separated format of 
  // energy | efficiency | efficiency uncertainty
  // Your multiple-isotope efficiencies should be scaled before you add them
  // to this file, if applicable.
  // Order doesn't matter, EffPeakSet sorts them in order by energy already.
  EffPeakSet createEffPeakSet(string infilename);
 
  static Double_t totalFit(Double_t *x, Double_t *par);
 
  static Double_t fullGeFit(Double_t *x, Double_t *par);
 
  vector<double> calcEffCurve(EffPeakSet const& a_ps, string a_fitType) const;
 
  void calcAddBackFactor();
  // EffPeakSet U235Efficiency(int a_channel, 
  //                           double a_activity,
  //                           double a_activityErr, 
  //                           double a_countTime);
 
  map<int, TH1*> getHists();
  map<int,vector<TH1*>> getCloverHists();
  TH1* getHist(DetType a_detType,
                int a_detID,
                int a_elemID=0,
                bool a_addBack=false,
                bool a_bgoVeto=false);
  map<int, map<double, PeakInfo>> getPeakFits();
 
  vector<double> ChannelResolution(DetType a_detType,
                                  int a_detID,
                                  int a_elemID=0);
 
  int autoEu152EnergyCalibration(int a_channel, int a_guess=5000);
  int autoEu152EnergyCalibration(TH1* a_hist, int a_guess=5000);
 
  int autoEu152EnergyCalibration1408Guess(int a_channel, int a_guess);
  int autoEu152EnergyCalibration1408Guess(TH1* a_hist, int a_guess);
  void removePeakFit(DetType a_detType,
                    int a_detID,
                    int a_elemID,
                    double a_energy);
 
  void setUpMinimizer();
 
  double fitPeakConstrained(const double* a_params);
 
  std::vector<double> fitSTOFBandProjection(TH1F* a_histogram,
                                            double a_minLight,
                                            double a_maxLight,
                                            vector<double> a_peakCentroid,
                                            int a_polOrder
                                           );
 
/* 
  void findPeaks(int a_channel, double a_minHeight, int a_numPeaks);
 
  void FitPeakBasic(int a_channel, double a_peakCentroid, double a_low= 5 , double a_high=5);;
  void FitDoublePeak(int a_channel, double a_peakCentroid, double a_peakSeparation, 
                                      double a_relativeHeight, double a_fitWindow);
 
  map<string, vector<vector<double>>> m_zeroCrossings;
  void writeToConfig();
*/
 
map<DetType, map<int, map<int,TH1*>>> m_histsD;
private:
  double m_countTime=0;
  double m_bkgCountTime=0;
 
  map<int, TH1*> m_hists;
  map<int, TH1*> m_bkgHists;
  map<int, vector<TH1*>> m_cloverHists;
  map<int, vector<TH1*>> m_bkgCloverHists;
  // map<DetType, map<int, map<int,TH1*>>> m_histsD;
  map<DetType, map<int, map<int,TH1*>>> m_histsBkgD;
  map<int, map<double, PeakInfo>> m_peakFits;
  map<DetType, map<int, map<int, map<double, PeakInfo>>>> m_peakFitsD;
 
  TH1F* m_currHist;
  TF1* m_currFit;
  double m_currSig;
  int m_numParams;
  vector<int> m_sigParamIndex;
  ROOT::Math::Minimizer* m_minimum;
  ROOT::Math::Functor m_chiSquareFunc;
  vector<double> m_resParams;
  bool m_fixRes;
 
  bool m_fixMean;
  vector<int> m_meanParamIndex;
  int m_verbosity;
};
#endif