diff --git a/Projects/e793s/autosim.sh b/Projects/e793s/autosim.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e077c83d1850d1f7580b3be0c2e146e327485397
--- /dev/null
+++ b/Projects/e793s/autosim.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+#./sim.sh 47Kdp 0000
+#./sim.sh 47Kdp 0143
+#./sim.sh 47Kdp 1981
+#./sim.sh 47Kdp 1410
+#./sim.sh 47Kdp 0279
+#./sim.sh 47Kdp 0728
+#./sim.sh 47Kdp 0968
+#./sim.sh 47Kdp 2412
+#./sim.sh 47Kdp 2910
+
+#./sim.sh 47Kdt 1945
+#./sim.sh 47Kdt 3344
+#./sim.sh 47Kdt 2730
+
+./sim.sh 47Kdt 0000
+./sim.sh 47Kdt 1000
+./sim.sh 47Kdt 2000
+./sim.sh 47Kdt 3000
+./sim.sh 47Kdt 4000
+./sim.sh 47Kdt 5000
+./sim.sh 47Kdt 6000
diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
index 4853419b10b71e32517e5980618e36aa6e0ea0b1..3949ca6d5db02abc52ccbad7bf1f95351c325bea 100644
--- a/Projects/e793s/macro/CS2.h
+++ b/Projects/e793s/macro/CS2.h
@@ -14,9 +14,11 @@ TGraph* currentThry;
TGraphErrors* staticExp;
int indexE;
double globalS, globalSerr;
+int numAngleBins = 20;
+vector<int> removePoint;
/* Output volume toggle */
-bool loud = 1;
+bool loud = 0;//1;
/* Scale method toggle */
bool scaleTogether = 1;
@@ -114,9 +116,11 @@ void CS(double Energy, double Spin, double spdf, double angmom){
//double ElasticNorm = 2.363, ElasticNormErr = 0.000; // DeuteronNorm in elastics, reconstructed from determined thickness
//double ElasticNorm = 3.7, ElasticNormErr = 0.000; // Estimated 'goal' normalization
- double ElasticNorm = 0.000220, ElasticNormErr = 0.000; //14Oct22
+ //double ElasticNorm = 0.000220, ElasticNormErr = 0.000; //14Oct22
+ double ElasticNorm = 0.000234, ElasticNormErr = 0.000; //18Oct22
inputdate = "18Oct22";
+ // Extract spin orbit name
orbitalname.clear();
orbital.clear();
if(spdf==1){
@@ -146,13 +150,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){
orbital="???";
}
- /* Reduce by factor of 10,000 */
- //ElasticNorm /= 10000.;
- ////ElasticNorm /= 1000.;
- //ElasticNormErr /= 10000.;
- ////ElasticNormErr /= 1000.;
+ // Number of nodes
double nodes;
-
if(spdf==1){
nodes=1;
}
@@ -193,16 +192,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){
}
/* Solid Angle (from simulation) */
- /* ADD OPTION TO CHANGE SOLID ANGLE FILE DEPENDING ON PEAK!!!!*/
- //auto file = new TFile("../SolidAngle_HistFile_47Kdp_26May22_v2_0143.root");
- //auto file = new TFile("../SolidAngle_HistFile_19Jul22_47Kdp_0p000.root");
- //auto file = new TFile("../SolidAngle_HistFile_19Jul22_47Kdp_1p981.root");
- //auto file = new TFile("../SolidAngle_HistFile_19Jul22_47Kdp_4p393.root");
- //auto file = new TFile("../SolidAngle_HistFile_30Jul22_47Kdp_0p000_ThetaBin0p5.root");
//auto file = new TFile("SolidAngle_HistFiles/SolidAngle_HistFile_10Aug22_TrueStripRemoval.root");
-
string backupFileName = "SolidAngle_HistFiles/SolidAngle_HistFile_10Aug22_TrueStripRemoval.root";
-
string saFileName = "SolidAngle_HistFiles/SAHF_";
saFileName.append(inputdate);
saFileName.append("_47Kdp_");
@@ -224,14 +215,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){
cout << BOLDRED << "FAILED TO OPEN MAIN OR BACKUP SOLID ANGLE FILE" << endl;
cout << RED << "Check SolidAngle file exists..." << RESET << endl;
}
- //cout << "MADE IT OUT" << endl;
-
-
//auto file = new TFile("../SolidAngle_HistFile_New.root");
- /* ADD OPTION TO CHANGE SOLID ANGLE FILE DEPENDING ON PEAK!!!!*/
TH1F* SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_Lab_MG");
- //cout << BLUE << "MADE IT HERE" << endl;
TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
SolidAngle->Draw("HIST");
SolidAngle->GetXaxis()->SetRangeUser(100.,160.);
@@ -310,6 +296,12 @@ void CS(double Energy, double Spin, double spdf, double angmom){
<< setprecision(3)
<< endl;
}
+
+ // Exclude bad angles
+ if(AoSA[i]<1E-6 ||AoSAerr[i]>1e2){
+ removePoint.push_back(i);
+ }
+
}
//delete c_SolidAngle;
@@ -344,6 +336,15 @@ void CS(double Energy, double Spin, double spdf, double angmom){
gdSdO->Draw();
c_dSdO->Update();
+// for(int i=0; i<numAngleBins; i++){
+// if(binary_search(removePoint.begin(), removePoint.end(), i)){
+// cout << BLUE << "removing point " << i << RESET << endl;
+// gAoSA->RemovePoint(i);
+// gdSdO->RemovePoint(i);
+// }
+// }
+
+
/* TWOFNR diff. cross section, in mb/msr */
TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
c_TWOFNR->SetLogy();
@@ -431,7 +432,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){
gdSdO->SetTitle(textstring.c_str());
gdSdO->GetYaxis()->SetTitleOffset(1.3);
gdSdO->GetYaxis()->SetTitleSize(0.042);
- gdSdO->GetXaxis()->SetRangeUser(103.,157.);
+ gdSdO->GetXaxis()->SetRangeUser(103.,157.);
+ //gdSdO->GetXaxis()->SetRangeUser(105.,155.);
+ gdSdO->GetXaxis()->SetNdivisions(512,kTRUE);
gdSdO->Draw("AP");
Final->Draw("SAME");
@@ -447,6 +450,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){
TLine* markzero = new TLine(103.,0.,157.,0.);
gResid->SetTitle("");
gResid->GetXaxis()->SetRangeUser(103.,157.);
+ //gResid->GetXaxis()->SetRangeUser(105.,155.);
+ gResid->GetXaxis()->SetNdivisions(512,kTRUE);
gResid->GetYaxis()->SetTitle("Residuals");
gResid->GetYaxis()->SetTitleSize(0.15);
gResid->GetYaxis()->SetTitleOffset(0.36);
@@ -486,7 +491,6 @@ vector<vector<double>> GetExpDiffCross(double Energy){
vector<vector<double>> AllPeaks_OneGate;
vector<vector<double>> OnePeak_AllGates;
/****CHANGE ANGLE GATING****/
- int numAngleBins = 20;
double widthAngleBins = 2.5;
double firstAngle = 105.;
/***************************/
@@ -536,7 +540,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
// TEMPORARY!!! REMOVE LAST THREE BINS ON HIGH ENERGY STATES!!!
- if(means[indexE] > 3.0){numAngleBins-=3;}
+ //if(means[indexE] > 3.0){numAngleBins-=3;}
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
diff --git a/Projects/e793s/macro/CS2_dt.h b/Projects/e793s/macro/CS2_dt.h
index cf4e0ae2346b0c7b337576aea965891508c03190..a479ea275bdea7a908d86d33d8c4736453ca1ba7 100644
--- a/Projects/e793s/macro/CS2_dt.h
+++ b/Projects/e793s/macro/CS2_dt.h
@@ -14,6 +14,12 @@ TGraph* currentThry;
TGraphErrors* staticExp;
int indexE;
double globalS, globalSerr;
+int numAngleBins = 18;//10;
+int numAngleBinsBase = 18;//10;
+double widthAngleBins = 1.;//2.;
+double firstAngle = 2.;//3.;
+
+vector<int> removePoint;
/* Output volume toggle */
bool loud = 1;
@@ -88,6 +94,8 @@ void CS_Diagnosis(){
void CS(){
/* Overload function */
cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
+ cout << "---- 1.945, p1/2 = CS(1.945, 1, 0, 0.5) "<< endl;
+ cout << "---- 1.945, d3/2 = CS(1.945, 1, 2, 1.5) "<< endl;
}
////////////////////////////////////////////////////////////////////////////////
@@ -96,7 +104,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){
// J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
double J0 = 0.5;
- double ElasticNorm = 0.000220, ElasticNormErr = 0.000; //14Oct22
+ //double ElasticNorm = 0.000220, ElasticNormErr = 0.000; //14Oct22
+ double ElasticNorm = 0.000234, ElasticNormErr = 0.000016; //18Oct22
inputdate = "18Oct22";
// Extract spin orbit name
@@ -154,6 +163,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
globalS=0.;
globalSerr=0.;
peakFitList->Clear();
+ numAngleBins=numAngleBinsBase;
/* Retrieve array index of the entered peak energy */
/* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
@@ -165,7 +175,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
indexE = i;
found = 1;
stringstream ss;
- ss << setfill('0') << setw(4) << (int)(means[i]*1000.);
+ ss << setfill('0') << setw(4) << (int)(means_dt[i]*1000.);
statename = ss.str();
}
}
@@ -249,11 +259,16 @@ void CS(double Energy, double Spin, double spdf, double angmom){
SA = SA*1000.; //sr->msr
SAerr = SAerr*1000.; //sr->msr
- /* Area over Solid angle ONLY */
- AoSA.push_back(areaArray[i][1]/SA);
- AoSAerr.push_back((areaArray[i][1]/SA)
- * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
- + pow(SAerr/SA,2)));
+ /* Area over Solid angle (exclude errors) */
+ if(SA<=0.){
+ SA = -20.; SAerr = -20.;
+ AoSA.push_back(-20.); AoSAerr.push_back(-20.);
+ } else {
+ AoSA.push_back(areaArray[i][1]/SA);
+ AoSAerr.push_back((areaArray[i][1]/SA)
+ * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
+ + pow(SAerr/SA,2)));
+ }
/* Differential Cross Section */
/* NOTE: DON'T INCLUDE NORM ERROR IN ERR PROPAGATION AS IT IS SYSTEMATIC! */
@@ -279,9 +294,15 @@ void CS(double Energy, double Spin, double spdf, double angmom){
<< setprecision(3)
<< endl;
}
+
+ // Exclude bad angles
+ if(AoSA[i]<1E-6 ||AoSAerr[i]>1e2){
+ removePoint.push_back(i);
+ }
+
}
//delete c_SolidAngle;
-
+
/* Graph of Area/Solid Angle*/
TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
c_AoSA->SetLogy();
@@ -313,6 +334,16 @@ void CS(double Energy, double Spin, double spdf, double angmom){
gdSdO->Draw();
c_dSdO->Update();
+// for(int i=0; i<numAngleBins; i++){
+// if(binary_search(removePoint.begin(), removePoint.end(), i)){
+// cout << BLUE << "removing point " << i << RESET << endl;
+// gAoSA->RemovePoint(i);
+// gdSdO->RemovePoint(i);
+// }
+// }
+
+
+
/* TWOFNR diff. cross section, in mb/msr */
TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
c_TWOFNR->SetLogy();
@@ -324,7 +355,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
/** TEMP **/
cout << "UNSCALED THEORY DIFF CROSS SECTION EVALUATED AT DATA POINTS:::" << endl;
cout << setprecision(6);
- for(int i=0; i<20; i++){
+ for(int i=0; i<numAngleBins; i++){
cout << anglecentres.at(i) << "\t" << TheoryDiffCross->Eval(anglecentres.at(i)) << endl;
}
@@ -375,6 +406,11 @@ void CS(double Energy, double Spin, double spdf, double angmom){
pad2->Draw();
pad1->cd();
+
+ double graphMin = (double) floor(gdSdO->GetPointX(0));
+ double graphMax = (double) ceil(gdSdO->GetPointX(numAngleBins-1));
+ cout << RED << graphMin << " " << graphMax << RESET << endl;
+
TGraph* Final = FindNormalisation(TheoryDiffCross,gdSdO);
gdSdO->SetLineColor(kRed);
gdSdO->SetMarkerColor(kRed);
@@ -400,7 +436,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
gdSdO->SetTitle(textstring.c_str());
gdSdO->GetYaxis()->SetTitleOffset(1.3);
gdSdO->GetYaxis()->SetTitleSize(0.042);
- gdSdO->GetXaxis()->SetRangeUser(103.,157.);
+ gdSdO->GetXaxis()->SetRangeUser(graphMin,graphMax);
gdSdO->Draw("AP");
Final->Draw("SAME");
@@ -413,9 +449,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){
gResid->SetPoint(n,x,residual);
gResid->SetPointError(n,0,gdSdO->GetErrorY(n));
}
- TLine* markzero = new TLine(103.,0.,157.,0.);
+ TLine* markzero = new TLine(graphMin,0.,graphMax,0.);
gResid->SetTitle("");
- gResid->GetXaxis()->SetRangeUser(103.,157.);
+ gResid->GetXaxis()->SetRangeUser(graphMin,graphMax);
gResid->GetYaxis()->SetTitle("Residuals");
gResid->GetYaxis()->SetTitleSize(0.15);
gResid->GetYaxis()->SetTitleOffset(0.36);
@@ -454,10 +490,6 @@ vector<vector<double>> GetExpDiffCross(double Energy){
cout << "========================================================" << endl;
vector<vector<double>> AllPeaks_OneGate;
vector<vector<double>> OnePeak_AllGates;
- /****CHANGE ANGLE GATING****/
- int numAngleBins = 7;//10;
- double widthAngleBins = 5.;//2.;
- double firstAngle = 0.;//3.;
/***************************/
double x[numAngleBins], y[numAngleBins];
//TList* list = new TList();
@@ -504,8 +536,23 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- // TEMPORARY!!! REMOVE LAST THREE BINS ON HIGH ENERGY STATES!!!
- if(means_dt[indexE] > 3.0){numAngleBins-=3;}
+ // TEMPORARY FIX FOR EX:THETACM PROBLEM!!
+ // Assuming angular bins are 2-20 in 1 degree sections...
+ if(means_dt[indexE] > 4.0){
+ numAngleBins-= 4; //Up to 16 deg
+ }
+ else if(means_dt[indexE] > 3.0){
+ numAngleBins-= 5; //Up to 15 deg
+ }
+ else if(means_dt[indexE] > 2.0){
+ numAngleBins-= 6; //Up to 14 deg
+ }
+ else if(means_dt[indexE] > 1.0){
+ numAngleBins-= 8; //Up to 12 deg
+ }
+ else if(means_dt[indexE] > 0.0){
+ numAngleBins-=12; //Up to 8 deg
+ }
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
@@ -602,17 +649,15 @@ vector<vector<double>> GetExpDiffCross(double Energy){
delete c_peakFits;
}
- /* Write PS-subtracted spectrum to file */
- //TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
- //list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
-
return OnePeak_AllGates;
}
////////////////////////////////////////////////////////////////////////////////
TH1F* PullThetaCMHist(int i, double minTheta, double gatesize){
//TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p2.root","READ");
- TFile* file = new TFile("GateThetaCMHistograms_21Oct22_47Kdt.root","READ");
+ //TFile* file = new TFile("GateThetaCMHistograms_21Oct22_47Kdt.root","READ");
+ //TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p1_10angles_0to20.root","READ");
+ TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p1.root","READ");
string histname = "cThetaCMGate_"
+ to_string((int) (minTheta+(i*gatesize))) + "-"
@@ -624,24 +669,6 @@ TH1F* PullThetaCMHist(int i, double minTheta, double gatesize){
return hist;
}
-////////////////////////////////////////////////////////////////////////////////
-/*
-TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize){
- //TFile* file = new TFile("GatePhaseSpaceThetaCMHistograms_ReadMe.root","READ");
- //TFile* file = new TFile("GatePhaseSpaceThetaCMHistograms_2p5degAngles_26May22v2.root","READ");
- //TFile* file = new TFile("GatePhaseSpaceThetaCMHistograms_11Jul22.root","READ");
- TFile* file = new TFile("GatePhaseSpaceThetaCMHistograms_29Aug22_TrueStripRemoval_0p05.root","READ");
- string histname = "cPSpaceThetaCMGate_"
- + to_string((int) (minTheta+(i*gatesize))) + "-"
- + to_string((int) (minTheta+((i+1)*gatesize)));
- cout << "Loading " << histname << endl;
- TList *list = (TList*)file->Get("GatePhaseSpaceThetaCMHistograms");
- TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
- file->Close();
- return hist;
-}
-*/
-
////////////////////////////////////////////////////////////////////////////////
void Scale(TGraph* g , TGraphErrors* ex){
double scale;
@@ -809,7 +836,7 @@ double Chi2(TGraph* theory, TGraphErrors* exper){
if(exper->GetPointY(i)>1.0e-8){ //0){
//chi=(exper->Eval(anglecentres[i])-theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
chi=(exper->GetPointY(i) - theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
- //cout << "COMPARE::::: " << exper->Eval(anglecentres[i]) << " TO " << exper->GetPointY(i) << endl;
+ cout << "COMPARE::::: Thr " << theory->Eval(anglecentres[i]) << " TO Exp " << exper->GetPointY(i) << endl;
Chi2 +=chi*chi;
}
}
@@ -859,10 +886,9 @@ TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){
parameter[i] = 0.8;
char name[4];
sprintf(name,"S%d",i+1);
- min->SetLimitedVariable(i,name,parameter[i],0.01,0,10);
+ min->SetLimitedVariable(i,name,parameter[i],0.01,0,20);
}
-
///// TO IMPROVE: FIND WAY OF OBTAINING NDF AND PRINT CHI2/NDF /////
// Minimise
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index b38eeb6ef85174db1f2ae27717120445f5a92ea2..42b6a3ccd6cfe4d6ac0f5476bbf7cc417a2fba2d 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -24,6 +24,14 @@ array<double,numPeaks> means = { 0.000,
5.82
};
+const int numPeaks_dt = 3;
+array<double,numPeaks_dt> means_dt = { 0.000,
+ 1.945,
+ 3.344,
+ };
+
+
+
array<double,27> knowngammas = { 0.143,
0.279,
0.449,
@@ -88,6 +96,21 @@ Double_t f_full(Double_t *x, Double_t *par) {
return result;
}
+Double_t f_full_dt(Double_t *x, Double_t *par) {
+ float xx = x[0];
+ double result, norm;
+ // Flat background
+// result = par[0];
+ result = 0;
+ // Add N peaks
+ for(int pk=0; pk<numPeaks_dt; pk++){
+ result += (par[3+(pk*3)]/(par[1+(pk*3)]*sqrt(2*pi)))
+ //* exp(-0.5*pow((xx-par[2+(pk*3)])/par[1+(pk*3)],2));
+ * exp(-0.5*pow((xx-par[2+(pk*3)]-par[0])/par[1+(pk*3)],2)); //added par 0 as shift in energy
+ }
+ return result;
+}
+
vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){
double minFit=-1.0, maxFit=8.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
@@ -159,8 +182,9 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){
//full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
//full->FixParameter(0,0.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
full->FixParameter(9,0.); /* FIX 0.279 AREA TO ZERO */
- full->SetParLimits(0,-0.5,+0.5); /* FOR WHEN PAR[0] IS VARIABLE ENERGY CENTRIOD SLIDER */
+ //full->SetParLimits(0,-0.5,+0.5); /* FOR WHEN PAR[0] IS VARIABLE ENERGY CENTRIOD SLIDER */
//full->FixParameter(0,slideshift); /* FOR WHEN PAR[0] IS FIXED ENERGY CENTRIOD SLIDER */
+ full->FixParameter(0,0.0); /* FOR WHEN FIXING ENERGY SLIDER TO ZERO */
full->SetParameter(52,0.39); /* SET 5.8MeV SIGMA */
full->SetParLimits(52,0.34,0.44); /* SET 5.8MeV SIGMA */
@@ -219,6 +243,130 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){
return allpeaks;
}
+vector<vector<double>> FitKnownPeaks_dt_RtrnArry(TH1F* hist, double slideshift){
+ double minFit=-1.0, maxFit=8.0;
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+ double sigma = 0.18;
+
+ hist->Sumw2();
+
+ /* Construct flat BG to subtract */
+ /**
+ cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
+ cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
+ cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
+ double ConstBG = 36.0; double ErrBG = 1.0;
+ int xbins = hist->GetXaxis()->GetNbins();
+ double xmin = hist->GetXaxis()->GetXmin();
+ double xmax = hist->GetXaxis()->GetXmax();
+ TH1F *FlatBG = new TH1F("FlatBG","FlatBG", xbins, xmin, xmax);
+ for(int i=0; i<xbins;i++){
+ FlatBG->SetBinContent(i,ConstBG);
+ FlatBG->SetBinError(i,ErrBG);
+ }
+ hist->Add(FlatBG,-1);
+ **/
+
+ //Build individual peak fit functions
+ string nameBase = "Peak ";
+ string function = "([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2))";
+ TF1 **allPeaks = new TF1*[numPeaks_dt];
+ for(int i=0; i<numPeaks_dt; i++) {
+ string nameHere = nameBase;
+ nameHere +=to_string(i);
+
+ allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit);
+ //allPeaks[i] = new TF1(nameHere.c_str(), f_peak, -1, 5);
+ allPeaks[i]->SetLineColor(kBlack);
+ allPeaks[i]->SetLineStyle(7);
+ allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth");
+ }
+
+ //Subtract flat background equal to smallest bin in range
+
+ hist->GetXaxis()->SetRange(hist->FindBin(-0.9), hist->FindBin(-0.4));
+ double bgmin = hist->GetBinContent(hist->GetMinimumBin());
+ hist->GetXaxis()->UnZoom();
+ cout << "Subtracting background of " << bgmin << endl;
+ for(int b=1; b<hist->GetNbinsX() ; b++){
+ hist->SetBinContent(b,hist->GetBinContent(b)-bgmin);
+ }
+
+
+ //Build background function
+ TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+ bg->SetParNames("Background");
+
+ //Build IMPROVED total function
+ TF1 *full = new TF1("fitAllPeaks", f_full_dt, minFit, maxFit, (int) 1+(3*numPeaks_dt));
+ full->SetLineColor(kRed);
+ const int numParams = (numPeaks_dt*3)+1;
+ for(int i=0; i<numPeaks_dt; i++) {
+ full->FixParameter((i*3)+1,sigma);
+ full->FixParameter((i*3)+2,means_dt.at(i));
+ full->SetParameter((i*3)+3,1e1);
+ full->SetParLimits((i*3)+3,0.0,1e5);
+ //full->SetParLimits((i*3)+3,10.0,1e5);
+ }
+ //full->SetParLimits(0,0.,40.); /* FOR TOTAL SPECTRUM FITTING */
+ //full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
+ //full->FixParameter(0,0.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
+
+ //Fit full function to histogram
+ hist->Fit(full, "RWQB", "", minFit, maxFit);
+ hist->Draw();
+
+ //Extract fitted variables, assign them to individual fits, and draw them
+ const Double_t* finalPar = full->GetParameters();
+ const Double_t* finalErr = full->GetParErrors();
+ for (int i=0; i<numPeaks_dt; i++){
+ allPeaks[i]->SetParameters(sigma, means_dt.at(i), finalPar[3+(i*3)]);
+ }
+ bg->SetParameter(0,finalPar[0]);
+ bg->Draw("SAME");
+ full->Draw("SAME");
+
+ for (int i=0; i<numPeaks_dt; 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;
+
+ vector<vector<double>> allpeaks;
+ for(int i=0; i<numPeaks_dt; i++){
+ double A = finalPar[(i*3)+3]/binWidth;
+ double deltaA = A * (finalErr[(i*3)+3]/finalPar[(i*3)+3]);
+
+ cout << fixed << setprecision(3)
+ << " #" << i << " "
+ << finalPar[(i*3)+2] << "\t" << setprecision(0)
+ << A << "\t+- "
+ << deltaA << setprecision(3)
+ << endl;
+
+ vector<double> onepeak; //energy, area and error for one peak
+ onepeak.push_back(finalPar[(i*3)+2]);
+ onepeak.push_back(A);
+ onepeak.push_back(deltaA);
+ allpeaks.push_back(onepeak);
+ }
+ cout << " BG " << full->GetParameter(0)
+ << " +- " << full->GetParError(0) << endl;
+
+ return allpeaks;
+}
+
void FitKnownPeaks(TH1F* hist){
//Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0);
diff --git a/Projects/e793s/macro/Plots_47Kdt.C b/Projects/e793s/macro/Plots_47Kdt.C
index 289d7de22c6c689d120d2d06070769a3b78369e6..74e7d97e9a489f4f3f769289951789ff45bef11f 100644
--- a/Projects/e793s/macro/Plots_47Kdt.C
+++ b/Projects/e793s/macro/Plots_47Kdt.C
@@ -115,7 +115,11 @@ void Plots_47Kdt(){
cout << "==============================================" << endl;
cout << "=============== (d,t) reaction ===============" << endl;
cout << "==============================================" << endl;
- cout << " - (d,t) selection: 'cutTritons' " << endl;
-
+ cout << " (d,t) selection: 'cutTritons' && 'cutTime' " << endl;
+ cout << "==============================================" << endl;
+ cout << ""<< endl;
+ CS();
+ cout << ""<< endl;
+ cout << "==============================================" << endl;
}
diff --git a/Projects/e793s/sim.sh b/Projects/e793s/sim.sh
index a3f089c2a778137055ad6fafac30a350d031fcac..ea06cb9e117e453520867988559921065f2db363 100755
--- a/Projects/e793s/sim.sh
+++ b/Projects/e793s/sim.sh
@@ -1,32 +1,68 @@
#!/bin/bash
rm RunToTreat_AutoGenerated.txt
+rm ./Reaction/Auto.reaction
-#====================================================
-#rfile='./Reaction/47Kdp_11Jul22_Sim0p000.reaction'
-#rfile='./Reaction/47Kdp_11Jul22_Sim1p981.reaction'
-#rfile='./Reaction/47Kdp_11Jul22_Sim4p045.reaction'
-#rfile='./Reaction/47Kdp_11Jul22_Sim4p393.reaction'
-#------
-#rfile='./Reaction/47Kdd_11Jul22.reaction'
-#rfile='./Reaction/47Kdd_11Jul22_Sim.reaction'
-#------
-#rfile='./Reaction/47Kpp_11Jul22.reaction'
-#----------------------------------------------------
-#rfile='./Reaction/47Kdp_18Oct22.reaction'
-rfile='./Reaction/47Kdt_18Oct22.reaction'
-#====================================================
-#====================================================
-#dfile='Detector/mugast_08Nov.detector'
-#----------------------------------------------------
-#dfile='Detector/mugast_11Jul22.detector'
-#dfile='Detector/mugast_testingMG1dip.detector'
-#----------------------------------------------------
-#dfile='Detector/mugast_01Sep22.detector'
-#----------------------------------------------------
-dfile='Detector/mugast_18Oct22.detector'
-#dfile='Detector/mugast_18Oct22_Tx2.detector'
-#dfile='Detector/mugast_18Oct22_Tx0p5.detector'
+XYZTdate="18Oct22"
+sp=" "
+
+rfile="./Reaction/Auto.reaction"
+E=$(echo "scale=3; $2/1000" | bc -l)
+
+#Boilerplate
+echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./Reaction/Auto.reaction
+echo -e "Beam" >> ./Reaction/Auto.reaction
+echo -e "$sp""Particle=""$sp""47K" >> ./Reaction/Auto.reaction
+echo -e "$sp""ExcitationEnergy=""$sp""0""$sp""MeV" >> ./Reaction/Auto.reaction
+echo -e "$sp""Energy=""$sp""354.75""$sp""MeV" >> ./Reaction/Auto.reaction
+echo -e "$sp""SigmaEnergy=""$sp""0.""$sp"" MeV" >> ./Reaction/Auto.reaction
+echo -e "$sp""SigmaThetaX=""$sp""0.0""$sp""deg" >> ./Reaction/Auto.reaction
+echo -e "$sp""SigmaPhiY=""$sp""0.0""$sp""deg" >> ./Reaction/Auto.reaction
+echo -e "$sp""SigmaX=""$sp""2.0""$sp""millimeter" >> ./Reaction/Auto.reaction
+echo -e "$sp""SigmaY=""$sp""2.0""$sp""millimeter" >> ./Reaction/Auto.reaction
+echo -e "$sp""MeanThetaX=""$sp""0""$sp""deg" >> ./Reaction/Auto.reaction
+echo -e "$sp""MeanPhiY=""$sp""0""$sp""deg" >> ./Reaction/Auto.reaction
+
+#Define XY of given XYZT date
+echo -e "$sp""MeanX=""$sp""-4.1561""$sp""millimeter" >> ./Reaction/Auto.reaction
+echo -e "$sp""MeanY=""$sp""+0.4701""$sp""millimeter" >> ./Reaction/Auto.reaction
+echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./Reaction/Auto.reaction
+echo -e "TwoBodyReaction" >> ./Reaction/Auto.reaction
+echo -e "$sp""Beam=""$sp""47K" >> ./Reaction/Auto.reaction
+
+#Define reaction
+if [[ "$1" == "47Kdp" ]]; then
+ echo "Running 47Kdp"
+ echo -e "$sp""Target=""$sp""2H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Light=""$sp""1H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Heavy=""$sp""48K" >> ./Reaction/Auto.reaction
+elif [[ "$1" == "47Kdt" ]]; then
+ echo "Running 47Kdt"
+ echo -e "$sp""Target=""$sp""2H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Light=""$sp""3H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Heavy=""$sp""46K" >> ./Reaction/Auto.reaction
+elif [[ "$1" == "47Kdd" ]]; then
+ echo "Running 47Kdd"
+ echo -e "$sp""Target=""$sp""2H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Light=""$sp""2H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Heavy=""$sp""47K" >> ./Reaction/Auto.reaction
+elif [[ "$1" == "47Kpp" ]]; then
+ echo "Running 47Kpp"
+ echo -e "$sp""Target=""$sp""1H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Light=""$sp""1H" >> ./Reaction/Auto.reaction
+ echo -e "$sp""Heavy=""$sp""47K" >> ./Reaction/Auto.reaction
+fi
+
+#Define heavy recoil excitation & finish
+echo -e "$sp""ExcitationEnergyLight=""$sp""0.0""$sp""MeV" >> ./Reaction/Auto.reaction
+echo -e "$sp""ExcitationEnergyHeavy=""$sp""$E""$sp""MeV" >> ./Reaction/Auto.reaction
+echo -e "$sp""CrossSectionPath=""$sp""flat.txt""$sp""CSR" >> ./Reaction/Auto.reaction
+echo -e "$sp""ShootLight=""$sp""1" >> ./Reaction/Auto.reaction
+echo -e "$sp""ShootHeavy=""$sp""0" >> ./Reaction/Auto.reaction
+
+
+#rfile="./Reaction/18Oct22_47Kdp_"$2".reaction"
+dfile="./Detector/mugast_"$XYZTdate".detector"
#====================================================
cd ~/Programs/nptool/Projects/e793s;
@@ -34,31 +70,33 @@ cmake ./;
make -j6;
directory='/home/charlottepaxman/Programs/nptool/Outputs/Simulation/'
-dotroot='.root'
-dash='-'
-for x in 1 #2 3 4 5 6 7 8 9 10
+nameconstruct=$XYZTdate"_"$1"_"$2
+echo "Name: $nameconstruct"
+echo "Reaction file: $rfile"
+echo "Detector file: $dfile"
+
+for x in 1 #2 3 4 #5 6 7 8 9 10
do
- outname=$1$dash$x
- #npsimulation -D $dfile -E $rfile -B ./runsimulation.mac -O $outname;
- npsimulation -D $dfile -E $rfile -B ./runsimulation_small.mac -O $outname;
- filename=$directory$1$dash$x$dotroot
+ outname=$nameconstruct"-"$x
+ npsimulation -D $dfile -E $rfile -B ./runsimulation.mac -O $outname;
+ #npsimulation -D $dfile -E $rfile -B ./runsimulation_small.mac -O $outname;
+ #filename=$directory$nameconstruct"-"$x".root"
done
-space=" "
-sim='Sim_'
-outfile=$sim$1
+outfile="Sim_"$nameconstruct
echo "TTreeName" >> RunToTreat_AutoGenerated.txt
echo " SimulatedTree" >> RunToTreat_AutoGenerated.txt
echo "RootFileName" >> RunToTreat_AutoGenerated.txt
-for x in 1 #2 3 4 5 6 7 8 9 10
+for x in 1 #2 3 4 #5 6 7 8 9 10
do
- echo "$space""$directory""$1""$dash""$x""$dotroot" >> RunToTreat_AutoGenerated.txt
+ echo "$sp""$directory""$nameconstruct""-""$x"".root" >> RunToTreat_AutoGenerated.txt
done
npanalysis --definition Sim -R RunToTreat_AutoGenerated.txt -E $rfile -D $dfile -O $outfile;
-#npanalysis --definition Sim --definition ExcludeThePoor -R RunToTreat_AutoGenerated.txt -E $rfile -D $dfile -O $outfile;
-
+#####npanalysis --definition Sim --definition ExcludeThePoor -R RunToTreat_AutoGenerated.txt -E $rfile -D $dfile -O $outfile;
+mvname="./macro/SolidAngle_HistFiles/SAHF_"$nameconstruct".root"
+mv SolidAngle_HistFile_New.root $mvname