diff --git a/Projects/Sofia/macro/ShowResultSimu.C b/Projects/Sofia/macro/ShowResultSimu.C
new file mode 100644
index 0000000000000000000000000000000000000000..45e9d161255f7d8304db60749c24e69036b7145d
--- /dev/null
+++ b/Projects/Sofia/macro/ShowResultSimu.C
@@ -0,0 +1,46 @@
+TChain* chain=NULL;
+
+//////////////////////////////////////////////////////////////////
+void LoadRootFile(){
+ chain = new TChain("SimulatedTree");
+
+ chain->Add("../../../Outputs/Simulation/sofia_simu_1.root");
+ chain->Add("../../../Outputs/Simulation/sofia_simu_2.root");
+ chain->Add("../../../Outputs/Simulation/sofia_simu_3.root");
+ chain->Add("../../../Outputs/Simulation/sofia_simu_4.root");
+}
+
+//////////////////////////////////////////////////////////////////
+void ShowResultSimu()
+{
+ LoadRootFile();
+
+ chain->Draw("fDetected_Position_Y:fDetected_Position_X>>hpos(500,-1700,-700,500,-350,350)","fTOF_Energy@.size()==2","colz");
+ TH2F* hpos = (TH2F*)gDirectory->FindObjectAny("hpos");
+
+ chain->Draw("fFC_Fragment_Z>>h1(35,30,65)");
+ TH1F* h1 = (TH1F*)gDirectory->FindObjectAny("h1");
+
+ chain->Draw("fDetected_Z>>h2(35,30,65)","fTOF_Energy@.size()==2");
+ TH1F* h2 = (TH1F*)gDirectory->FindObjectAny("h2");
+
+ h1->Sumw2();
+ h2->Sumw2();
+
+ TH1F* hratio = (TH1F*)h2->Clone();
+ hratio->Divide(h1);
+
+ TCanvas* c1 = new TCanvas("c1","c1",1200,1200);
+ c1->Divide(2,2);
+
+ c1->cd(1);
+ hpos->Draw("colz");
+ c1->cd(2);
+ h1->Draw();
+ c1->cd(3);
+ h2->Draw();
+ c1->cd(4);
+ hratio->Draw("E1");
+
+
+}
diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
index 89cc9dffb3f025d9c57ea7cc3d40708517ac0422..a8dd2197b03c0db803dd3f095076846d7c64bdf4 100644
--- a/Projects/e793s/macro/CS2.h
+++ b/Projects/e793s/macro/CS2.h
@@ -6,12 +6,14 @@ 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* theory, TGraphErrors* experiment);
+TList* peakFitList = new TList();
/* Global variables */
vector<Double_t> anglecentres, anglewidth;
TGraph* currentThry;
TGraphErrors* staticExp;
int indexE;
+double globalS, globalSerr;
/* Output volume toggle */
bool loud = 0;
@@ -45,7 +47,11 @@ void CS(double Energy, double Spin, double spdf, double angmom){
// p3/5 -> spdf = 1, angmom = 1.5
// J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
double J0 = 0.5;
- double ElasticNorm = 5.8, ElasticNormErr = 1.3; // DeuteronNorm in elastics, 5.6 +- 1.3
+ double ElasticNorm = 5.8, ElasticNormErr = 1.3; // DeuteronNorm in elastics, 5.8 +- 1.3
+ /* Reduce by factor of 10,000 */
+ /* WHY DOES THIS WORK???????*/
+ ElasticNorm /= 10000.;
+ ElasticNormErr /= 10000.;
double nodes;
if(spdf==1){
@@ -63,6 +69,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){
indexE = 100;
anglecentres.clear();
anglewidth.clear();
+ globalS=0.;
+ globalSerr=0.;
+ peakFitList->Clear();
/* Retrieve array index of the entered peak energy */
/* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
@@ -149,7 +158,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
<< " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"<< endl;
}
}
- //delete c_SolidAngle;
+ delete c_SolidAngle;
/* Graph of Area/Solid Angle*/
TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
@@ -212,8 +221,22 @@ void CS(double Energy, double Spin, double spdf, double angmom){
gdSdO->SetLineColor(kRed);
gdSdO->SetMarkerColor(kRed);
gdSdO->SetMarkerStyle(21);
+ /* Construct title string */
+ /**/ ostringstream textstream;
+ /**/ textstream << std::fixed << setprecision(3);
+ /**/ textstream << "Peak " << means[indexE];
+ /**/ string tempstr = textstream.str();
+ /**/ textstream << ": S = " << globalS
+ /**/ << " +- " << globalSerr;
+ /**/ string textstring = textstream.str();
+ gdSdO->SetTitle(textstring.c_str());
gdSdO->Draw("AP");
Final->Draw("SAME");
+ string savestring = "./CS2_Figures/"+tempstr+".root";
+ c_Chi2Min->SaveAs(savestring.c_str());
+
+ //delete c_AoSA;
+ //delete c_dSdO;
}
////////////////////////////////////////////////////////////////////////////////
@@ -255,13 +278,18 @@ vector<vector<double>> GetExpDiffCross(double Energy){
double max = min + bin;
cout << "===================================" << endl;
cout << "min: " << min << " max: " << max << endl;
+
+ stringstream tmp; tmp << fixed << setprecision(0);
+ tmp << "c_peakFits_" << min << "_" << max;
+ string tmp2 = tmp.str();
+ TCanvas* c_peakFits = new TCanvas("c_peakFits",tmp2.c_str(),1000,1000);
/* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
/* To change angle gates, run GateThetaLab_MultiWrite() */
TH1F* gate = PullThetaLabHist(i,105.,5.);
TH1F* pspace = PullPhaseSpaceHist(i,105.,5.);
- /* Scale Phase Space... */
+ /* Scale the Phase Space at this angle... */
/* ... for all angles together */
if(scaleTogether){
gate->Add(pspace,-trackScale);
@@ -270,6 +298,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){
else {
if(pspace->Integral() > 50.){ // Non-garbage histogram
pspace->Scale(0.01);
+ trackScale=0.01;
int numbins = gate->GetNbinsX();
for(int b=0; b<numbins; b++){
if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
@@ -277,11 +306,13 @@ vector<vector<double>> GetExpDiffCross(double Energy){
}
while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
pspace->Scale(0.9999);
+ trackScale*=0.9999;
}
if(loud){cout << " TO " << pspace->GetBinContent(b) <<
" > " << gate->GetBinContent(b) << endl;
}
}
+ cout << " !!! SCALE FOR THIS ANGLE = " << trackScale << endl;
gate->Add(pspace,-1);
}
}
@@ -292,7 +323,11 @@ vector<vector<double>> GetExpDiffCross(double Energy){
AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate);
/* Write PS-subtracted spectrum to list */
- list->Add(gate);
+ //list->Add(gate);
+ //list->Add(c_peakFits);
+ string filename = "./CS2_Figures/";
+ filename += tmp2 + ".root";
+ c_peakFits->SaveAs(filename.c_str());
/* Check correct OneGate vector is selected */
cout << "area of " << means[indexE] << " = "
@@ -306,11 +341,12 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Push relevant OneGate vector to end of AllGates vector */
OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
+ delete c_peakFits;
}
/* Write PS-subtracted spectrum to file */
- TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
- list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
+ //TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
+ //list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
return OnePeak_AllGates;
}
@@ -362,7 +398,7 @@ void Scale(TGraph* g , TGraphErrors* ex){
//scaleTWOFNR = mean;
cout << "SCALED THEORY BY " << mean << endl;
cout << " therefore S = " << 1/mean << " ??" << endl;
-
+
double* x = g->GetX();
double* y = g->GetY();
for(unsigned int i = 0 ; i < g->GetN() ; i++)
@@ -547,6 +583,10 @@ TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){
cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
}
+ /* Store S value in global variable, to access for drawing on plots */
+ globalS = xs[0];
+ globalSerr = err[0];
+
// Return the Fitted CS
TGraph* g = new TGraph();
double* X = theory->GetX();
diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C
index a7d5258fd3971e624b963a07a451ed37e3ea7d93..dffa20c2c1b9d737ebb2c44fa3bb321dc0323acb 100755
--- a/Projects/e793s/macro/DrawPlots.C
+++ b/Projects/e793s/macro/DrawPlots.C
@@ -560,8 +560,18 @@ void DrawPlots(){
cout << "\t- AGATA_efficiency(double Energy_kev) "<< endl;
cout << "\t- CorrectForAGATAEffic(TH1F* hist) "<< endl;
cout << "\t- CS(stateEnergy, stateSpin, orbital_l, orbital_j, nodes) "<< endl;
- cout << "\t---- 0.143, p3/2 = CS(0.143, 2, 1, 1.5, ?) "<< endl;
- cout << "\t---- 0.968, p1/2 = CS(0.968, 0, 1, 0.5, ?) "<< endl;
+ cout << "\t---- 0.143, p3/2 = CS(0.143, 2, 1, 1.5) "<< endl;
+ cout << "\t---- 0.279, p3/2 = CS(0.279, 2, 1, 1.5) "<< endl;
+ cout << "\t---- 0.728, f7/2 = CS(0.728, 3, 3, 3.5) "<< endl;
+ cout << "\t---- 0.968, p1/2 = CS(0.968, 0, 1, 0.5) "<< endl;
+ cout << "\t---- 1.410, p3/2 = CS(1.410, 2, 1, 1.5) "<< endl;
+ cout << "\t---- 1.981, p3/2 = CS(1.981, 2, 1, 1.5) "<< endl;
+ cout << "\t---- 2.410, p3/2 = CS(2.410, 2, 1, 1.5) "<< endl;
+ cout << "\t---- 3.2 , f7/2 = CS(3.2 , 3, 3, 3.5) "<< endl;
+ cout << "\t---- 3.6 , f7/2 = CS(3.6 , 3, 3, 3.5) "<< endl;
+ cout << "\t---- 3.8 , f7/2 = CS(3.8 , 3, 3, 3.5) "<< endl;
+ cout << "\t---- 4.1 , f7/2 = CS(4.1 , 3, 3, 3.5) "<< endl;
+ cout << "\t---- 4.4 , f7/2 = CS(4.4 , 3, 3, 3.5) "<< endl;
cout << ""<< endl;
cout << "==========================================" << endl;
diff --git a/Projects/e793s/macro/GausFit.h b/Projects/e793s/macro/GausFit.h
index 0e37d1ec016357ccfa9ab19298c0d5b9bf0b4ef0..22bb207c0388c9842ccda0afda3b2ee09627e9c0 100755
--- a/Projects/e793s/macro/GausFit.h
+++ b/Projects/e793s/macro/GausFit.h
@@ -14,7 +14,7 @@ Double_t pi = 3.14159265358979323846;
////////////////////////////////////////////////////////
void DoubleGaus(TH1F* hist){
- bool repeat=true;
+ bool repeat=true, bgbool = true;
int repeatInt;
double minFit, maxFit, mean1, mean2;
@@ -31,6 +31,8 @@ void DoubleGaus(TH1F* hist){
cin >> mean1;
cout << " Peak 2 = ";
cin >> mean2;
+ cout << " Background, yes or no?" << endl;
+ cin >> bgbool;
//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))",
@@ -64,13 +66,18 @@ void DoubleGaus(TH1F* hist){
g1->SetParameter(1, mean1);
g1->SetParLimits(1, mean1-0.5, mean1+0.5);
g1->SetParameter(2, 0.13);
+ //g1->SetParameter(2, 0.14);
g1->SetParLimits(2, 0.05, 0.30);
+ //g1->SetParLimits(2, 0.13, 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->SetParameter(2, 0.14);
g2->SetParLimits(2, 0.05, 0.30);
+ //g2->SetParLimits(2, 0.13, 0.30);
+
hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit);
hist->Fit(g2, "WWR", "", mean2-5, maxFit);//minFit, maxFit);
@@ -81,6 +88,12 @@ void DoubleGaus(TH1F* hist){
bg->GetParameters(&par[6]);
f1->SetParameters(par);
+ /* JUST FOR FITTING 0.36-GATED 5.3MeV PEAK! */
+ // f1->SetParLimits(2, 0.137, 0.40);
+ // f1->SetParLimits(5, 0.137, 0.40);
+
+ if(bgbool==false){bg->FixParameter(0,0.); f1->FixParameter(6,0.);}
+
hist->Fit(f1, "WWR", "", minFit, maxFit);
hist->Draw();
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index c70c4d6d06a68935caa978f520bb7eca77588790..b6d6be314a6e141732b8457c01a6c38d54942145 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -104,20 +104,27 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
for(int i=0; i<numPeaks; i++) {
full->FixParameter((i*3)+1,sigma);
full->FixParameter((i*3)+2,means.at(i));
+ // Set max 279 counts to 100
+ //full->SetParameter((i*3)+3,0.);//1e1);
+ //if(i==2){full->SetParLimits((i*3)+3,0.0,1e2);}
+ //else{full->SetParLimits((i*3)+3,0.0,1e5);}
full->SetParameter((i*3)+3,1e1);
full->SetParLimits((i*3)+3,0.0,1e5);
}
//full->SetParameter(0,30.);
- //full->SetParLimits(0,0.,40.);
- full->FixParameter(0,0.);
+ //full->SetParLimits(0,0.,40.); /* FOR TOTAL SPECTRUM FITTING */
+ full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
+ //full->FixParameter(0,0.);
+ full->FixParameter(9,0.); //??
// Specific limits
// Set max 279 counts to 100
- full->SetParLimits(9,0.0,1e2); // Doesnt work???
+ //full->FixParameter(9,0.0);
+ //full->SetParLimits(9,0.0,1e2); // Doesnt work???
allPeaks[14]->SetLineColor(kOrange);
//Fit full function to histogram
- hist->Fit(full, "RQ", "", minFit, maxFit);
+ hist->Fit(full, "RWQB", "", minFit, maxFit);
hist->Draw();
//Extract fitted variables, assign them to individual fits, and draw them
@@ -165,6 +172,7 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
}
cout << " BG " << full->GetParameter(0)
<< " +- " << full->GetParError(0) << endl;
+
return allpeaks;
}