Skip to content
Snippets Groups Projects
Commit 6de9ec75 authored by Charlie Paxman's avatar Charlie Paxman
Browse files

* e793s - Improved peak fitting 

+ easier to add peaks
* moved BG to parameter 0
parent c8095cd0
No related branches found
No related tags found
No related merge requests found
Pipeline #164695 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/e793s/macro/DrawPlots.C
+ 4
− 0
View file @ 6de9ec75
#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;
......
Projects/e793s/macro/DrawPlots.h
+ 33
− 28
View file @ 6de9ec75
......@@ -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");
......
Projects/e793s/macro/KnownPeakFitter.h
+ 75
− 134
View file @ 6de9ec75
......@@ -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;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment