From 2094bde7d88b6bdb575bedd377faf01992a271eb Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Thu, 3 Feb 2022 11:58:56 +0000
Subject: [PATCH] *e793s - Adding diff. cross sect. analysis
---
Projects/e793s/macro/CS2.h | 313 +++++++++
Projects/e793s/macro/DrawPlots.C | 8 +-
Projects/e793s/macro/DrawPlots.h | 178 ++++-
.../e793s/macro/GateThetaLabHistograms.root | Bin 0 -> 6083 bytes
Projects/e793s/macro/GausFit.h | 658 ++++++++++++++++++
Projects/e793s/macro/KnownPeakFitter.h | 194 +++++-
Projects/e793s/macro/ThreeBodyBreakup.h | 39 ++
7 files changed, 1369 insertions(+), 21 deletions(-)
create mode 100644 Projects/e793s/macro/CS2.h
create mode 100644 Projects/e793s/macro/GateThetaLabHistograms.root
create mode 100755 Projects/e793s/macro/GausFit.h
mode change 100644 => 100755 Projects/e793s/macro/KnownPeakFitter.h
create mode 100644 Projects/e793s/macro/ThreeBodyBreakup.h
diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
new file mode 100644
index 000000000..7c1a0c00b
--- /dev/null
+++ b/Projects/e793s/macro/CS2.h
@@ -0,0 +1,313 @@
+//TGraph* ExpDiffCross(double Energy);
+double* ExpDiffCross(double Energy);
+TH1F* GateThetaLab_RtrnHist(double minTheta, double maxTheta, double binsize);
+TH1F* PullThetaLabHist(int i, double minTheta, double gatesize);
+void Scale(TGraph* g , TGraphErrors* ex);
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
+double ToMininize(const double* parameter);
+//TGraph* FindNormalisation(TGraph* cs1, TH1* hexp);
+TGraph* FindNormalisation(TGraph* cs1, double* hexp);
+
+TGraph* g1;
+//TH1* current;
+double* current;
+TH1F* SolidAngle;
+//auto diffcrossexp = new TGraph("diffcrossexp","diffcrossexp",180,0,180);
+
+////////////////////////////////////////////////////////////////////////////////
+void CS(double Energy, double Spin, double spdf, double angmom, double nodes){
+ auto file = new TFile("../SolidAngle_HistFile_06Dec_47Kdp.root");
+ SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_Lab_MG");
+
+ // p3/5 -> spdf = 1, angmom = 1.5
+
+ // J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
+ double J0 = 0.5;
+
+ // Solid Angle
+ TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
+ SolidAngle->Draw();
+
+ //
+ TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
+
+ // some work around here.
+ // TWOFNR is calling the ADWA calculation
+ auto diffcross = TWOFNR(Energy, J0, Spin, nodes, spdf, angmom);
+ diffcross->Draw();
+ // need to define hexp, best to have to have some function that take Ex and Edc
+ // min and max, make the angular distrib, divide by solid angle and return the
+ // normalised CS.
+ //
+ //
+
+ //TCanvas* c_ExpDiffCross = new TCanvas("c_ExpDiffCross","c_ExpDiffCross",1000,1000);
+ auto temp = ExpDiffCross(0.143);
+ //TCanvas* c_ExpDiffCross = new TCanvas("c_ExpDiffCross","c_ExpDiffCross",1000,1000);
+ //temp->Draw();
+
+// TCanvas* c_Final = new TCanvas("c_Final","c_Final",1000,1000);
+// temp->Divide(SolidAngle);
+// temp->Draw();
+
+ //GateThetaLab(105,110,0.1);FitKnownPeaks(Ex_ThetaLabGate)
+ //
+ FindNormalisation(diffcross,temp);
+}
+
+//TGraph* ExpDiffCross(double Energy){
+double* ExpDiffCross(double Energy){
+ // Implement some kind of energy selection later, for now, hard code that we are selecting 0.143, array index 1
+ int numbins = 10;
+ double x[numbins], y[numbins];
+ //for(int i=0; i<10;i++){
+ for(int i=0; i<numbins;i++){
+ double* peakAreas;
+ double bin = 5.;
+ double min = 105. + (i*bin);
+ double max = min + bin;
+ cout << "min: " << min << " max: " << max << endl;
+
+ TH1F* gate = PullThetaLabHist(i,105.,5.);
+ peakAreas = FitKnownPeaks_RtrnArry(gate);
+ cout << "area of 0.143 = " << peakAreas[1] << endl;
+// delete gate;
+ //diffcrossexp->Fill(min+(bin/2.),peakAreas[1]);
+
+// x[i] = min+(bin/2.);
+ y[i] = peakAreas[1];
+ }
+
+// TGraph* diffcrossexp = new TGraph(numbins,x,y);
+
+// return diffcrossexp;
+ return y;
+}
+
+/*
+TH1F* GateThetaLab_RtrnHist(double minTheta, double maxTheta, double binsize){
+ string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > "
+ + to_string(minTheta)
+ + " && ThetaLab < "
+ + to_string(maxTheta);
+
+ 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)";
+
+ TCanvas *cEx_ThetaLabGate = new TCanvas("cEx_ThetaLabGate","cEx_ThetaLabGate",1000,1000);
+ chain->Draw(draw.c_str(),gating.c_str(),"colz");
+ TH1F* Ex_ThetaLabGate = (TH1F*) gDirectory->Get("Ex_ThetaLabGate");
+ Ex_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
+ Ex_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
+ Ex_ThetaLabGate->SetTitle(title.c_str());
+
+ return Ex_ThetaLabGate;
+}
+*/
+
+TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
+ TFile* file = new TFile("GateThetaLabHistograms.root","READ");
+
+ string histname = "cThetaLabGate_" + to_string(minTheta+(i*gatesize)) + "-" + to_string(minTheta+((i+1)*gatesize));
+ TList *list = (TList*)file->Get("GateThetaLabHistograms");
+ TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+ return hist;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+/*
+ double Chi2(TGraph* g , TH1* h){
+ double Chi2 = 0 ;
+ double chi;
+
+ for(int i = 1 ; i < h->GetNbinsX() ; i++){
+ if(h->GetBinContent(i)>0) {
+ chi = (h->GetBinContent(i) - g->Eval(h->GetBinCenter(i)) ) / (h->GetBinError(i));
+ Chi2 +=chi*chi;
+ }
+ }
+
+ return Chi2;
+}
+*/
+
+double Chi2(TGraph* g , double* h){
+ double Chi2 = 0 ;
+ double chi;
+
+ for(int i = 1 ; i < 10 ; i++){
+ if(h[i]>0) {
+ chi = (h[i] - g->Eval(h[i]) ); // DIVIDE BY ERROR!!! /(h->GetBinError(i));
+ Chi2 +=chi*chi;
+ }
+ }
+
+ return Chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double ToMininize(const double* parameter){
+static int f = 0 ;
+ TGraph* g = new TGraph();
+ double* X = g1->GetX();
+ double* Y = g1->GetY();
+ for(int i = 0 ; i < g1->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i]);
+ }
+ double chi2 = Chi2(g,current);
+ return chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Scale(TGraph* g , TGraphErrors* ex){
+ double scale;
+ double mean = 0 ;
+ double* eX = ex->GetX();
+ double* eY = ex->GetY();
+ double totalW = 0;
+ double W = 0;
+ for(int i = 0 ; i < ex->GetN() ; i++){
+ if(eY[i]>1 && eY[i] <1e4){
+ // Incremental Error weighted average
+ W = 1./ex->GetErrorY(i);
+ scale = eY[i]/g->Eval(eX[i]);
+ totalW +=W;
+ mean = mean + W*(scale - mean)/(totalW);
+ }
+ }
+
+ double* x = g->GetX();
+ double* y = g->GetY();
+ for(unsigned int i = 0 ; i < g->GetN() ; i++)
+ g->SetPoint(i,x[i],y[i]*mean);
+}
+////////////////////////////////////////////////////////////////////////////////
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
+ cout << " ==================================== " << endl;
+ char origDirchar[200];
+ getcwd(origDirchar,200);
+ 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
+ system("pwd");
+ cout << " ==================================== " << endl;
+
+ double BeamEnergy = 7.7;
+ double QValue = 2.274 - E;
+
+ std::ofstream Front_Input("in.front");
+ Front_Input << "jjj" << std::endl;
+ Front_Input << "pipo" << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << BeamEnergy << std::endl;
+ Front_Input << 47 << " " << 19 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << "0 0 0" << std::endl;
+ Front_Input << l << " " << j << std::endl;
+ Front_Input << n << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << QValue << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << J0 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 5 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << J << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 6 << std::endl;
+ Front_Input << 4 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1.25 << " " << 0.65 << std::endl;
+ Front_Input << 6.0 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << 0 << std::endl;
+
+ Front_Input.close();
+
+ cout << "Filled front input file." << endl;
+ cout << "Executing front20..." << endl;
+ system("(exec ~/Programs/TWOFNR/front20 < in.front > /dev/null)");
+ cout << "Executing twofnr20..." << endl;
+ system("(exec ~/Programs/TWOFNR/twofnr20 < in.twofnr > /dev/null)");
+
+ cout << "twofnr20 complete." << endl;
+
+
+ TGraph* CS = new TGraph("24.jjj");
+
+ cout << " ==================================== " << endl;
+ cout << "Current directory " << twofnrDir << endl;
+
+ cout << "Moving to directory " << origDir << endl;
+// string line = "cd ";
+// line += twofnrDir;
+// system(line.c_str());
+ chdir(origDir.c_str());
+ //Check
+ system("pwd");
+ cout << " ==================================== " << endl;
+
+
+ return CS;
+ //return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//TGraph* FindNormalisation(TGraph* cs1, TH1* hexp){
+TGraph* FindNormalisation(TGraph* cs1, double* hexp){
+ // Set global variable
+ g1 = cs1;
+ current = hexp;
+
+ const char* minName ="Minuit";
+ const char* algoName="Migrad";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+ // Number of parameter
+ int mysize = 2;
+
+ // create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ parameter[i] = 0.5;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ min->SetLimitedVariable(i,name,parameter[i],0.1,0,1000);
+ }
+
+ // do the minimization
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] <<"("<<err[i] << ") ";
+ }
+ cout << endl;
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = cs1->GetX();
+ double* Y = cs1->GetY();
+ for(int i = 0 ; i < cs1->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],xs[0]*Y[i]);
+ }
+ return g;
+}
+
diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C
index 548fc50e9..314a6b4dd 100755
--- a/Projects/e793s/macro/DrawPlots.C
+++ b/Projects/e793s/macro/DrawPlots.C
@@ -2,6 +2,8 @@
#include "KnownPeakFitter.h"
#include "DrawPlots.h"
+#include "CS2.h"
+#include "ThreeBodyBreakup.h"
/* USE THIS SPACE TO TEST NEW FEATURES */
void thickness(){
@@ -514,7 +516,7 @@ void ForPoster_DiffCrossSec(){
/* MAIN FUNCTION */
void DrawPlots(){
- gStyle->SetOptStat("nei");
+ gStyle->SetOptStat("nemMrRi");
LoadChainNP();
cout << "==========================================" << endl;
@@ -525,7 +527,10 @@ cout << " 2D Matrices " << endl;
cout << ""<< endl;
cout << " Ungated histograms " << endl;
cout << "\t- Draw_1DParticle() "<< endl;
+ cout << "\t- Draw_1DParticle_MUST2() "<< endl;
cout << "\t- Draw_1DGamma() "<< endl;
+ cout << "\t- Draw_1DGamma_MG() "<< endl;
+ cout << "\t- Draw_1DGamma_MM() "<< endl;
cout << ""<< endl;
cout << " Gated histograms " << endl;
cout << "\t- GateParticle_SeeGamma(particle, width) "<< endl;
@@ -550,6 +555,7 @@ cout << " 2D Matrices " << 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 << ""<< endl;
cout << "==========================================" << endl;
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 6d02e9ed6..055c48c4d 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -87,12 +87,21 @@ void DrawParticleStates(TCanvas* canvas){
TLine *l2907 = new TLine(2.907, 0.0, 2.907, max);
l2907->SetLineStyle(kDotted);
l2907->Draw("same");
- TLine *l3600 = new TLine(3.600, 0.0, 3.600, max);
- l3600->SetLineStyle(kDotted);
- l3600->Draw("same");
+ TLine *l3200 = new TLine(3.2, 0.0, 3.2, max);
+ l3200->SetLineStyle(kDotted);
+ l3200->Draw("same");
+ TLine *l3605 = new TLine(3.605, 0.0, 3.605, max);
+ l3605->SetLineStyle(kDotted);
+ l3605->Draw("same");
TLine *l3800 = new TLine(3.792, 0.0, 3.792, max);
l3800->SetLineStyle(kDotted);
l3800->Draw("same");
+ TLine *l4100 = new TLine(4.1, 0.0, 4.1, max);
+ l4100->SetLineStyle(kDotted);
+ l4100->Draw("same");
+ TLine *l4400 = new TLine(4.4, 0.0, 4.4, max);
+ l4400->SetLineStyle(kDotted);
+ l4400->Draw("same");
}
void plot_kine(NPL::Reaction r, double Ex,Color_t c,int w, int s){
@@ -132,9 +141,27 @@ void Draw_1DGamma(){
Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
}
+void Draw_1DGamma_MG(){
+ TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
+ gStyle->SetOptStat(0);
+ chain->Draw("AddBack_EDC>>Eg(5000,0,5)","abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8");
+ TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
+ Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
+}
+
+void Draw_1DGamma_MM(){
+ TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
+ gStyle->SetOptStat(0);
+ chain->Draw("AddBack_EDC>>Eg(5000,0,5)","abs(T_MUGAST_VAMOS-2777)<600 && MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<5");
+ TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
+ Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
+}
+
void Draw_1DParticle(){
TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000);
- chain->Draw("Ex>>Ep(120,-1,5)",
+ chain->Draw("Ex>>Ep(160,-1,7)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"");
TH1F* Ep = (TH1F*) gDirectory->Get("Ep");
@@ -145,6 +172,17 @@ void Draw_1DParticle(){
DrawParticleStates(cEx);
}
+void Draw_1DParticleMUST2(){
+ TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000);
+ chain->Draw("Ex>>Ep(160,-1,7)",
+ "abs(T_MUGAST_VAMOS-2777)<600 && MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<4",
+ "");
+ TH1F* Ep = (TH1F*) gDirectory->Get("Ep");
+// Ep->SetTitle("Ex");
+ Ep->GetXaxis()->SetTitle("Ex [MeV]");
+ Ep->GetYaxis()->SetTitle("Counts / 0.05 MeV");
+}
+
void Draw_2DParticleGamma(){
TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000);
chain->Draw("AddBack_EDC:Ex>>ExEg(140,-1,7,2500,0,5)",
@@ -176,8 +214,40 @@ void Draw_2DGammaGamma(){
XeqY->Draw("same");
}
+/*
+void gg(){
+ // Example filling of GG matrix
+ unsigned int size = AddBack_EDC.size();
+ double e1,e2;
+ auto h=new TH2("gg","gg",1000,0,10,1000,0,10);
+ for(unsigned int i = 0 ; i < size-1 ; i++){
+ e1=AddBack_EDC[i] ; e2 = AddBack_EDC[i+1];
+ // Folding of the matrix, always fill big first
+ if (e1>e2) {h->Fill(e1,e2);}
+ else {h->Fill(e2,e1);}
+ }
+}
+*/
+
+void Draw_2DGammaGamma_TimeGated(){
+ TCanvas *cEgEg = new TCanvas("cEgEg","cEgEg",1000,1000);
+ chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(999,0.005,5,999,0.005,5)","abs(T_MUGAST_VAMOS-2777)<600","colz");
+ TH1F* EgEg = (TH1F*) gDirectory->Get("EgEg");
+ chain->Draw("AddBack_EDC2:AddBack_EDC>>EgEg2(999,0.005,5,999,0.005,5)","abs(T_MUGAST_VAMOS-2777)<600","colz");
+ TH1F* EgEg2 = (TH1F*) gDirectory->Get("EgEg2");
+ EgEg->Add(EgEg2,1);
+ EgEg->SetTitle("Egamma-Egamma");
+ EgEg->GetXaxis()->SetTitle("Eg [MeV]");
+ EgEg->GetYaxis()->SetTitle("Eg [MeV]");
+ EgEg->Draw("colz");
+ TLine *XeqY = new TLine(0,0,5,5);
+ XeqY->SetLineColor(kRed);
+ XeqY->SetLineStyle(kDashed);
+ XeqY->Draw("same");
+}
+
void GateGamma_SeeParticle(double gamma, double width, double binsize){
- string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-"
+ string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
//string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber!=3 && abs(AddBack_EDC-"
+ to_string(gamma)
+ ")<"
@@ -207,11 +277,11 @@ cout << " NO MG3 IN THIS ONE!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
}
void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){
- string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-"
+ string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ to_string(gamma)
+ ")<"
+ to_string(width);
- string bggate = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-"
+ string bggate = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ to_string(bg)
+ ")<"
+ to_string(width);
@@ -241,7 +311,48 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){
DrawParticleStates(cEx_Gate);
}
+void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg, double widthbg){
+ string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ + to_string(gamma)
+ + ")<"
+ + to_string(width);
+ string bggate = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ + to_string(bg)
+ + ")<"
+ + to_string(widthbg);
+
+ double ratio = width/widthbg;
+
+ string title = "Gate: "+to_string(gamma-width)+" to "+to_string(gamma+width)+"."
+ + " 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(120,-1,5)",gating.c_str(),"");
+ TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate");
+ ExGate->GetXaxis()->SetTitle("Ex [MeV]");
+ ExGate->GetYaxis()->SetTitle("Counts / 0.10 MeV");
+ //ExGate->GetYaxis()->SetTitle("Counts / 0.05 MeV");
+ ExGate->SetLineColor(kGreen);
+ ExGate->SetFillColor(kGreen);
+ ExGate->SetFillStyle(3154);
+ ExGate->SetTitle(title.c_str());
+
+ chain->Draw("Ex>>ExBG(80,-1,7)",bggate.c_str(),"same");
+ //chain->Draw("Ex>>ExBG(120,-1,5)",bggate.c_str(),"same");
+ TH1F* ExBG = (TH1F*) gDirectory->Get("ExBG");
+ ExBG->Scale(ratio);
+ ExBG->SetLineColor(kRed);
+ ExBG->SetFillColor(kRed);
+ ExBG->SetFillStyle(3345);
+ //ExBG->Draw("BSAME");
+
+ DrawParticleStates(cEx_Gate);
+}
+
void GateParticle_SeeGamma(double particle, double width){
+ gStyle->SetOptStat("nemMrRi");
+
string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber<8 && Mugast.TelescopeNumber>0 && abs(Ex-"
+ to_string(particle)
+ ")<"
@@ -261,9 +372,25 @@ void GateParticle_SeeGamma(double particle, double width){
cEg_Gate->Update();
TLine *limit = new TLine(particle, 0.0, particle, cEg_Gate->GetUymax());
limit->SetLineColor(kRed);
+ TLine *sub0143 = new TLine(particle-0.143, 0.0, particle-0.143, cEg_Gate->GetUymax());
+ sub0143->SetLineColor(kBlack); sub0143->SetLineStyle(kDotted);
+ TLine *sub0279 = new TLine(particle-0.279, 0.0, particle-0.279, cEg_Gate->GetUymax());
+ sub0279->SetLineColor(kBlack); sub0279->SetLineStyle(kDotted);
+ TLine *sub0728 = new TLine(particle-0.728, 0.0, particle-0.728, cEg_Gate->GetUymax());
+ sub0728->SetLineColor(kBlack); sub0728->SetLineStyle(kDotted);
+
+ TLine *g0143 = new TLine(0.143, 0.0, 0.143, cEg_Gate->GetUymax());
+ g0143->SetLineColor(kBlack); g0143->SetLineStyle(kDotted);
+ TLine *g0279 = new TLine(0.279, 0.0, 0.279, cEg_Gate->GetUymax());
+ g0279->SetLineColor(kBlack); g0279->SetLineStyle(kDotted);
+ TLine *g0449 = new TLine(0.449, 0.0, 0.449, cEg_Gate->GetUymax());
+ g0449->SetLineColor(kBlack); g0449->SetLineStyle(kDotted);
EgGate->Draw();
limit->Draw();
+ sub0143->Draw();
+ sub0279->Draw();
+ sub0728->Draw();
}
void GateParticle_SeeGamma_WithBG(double particle, double width, double bg){
@@ -859,6 +986,8 @@ void ELabThetaLab(){
plot_kine(Tidp, 0.000, kBlack, 2, 1);
plot_kine(Tidp, 5.652, kBlack, 2, 6); //strongest populated state according to PDBarnes(1965)
+
+ plot_kine(K12C12C, 0.000, kBlack, 2, 1);
}
void XYMust2(){
@@ -1028,7 +1157,7 @@ void ElasticsGate(double EMin, double EMax){
}
void GateThetaCM(double minTheta, double maxTheta, double binsize){
- string gating = "abs(T_MUGAST_VAMOS-2777)<600 && ThetaCM > "
+ string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaCM > "
+ to_string(minTheta)
+ " && ThetaCM < "
+ to_string(maxTheta);
@@ -1048,7 +1177,7 @@ void GateThetaCM(double minTheta, double maxTheta, double binsize){
}
void GateThetaLab(double minTheta, double maxTheta, double binsize){
- string gating = "abs(T_MUGAST_VAMOS-2777)<600 && ThetaLab > "
+ string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > "
+ to_string(minTheta)
+ " && ThetaLab < "
+ to_string(maxTheta);
@@ -1067,4 +1196,33 @@ void GateThetaLab(double minTheta, double maxTheta, double binsize){
DrawParticleStates(cEx_ThetaLabGate);
}
-
+void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize){
+ string core = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > ";
+ string ytitle = "Counts / " + to_string(binsize) + " MeV";
+ double gatesize = (finishTheta-startTheta)/numGates;
+ TList* list = new TList();
+
+ for (int i=0; i<numGates; i++){
+ double minTheta = startTheta + (i * gatesize);
+ string title = to_string(minTheta)+" < ThetaLab < "+to_string(minTheta+gatesize);
+ string gating = core
+ + to_string(minTheta)
+ + " && 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)";
+
+ TCanvas *cEx_ThetaLabGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000);
+ chain->Draw(draw.c_str(),gating.c_str(),"colz");
+ TH1F* Ex_ThetaLabGate = (TH1F*) gDirectory->Get(histname.c_str());
+ Ex_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
+ Ex_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
+ Ex_ThetaLabGate->SetTitle(title.c_str());
+ list->Add(Ex_ThetaLabGate);
+ delete cEx_ThetaLabGate;
+ }
+
+ TFile* file = new TFile("GateThetaLabHistograms.root","RECREATE");
+ list->Write("GateThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
diff --git a/Projects/e793s/macro/GateThetaLabHistograms.root b/Projects/e793s/macro/GateThetaLabHistograms.root
new file mode 100644
index 0000000000000000000000000000000000000000..a527f763565d61a3a6c743846c0ce4272a28b2c2
GIT binary patch
literal 6083
zcmb`LWl$X5o5cqWWN>#2g9i->7CgZrA?N_XW(H?)hv0$W65L5}3+{stkl+r%HMrX(
zu)JCOvbF#Huw8Y!yXy9R`qy=Db)Dxp!Qn0dz$pp<06+i$%I*jGJlrch2<<^&%!d;c
z0C>C#0N@b;&}K68dV*((_Y0YTXhj|EM-QX_ry~jQPa%EV&|FKv%fFQ$-Ua|*X}z|x
zh0-dRx<Iuop)RJXre=y(&Mt5ZCsR9TUZ;Q9{?|bO);|>=^zaRUhzEfm^tUnqKymV~
zaB;W)p#sJJrBeP&b^FhJEG^Xs2#WuiibNjnYG&)nWNT$_1BE<XoL#`GIElD`SP;Mi
z?I7)gt@b80(pW{nT@SYC;MTQqSg$VqsR2wkjDGZT8{J4pPbSY59lahX_EYR>!I^G!
zq{(8~@Q}zFPSaW!j+`!aJXpL@+T5Wt1OiLzZThZC^K&w>m-@~`N^->c(EBhZ^2>SU
z-B^Op@Y$Wm665yyP1bfsI@tN*xzLy)M1-wB+&o^Gl;A$u9O~R{CI>~NnyB@HNo)Fk
z23DR$Pc79s$tlb<C%Nt=m+kE~)-a2@t82#GhK54gQ#fk-($|Kw%{oJ@%nS0?3<QNY
zBh;|cZ~Vtum202kMbdGX38X{ZVb)`8KY9*b7u}Z`4JPRk=2pn^(q9!qpZg5Y8C}y{
z&wy}^v=6^ZQq}IoZWC2mp<Xa*b1)&H^H$Y_tap2jFsz1GA|JBK*x!|mT3$F`tU#^L
z^d>mdoUg6h1`%&BhLma60^8EZ>-iD1ocD-s8t2NxwKGtQ#~ZiK>|eiuNk*4plpB=}
z^9<?a{jmM9G)<J>RG0J7b6&OHicFsY<LsuVsB<nyp<(N<*N~_to=>%wz*h{v$A%t%
zC=n=TeM4-^s}btOZu_FoMyRhc^yZD6FZaQ<5&~%~=}M!5sP)7~5#YUB0*nAY9B5G(
zU}5}_SwsY$sw;Ne_QxbU?xufvO5oa&vCDpTA?SwY!zdVV(()~0#X1GIClx{FhuP4X
z?9Y8s87z8?zasl*lWC@oN=KUeCxPIM+$h@kJ%`Vuq>l2AKc^N|giTv`{bE6AFNRci
zkS>N{K2L^Nj0-p9WePS=U3uLM+e6K52d_BLI!$Btx#|j&7KL0UoHYu$vK#)`NQ!Mh
z_R7R4v-kTqycOagHKwB4QNdL9#f7mE$e=E9g%stTiETT~xC}+z+gC}6ZB+fUHMyTh
zF9~D_C#C7!gph_rxEcQJxxn}ikI#P6`dc<1DPo1OtfU6?;<ncvnE3A^9QZc(R`(o8
z!)hBa`w)cbp(B93Uw6r->C#(`Skv8uKAz0eK_DjYL}Y-X!KEOkZ-;TOdFx@qo+@-_
z8PceqA6H2V@$J3{*{9`fQf-fmpiky|;S)J{b^LlR{To7D$NIjtb?rt$G6j0ywx;_1
z^bvzgwX05MgoYdB)R%{!MK}EkSuaekbWVgxD=twv((f~R>c|4HUsle8@KVGzZu^-u
zTUr#<kdR3iF5VMx(-!nYZtKD(H}KZEgM9Bv9(9Dqr#cA3z-?CTKx<7dx!*`*sfmrr
zZ)3A7e+ZIx1$eA&0>68CjN^mFWy%s(Ti6X{8c%eh_aT-9vN-#+X{=!OFKrUDqViw5
z3bsz~JPkpa%S=;!_IV=n?+Z2X9}D&KSb%FefD!`lySy)%ygQlg8LW@KxJ>BD!=SBc
z5T6q(9qwfp7rb$+(Tc8$lp_rDi47kz6JLFns}g$7j$yjLj0ItXONG^z95a_`r&;lw
zc|e5V47Q=BCcH6CwK3KA!Ykiaja3%J9RO}Jz)68+N1n6Y--5~g4a9oO3rM}V1)A9(
zbQ-2d!I=CyP_9kd?&d2_b3uN<{se-|o5BNcXMP-k;TE_b=baNn&@FlVZeTp=le0$P
z@q7IQzIU;x4R})MTgb9v%?17ovTFDwm~FW~1d2rXB1?P-Nk&NsDv2C=FY(V#3sKpB
zMhOnQJ(r{;z4`~H^<8--=i&yEDHSDvxv9)r#Sy{pe1n~Y=~(1n&eF!7n4sAm1uCle
zfPQ#$B6H}|{7-H(m37XrbAewnOUL1ZL$hAHER57Gj~K91`GP*Xu6tTre(&VJfaR1?
zlQr$B9J<!FJ>u&mxTf{NtN!5mqrdff*AF@(_z~JZ;UD7%x2QtqlIj^SSScro;Bpx7
zdvd!f_I}gfx-E_y2HsyrVTSt{S*wvH8uA)%n;|#9<1^;ZU1_(!7PVbY)I@o$T})Lj
z&YEfX(pCO8Jj++FA1-k|8HiSGH_+q4K--WlT}=l`c;gIN9lJwh^-zWA3DNrN_iUqZ
zrOzt#e*GP#!7vY!`5b!rOT!R|%8^_w<)+lzWcUO<8efOr4XzE2yjL4PEle6OjVU#o
zQEA|NGsMD^dWOd!s15sCl+$Z$QC9!GM_eM(^%Wfb2X_qB=Tbuhz_9Y1)4?Q`=sCme
zvEM8p9RnjnYgiN8b8=wa)aBGQunca=^RQS0(9<YGJ$DIIwDwK=>sTd`^RpJKZ}h4e
z)LA~-m4Ck1w=BOulm3`cUbAcHQ~f02h_8#NM;y=iF-dbVo~b@oxdTP`z98F6VVm)~
z>Fc9+#%e8)*5^b`y)s`KNd3@hO^Bp!*e*C#lD@{6T`zj!wTR=O?v$2siExFoA<jlE
z^xV$hp9DMJigfPgoBE%eTBgWnR>O`M${;h(`cq!7&}i-}_{nJK$J(r?X*Fk>X2K?Q
z6QRYkOtuMbqGcUQUuoYA1D=1zo9S5Xs)4gAZtNvyaBwskwhxEcMSgcL8@!rAM$*;8
zWr|-)oV{tj|5WB%>Uf(>>5lmIUq3%ojgR71mjfgbMuV#|RHE34pmIy_EYwTLVx+pq
zyFlSA@H|^}m%uZP>0GK|BC-A``P)Ng8osz&RYC3f#4ae{gwTw-lCK4BdW8|M0zI`E
z$!GZ9aFDCUnvySxr4s;|ZDX2y&Jge9_tH>-0>XQH?gjp4V^zb<^6HV=TYV50{!KSF
z)(3j-3%rdk0jKOAMH_pV`pf(cvA$10$pbH(zOA2bhQ^E~Z(~R$wNaxM*p&~sml^Dh
zB`Rmy-!7bX!`pA1*$phFnyw)VY{u^)MpU~%Goon_<XsMegOFLsk44!yFql+-o$%fE
zqVxxQ!E61&&rEafv=VrB{y%?zmw6xiNE$M&vG@pbC`C*qrA{k2-Cm~$;pH2+6#`>(
zJM9Otj#e%gjE^vA%};P-2fr_Z9idH|)ZFa#n{-}`g^bJ2DK%Xa%?+-WCEB&~PUoH$
zRP2RqzLnVyu~=4&?VsTkef!3#$tM$tq8e#8*KX$n;bM20E+=WC)+!5ir+9(;ap{$w
zaGR}+|JrGMnq5MV^@q47LBF*H^#+P0^~dgm6sf#1E56Vdf9lb;rtnUx#~)Rq4)?TS
z{s<?hptcU(*p(10zcc9Z8Qk_4J#P^mZcOM98LF+_i-2=el!X|ERoGY-su_o+>VDxy
z>FSf~cbn3qVMUvTXH`}R){WCLfC*3cM!$FuR4uu2-DBAAeJ%RyTOVv5tIJ=TDE*Jk
z`^U}xr_IBA=i&r4wSzh-*~8%fv$;b80WPc$=57kS0H@QXzL&V14?ma}apYAq@WN?h
zRUe;K&$U+mKI>Q(T;Qk;r`gV$Tcq}mpQE;bJ+Jr9MFoP3Ng8x=N_^M&R-%*@anKV#
zvYZ4GDdMnttY8Ir{36x7_Fc}PAF-XeoX1zBv)TIbBTE!<W;E<hn||KULS7U!iYJ+i
zpE7x}YE-*sOE<4GMWBIFm-e{Dn71dvH*ahm4R>}w;d!mKGa_|6uX}9&=ye<0(Qu-@
zYH)5{*U@6mu!KDOx=9I`05<-~*3T<S%u!OoNJJ_rkoc_zan1au2hNklYA-8TK@qxV
zAS1iu#(V?g^{=@7po)R;xdab?Va*H-?92+Ok<sLVjvHuzA1lzvA!RrhQ>}@7A;t(<
zC&!Z3@r@FhR{oJ17m%>yl_{1bBc*^n%!eC2xaK1M>c!>0A3oU`z8yM-yh}z2sgW%F
zQ={A02{}W2f_NCg_dP$nd&%2e7$-FVR?wPio^;&aSj~}l+TQ_C*bzxXCmwU#&Nx~s
znJX;3wAQqC%C=barJ4&dADT8kd*Z&-z&&_jo=^m4H%Dtz!ycZ&yppD9SW#_utS#N>
zd%SkS#r3JYk)rr3QyI}1Q}iWT>NVkBm+EaUzq7_CESNoCR&`#8kUiN%CBpC4A_d$S
zLfN72m7kjFA%`tO<$z>uR?ZZm3<*`3@C+|{vmxrPid7fikS8&Nz`cDLYBl!L`0h7y
z4vGJi?-(pu!X9xS&f{nCIPfbg(-&htlVfkHAHe%zFwnTqU!1Cn1E>*vrzOVT=|pC6
zV8D(((#_zCagX+M$`)a302<NRJa0-%kq8*ew!Xj2y5T@o$$WA+cySS&XpGBASWw0E
z6*r0IEafV&pnZ)?Cfm48ONaB5Tq#{CNm&w?w+<O}W~ou~ekMlB($IN1r<io@mj@!t
zcG5k#09Er0YO17x?)powgS%C~@&&Ls!Qw^umcR<Z?R@~T>t47(O+R1-^mZeNlNghp
zG{0%#<EB<X{#r;G#mW^f-2r%5`slZ(pel<>NAzYVY*E*CbEX6&4iQ<XDdnFaR#kX?
z0~ja+&)~Rp%~tv7?1uZ;hx@n{TBev7II={KG`6c9U1mr`;o?V1MV;_ugx3W-8<EBd
z-=ik11XL=!on?mttCA}{R?)<%<rsIs4YIXF!L}?5Jt<MT65^d<*CaJd&XH&Kp!bUW
zN!QiBcw2Yf7r`j+>6N?MJx`Y_!Y9j{R%&@4rb$s34X{HYqm)?Np+DS^mBZGe(RI(m
zvy84kGVPiaYuQi3a@3fJWw0$1oJ<&e^Rte3Tvx}I5~(aL`D}K8Iv%7lgPk9j?!<AE
z#b77cdpl=ChHFNhh9#D~P`O;>c(_2pvJ#A3@`8QFz(jLFq^kgXn9UQ@OF1no0<yyE
z*w8dup)t|Dm9LRNx|(Lv;ZSD7`8b+V&z5^f)reK*9ZUFB4$#)F=Yu@`$)A)l{&0p8
z_QB>-eLPlt^!XTKy555XT=y>|J00<#3i}Fv@g4cV9|wu6<!SII8C^-GSqSm_>uyhx
zR=$Z_rF?Zk`1PHcDtp1t7jBTj@)9h0nob`}?H3Mqgp$xZba+Ae>U@1_^;8egh`3e1
z^SeM>A;(<jd5(*1MM%O^-<{rejl0izD|u*>XoI&1ZOZ2t2u@oXWh?ab6t?lwp1fOb
zs>i$+#`as|lB$D-UKkjMImYqYy3nbSuR@w>2a#l}11ndHilQ<_JiQDRdGCl_{f`K^
zaTf{KzLJVKsckF&s!3!Z&2@RH-<8&By6b`6E2Sr(8KL>K-oQy6RiMF4E!|ukHz>Q0
z#q7qaQpWU+n(lpF$X>UV`A?F=&`@sNpBWvSfgzD?G~+D}zR&C?!Fg>3FML%)l{SDh
zmn3(iRiI*|s6qWpKCxOFzC=SQG5+VS$5rvhlL05iLA#?(`9ci`M1xH|bX&SJ!Z>td
zDe`Z5p=nA^S`r9(>|rJA^oX05Oi`P;sL(Ak_YyOO%~7Ms-YY2;c?9wsjGK)8I9^Gs
zVqyigGZWCn#7}v*z)y*_wPq3MsSOlh!FUE*IGHOW^3GJ^&T!1=k7j(y|6EK=jHG`^
zW-bZDt+k28Esq(dwmWx*rZPJEsCoJ+{uiP@N)?_Fc%%!jTy|1<e!GPxQ*Fm#8U<!}
zV-Sclt>|J&%<(qqLDQYbk<$(NxizK4@Y7(TtX&MFy=kM>{qg4@Q=%_64;$?pMpGgU
zD7;9=h=v#Y+h*|`_;dM1{J}6^&yP_x=~rn*>F7|PLiP`*TJ`{%4<aAQ7fYClg2)<h
zA|zgM{Mh0YRVg%0cHcB1B~jRVuJA;n++RDmuT*wDjwfn?geDnzRV%;1^S-pz@>$Fk
z`lVT&&SIgr>%ysZy=~kpzA$rMUpwW<1L<Zz^wjfgYA3mhOHqrUDzT_w37n5N<EKA$
zc`9T%6d3HIdm6yAN4VJ$sa0~=!)GkBe4Et1EXV6P6p|S|&m42O>}z-Qxx=D2A5rBl
z)FackMCoWD?ScbjIrcMFl;|=Z0vOv^Xm?H&jSHOxhBV9Bp_;Q5t@_}@^dY8h9Hlxz
zn_#D6=@_A!ehCI}6g7rt$20YnS^vIT-dh=<AA+hQnZ*-rtTG|sEIw8;6Vx`7T5+Rk
z`zcH0G{a9W;LMOWp{V+#Mr<XAzlaBsXB|iHWOAuZXWKSyyQmpPEASRR+H#9d-KerQ
zmvMjS%w`j9{xL)Te9z`<O_;qkpQQ&M60^7l=5bsi%^nhCSX`&jebADq-%TDeuF#mL
zT+PK}LWXSaxnyN94MswXS2^b8iS?XU@C6+jSn~Qs*U_}z|Iybn&<JI&T?2-*ux3$_
zQ&n`!ve{Eh_7HS4LCEL$GicDSPuAH5gkbm{@BKgIB%kt9H=+u7Hhar-@I+dppV@mm
z#Q|kWf|$QoZRtzpV3cg>?jKu@;v6-qCstIKPaUZ0Bx3|MyTL0})W;W@MU~#oeU=k$
z_Lw(j-PmCxzYFW^SALAJ8&*E-2I$<~r<VnD(rzSxr5q4+KcU=-^{?+Hpm`h>P?R8b
zqxV}#dTFeyxP|B&M!hJn>!kTn5PjN+O|bBKW}uL+>mQ@A>(551$%9|lL;2YO=mhXQ
z*B9f~mLX>4J?Jw7{NeFs;>l<HiFK_l_5C#S069)i8@)A9wVs;kPdaC74-)I)H?bri
zvdzo0F_h9zfaGer_41R!!b5bl>NkA`>f&-Hc8v#hManen2F(LE*@4d3Kj!CFW~blS
zTahK>k#Z&$8S+Z(C&881Tb|JRwW14$rMGik1hHWt)9jf3L{1frw||uxOO>z0**Vjv
zosJ<4(DvrSW$Ht}=g>^Y1`>bEZS~J<mUl8-i8fo1MeABUZJ=+4E)j`4D7HtmL-)9_
z0xU?T>e7wbFx+#?K*>B;1K<MDl244;)n+IG8HEE3LIQ5@n?@Sn+nw3#vBGQ}^L1p*
zOHZG~#UU6vpisR<3ft@X<h5B6<4Yn}+|CKW2U>IGmrYoBOW+T3mj%neUh_{+DDGkI
z@L`hszw)MeGnLcD=zT(5VW&zjX5ttU(CiT;+u#v2Rf!WDm8jRlY?+QjI(xw$`B84?
z)~H5q2P0SY=>`Y};;t;_sU*VF{kdQt+|MIcYt1oX5nmo1eT|%C4ypC*VMK;V^$EOf
zF(5Te6oxci%eln!D*A|^%sV9_7F;X^q`<P14qXPVflGxUx`@Zig3~!u+ZsW*~
z-gRi1M2<7Wm&Z-tC<8XD+N^rcgM8wNEP&GGt#ZSM=j^N?tmqMtk^0#Jk4Cy3*~Q6f
zc=+0}tG|XpkIdV*lX<GxZz{t$2cQUB<-r#tvv)qN{24BJR31eo*Lt1cF-d(|Yi5JD
zzaPEoYFX8$sJ;V{&#SH4aa1p1#)vlNclJFs%QfLD<B2mpl`3#>aS>_hPWhap=HhEn
zaw1xLS=X)|#d;MW{jTpZLTKu?+<)u?Ac$?H;lcSw{0*bk|Mmla!{~>o_-{Y(FT?1E
m#~#{;|F7Kmp_;;9l*j*#@;7~dNTB~U$^!^xx4Up6;NJkJO<lYI
literal 0
HcmV?d00001
diff --git a/Projects/e793s/macro/GausFit.h b/Projects/e793s/macro/GausFit.h
new file mode 100755
index 000000000..0e37d1ec0
--- /dev/null
+++ b/Projects/e793s/macro/GausFit.h
@@ -0,0 +1,658 @@
+////////////////////////////////////////////////////////
+/* FUNCTIONS TO FIT ONE, TWO AND THREE GAUSSIAN PEAKS */
+/* SingleGaus(), DoubleGaus(), TripleGaus() */
+////////////////////////////////////////////////////////
+
+#include "TMath.h"
+#include "math.h"
+#include <cmath>
+#include "stdlib.h"
+
+Double_t pi = 3.14159265358979323846;
+
+////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////
+
+void DoubleGaus(TH1F* hist){
+ bool repeat=true;
+ int repeatInt;
+ double minFit, maxFit, mean1, mean2;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+
+ while (repeat){
+ cout << "====================================================================" << endl;
+ cout << " Input range to fit:" << endl;
+ cout << " Min = ";
+ cin >> minFit;
+ cout << " Max = ";
+ cin >> maxFit;
+ cout << " Peak 1 = ";
+ cin >> mean1;
+ cout << " Peak 2 = ";
+ cin >> mean2;
+
+ //TF1 *g1 = new TF1 ("m1", equation, minFit, maxFit);
+ TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+ //TF1 *g2 = new TF1 ("m1", equation, minFit, maxFit);
+ TF1 *g2 = new TF1 ("m2", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g2->SetLineColor(kBlue);
+ g2->SetLineStyle(2);
+ TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+
+ g1->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ g2->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ bg->SetParNames("Background");
+
+ TF1 *f1 = new TF1("double_gaus",
+ "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2)) + ([3]/([5]*sqrt(2*pi)))*exp(-0.5*pow((x-[4])/[5],2)) + [6]",
+ minFit, maxFit);
+
+ f1->SetParNames("Area*BinWidth 1", "Mean 1", "Sigma 1",
+ "Area*BinWidth 2", "Mean 2", "Sigma 2",
+ "Background");
+ f1->SetLineColor(kBlack);
+
+ g1->SetParameter(0, 100);
+ g1->SetParLimits(0, 0.0, 500.0);
+ g1->SetParameter(1, mean1);
+ g1->SetParLimits(1, mean1-0.5, mean1+0.5);
+ g1->SetParameter(2, 0.13);
+ g1->SetParLimits(2, 0.05, 0.30);
+ g2->SetParameter(0, 100);
+ g2->SetParLimits(0, 0.0, 500.0);
+ g2->SetParameter(1, mean2);
+ g2->SetParLimits(1, mean2-1.0, mean2+1.0);
+ g2->SetParameter(2, 0.13);
+ g2->SetParLimits(2, 0.05, 0.30);
+
+ hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit);
+ hist->Fit(g2, "WWR", "", mean2-5, maxFit);//minFit, maxFit);
+
+ Double_t par[7];
+ g1->GetParameters(&par[0]);
+ g2->GetParameters(&par[3]);
+ bg->GetParameters(&par[6]);
+ f1->SetParameters(par);
+
+ hist->Fit(f1, "WWR", "", minFit, maxFit);
+ hist->Draw();
+
+ Double_t finalPar[7];
+ Double_t finalErr[7];
+ f1->GetParameters(&finalPar[0]);
+ for (int i=0; i<7; i++){finalErr[i] = f1->GetParError(i);}
+ g1->SetParameters(finalPar[0], finalPar[1], finalPar[2]);
+ g2->SetParameters(finalPar[3], finalPar[4], finalPar[5]);
+ bg->SetParameter(0,finalPar[6]);
+
+ g1->Draw("SAME");
+ g2->Draw("SAME");
+ bg->Draw("SAME");
+ f1->Draw("SAME");
+
+
+/* Error propogation:
+ * (Abin) +- deltaAbin, B+-0 (no uncertainty)
+ * A = Abin/B
+ * deltaA/A = deltaAbin/Abin
+ * deltaA = A x deltaAbin/Abin
+ */
+/*
+ cout << "Area Red : " << finalPar[0]/binWidth
+ << " +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+ << "\nArea Blue : " << finalPar[3]/binWidth
+ << " +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3])
+ << endl;
+*/
+
+ cout << fixed << setprecision(5);
+
+ cout << "\033[91m Mean: \t" << finalPar[1]
+ << "\t +- " << finalErr[1]
+ << endl;
+ cout << "\033[91m Sigm: \t" << finalPar[2]
+ << "\t +- " << finalErr[2]
+ << endl;
+ cout << "\033[91m Area: \t" << finalPar[0]/binWidth
+ << "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+ << endl;
+
+ cout << "\033[36m Mean: \t" << finalPar[4]
+ << "\t +- " << finalErr[4]
+ << endl;
+ cout << "\033[36m Sigm: \t" << finalPar[5]
+ << "\t +- " << finalErr[5]
+ << endl;
+ cout << "\033[36m Area: \t" << finalPar[3]/binWidth
+ << "\t +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3])
+ << endl;
+
+
+ TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
+ legend->AddEntry(f1,"Total fit","l");
+ legend->AddEntry(g1,"Peak 1","l");
+ legend->AddEntry(g2,"Peak 2","l");
+ legend->AddEntry(bg,"Background","l");
+ legend->Draw();
+
+ //cGate->Draw("SAME");
+ gPad->Modified();
+ gPad->Update();
+
+ cout << "\33[37m Refit? " << endl;
+ cin >> repeatInt;
+ if(repeatInt!=1){ repeat=false; }
+ }
+
+}
+
+
+////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////
+
+/*
+void SingleGaus(TH1F* hist){
+ bool repeat=true;
+ int repeatInt;
+ double minFit, maxFit, mean;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+
+ while (repeat){
+ cout << "====================================================================" << endl;
+ cout << " Input range to fit:" << endl;
+ cout << " Min = ";
+ cin >> minFit;
+ cout << " Max = ";
+ cin >> maxFit;
+ cout << " Peak = ";
+ cin >> mean;
+
+ TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+
+ TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+
+ g1->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ bg->SetParNames("Background");
+
+ TF1 *f1 = new TF1("single_gaus", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))+[3]",
+ minFit, maxFit);
+ f1->SetParNames("Area*BinWidth 1", "Mean 1", "Sigma 1",
+ "Background");
+ f1->SetLineColor(kBlack);
+
+
+ g1->SetParameter(0, 100);
+ g1->SetParLimits(0, 0.0, 500.0);
+ g1->SetParameter(1, mean);
+ g1->SetParLimits(1, mean-1.0, mean+1.0);
+ g1->SetParameter(2, 0.15);
+ g1->SetParLimits(2, 0.05, 1.00);
+// bg->FixParameter(0, 0);
+ bg->SetParameter(0, 0);
+ bg->SetParLimits(0, 0., 50.);
+
+ hist->Fit(g1, "WWR", "", minFit, maxFit);//maxFit);
+
+ Double_t par[4];
+ g1->GetParameters(&par[0]);
+ bg->GetParameters(&par[3]);
+ f1->SetParameters(par);
+
+ hist->Fit(f1, "WWR", "", minFit, maxFit);
+ hist->Draw();
+
+ Double_t finalPar[4];
+ Double_t finalErr[4];
+ f1->GetParameters(&finalPar[0]);
+ for (int i=0; i<4; i++){finalErr[i] = f1->GetParError(i);}
+ g1->SetParameters(finalPar[0], finalPar[1], finalPar[2]);
+ bg->SetParameter(0,finalPar[3]);
+
+ g1->Draw("SAME");
+ bg->Draw("SAME");
+ f1->Draw("SAME");
+
+ cout << fixed << setprecision(5);
+
+ cout << "\033[91m Mean: \t" << finalPar[1]
+ << "\t +- " << finalErr[1]
+ << endl;
+ cout << "\033[91m Sigm: \t" << finalPar[2]
+ << "\t +- " << finalErr[2]
+ << endl;
+ cout << "\033[91m Area: \t" << finalPar[0]/binWidth
+ << "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+ << endl;
+
+ TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
+ legend->AddEntry(f1,"Total fit","l");
+ legend->AddEntry(g1,"Peak","l");
+ legend->AddEntry(bg,"Background","l");
+ legend->Draw();
+
+ //cGate->Draw("SAME");
+ gPad->Modified();
+ gPad->Update();
+
+ cout << "\033[37m Refit? " << endl;
+ cin >> repeatInt;
+ if(repeatInt!=1){ repeat=false; }
+ }
+
+}
+*/
+////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////
+
+void TripleGaus(TH1F* hist){
+ bool repeat=true;
+ int repeatInt;
+ double minFit, maxFit, mean1, mean2, mean3;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+
+ while (repeat){
+ cout << "====================================================================" << endl;
+ cout << " Input range to fit:" << endl;
+ cout << " Min = ";
+ cin >> minFit;
+ cout << " Max = ";
+ cin >> maxFit;
+ cout << " Peak 1 = ";
+ cin >> mean1;
+ cout << " Peak 2 = ";
+ cin >> mean2;
+ cout << " Peak 3 = ";
+ cin >> mean3;
+
+ TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+ TF1 *g2 = new TF1 ("m2", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g2->SetLineColor(kBlue);
+ g2->SetLineStyle(2);
+ TF1 *g3 = new TF1 ("m3", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g3->SetLineColor(kViolet);
+ g3->SetLineStyle(2);
+
+ TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+
+ g1->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ g2->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ g3->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ bg->SetParNames("Background");
+
+ TF1 *f1 = new TF1("double_gaus", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2)) + ([3]/([5]*sqrt(2*pi)))*exp(-0.5*pow((x-[4])/[5],2)) + ([6]/([8]*sqrt(2*pi)))*exp(-0.5*pow((x-[7])/[8],2)) + [9]" , minFit, maxFit);
+
+ f1->SetParNames("Area*BinWidth 1", "Mean 1", "Sigma 1",
+ "Area*BinWidth 2", "Mean 2", "Sigma 2",
+ "Area*BinWidth 3", "Mean 3", "Sigma 3",
+ "Background");
+ f1->SetLineColor(kBlack);
+
+/*
+ TF1 *g1 = new TF1 ("m1", "gaus", minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+ TF1 *g2 = new TF1 ("m2", "gaus", minFit, maxFit);
+ g2->SetLineColor(kBlue);
+ g2->SetLineStyle(2);
+ TF1 *bg = new TF1 ("bg", "[0]", minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+ TF1 *f1 = new TF1("double_gaus", "gaus(0) + gaus(3) + [6]", minFit, maxFit);
+ f1->SetParNames("Constant 1", "Mean 1", "Sigma 1",
+ "Constant 2", "Mean 2", "Sigma 2",
+ "Background");
+ f1->SetLineColor(kBlack);
+*/
+
+ g1->SetParameter(0, 100);
+ g1->SetParLimits(0, 0.0, 500.0);
+ g1->SetParameter(1, mean1);
+ g1->SetParLimits(1, mean1-1.0, mean1+1.0);
+ g1->SetParameter(2, 0.13);
+ g1->SetParLimits(2, 0.05, 0.30);
+ g2->SetParameter(0, 100);
+ g2->SetParLimits(0, 0.0, 500.0);
+ g2->SetParameter(1, mean2);
+ g2->SetParLimits(1, mean2-1.0, mean2+1.0);
+ g2->SetParameter(2, 0.13);
+ g2->SetParLimits(2, 0.05, 0.30);
+ g3->SetParameter(0, 100);
+ g3->SetParLimits(0, 0.0, 500.0);
+ g3->SetParameter(1, mean3);
+ g3->SetParLimits(1, mean3-1.0, mean3+1.0);
+ g3->SetParameter(2, 0.13);
+ g3->SetParLimits(2, 0.05, 0.30);
+
+
+ hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit);
+ hist->Fit(g2, "WWR", "", mean2-5, maxFit);//minFit, maxFit);
+
+ Double_t par[10];
+ g1->GetParameters(&par[0]);
+ g2->GetParameters(&par[3]);
+ g3->GetParameters(&par[6]);
+ bg->GetParameters(&par[9]);
+ f1->SetParameters(par);
+
+ hist->Fit(f1, "WWR", "", minFit, maxFit);
+ //hist->GetXaxis()->SetTitle("ThetaLab (degrees)");
+ //hist->GetYaxis()->SetTitle("Counts (per 0.5 degrees)");
+ //hist->SetTitle(titleString.c_str());
+ hist->Draw();
+
+ Double_t finalPar[10];
+ Double_t finalErr[10];
+ f1->GetParameters(&finalPar[0]);
+ for (int i=0; i<10; i++){finalErr[i] = f1->GetParError(i);}
+ g1->SetParameters(finalPar[0], finalPar[1], finalPar[2]);
+ g2->SetParameters(finalPar[3], finalPar[4], finalPar[5]);
+ g3->SetParameters(finalPar[6], finalPar[7], finalPar[8]);
+ bg->SetParameter(0,finalPar[9]);
+
+ g1->Draw("SAME");
+ g2->Draw("SAME");
+ g3->Draw("SAME");
+ bg->Draw("SAME");
+ f1->Draw("SAME");
+
+// cout << "Area Red : " << finalPar[0]/binWidth
+// << " +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+// << "\nArea Blue : " << finalPar[3]/binWidth
+// << " +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3])
+// << "\nArea Violet : " << finalPar[6]/binWidth
+// << " +- " << (finalPar[6]/binWidth) * (finalErr[6]/finalPar[6])
+// << endl;
+
+ cout << fixed << setprecision(5);
+
+ cout << "\033[91m Mean: \t" << finalPar[1]
+ << "\t +- " << finalErr[1]
+ << endl;
+ cout << "\033[91m Sigm: \t" << finalPar[2]
+ << "\t +- " << finalErr[2]
+ << endl;
+ cout << "\033[91m Area: \t" << finalPar[0]/binWidth
+ << "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+ << endl;
+
+ cout << "\033[36m Mean: \t" << finalPar[4]
+ << "\t +- " << finalErr[4]
+ << endl;
+ cout << "\033[36m Sigm: \t" << finalPar[5]
+ << "\t +- " << finalErr[5]
+ << endl;
+ cout << "\033[36m Area: \t" << finalPar[3]/binWidth
+ << "\t +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3])
+ << endl;
+
+ cout << "\033[95m Mean: \t" << finalPar[7]
+ << "\t +- " << finalErr[7]
+ << endl;
+ cout << "\033[95m Sigm: \t" << finalPar[8]
+ << "\t +- " << finalErr[8]
+ << endl;
+ cout << "\033[95m Area: \t" << finalPar[6]/binWidth
+ << "\t +- " << (finalPar[6]/binWidth) * (finalErr[6]/finalPar[6])
+ << endl;
+
+
+ TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
+ legend->AddEntry(f1,"Total fit","l");
+ legend->AddEntry(g1,"Peak 1","l");
+ legend->AddEntry(g2,"Peak 2","l");
+ legend->AddEntry(g3,"Peak 3","l");
+ legend->AddEntry(bg,"Background","l");
+ legend->Draw();
+
+ //cGate->Draw("SAME");
+ gPad->Modified();
+ gPad->Update();
+
+ cout << "\33[37m Refit? " << endl;
+ cin >> repeatInt;
+ if(repeatInt!=1){ repeat=false; }
+ }
+
+}
+
+
+////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////
+
+void SingleGausNoBG(TH1F* hist){
+ bool repeat=true;
+ int repeatInt;
+ double minFit, maxFit, mean;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+
+ while (repeat){
+ cout << "====================================================================" << endl;
+ cout << " Input range to fit:" << endl;
+ cout << " Min = ";
+ cin >> minFit;
+ cout << " Max = ";
+ cin >> maxFit;
+ cout << " Peak = ";
+ cin >> mean;
+
+ TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+
+ g1->SetParNames("Area*BinWidth", "Mean", "Sigma");
+
+ TF1 *f1 = new TF1("single_gaus", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ f1->SetParNames("Area*BinWidth 1", "Mean 1", "Sigma 1");
+ f1->SetLineColor(kBlack);
+
+/*
+ TF1 *g1 = new TF1 ("m1", "gaus", minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+ TF1 *g2 = new TF1 ("m2", "gaus", minFit, maxFit);
+ g2->SetLineColor(kBlue);
+ g2->SetLineStyle(2);
+ TF1 *bg = new TF1 ("bg", "[0]", minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+ TF1 *f1 = new TF1("double_gaus", "gaus(0) + gaus(3) + [6]", minFit, maxFit);
+ f1->SetParNames("Constant 1", "Mean 1", "Sigma 1",
+ "Constant 2", "Mean 2", "Sigma 2",
+ "Background");
+ f1->SetLineColor(kBlack);
+*/
+
+ g1->SetParameter(0, 100);
+ g1->SetParLimits(0, 0.0, 500.0);
+ g1->SetParameter(1, mean);
+ g1->SetParLimits(1, mean-1.0, mean+1.0);
+ g1->SetParameter(2, 0.15);
+ g1->SetParLimits(2, 0.05, 1.00);
+
+ hist->Fit(g1, "WWR", "", minFit, maxFit);//maxFit);
+
+ Double_t par[4];
+ g1->GetParameters(&par[0]);
+ f1->SetParameters(par);
+
+ hist->Fit(f1, "WWR", "", minFit, maxFit);
+ hist->Draw();
+
+ Double_t finalPar[4];
+ Double_t finalErr[4];
+ f1->GetParameters(&finalPar[0]);
+ for (int i=0; i<4; i++){finalErr[i] = f1->GetParError(i);}
+ g1->SetParameters(finalPar[0], finalPar[1], finalPar[2]);
+
+ g1->Draw("SAME");
+ f1->Draw("SAME");
+
+ cout << fixed << setprecision(5);
+
+ cout << "\033[91m Mean: \t" << finalPar[1]
+ << "\t +- " << finalErr[1]
+ << endl;
+ cout << "\033[91m Sigm: \t" << finalPar[2]
+ << "\t +- " << finalErr[2]
+ << endl;
+ cout << "\033[91m Area: \t" << finalPar[0]/binWidth
+ << "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+ << endl;
+
+ TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
+ legend->AddEntry(f1,"Total fit","l");
+ legend->AddEntry(g1,"Peak","l");
+ legend->Draw();
+
+ //cGate->Draw("SAME");
+ gPad->Modified();
+ gPad->Update();
+
+ cout << "\033[37m Refit? " << endl;
+ cin >> repeatInt;
+ if(repeatInt!=1){ repeat=false; }
+ }
+
+}
+
+
+
+////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////
+
+void SingleGaus(TH1F* hist, bool isGamma){
+ bool repeat=true;
+ int repeatInt;
+ double minFit, maxFit, mean;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+
+ while (repeat){
+ cout << "====================================================================" << endl;
+ cout << " Input range to fit:" << endl;
+ cout << " Min = ";
+ cin >> minFit;
+ cout << " Max = ";
+ cin >> maxFit;
+ cout << " Peak = ";
+ cin >> mean;
+
+ TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
+ minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+
+ TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+
+ g1->SetParNames("Area*BinWidth", "Mean", "Sigma");
+ bg->SetParNames("Background");
+
+ TF1 *f1 = new TF1("single_gaus", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))+[3]",
+ minFit, maxFit);
+ f1->SetParNames("Area*BinWidth 1", "Mean 1", "Sigma 1",
+ "Background");
+ f1->SetLineColor(kBlack);
+
+ double areaSet, areaMax, meanRange;
+ double sigSet, sigMin, sigMax;
+ double bgMax;
+
+ if(isGamma){
+ areaSet = 1000.; areaMax = 10000.; meanRange = 0.005;
+ sigSet = 0.001; sigMin = 0.0005; sigMax = 0.005;
+ bgMax = 1000.;
+ } else {
+ areaSet = 100.; areaMax = 1000.; meanRange = 1.0;
+ sigSet = 0.15; sigMin = 0.05; sigMax = 0.3;
+ bgMax = 50.;
+ }
+
+ g1->SetParameter(0, areaSet);
+ g1->SetParLimits(0, 0.0, areaMax);
+ g1->SetParameter(1, mean);
+ g1->SetParLimits(1, mean-meanRange, mean+meanRange);
+ g1->SetParameter(2, sigSet);
+ g1->SetParLimits(2, sigMin, sigMax);
+// bg->FixParameter(0, 0);
+ bg->SetParameter(0, 0);
+ bg->SetParLimits(0, 0., bgMax);
+
+ hist->Fit(g1, "WWR", "", minFit, maxFit);//maxFit);
+
+ Double_t par[4];
+ g1->GetParameters(&par[0]);
+ bg->GetParameters(&par[3]);
+ f1->SetParameters(par);
+
+ hist->Fit(f1, "WWR", "", minFit, maxFit);
+ hist->Draw();
+
+ Double_t finalPar[4];
+ Double_t finalErr[4];
+ f1->GetParameters(&finalPar[0]);
+ for (int i=0; i<4; i++){finalErr[i] = f1->GetParError(i);}
+ g1->SetParameters(finalPar[0], finalPar[1], finalPar[2]);
+ bg->SetParameter(0,finalPar[3]);
+
+ g1->Draw("SAME");
+ bg->Draw("SAME");
+ f1->Draw("SAME");
+
+ cout << fixed << setprecision(5);
+
+ cout << "\033[91m Mean: \t" << finalPar[1]
+ << "\t +- " << finalErr[1]
+ << endl;
+ cout << "\033[91m Sigm: \t" << finalPar[2]
+ << "\t +- " << finalErr[2]
+ << endl;
+ cout << "\033[91m Area: \t" << finalPar[0]/binWidth
+ << "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
+ << endl;
+
+ TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
+ legend->AddEntry(f1,"Total fit","l");
+ legend->AddEntry(g1,"Peak","l");
+ legend->AddEntry(bg,"Background","l");
+ legend->Draw();
+
+ //cGate->Draw("SAME");
+ gPad->Modified();
+ gPad->Update();
+
+ cout << "\033[37m Refit? " << endl;
+ cin >> repeatInt;
+ if(repeatInt!=1){ repeat=false; }
+ }
+
+}
+
+
+
+
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
old mode 100644
new mode 100755
index b3e0b0c20..9b7918a9e
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -21,14 +21,19 @@ void FitKnownPeaks(TH1F* hist){
* 1.981
* 2.410
* 2.910
- * 3.600
+ * 3.605
* 3.792
* ...
+ * NEW:
+ * 3.2
+ * 4.1
+ * 4.4
+ *
*/
//Assign peak centroids
- const int numPeaks = 11;
- array<double,11> means = { 0.000,
+ const int numPeaks = 14;
+ array<double,14> means = { 0.000,
0.143,
0.279,
0.728,
@@ -37,8 +42,11 @@ void FitKnownPeaks(TH1F* hist){
1.981,
2.410,
2.910,
- 3.600,
- 3.792
+ 3.2,
+ 3.605,
+ 3.792,
+ 4.1,
+ 4.4
};
//Build individual peak fit functions
@@ -72,7 +80,10 @@ void FitKnownPeaks(TH1F* hist){
"+ ([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)) "
- "+ [23]" , minFit, maxFit);
+ "+ ([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
@@ -89,7 +100,10 @@ void FitKnownPeaks(TH1F* hist){
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,"Background");
+ 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;
@@ -99,9 +113,9 @@ void FitKnownPeaks(TH1F* hist){
full->SetParameter(2*i+2,1.0);
full->SetParLimits(2*i+2,0.0,1e5);
}
- full->FixParameter(numParams-1,0.0);
- //full->SetParLimits(numParams-1,0.0,1e1);
- //full->SetParLimits(numParams-1,0.0,0.0);
+ //full->FixParameter(numParams-1,0.0);
+ full->SetParameter(numParams-1,1.0);
+ full->SetParLimits(numParams-1,0.0,1e1);
//Fit full function to histogram
hist->Fit(full, "WWR", "", minFit, maxFit);
@@ -146,3 +160,163 @@ void FitKnownPeaks(TH1F* hist){
}
+double* FitKnownPeaks_RtrnArry(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,
+ 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];
+ 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 to histogram
+ 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;
+
+ double areas[numPeaks];
+ 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;
+
+ areas[i] = finalPar[2*i+2]/binWidth;
+ }
+
+ return areas;
+}
+
diff --git a/Projects/e793s/macro/ThreeBodyBreakup.h b/Projects/e793s/macro/ThreeBodyBreakup.h
new file mode 100644
index 000000000..62ca664d4
--- /dev/null
+++ b/Projects/e793s/macro/ThreeBodyBreakup.h
@@ -0,0 +1,39 @@
+void ThreeBodyBreakup(){
+
+ TGenPhaseSpace PS1p;
+
+ //47K(d,p)47K+p+n
+ NPL::Reaction R("47K(d,p)48K@355");
+ NPL::Particle K47("47K");
+ NPL::Particle p("1H");
+ NPL::Particle d("2H");
+ NPL::Particle n("n");
+
+ TLorentzVector target(0,0,0,d.Mass());
+ K47.SetKineticEnergy(354.75);
+ double E = K47.GetGamma()*K47.Mass();
+ double pc=sqrt(E*E-K47.Mass()*K47.Mass());
+ TLorentzVector beam(0,0,pc,K47.GetGamma()*K47.Mass());
+ TLorentzVector W = target+beam;
+
+ 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);
+ PS1p.SetDecay(W,3,masses1p);
+
+ for(unsigned int i = 0 ; i < 1000000 ; i++){
+ // 1p case
+ double weight = PS1p.Generate();
+ auto pd = PS1p.GetDecay(0);
+ double theta = pd->Vect().Angle(TVector3(0,0,1));
+ double beta = pd->Beta();
+ p.SetBeta(beta);
+ double ELab=p.GetEnergy();
+ double ex=R.ReconstructRelativistic(ELab,theta);
+ //cout << ELab << " " << ex << endl;
+ Ex1p->Fill(ex,weight);
+ }
+
+ Ex1p->Draw("");
+}
+
--
GitLab