diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c1a0c00b31db568ecf7cfcf3e506481f971bb93
--- /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 548fc50e9cc79a20dff275891c9af6a79f37d2d3..314a6b4ddce39096232221b85f291db4f1fa0bcf 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 6d02e9ed61989d7e009b41867b54046e953edbb9..055c48c4dc2b152f836d2c3f044fe90f3c77100b 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
Binary files /dev/null and b/Projects/e793s/macro/GateThetaLabHistograms.root differ
diff --git a/Projects/e793s/macro/GausFit.h b/Projects/e793s/macro/GausFit.h
new file mode 100755
index 0000000000000000000000000000000000000000..0e37d1ec016357ccfa9ab19298c0d5b9bf0b4ef0
--- /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 b3e0b0c2066356804f98e9b9821e0f084e322c86..9b7918a9ed61fc22adfa7069496cec7934c3157a
--- 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 0000000000000000000000000000000000000000..62ca664d4335edf994b0dff855d0004a63404378
--- /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("");
+}
+