From b1aa18174269e71764c407ce0674507e1ddfb57d Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Wed, 9 Feb 2022 11:48:28 +0000
Subject: [PATCH] * e793s - CS() expanded to any given state
---
Projects/e793s/macro/CS2.h | 102 ++++++---
Projects/e793s/macro/DrawPlots.C | 2 +-
Projects/e793s/macro/DrawPlots.h | 38 +++-
Projects/e793s/macro/KnownPeakFitter.h | 243 ++++------------------
Projects/e793s/macro/minimizer.cxx | 276 ++++++++++++++++++++-----
5 files changed, 365 insertions(+), 296 deletions(-)
mode change 100644 => 100755 Projects/e793s/macro/minimizer.cxx
diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
index 2fabc6730..e3f0e915b 100644
--- a/Projects/e793s/macro/CS2.h
+++ b/Projects/e793s/macro/CS2.h
@@ -9,35 +9,58 @@ vector<Double_t> anglecentres, anglewidth;
TGraph* currentThry;
TGraphErrors* staticExp;
-//BASIC RUN: CS(0.143,2,1,1.5,1)
+int indexE;
+bool loud = 0;
////////////////////////////////////////////////////////////////////////////////
void CS(double Energy, double Spin, double spdf, double angmom, double nodes){
+
+ /* BASIC RUN: CS(0.143,2,1,1.5,1) */
+
// p3/5 -> spdf = 1, angmom = 1.5
// J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
double J0 = 0.5;
+ /* Retrieve array index of the entered peak energy */
+ /* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
+ bool found = 0;
+ for(int i=0;i<numPeaks;i++){
+ if(abs(Energy-means[i])<0.01){
+ cout << "========================================================" << endl;
+ cout << "Identified as state #" << i << ", E = " << means[i] << endl;
+ indexE = i;
+ found = 1;
+ }
+ }
+ if(!found){
+ cout << "========================================================" << endl;
+ cout << "NO STATE AT THAT ENERGY INDENTIFIED!! CHECK KNOWN PEAKS!!" << endl;
+ return;
+ }
+
/* Solid Angle (from simulation) */
auto file = new TFile("../SolidAngle_HistFile_06Dec_47Kdp.root");
TH1F* SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_Lab_MG");
TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
SolidAngle->Draw();
+ // canvas deleted after Area/SA calculation
/* Area of experimental peaks */
TCanvas* c_PeakArea = new TCanvas("c_PeakArea","c_PeakArea",1000,1000);
- vector<vector<double>> areaArray = ExpDiffCross(0.143);
- //cout << std::setprecision(3) << std::fixed;
+ vector<vector<double>> areaArray = ExpDiffCross(means[indexE]);
delete c_PeakArea;
- // Array: i, peakenergy, peakarea, areaerror, anglemin, anglemax
-// for(int i=0; i<areaArray.size();i++){
-// cout << i << " "
-// << areaArray[i][0] << " "
-// << areaArray[i][1] << " "
-// << areaArray[i][2] << " "
-// << areaArray[i][3] << " "
-// << areaArray[i][4] << endl;
-// }
+ // Array: peakenergy, peakarea, areaerror, anglemin, anglemax
+ if(loud){
+ for(int i=0; i<areaArray.size();i++){
+ cout << i << " "
+ << areaArray[i][0] << " "
+ << areaArray[i][1] << " "
+ << areaArray[i][2] << " "
+ << areaArray[i][3] << " "
+ << areaArray[i][4] << endl;
+ }
+ }
/* Area over Solid Angle */
vector<Double_t> AoSA, AoSAerr;
@@ -52,15 +75,19 @@ void CS(double Energy, double Spin, double spdf, double angmom, double nodes){
double SA = 1000. * SolidAngle->IntegralAndError(binmin,binmax,SAerr);
//SAerr = SAerr*1000.; //????
- AoSA.push_back(areaArray[i][1] / SA);
+ AoSA.push_back(areaArray[i][1]/SA);
AoSAerr.push_back((areaArray[i][1]/SA)
* sqrt( pow(areaArray[i][2]/areaArray[i][1],2) + pow(SAerr/SA,2) ) );
- cout << "Angle " << areaArray[i][3] << " to " << areaArray[i][4]
- << ", centre " << anglecentres[i]
- << ": SA = " << SA << " +- " << SAerr
- << " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"<< endl;
+ if(loud){
+ cout << "Angle " << areaArray[i][3] << " to " << areaArray[i][4]
+ << ", centre " << anglecentres[i]
+ << ": Area = " << areaArray[i][1] << " +- " << areaArray[i][2]
+ << " SA = " << SA << " +- " << SAerr
+ << " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"<< endl;
+ }
}
+ delete c_SolidAngle;
TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
TGraphErrors* gAoSA = new TGraphErrors(
@@ -75,7 +102,7 @@ void CS(double Energy, double Spin, double spdf, double angmom, double nodes){
/* TWOFNR diff. cross section */
TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
- TGraph* DiffCross = TWOFNR(Energy, J0, Spin, nodes, spdf, angmom);
+ TGraph* DiffCross = TWOFNR(means[indexE], J0, Spin, nodes, spdf, angmom);
DiffCross->Draw();
/* Scaled and compared on same plot */
@@ -103,12 +130,9 @@ void CS(double Energy, double Spin, double spdf, double angmom, double nodes){
Final->Draw("SAME");
}
-
-
////////////////////////////////////////////////////////////////////////////////
vector<vector<double>> ExpDiffCross(double Energy){
cout << "========================================================" << endl;
- // Implement some kind of energy selection later, for now, hard code that we are selecting 0.143, array index 1
vector<vector<double>> AllPeaks_OneGate;
vector<vector<double>> OnePeak_AllGates;
int numbins = 10;
@@ -120,15 +144,27 @@ vector<vector<double>> ExpDiffCross(double Energy){
cout << "===================================" << endl;
cout << "min: " << min << " max: " << max << endl;
+ /* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
+ /* To change angle gates, run GateThetaLab_MultiWrite() */
TH1F* gate = PullThetaLabHist(i,105.,5.);
+
+ /* Retrieve array containing all fits, for one angle gate. *
+ * Specific peak of interest selected from the vector by *
+ * global variable indexE */
AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate);
- //cout << "area of " << peakAreas[1][0] << " = " << peakAreas[1][1]
- cout << "area of 0.143 = " << AllPeaks_OneGate[1][1]
- << " +- " << AllPeaks_OneGate[1][2]
+
+ /* Check correct OneGate vector is selected */
+ cout << "area of " << means[indexE] << " = "
+ << AllPeaks_OneGate[indexE][1]
+ << " +- " << AllPeaks_OneGate[indexE][2]
<< endl;
- AllPeaks_OneGate[1].push_back(min);
- AllPeaks_OneGate[1].push_back(max);
- OnePeak_AllGates.push_back(AllPeaks_OneGate[1]);
+
+ /* Add min and max angle to end of relevant OneGate vector */
+ AllPeaks_OneGate[indexE].push_back(min);
+ AllPeaks_OneGate[indexE].push_back(max);
+
+ /* Push relevant OneGate vector to end of AllGates vector */
+ OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
}
return OnePeak_AllGates;
}
@@ -174,6 +210,8 @@ void Scale(TGraph* g , TGraphErrors* ex){
////////////////////////////////////////////////////////////////////////////////
TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
+ /* This function mved between directories in order to run TWOFNR in proper *
+ * location. This is, weirdly, the least tempremental way of doing this. */
cout << "========================================================" << endl;
char origDirchar[200];
@@ -181,7 +219,6 @@ TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
string origDir{origDirchar};
string twofnrDir = "/home/charlottepaxman/Programs/TWOFNR";
cout << "Current directory " << origDir << endl;
-
cout << "Moving to directory " << twofnrDir << endl;
chdir(twofnrDir.c_str());
//Check
@@ -228,23 +265,20 @@ TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
cout << "Filled front input file." << endl;
cout << "Executing front20..." << endl;
system("(exec ~/Programs/TWOFNR/front20 < in.front > /dev/null)");
+ cout << "-> front20 complete!" << endl;
cout << "Executing twofnr20..." << endl;
system("(exec ~/Programs/TWOFNR/twofnr20 < in.twofnr > /dev/null)");
-
- cout << "twofnr20 complete." << endl;
-
+ cout << "-> twofnr20 complete!" << endl;
TGraph* CS = new TGraph("24.jjj");
cout << "===================================" << endl;
cout << "Current directory " << twofnrDir << endl;
-
cout << "Moving to directory " << origDir << endl;
chdir(origDir.c_str());
system("pwd");
cout << "========================================================" << endl;
-
return CS;
}
@@ -259,7 +293,7 @@ double Chi2(TGraph* theory , TGraphErrors* exper){
Chi2 +=chi*chi;
}
}
- cout << "Chi2 = " << Chi2 << endl;
+ if(loud){cout << "Chi2 = " << Chi2 << endl;}
return Chi2;
}
diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C
index 314a6b4dd..9ec479d5d 100755
--- a/Projects/e793s/macro/DrawPlots.C
+++ b/Projects/e793s/macro/DrawPlots.C
@@ -554,8 +554,8 @@ cout << " 2D Matrices " << endl;
cout << ""<< endl;
cout << " Analysis functions" << endl;
cout << "\t- FitKnownPeaks(histogram) "<< endl;
- cout << "\t\t-- Fits Ex peaks to an excitation spectrum "<< endl;
cout << "\t- AGATA_efficiency(double Energy_kev) "<< endl;
+ cout << "\t- CorrectForAGATAEffic(TH1F* hist) "<< endl;
cout << ""<< endl;
cout << "==========================================" << endl;
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 9c6462d00..55122547f 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -21,6 +21,11 @@ NPL::Reaction Ti12C12C("47Ti(12C,12C)47Ti@355");
void KnownLines_Ex(bool isVertical, double rangemin, double rangemax, Style_t lType, Color_t lColour);
/* BASE FUNCTIONS */
+TF1* f_efficAGATA(){
+ TF1 *f_E = new TF1("fit_1","TMath::Exp([0]+[1]*TMath::Log(x)+[2]*pow(TMath::Log(x),2.0)+[3]*pow(TMath::Log(x),3.0)+[4]*pow(TMath::Log(x),4.0))",10,6000);
+ f_E->SetParameters(-6.34543e+01, +4.24746e+01, -1.00304e+01, +1.03468e+00, -3.97076e-02);
+ return f_E;
+}
TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventList){
TChain* chain = new TChain(TreeName.c_str());
@@ -54,6 +59,24 @@ void LoadChainNP(){
chain = Chain("PhysicsTree",files,true);
}
+void CorrectForAGATAEffic(TH1F* hist){
+ int nbins = hist->GetNbinsX();
+ TF1* effic_keV = f_efficAGATA();
+
+ for(int i=1;i<nbins+1;i++){
+ double x = hist->GetBinCenter(i);
+ double val = effic_keV->Eval(x*1000.);
+// cout << x << " " << val << endl;
+
+ double bincontent = hist->GetBinContent(i);
+ bincontent = bincontent/(val/100.);
+ hist->SetBinContent(i,bincontent);
+ }
+ hist->SetFillStyle(3244);
+ hist->SetFillColor(kBlue);
+ hist->Draw();
+}
+
void DrawParticleStates(TCanvas* canvas){
canvas->cd();
@@ -1134,17 +1157,10 @@ void AGATA_efficiency(){
}
void AGATA_efficiency(double Energy_keV){
- TF1 *fit_1 = new TF1("fit_1","TMath::Exp([0]+[1]*TMath::Log(x)+[2]*pow(TMath::Log(x),2.0)+[3]*pow(TMath::Log(x),3.0)+[4]*pow(TMath::Log(x),4.0))",10,5000);
- fit_1->SetParameters(-6.34543e+01,
- +4.24746e+01,
- -1.00304e+01,
- +1.03468e+00,
- -3.97076e-02);
- fit_1->Draw();
- cout << "At E = "
- << Energy_keV
- << " keV, AGATA efficiency = "
- << fit_1->Eval(Energy_keV)
+ TF1* func = f_efficAGATA();
+ func->Draw();
+ cout << "At E = " << Energy_keV
+ << " keV, AGATA efficiency = " << func->Eval(Energy_keV)
<< " %" << endl;
}
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index a249438cf..2c24c4395 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -3,37 +3,8 @@
#include <cmath>
#include "stdlib.h"
-void FitKnownPeaks(TH1F* hist){
- double minFit=-1.0, maxFit=5.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))";
-
- /* 11 KNOWN PEAKS as of 12/10/21
- * 0.000
- * 0.143
- * 0.279
- * 0.728
- * 0.968
- * 1.410
- * 1.981
- * 2.410
- * 2.910
- * 3.605
- * 3.792
- * ...
- * NEW:
- * 3.2
- * 4.1
- * 4.4
- *
- */
-
- //Assign peak centroids
- const int numPeaks = 14;
- array<double,14> means = { 0.000,
+ const int numPeaks = 14;
+ array<double,14> means = { 0.000,
0.143,
0.279,
0.728,
@@ -49,122 +20,6 @@ void FitKnownPeaks(TH1F* hist){
4.4
};
- //Build individual peak fit functions
- 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]->SetLineColor(kBlack);
- allPeaks[i]->SetLineStyle(7);
- allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth");
- }
-
- //Build background function
- TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
- bg->SetLineColor(kGreen);
- bg->SetLineStyle(9);
- bg->SetParNames("Background");
-
- //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);
-
- //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");
-
- //Fix sigma & centroid, only allow area to vary
- const int numParams = (numPeaks*2)+2;
- full->FixParameter(0, sigma);
- 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(numParams-1,0.0);
- full->SetParameter(numParams-1,1.0);
- full->SetParLimits(numParams-1,0.0,1e1);
-
- //Fit full function tohistogram
- hist->Fit(full, "WWR", "", minFit, maxFit);
- hist->Draw();
-
- //Extract fitted variables, assign them to individual fits, and draw them
- Double_t finalPar[numParams];
- Double_t finalErr[numParams];
- full->GetParameters(&finalPar[0]);
- for (int i=0; i<numPeaks; i++){
- finalErr[2*i+1] = full->GetParError(2*i+1);
- finalErr[2*i+2] = full->GetParError(2*i+2);
-
- allPeaks[i]->SetParameters(sigma, means.at(i), finalPar[2*i+2]);
- }
- bg->SetParameter(0,finalPar[numParams-1]);
- bg->Draw("SAME");
- full->Draw("SAME");
-
- for (int i=0; i<numPeaks; i++){
- allPeaks[i]->Draw("SAME");
- }
-
-
- /* Error propogation:
- * (Abin) +- deltaAbin, B+-0 (no uncertainty)
- * A = Abin/B
- * deltaA/A = deltaAbin/Abin
- * deltaA = A x deltaAbin/Abin
- */
-
- //Write to screen
- cout << "===========================" << endl;
- cout << "== PEAK =========== AREA ==" << endl;
-
- for(int i=0; i<numPeaks; i++){
- cout << fixed << setprecision(3) << finalPar[2*i+1]
- << "\t" << setprecision(0)<< finalPar[2*i+2]/binWidth
- << "\t+- " << (finalPar[2*i+2]/binWidth)*(finalErr[2*i+2]/finalPar[2*i+2])
- << endl;
- }
-
-}
-
-
-
-
-
-
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];
@@ -176,6 +31,12 @@ Double_t f_bg(Double_t *x, Double_t *par){
return f;
}
+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));
+ return f;
+}
+
Double_t f_full(Double_t *x, Double_t *par){
Float_t xx = x[0];
Double_t f;
@@ -204,8 +65,6 @@ Double_t f_full(Double_t *x, Double_t *par){
return f;
}
-
-
vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
double minFit=-1.0, maxFit=5.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
@@ -214,42 +73,23 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
string nameBase = "Peak ";
string function = "([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2))";
- /* 11 KNOWN PEAKS as of 12/10/21
- * 0.000
- * 0.143
- * 0.279
- * 0.728
- * 0.968
- * 1.410
- * 1.981
- * 2.410
- * 2.910
- * 3.605
- * 3.792
- * ...
- * NEW:
- * 3.2
- * 4.1
- * 4.4
- *
- */
//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
- };
+// 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
+// };
//Build individual peak fit functions
TF1 **allPeaks = new TF1*[numPeaks];
@@ -258,6 +98,7 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
nameHere +=to_string(i);
allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), -1, 5);
+ //allPeaks[i] = new TF1(nameHere.c_str(), f_peak, -1, 5);
allPeaks[i]->SetLineColor(kBlack);
allPeaks[i]->SetLineStyle(7);
allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth");
@@ -332,17 +173,13 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
full->SetParLimits(numParams-1,0.0,1e1);
//Fit full function to histogram
- hist->Fit(full, "WWRQ", "", minFit, maxFit);
+ hist->Fit(full, "RQ", "", minFit, maxFit);
hist->Draw();
-
+
//Extract fitted variables, assign them to individual fits, and draw them
- Double_t finalPar[numParams];
- Double_t finalErr[numParams];
- full->GetParameters(&finalPar[0]);
+ const Double_t* finalPar = full->GetParameters();
+ const Double_t* finalErr = full->GetParErrors();
for (int i=0; i<numPeaks; i++){
- finalErr[2*i+1] = full->GetParError(2*i+1);
- finalErr[2*i+2] = full->GetParError(2*i+2);
-
allPeaks[i]->SetParameters(sigma, means.at(i), finalPar[2*i+2]);
}
bg->SetParameter(0,finalPar[numParams-1]);
@@ -366,18 +203,26 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
vector<vector<double>> allpeaks;
for(int i=0; i<numPeaks; i++){
- cout << fixed << setprecision(3) << finalPar[2*i+1]
- << "\t" << setprecision(0)<< finalPar[2*i+2]/binWidth
- << "\t+- " << (finalPar[2*i+2]/binWidth)*(finalErr[2*i+2]/finalPar[2*i+2])
- << setprecision(3) << endl;
+ double A = finalPar[2*i+2]/binWidth;
+ double deltaA = A * (finalErr[2*i+2]/finalPar[2*i+2]);
+
+ cout << fixed << setprecision(3)
+ << " #" << i << " "
+ << finalPar[2*i+1] << "\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[2*i+2]/binWidth);
- onepeak.push_back((finalPar[2*i+2]/binWidth)*(finalErr[2*i+2]/finalPar[2*i+2]));
+ onepeak.push_back(A);
+ onepeak.push_back(deltaA);
allpeaks.push_back(onepeak);
}
return allpeaks;
}
-
+void FitKnownPeaks(TH1F* hist){
+ //Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
+ vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist);
+}
diff --git a/Projects/e793s/macro/minimizer.cxx b/Projects/e793s/macro/minimizer.cxx
old mode 100644
new mode 100755
index 6fb732460..2056cb4d3
--- a/Projects/e793s/macro/minimizer.cxx
+++ b/Projects/e793s/macro/minimizer.cxx
@@ -1,46 +1,200 @@
+#include "Math/Minimizer.h"
+#include "Math/Factory.h"
+#include "Math/Functor.h"
+#include "TRandom2.h"
+#include "TError.h"
+#include <iostream>
-double ref = 0.143; // the energy of the selected states
+double refE = 0.143; // the energy of the selected states
vector<TVector3*> pos;
vector<double> energy;
+vector<int> detnum;
NPL::EnergyLoss CD2("proton_CD2.G4table","G4Table",100);
NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
+using namespace std;
+
+bool flagDraw = 0;
+static auto h = new TH1D("h","h", 80,-1.,1.);
+static auto h1 = new TH1D("h1","h1", 40,-1.,1.);
+static auto h2 = new TH1D("h2","h2", 40,-1.,1.);
+static auto h3 = new TH1D("h3","h3", 40,-1.,1.);
+static auto h4 = new TH1D("h4","h4", 40,-1.,1.);
+static auto h5 = new TH1D("h5","h5", 40,-1.,1.);
+static auto h7 = new TH1D("h7","h7", 40,-1.,1.);
+
+void InitiliseCanvas(){
+ TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
+ gStyle->SetOptStat(0);
+ canv->Divide(2,1,0.005,0.005,0);
+ canv->cd(1)->SetLeftMargin(0.15);
+ canv->cd(1)->SetBottomMargin(0.15);
+ gPad->SetTickx();
+ gPad->SetTicky();
+ canv->cd(2)->SetLeftMargin(0.15);
+ canv->cd(2)->SetBottomMargin(0.15);
+ gPad->SetTickx();
+ gPad->SetTicky();
+
+ canv->cd(1);
+ h1->SetMaximum(75.);
+ h1->GetXaxis()->SetTitle("Ex [MeV]");
+ h1->GetYaxis()->SetTitle("Counts");
+
+ // ----- MG1 -----
+ h1->SetStats(0);
+ h1->SetLineColor(kRed);
+ h1->SetFillStyle(3244);
+ h1->SetFillColor(kRed);
+ h1->Draw();
+ h1->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG2 -----
+ h2->SetStats(0);
+ h2->SetLineColor(kOrange);
+ h2->SetFillStyle(3244);
+ h2->SetFillColor(kOrange);
+ h2->Draw("same");
+ h2->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG3 -----
+ h3->SetStats(0);
+ h3->SetLineColor(kGreen);
+ h3->SetFillStyle(3344);
+ h3->SetFillColor(kGreen);
+ h3->Draw("same");
+ h3->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG4 -----
+ h4->SetStats(0);
+ h4->SetLineColor(kTeal);
+ h4->SetFillStyle(3444);
+ h4->SetFillColor(kTeal);
+ h4->Draw("same");
+ h4->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG5 -----
+ h5->SetStats(0);
+ h5->SetLineColor(kBlue);
+ h5->SetFillStyle(3544);
+ h5->SetFillColor(kBlue);
+ h5->Draw("same");
+ h5->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG7 -----
+ h7->SetStats(0);
+ h7->SetLineColor(kViolet);
+ h7->SetFillStyle(3644);
+ h7->SetFillColor(kViolet);
+ h7->Draw("same");
+ h7->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // Format legend
+ auto legend = new TLegend(0.15,0.7,0.35,0.9);
+ legend->AddEntry(h1,"MUGAST 1","f");
+ legend->AddEntry(h2,"MUGAST 2","f");
+ legend->AddEntry(h3,"MUGAST 3","f");
+ legend->AddEntry(h4,"MUGAST 4","f");
+ legend->AddEntry(h5,"MUGAST 5","f");
+ legend->AddEntry(h7,"MUGAST 7","f");
+ legend->Draw();
+
+ // ----- ALL -----
+ canv->cd(2);
+ h->SetStats(0);
+ h->GetXaxis()->SetTitle("Ex [MeV]");
+ h->GetYaxis()->SetTitle("Counts");
+ h->Draw();
+ h->Fit("gaus", "WQ");
+ gPad->Update();
+ //delete c1;
+ //c1 = 0;
+}
////////////////////////////////////////////////////////////////////////////////
double devE(const double* parameter){
- // My reaction at 8A MeV*47A=376 MeV
- static NPL::Reaction reaction("47K(d,p)48K@376");
- // Return the standard deviation in Brho
- unsigned int size = pos1.size();
+ //Beam energy: 7.7 [MeV/A] * 47 [A] = 361.9 [MeV]
+ static NPL::Reaction reaction("47K(d,p)48K@362");
+ //Beam spot offset
TVector3 offset(parameter[0],parameter[1],parameter[2]);
+ unsigned int size = pos.size();
double dE,Theta;
TVector3 dir;
- static auto h = new TH1D("h","h", 1000,-0,100);
+
+ //Initilize histogram
h->Reset();
+ h1->Reset();
+ h2->Reset();
+ h3->Reset();
+ h4->Reset();
+ h5->Reset();
+ h7->Reset();
+
+ //Loop over events
for(unsigned int i = 0 ; i < size ; i++){
- dir=*(pos[i])-o1;
+ //Particle path vector
+ dir=*(pos[i])-offset;
+
+ //Detected energy, and angle of particle leaving target
double Theta= dir.Angle(TVector3(0,0,1));
double Energy = energy[i];
- // do some energy loss correction here:
- // This need to take parameter[4]*0.5*micrometer as the target thickness
- Energy=CD2.EvaluateInitialEnergy(Energy,0,5*parameter[4]*micrometer,Theta);
-
- double Ex=reaction->ReconstructRelativistic(Energy,Theta);
- //cout << (pB-pM).Mag2() << " " << dR << endl;
- h->Fill(Ex-ref);
- }
- cout << h->GetMean() << " " << h->GetStdDev() << " " << parameter[4] << endl;
+
+ //NOTE!!! Not calucualting energy loss in Al???
+ //Energy loss in target
+ Energy=CD2.EvaluateInitialEnergy(
+ Energy, //energy after leaving target
+ 0.5*parameter[4]*micrometer, //pass through half target
+ Theta); //angle leaving target
+
+ //Final value of Ex
+ double Ex = reaction.ReconstructRelativistic(Energy,Theta);
+
+ //Fill histogram with ERROR in Ex!
+ //h->Fill(Ex-refE);
+ h->Fill(Ex);
+
+ switch(detnum[i]){
+ case 1:
+ h1->Fill(Ex);
+ break;
+ case 2:
+ h2->Fill(Ex);
+ break;
+ case 3:
+ h3->Fill(Ex);
+ break;
+ case 4:
+ h4->Fill(Ex);
+ break;
+ case 5:
+ h5->Fill(Ex);
+ break;
+ case 7:
+ h7->Fill(Ex);
+ break;
+ default:
+ cout << "ERROR! Invalid detnum: " << detnum[i] << " @" << i << endl;
+ return 1; // Exit code
+ }
+
+ }
+ //End loop over events
+
+ //Write vals to screen
+ cout << "Mean: " << h->GetMean() << "\t StdDev: " << h->GetStdDev() << "\t Thickness??: " << parameter[4] << endl;
h->Draw();
- gPad->Update();
- // adapt the metric as needed
- return (h->GetMean()+h->GetStdDev() );
+
+ if(flagDraw){ InitiliseCanvas(); }
+
+ //Adapt the metric as needed
+ return sqrt( pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
}
////////////////////////////////////////////////////////////////////////////////
void LoadFile(){
-
- ifstream file("myFile.txt");
+ // Open XYZE gamma-gated file
+ ifstream file("BeamSpot/XYZE_Full_02June.txt");
if(!file.is_open()){
cout << "fail to load file" << endl;
exit(1);
@@ -49,44 +203,64 @@ void LoadFile(){
cout << "Success opening file" << endl;
}
+ // Read in
+ int mg;
double x,y,z,e;
-
- while(file >> x >> y >> z >> e ){
- auto p1 = new TVector3(x,y,z);
- double e;
- pos.push_back(p1);
+ while(file >> mg >> x >> y >> z >> e ){
+ auto p = new TVector3(x,y,z);
+ detnum.push_back(mg);
+ pos.push_back(p);
energy.push_back(e);
+// cout << "MG" << mg << " X " << x << " Y " << y << " Z " << z << " E " << e << endl;
}
file.close();
}
////////////////////////////////////////////////////////////////////////////////
-void BDC(){
+void minimizer(){
+ flagDraw = 0; //stop from drawing canvas on every run, just minimized one
+
+ gErrorIgnoreLevel = kWarning;
- // Data
+ // Read data in
LoadFile();
- // Minimizer
- auto min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ // Start with beam (0,0,0) and 4.7um 0.5mg/c2 target
+ double parameter[4] = {0.0, 0.0, 0.0, 4.7};
+ devE(parameter);
+
// Function with 4 parameter XYZ and Target thickness
- auto func=ROOT::Math::Functor(&devE,4);
- min->SetFunction(func);
- min->SetPrintLevel(0);
- min->SetPrecision(1e-10);
- // Start value: Center beam (0,0,0) and 4.7um target 0.5mg/cm2 target
- double parameter[6]={0,0,0,4.7};
- devR(parameter);
-
- // min->SetLimitedVariable(0,"AM",parameter[0],1,-90,90);
- min->SetLimitedVariable(0,"X",parameter[0],0.1,-10,10);
- min->SetLimitedVariable(1,"Y",parameter[1],0.1,-10,10);
- min->SetLimitedVariable(2,"Z",parameter[2],0.1,-5,5);
- min->SetLimitedVariable(3,"T",parameter[3],0.1,parameter[2]-0.1*parameter[2],parameter[2]+0.1*parameter[2]);
- min->Minimize();
-
- const double* x = min->X();
- cout << "X =" << x[0]<<endl;
- cout << "Y =" << x[1]<<endl;
- cout << "Z =" << x[2]<<endl;
- cout << "Y =" << x[3]<<endl;
- cout << "Minimum: " << devR(x) << endl;
+ auto func = ROOT::Math::Functor(&devE,4);
+
+ // Minimizer
+ auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
+
+ minim->SetPrintLevel(0);
+ minim->SetPrecision(1e-10);
+
+ // Set minimizer function
+ minim->SetFunction(func);
+
+ // Assign variable limits
+ minim->SetLimitedVariable(0,"X",parameter[0],0.01,-10,10);
+ minim->SetLimitedVariable(1,"Y",parameter[1],0.01,-10,10);
+ minim->SetLimitedVariable(2,"Z",parameter[2],0.01,-5,5);
+ minim->SetLimitedVariable(3,"T",parameter[3],0.01,4.7-3.0,4.7+3.0);
+// parameter[2]-0.1*parameter[2],
+// parameter[2]+0.1*parameter[2]);
+
+ // Shrink it, babeyyy
+ minim->Minimize();
+
+ // Draw minimal value
+ flagDraw = 1;
+
+ // Pull values from minimizer
+ const double* x = minim->X();
+ cout << "========================================" << endl;
+ cout << "\t\tX =" << x[0] << endl;
+ cout << "\t\tY =" << x[1] << endl;
+ cout << "\t\tZ =" << x[2] << endl;
+ cout << "\t\tT =" << x[3] << endl;
+ cout << "Minimum: " << devE(x) << endl;
+ cout << "========================================" << endl;
}
--
GitLab