diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 055c48c4dc2b152f836d2c3f044fe90f3c77100b..9c6462d00afd2a229c1654249d1d18299e8a1a28 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -1,3 +1,4 @@
+#include <math.h>
#include "NPReaction.h"
#include <string>
#include <sstream>
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
old mode 100755
new mode 100644
index dddab6f902c392cc633b774646858350a6e1ac06..a249438cf6a4de806a4f1e0d5823b2b8f6a3c81e
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -117,7 +117,7 @@ void FitKnownPeaks(TH1F* hist){
full->SetParameter(numParams-1,1.0);
full->SetParLimits(numParams-1,0.0,1e1);
- //Fit full function to histogram
+ //Fit full function tohistogram
hist->Fit(full, "WWR", "", minFit, maxFit);
hist->Draw();
@@ -160,6 +160,52 @@ void FitKnownPeaks(TH1F* hist){
}
+
+
+
+
+
+Double_t f_bg(Double_t *x, Double_t *par){
+ // Flat bg [0] + semicircle [1]*sqrt(6.183^2 - (x-10.829)^2)
+ Float_t xx = x[0];
+ Double_t f;
+ Double_t a = TMath::Power(6.183,2);
+ Double_t b = TMath::Power(xx-10.829,2);
+ if(a > b){ f = par[0] + (par[1]*TMath::Sqrt(a-b)); }
+ else{ f = par[0]; }
+ return f;
+}
+
+Double_t f_full(Double_t *x, Double_t *par){
+ Float_t xx = x[0];
+ Double_t f;
+ Double_t peaks = (par[2]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[1])/par[0],2))
+ + (par[4]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[3])/par[0],2))
+ + (par[6]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[5])/par[0],2))
+ + (par[8]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[7])/par[0],2))
+ + (par[10]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[9])/par[0],2))
+ + (par[12]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[11])/par[0],2))
+ + (par[14]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[13])/par[0],2))
+ + (par[16]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[15])/par[0],2))
+ + (par[18]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[17])/par[0],2))
+ + (par[20]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[19])/par[0],2))
+ + (par[22]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[21])/par[0],2))
+ + (par[24]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[23])/par[0],2))
+ + (par[26]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[25])/par[0],2))
+ + (par[28]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[27])/par[0],2));
+
+ Double_t bg;
+ Double_t a = TMath::Power(6.183,2);
+ Double_t b = TMath::Power(xx-10.829,2);
+ if(a > b){ bg = par[29] + (par[30]*TMath::Sqrt(a-b)); }
+ else{ bg = par[29]; }
+
+ f = peaks + bg;
+ return f;
+}
+
+
+
vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
double minFit=-1.0, maxFit=5.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
@@ -223,6 +269,12 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
bg->SetLineStyle(9);
bg->SetParNames("Background");
+ //Build IMPROVED background function
+ //TF1 *bg = new TF1 ("bg",f_bg, minFit, maxFit);
+ //bg->SetLineColor(kGreen);
+ //bg->SetLineStyle(9);
+ //bg->SetParNames("Background","BreakupScale");
+
//Build total function
TF1 *full = new TF1("fitAllPeaks",
" ([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2)) "
@@ -242,6 +294,11 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
"+ [29]" , minFit, maxFit);
full->SetLineColor(kRed);
+ //Build IMPROVED total function
+ //TF1 *full = new TF1("fitAllPeaks", f_full, minFit, maxFit);
+ //full->SetLineColor(kRed);
+
+
//Annoyingly long parameter name assignment
//(SetParNames only works for up to 11 variables)
full->SetParName(0,"Sigma");
@@ -260,6 +317,7 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
full->SetParName(25,"Mean 13"); full->SetParName(26,"Area*BinWidth 13");
full->SetParName(27,"Mean 14"); full->SetParName(28,"Area*BinWidth 14");
full->SetParName(29,"Background");
+ full->SetParName(30,"BreakupScale");
//Fix sigma & centroid, only allow area to vary
const int numParams = (numPeaks*2)+2;
@@ -295,7 +353,6 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
allPeaks[i]->Draw("SAME");
}
-
/* Error propogation:
* (Abin) +- deltaAbin, B+-0 (no uncertainty)
* A = Abin/B
@@ -323,3 +380,4 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
return allpeaks;
}
+
diff --git a/Projects/e793s/macro/ThreeBodyBreakup.h b/Projects/e793s/macro/ThreeBodyBreakup.h
index 62ca664d4335edf994b0dff855d0004a63404378..5bdbe366e5519a29bc968ee0aecbc199d02dda8c 100644
--- a/Projects/e793s/macro/ThreeBodyBreakup.h
+++ b/Projects/e793s/macro/ThreeBodyBreakup.h
@@ -1,5 +1,15 @@
-void ThreeBodyBreakup(){
+Double_t f_semi(Double_t *x, Double_t *par){
+ Float_t xx = x[0];
+ Double_t f;
+ Double_t a = TMath::Power(par[0],2);
+ Double_t b = TMath::Power(xx-par[1],2);
+ if(a > b){ f = par[2]*TMath::Sqrt(a-b); }
+ else{ f = 0; }
+ return f;
+}
+
+void ThreeBodyBreakup(){
TGenPhaseSpace PS1p;
//47K(d,p)47K+p+n
@@ -18,10 +28,10 @@ void ThreeBodyBreakup(){
double masses1p[3]={p.Mass(),K47.Mass(),n.Mass()};
//cout << p.Mass() << " " << K47.Mass()<< " " << n.Mass() << endl;
- auto Ex1p=new TH1D("ex1p","ex1p",300,-10,40);
+ auto Ex1p=new TH1D("ex1p","ex1p",500,-10,40);
PS1p.SetDecay(W,3,masses1p);
- for(unsigned int i = 0 ; i < 1000000 ; i++){
+ for(unsigned int i = 0 ; i < 10000000 ; i++){
// 1p case
double weight = PS1p.Generate();
auto pd = PS1p.GetDecay(0);
@@ -34,6 +44,30 @@ void ThreeBodyBreakup(){
Ex1p->Fill(ex,weight);
}
- Ex1p->Draw("");
+ TCanvas* c_Ex1p = new TCanvas("c_Ex1p","c_Ex1p",1000,500);
+ Ex1p->Draw();
+
+// TCanvas* c_semi = new TCanvas("c_semi","c_semi",1000,500);
+ TF1 *f1 = new TF1("f_semi",f_semi,-10,40,3);
+ f1->SetParameters(6,10,100);
+ f1->SetParNames("Radius","Centre","VertScale");
+
+ TFitResultPtr result = Ex1p->Fit("f_semi","S+");
+
+ f1->Draw("same");
+
+ cout << "Chi2 = " << result->Chi2()
+ << " NDF = " << result->Ndf()
+ << " Chi2/NDF = " << result->Chi2()/result->Ndf()
+ << endl;
+
+ cout << "==========================================" << endl;
+ cout << " Three-Body Breakup, B(x);" << endl;
+ cout << " B(x) = "
+ << result->Value(2) << " * SQRT( ("
+ << result->Value(0) << ")^2 - (x - "
+ << result->Value(1) << ")^2 )"
+ << endl;
+
}