* e793s - fit function to 3BBreakup

 * using semicircle function, translated in x and scaled in y;
 * f(x) = s * sqrt( r^2 - (x - t)^2 )
 * attempted to add this to background function, but failed so far

Projects/e793s/macro/DrawPlots.h

+#include <math.h>
 #include "NPReaction.h"
 #include <string>
 #include <sstream>

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;
 }
 
+ 

Projects/e793s/macro/ThreeBodyBreakup.h

-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;
+
 }