diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C
index 290f6a48b998a9f389e95aea98ecff8d2392b94a..9bfeae5f0fcf362f3fb75bd67c4c75aaf64623df 100755
--- a/Projects/e793s/macro/DrawPlots.C
+++ b/Projects/e793s/macro/DrawPlots.C
@@ -1,9 +1,11 @@
#include "GausFit.h"
+//#include "KnownPeakFitter_With47K.h"
#include "KnownPeakFitter.h"
#include "DrawPlots.h"
#include "CS2.h"
#include "ThreeBodyBreakup.h"
+#include "ThreeBodyBreakup_FitPhaseSpace.h"
/* USE THIS SPACE TO TEST NEW FEATURES */
void thickness(){
@@ -523,10 +525,12 @@ void DrawPlots(){
cout << "========== AVAILABLE FUNCTIONS ===========" << endl;
cout << " 2D Matrices " << endl;
cout << "\t- Draw_2DParticleGamma() "<< endl;
+ cout << "\t- Load_2DParticleGamma() "<< endl;
cout << "\t- Draw_2DGammaGamma() "<< endl;
cout << ""<< endl;
cout << " Ungated histograms " << endl;
cout << "\t- Draw_1DParticle() "<< endl;
+ cout << "\t- Load_1DParticle() "<< endl;
cout << "\t- Draw_1DParticle_MUST2() "<< endl;
cout << "\t- Draw_1DGamma() "<< endl;
cout << "\t- Draw_1DGamma_MG() "<< endl;
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 55122547f19667b1ee35a2117aa8cd710b569ea2..ebd6eaf71bcea2123cf557dbb47dd108b60db351 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -40,22 +40,13 @@ TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventLi
void LoadChainNP(){
vector<string> files;
- //files.push_back("../../../Outputs/Analysis/OriginalValues_ptI.root");
- //files.push_back("../../../Outputs/Analysis/OriginalValues_ptII.root");
- //files.push_back("../../../Outputs/Analysis/OriginalValues_ptIII.root");
-
- //files.push_back("../../../Outputs/Analysis/47K_Full_07Sep_PartI.root");
- //files.push_back("../../../Outputs/Analysis/47K_Full_07Sep_PartII.root");
-
- //files.push_back("../../../Outputs/Analysis/47Kdp_11Oct_PartI.root");
- //files.push_back("../../../Outputs/Analysis/47Kdp_11Oct_PartII.root");
-
- //files.push_back("../../../Outputs/Analysis/47Kdp_13Oct_wildTry_PartI.root");
- //files.push_back("../../../Outputs/Analysis/47Kdp_13Oct_wildTry_PartII.root");
-
files.push_back("../../../Outputs/Analysis/47Kdp_08Nov_PartI.root");
files.push_back("../../../Outputs/Analysis/47Kdp_08Nov_PartII.root");
+// files.push_back("../../../Outputs/Analysis/47Kdp_17Feb_AGATA_RotateBYpi.root");
+
+// files.push_back("../../../Outputs/Analysis/47Kdp_19Feb_AGATA_RotateBXYZ.root");
+
chain = Chain("PhysicsTree",files,true);
}
@@ -183,9 +174,16 @@ void Draw_1DGamma_MM(){
Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
}
+void Load_1DParticle(){
+ TH1F *hEx = new TH1F("hEx","Loaded 1D Particle Spectrum",200,-1,9);
+ TFile *file = new TFile("LoadHistograms/Load_1DParticle.root","READ");
+ hEx = (TH1F*)file->Get("Ep");
+ hEx->Draw();
+}
+
void Draw_1DParticle(){
TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000);
- chain->Draw("Ex>>Ep(160,-1,7)",
+ chain->Draw("Ex>>Ep(200,-1,9)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"");
TH1F* Ep = (TH1F*) gDirectory->Get("Ep");
@@ -198,7 +196,7 @@ void Draw_1DParticle(){
void Draw_1DParticleMUST2(){
TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000);
- chain->Draw("Ex>>Ep(160,-1,7)",
+ chain->Draw("Ex>>Ep(200,-1,9)",
"abs(T_MUGAST_VAMOS-2777)<600 && MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<4",
"");
TH1F* Ep = (TH1F*) gDirectory->Get("Ep");
@@ -207,9 +205,16 @@ void Draw_1DParticleMUST2(){
Ep->GetYaxis()->SetTitle("Counts / 0.05 MeV");
}
+void Load_2DParticleGamma(){
+ TH2F *hExEg = new TH2F("hExEg","Loaded 2D Particle-Gamma",200,-1,9,2500,0,5);
+ TFile *file = new TFile("LoadHistograms/Load_2DParticleGamma.root","READ");
+ hExEg = (TH2F*)file->Get("ExEg");
+ hExEg->Draw();
+}
+
void Draw_2DParticleGamma(){
TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000);
- chain->Draw("AddBack_EDC:Ex>>ExEg(140,-1,7,2500,0,5)",
+ chain->Draw("AddBack_EDC:Ex>>ExEg(200,-1,9,2500,0,5)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"colz");
TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg");
@@ -287,7 +292,7 @@ cout << " NO MG3 IN THIS ONE!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
string title = to_string(gamma-width)+" < Eg < "+to_string(gamma+width);
string ytitle = "Counts / " + to_string(binsize) + " MeV";
- string draw = "Ex>>ExGate(" + to_string(8.0/binsize) + ",-1,7)";
+ string draw = "Ex>>ExGate(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000);
//chain->Draw("Ex>>ExGate(60,-1,5)",gating.c_str(),"colz");
@@ -314,7 +319,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){
+ " BG: "+to_string(bg-width)+" to "+to_string(bg+width)+".";
TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000);
- chain->Draw("Ex>>ExGate(80,-1,7)",gating.c_str(),"");
+ chain->Draw("Ex>>ExGate(100,-1,9)",gating.c_str(),"");
//chain->Draw("Ex>>ExGate(120,-1,5)",gating.c_str(),"");
TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate");
ExGate->GetXaxis()->SetTitle("Ex [MeV]");
@@ -325,7 +330,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){
ExGate->SetFillStyle(3154);
ExGate->SetTitle(title.c_str());
- chain->Draw("Ex>>ExBG(80,-1,7)",bggate.c_str(),"same");
+ chain->Draw("Ex>>ExBG(100,-1,9)",bggate.c_str(),"same");
//chain->Draw("Ex>>ExBG(120,-1,5)",bggate.c_str(),"same");
TH1F* ExBG = (TH1F*) gDirectory->Get("ExBG");
ExBG->SetLineColor(kRed);
@@ -351,7 +356,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg, double
+ " BG: "+to_string(bg-width)+" to "+to_string(bg+width)+".";
TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000);
- chain->Draw("Ex>>ExGate(80,-1,7)",gating.c_str(),"");
+ chain->Draw("Ex>>ExGate(100,-1,9)",gating.c_str(),"");
//chain->Draw("Ex>>ExGate(120,-1,5)",gating.c_str(),"");
TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate");
ExGate->GetXaxis()->SetTitle("Ex [MeV]");
@@ -362,7 +367,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg, double
ExGate->SetFillStyle(3154);
ExGate->SetTitle(title.c_str());
- chain->Draw("Ex>>ExBG(80,-1,7)",bggate.c_str(),"same");
+ chain->Draw("Ex>>ExBG(100,-1,9)",bggate.c_str(),"same");
//chain->Draw("Ex>>ExBG(120,-1,5)",bggate.c_str(),"same");
TH1F* ExBG = (TH1F*) gDirectory->Get("ExBG");
ExBG->Scale(ratio);
@@ -621,7 +626,7 @@ void CompareSimExp(){
TCanvas *cSimExp = new TCanvas("cSimExp","cSimExp",1000,1000);
gStyle->SetOptStat(0);
- chain->Draw("Ex>>hexp(70,-1,6)","abs(T_MUGAST_VAMOS-2777)<600","");
+ chain->Draw("Ex>>hexp(100,-1,9)","abs(T_MUGAST_VAMOS-2777)<600","");
TH1F* hexp = (TH1F*) gDirectory->Get("hexp");
hexp->SetTitle("Comparing Simulation to Experiment");
hexp->GetXaxis()->SetTitle("Ex [MeV]");
@@ -632,7 +637,7 @@ void CompareSimExp(){
//TFile* simfile = new TFile("../../../Outputs/Analysis/SimTest_Jun22_TWOFNR_p32.root",
"READ");
TTree* simtree = (TTree*) simfile->FindObjectAny("PhysicsTree");
- simtree->Draw("Ex>>hsimMGp32(70,-1,6)",
+ simtree->Draw("Ex>>hsimMGp32(100,-1,9)",
"Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same");
TH1F* hsimMGp32 = (TH1F*) gDirectory->Get("hsimMGp32");
hsimMGp32->SetLineColor(kBlue);
@@ -890,7 +895,7 @@ void ExPhiLab_ForPoster(){
void ExThetaLab(){
TCanvas *diagnoseTheta = new TCanvas("diagnoseTheta","diagnoseTheta",1000,1000);
chain->Draw(
- "Ex:ThetaLab>>thetaHist(60,100,160,120,-1,5)",
+ "Ex:ThetaLab>>thetaHist(60,100,160,100,-1,9)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"colz");
TH1F* thetaHist = (TH1F*) gDirectory->Get("thetaHist");
@@ -939,7 +944,7 @@ void ExThetaLab(double gamma, double width){
string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ to_string(gamma) + ") < " + to_string(width);
- chain->Draw("Ex:ThetaLab>>thetaHist(60,100,160,60,-1,5)", gating.c_str(), "colz");
+ chain->Draw("Ex:ThetaLab>>thetaHist(60,100,160,100,-1,9)", gating.c_str(), "colz");
TH1F* thetaHist = (TH1F*) gDirectory->Get("thetaHist");
thetaHist->GetXaxis()->SetTitle("Theta (degrees)");
thetaHist->GetYaxis()->SetTitle("Ex [MeV]");
@@ -1181,7 +1186,7 @@ void GateThetaCM(double minTheta, double maxTheta, double binsize){
string title = to_string(minTheta)+" < ThetaCM < "+to_string(maxTheta);
string ytitle = "Counts / " + to_string(binsize) + " MeV";
- string draw = "Ex>>Ex_ThetaCMGate(" + to_string(8.0/binsize) + ",-1,7)";
+ string draw = "Ex>>Ex_ThetaCMGate(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaCMGate = new TCanvas("cEx_ThetaCMGate","cEx_ThetaCMGate",1000,1000);
chain->Draw(draw.c_str(),gating.c_str(),"colz");
@@ -1201,7 +1206,7 @@ void GateThetaLab(double minTheta, double maxTheta, double binsize){
string title = to_string(minTheta)+" < ThetaLab < "+to_string(maxTheta);
string ytitle = "Counts / " + to_string(binsize) + " MeV";
- string draw = "Ex>>Ex_ThetaLabGate(" + to_string(8.0/binsize) + ",-1,7)";
+ string draw = "Ex>>Ex_ThetaLabGate(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaLabGate = new TCanvas("cEx_ThetaLabGate","cEx_ThetaLabGate",1000,1000);
chain->Draw(draw.c_str(),gating.c_str(),"colz");
@@ -1227,7 +1232,7 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
+ " && ThetaLab < "
+ to_string(minTheta+gatesize);
string histname = "cThetaLabGate_" + to_string(minTheta) + "-" + to_string(minTheta+gatesize);
- string draw = "Ex>>" + histname + "(" + to_string(8.0/binsize) + ",-1,7)";
+ string draw = "Ex>>" + histname + "(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaLabGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000);
chain->Draw(draw.c_str(),gating.c_str(),"colz");
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index c50f71d2cb5ff7300ff804477ab392d33dd0a9b3..3bacf35cd1125efa8167ade966514355f428d456 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -3,23 +3,24 @@
#include <cmath>
#include "stdlib.h"
- const int numPeaks = 14;
- array<double,14> means = { 0.000,
- 0.143,
- 0.279,
- 0.728,
- 0.968,
- 1.410,
- 1.981,
- 2.410,
- 2.910,
- 3.2,
- 3.605,
- 3.792,
- 4.1,
- 4.4
+const int numPeaks = 15;
+array<double,numPeaks> means = { 0.000,
+ 0.143,
+ 0.279,
+ 0.728,
+ 0.968,
+ 1.410,
+ 1.981,
+ 2.410,
+ 2.910,
+ 3.2,
+ 3.605,
+ 3.792,
+ 4.1,
+ 4.4,
+ 5.24
};
-
+/*
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];
@@ -32,72 +33,57 @@ Double_t f_bg(Double_t *x, Double_t *par){
}
Double_t f_peak(Double_t *x, Double_t *par){
- Float_t xx = x[0];
- Double_t f = (par[2]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[1])/par[0],2));
+ float xx = x[0];
+ double f = (par[2]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[1])/par[0],2));
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;
+*/
+
+Double_t f_full(Double_t *x, Double_t *par) {
+ float xx = x[0];
+ double result, norm;
+ // Flat background
+ result = par[0];
+ // Add N peaks
+ for(int pk=0; pk<numPeaks; pk++){
+ result += (par[3+(pk*3)]/(par[1+(pk*3)]*sqrt(2*pi)))
+ * exp(-0.5*pow((xx-par[2+(pk*3)])/par[1+(pk*3)],2));
+ }
+ return result;
}
vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
- double minFit=-1.0, maxFit=5.0;
+ double minFit=-1.0, maxFit=8.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
double sigma = 0.14;
- string nameBase = "Peak ";
- string function = "([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2))";
-
- //Assign peak centroids
-// const int numPeaks = 14;
-// array<double,14> means = { 0.000,
-// 0.143,
-// 0.279,
-// 0.728,
-// 0.968,
-// 1.410,
-// 1.981,
-// 2.410,
-// 2.910,
-// 3.2,
-// 3.605,
-// 3.792,
-// 4.1,
-// 4.4
-// };
+ hist->Sumw2();
+
+ /* Construct flat BG to subtract */
+ cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
+ cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
+ cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
+ double ConstBG = 36.0; double ErrBG = 1.0;
+ int xbins = hist->GetXaxis()->GetNbins();
+ double xmin = hist->GetXaxis()->GetXmin();
+ double xmax = hist->GetXaxis()->GetXmax();
+ TH1F *FlatBG = new TH1F("FlatBG","FlatBG", xbins, xmin, xmax);
+ for(int i=0; i<xbins;i++){
+ FlatBG->SetBinContent(i,ConstBG);
+ FlatBG->SetBinError(i,ErrBG);
+ }
+ hist->Add(FlatBG,-1);
+ /**/
//Build individual peak fit functions
+ string nameBase = "Peak ";
+ string function = "([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2))";
TF1 **allPeaks = new TF1*[numPeaks];
for(int i=0; i<numPeaks; i++) {
string nameHere = nameBase;
nameHere +=to_string(i);
- allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), -1, 5);
+ allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit);
//allPeaks[i] = new TF1(nameHere.c_str(), f_peak, -1, 5);
allPeaks[i]->SetLineColor(kBlack);
allPeaks[i]->SetLineStyle(7);
@@ -110,71 +96,24 @@ 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)) "
- "+ ([4]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[3])/[0],2)) "
- "+ ([6]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[5])/[0],2)) "
- "+ ([8]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[7])/[0],2)) "
- "+ ([10]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[9])/[0],2)) "
- "+ ([12]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[11])/[0],2)) "
- "+ ([14]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[13])/[0],2)) "
- "+ ([16]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[15])/[0],2)) "
- "+ ([18]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[17])/[0],2)) "
- "+ ([20]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[19])/[0],2)) "
- "+ ([22]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[21])/[0],2)) "
- "+ ([24]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[23])/[0],2)) "
- "+ ([26]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[25])/[0],2)) "
- "+ ([28]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[27])/[0],2)) "
- "+ [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");
- full->SetParName(1,"Mean 01"); full->SetParName(2,"Area*BinWidth 01");
- full->SetParName(3,"Mean 02"); full->SetParName(4,"Area*BinWidth 02");
- full->SetParName(5,"Mean 03"); full->SetParName(6,"Area*BinWidth 03");
- full->SetParName(7,"Mean 04"); full->SetParName(8,"Area*BinWidth 04");
- full->SetParName(9,"Mean 05"); full->SetParName(10,"Area*BinWidth 05");
- full->SetParName(11,"Mean 06"); full->SetParName(12,"Area*BinWidth 06");
- full->SetParName(13,"Mean 07"); full->SetParName(14,"Area*BinWidth 07");
- full->SetParName(15,"Mean 08"); full->SetParName(16,"Area*BinWidth 08");
- full->SetParName(17,"Mean 09"); full->SetParName(18,"Area*BinWidth 09");
- full->SetParName(19,"Mean 10"); full->SetParName(20,"Area*BinWidth 10");
- full->SetParName(21,"Mean 11"); full->SetParName(22,"Area*BinWidth 11");
- full->SetParName(23,"Mean 12"); full->SetParName(24,"Area*BinWidth 12");
- 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;
- full->FixParameter(0, sigma);
+ TF1 *full = new TF1("fitAllPeaks", f_full, minFit, maxFit, (int) 1+(3*numPeaks));
+ full->SetLineColor(kRed);
+ const int numParams = (numPeaks*3)+1;
for(int i=0; i<numPeaks; i++) {
- full->FixParameter(2*i+1,means.at(i));
- full->SetParameter(2*i+2,1.0);
- full->SetParLimits(2*i+2,0.0,1e5);
+ full->FixParameter((i*3)+1,sigma);
+ full->FixParameter((i*3)+2,means.at(i));
+ full->SetParameter((i*3)+3,1e1);
+ full->SetParLimits((i*3)+3,0.0,1e5);
}
+ //full->SetParameter(0,30.);
+ //full->SetParLimits(0,0.,40.);
+ full->FixParameter(0,0.);
- //TESTING!!!!!
- //full->FixParameter(6,0.);
-
- //full->FixParameter(numParams-1,0.0);
- full->SetParameter(numParams-1,1.0);
- full->SetParLimits(numParams-1,0.0,1e1);
+ // Specific limits
+ // Set max 279 counts to 100
+ full->SetParLimits(9,0.0,1e2); // Doesnt work???
+ allPeaks[14]->SetLineColor(kOrange);
//Fit full function to histogram
hist->Fit(full, "RQ", "", minFit, maxFit);
@@ -184,9 +123,9 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
const Double_t* finalPar = full->GetParameters();
const Double_t* finalErr = full->GetParErrors();
for (int i=0; i<numPeaks; i++){
- allPeaks[i]->SetParameters(sigma, means.at(i), finalPar[2*i+2]);
+ allPeaks[i]->SetParameters(sigma, means.at(i), finalPar[3+(i*3)]);
}
- bg->SetParameter(0,finalPar[numParams-1]);
+ bg->SetParameter(0,finalPar[0]);
bg->Draw("SAME");
full->Draw("SAME");
@@ -207,22 +146,24 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
vector<vector<double>> allpeaks;
for(int i=0; i<numPeaks; i++){
- double A = finalPar[2*i+2]/binWidth;
- double deltaA = A * (finalErr[2*i+2]/finalPar[2*i+2]);
+ double A = finalPar[(i*3)+3]/binWidth;
+ double deltaA = A * (finalErr[(i*3)+3]/finalPar[(i*3)+3]);
cout << fixed << setprecision(3)
<< " #" << i << " "
- << finalPar[2*i+1] << "\t" << setprecision(0)
+ << finalPar[(i*3)+2] << "\t" << setprecision(0)
<< A << "\t+- "
<< deltaA << setprecision(3)
<< endl;
vector<double> onepeak; //energy, area and error for one peak
- onepeak.push_back(finalPar[2*i+1]);
+ onepeak.push_back(finalPar[(i*3)+2]);
onepeak.push_back(A);
onepeak.push_back(deltaA);
allpeaks.push_back(onepeak);
}
+ cout << " BG " << full->GetParameter(0)
+ << " +- " << full->GetParError(0) << endl;
return allpeaks;
}