KVEvaporationSpectrum.cpp

#include "KVEvaporationSpectrum.h"
#include "TH1.h"
#include "TGraph.h"
 
ClassImp(KVEvaporationSpectrum)
 
 
 
KVEvaporationSpectrum::KVEvaporationSpectrum(const TString& name, int ncomponents, bool constrain_sigma)
   : spe_calc{ncomponents,constrain_sigma}, TF1(name, &spe_calc, 0., 500., ncomponents*(constrain_sigma ? 3 : 4))
{
   // \param[in] constrain_sigma if false, let the width of the gaussian be a free parameter
 
   if(ncomponents==1)
   {
      SetParNames("N", "B", "T");
      SetParLimits(0,1.,1.e+09);
      SetParLimits(1,0.,100.);
      SetParLimits(2,1.e-2,100.);
      if(!constrain_sigma)
      {
         SetParName(3, "#sigma");
         SetParLimits(3,.1,20.);
      }
   }
   else
   {
      int npar = constrain_sigma ? 3 : 4;
      for(int nc=0; nc<ncomponents; ++nc)
      {
         SetParName(nc*npar,Form("N_%d",nc+1));
         SetParName(nc*npar+1,Form("B_%d",nc+1));
         SetParName(nc*npar+2,Form("T_%d",nc+1));
         SetParLimits(nc*npar,1.,1.e+09);
         SetParLimits(nc*npar+1,0.,100.);
         SetParLimits(nc*npar+2,1.e-2,100.);
         if(!constrain_sigma)
         {
            SetParName(nc*npar+3, Form("#sigma_%d",nc+1));
            SetParLimits(nc*npar+3,.1,20.);
         }
      }
   }
   SetNpx(500);
   SetLineColor(kBlue+2);
}
 
 
 
 
void KVEvaporationSpectrum::SetInitialParameters(const TH1* spec)
{
   // Deduce initial parameter values from contents of histogram
   //
   //  - assuming standard evaporation spectrum
   //  - position of maximum => B+T
   //  - mean value of spectrum => B+2T
   //
   auto npar=get_npar();
   int ncomp = GetNComponents();
   SetParameter(0, spec->GetEntries()/ncomp);
   auto sp_max = spec->GetBinCenter(spec->GetMaximumBin());
   auto sp_mean = spec->GetMean();
   auto B = 2*sp_max - sp_mean;
   auto T = sp_mean - sp_max;
   SetParameter(1,B);
   SetParameter(2,T);
   if(npar>3)
      SetParameter(3, TMath::E()/TMath::Sqrt(TMath::TwoPi())*T);
   if(ncomp>1)
   {
      for(int nc=1; nc<ncomp; ++nc)
      {
         SetParameter(nc*npar, spec->GetEntries()/ncomp);
         SetParameter(nc*npar+1, B+nc*2);
         SetParameter(nc*npar+2, T+nc*5);
         if(npar>3)
            SetParameter(nc*npar+3, TMath::E()/TMath::Sqrt(TMath::TwoPi())*T);
      }
   }
}
 
 
 
 
TMultiGraph* KVEvaporationSpectrum::GetComponents(const std::vector<Color_t>& line_colors) const
{
   static std::vector<Style_t> line_styles = {2, 7, 9, 3, 4, 5, 6, 8, 10, 11};
   auto mg = new TMultiGraph;
 
   double tot_stat = 0;
   for(int nc=0; nc<GetNComponents(); ++nc)
      tot_stat += GetParameter(nc*get_npar());
 
   for(int nc=0; nc<GetNComponents(); ++nc)
   {
      auto g = new TGraph;
      g->SetLineWidth(2);
      g->SetLineStyle(line_styles[nc]);
      g->SetLineColor(line_colors[nc]);
 
      g->SetTitle(Form("B=%.1f MeV T=%.1f MeV (%.0f%%)", GetParameter(nc*get_npar()+1), GetParameter(nc*get_npar()+2), 100*GetParameter(nc*get_npar())/tot_stat));
 
      double xmin,xmax;
      GetRange(xmin,xmax);
      int npoints = GetNpx();
      auto dx = (xmax-xmin)/(npoints-1.);
      for(int i=0; i<npoints; ++i)
      {
         double x = xmin+i*dx;
         g->SetPoint(i, x, spe_calc.SC(&x,GetParameters()+nc*get_npar()));
      }
      mg->Add(g,"l");
   }
 
   return mg;
}
 
 
 
 
double KVEvaporationSpectrum::spectrum_calculator::operator()(double*x, double*par) const
{
   auto N=par[0];
   auto B=par[1];
   auto T=par[2];
   auto E=x[0];
   if(E>=(B+T)) return N*((E-B)/pow(T,2.))*TMath::Exp(-(E-B)/T);
   double sigma;
   if(con_sig)
   {
      sigma = norm_const*T;
      return N*TMath::Gaus(E,B+T,sigma,kTRUE);
   }
   sigma=par[3];
   return (norm_const*sigma/T)*N*TMath::Gaus(E,B+T,sigma,kTRUE);
}
 
 
 
 
double KVEvaporationSpectrum::multi_comp_spectrum_calculator::operator()(double*x, double*par) const
{
   // return sum of Nc components
   // par[0],[1],[2]: N,B,T of 1st component
   // par[3],[4],[5]: N,B,T of 2nd component
   // etc.
 
   double sum = 0;
   for(int i=0; i<Nc; ++i)
   {
      // temperatures of components must increase
      if(i && (par[i*3+2]<par[(i-1)*3+2]))
         return 0;
 
      sum += SC(x, &par[i*3]);
   }
   return sum;
}