diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
index 6a429bb2052dcd3373e1a7f1e01396cfa8880f72..89cc9dffb3f025d9c57ea7cc3d40708517ac0422 100644
--- a/Projects/e793s/macro/CS2.h
+++ b/Projects/e793s/macro/CS2.h
@@ -1,6 +1,7 @@
/* Predefine functions */
vector<vector<double>> GetExpDiffCross(double Energy);
TH1F* PullThetaLabHist(int i, double minTheta, double gatesize);
+TH1F* PullPhaseSpaceHist(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);
@@ -13,7 +14,10 @@ TGraphErrors* staticExp;
int indexE;
/* Output volume toggle */
-bool loud = 1;
+bool loud = 0;
+
+/* Scale method toggle */
+bool scaleTogether = 1;
////////////////////////////////////////////////////////////////////////////////
void CS(){
@@ -113,6 +117,14 @@ void CS(double Energy, double Spin, double spdf, double angmom){
anglecentres.push_back(((areaArray[i][4]-areaArray[i][3])/2.)+areaArray[i][3]);
anglewidth.push_back(areaArray[i][4]-areaArray[i][3]);
+ double tempsum=0, tempsumerr=0;
+ for(int x=binmin;x<binmax+1;x++){
+ tempsum += SolidAngle->GetBinContent(x);
+ tempsumerr += SolidAngle->GetBinError(x);
+ if(loud){cout << x << endl;}
+ }
+ if(loud){cout << "TEST CHECK " << tempsum << " +- " << tempsumerr << endl;}
+
double SAerr;
double SA = SolidAngle->IntegralAndError(binmin,binmax,SAerr);
SA = SA*1000.; //sr->msr
@@ -162,20 +174,20 @@ void CS(double Energy, double Spin, double spdf, double angmom){
0, &(expDCSerr[0]) );
gdSdO->SetTitle("Differential Cross Section");
gdSdO->GetXaxis()->SetTitle("ThetaLab [deg]");
- gdSdO->GetYaxis()->SetTitle("d#sigma/d#Theta [mb/msr]");
+ gdSdO->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
gdSdO->Draw();
/* TWOFNR diff. cross section, in mb/msr */
TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
c_TWOFNR->SetLogy();
TGraph* TheoryDiffCross = TWOFNR(means[indexE], J0, Spin, nodes, spdf, angmom);
- TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Theta [mb/msr]"); //msr set in func above
+ TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); //msr set in func above
TheoryDiffCross->GetXaxis()->SetTitle("ThetaLab [deg]");
TheoryDiffCross->Draw();
/* Convert AoSA into Differential Cross Section */
//cout << " SCALING BY NORMALIZATION = " << ElasticNorm << endl;
- //gAoSA->GetYaxis()->SetTitle("d#sigma/d#Theta [mb/msr]");
+ //gAoSA->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
/* Scaled and compared on same plot */
@@ -211,6 +223,32 @@ vector<vector<double>> GetExpDiffCross(double Energy){
vector<vector<double>> OnePeak_AllGates;
int numbins = 10;
double x[numbins], y[numbins];
+ TList* list = new TList();
+
+ /* Determine scaling factor for PhaseSpace */
+ double trackScale = 0.0;
+ if(scaleTogether){
+ TH1F* baseEx = PullThetaLabHist(0,105.,5.);
+ TH1F* basePS = PullPhaseSpaceHist(0,105.,5.);
+ for(int i=1; i<numbins;i++){
+ TH1F* addEx = PullThetaLabHist(i,105.,5.); baseEx->Add(addEx,1.);
+ TH1F* addPS = PullPhaseSpaceHist(i,105.,5.); basePS->Add(addPS,1.);
+ }
+ basePS->Scale(0.01);
+ trackScale = 0.01;
+ int numbinsScale = baseEx->GetNbinsX();
+ for(int b=0; b<numbinsScale; b++){
+ while(basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
+ basePS->Scale(0.99999);
+ trackScale *= 0.99999;
+ }
+ }
+ baseEx->Add(basePS,-1.);
+ baseEx->SetName("AllAngles");
+ list->Add(baseEx);
+ cout << " !!!!!!!!!!!!!!!FINAL SCALING = " << trackScale << endl;
+ }
+
for(int i=0; i<numbins;i++){
double bin = 5.;
double min = 105. + (i*bin);
@@ -221,12 +259,41 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* 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... */
+ /* ... for all angles together */
+ if(scaleTogether){
+ gate->Add(pspace,-trackScale);
+ }
+ /* ... or seperately for each angular bin */
+ else {
+ if(pspace->Integral() > 50.){ // Non-garbage histogram
+ pspace->Scale(0.01);
+ int numbins = gate->GetNbinsX();
+ for(int b=0; b<numbins; b++){
+ if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
+ " > " << gate->GetBinContent(b);
+ }
+ while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
+ pspace->Scale(0.9999);
+ }
+ if(loud){cout << " TO " << pspace->GetBinContent(b) <<
+ " > " << gate->GetBinContent(b) << endl;
+ }
+ }
+ gate->Add(pspace,-1);
+ }
+ }
/* Retrieve array containing all fits, for one angle gate. *
* Specific peak of interest selected from the vector by *
* global variable indexE */
AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate);
+ /* Write PS-subtracted spectrum to list */
+ list->Add(gate);
+
/* Check correct OneGate vector is selected */
cout << "area of " << means[indexE] << " = "
<< AllPeaks_OneGate[indexE][1]
@@ -240,17 +307,37 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Push relevant OneGate vector to end of AllGates vector */
OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
}
+
+ /* Write PS-subtracted spectrum to file */
+ TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
+ list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
+
return OnePeak_AllGates;
}
////////////////////////////////////////////////////////////////////////////////
TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
- TFile* file = new TFile("GateThetaLabHistograms.root","READ");
+ TFile* file = new TFile("GateThetaLabHistograms_ReadMe.root","READ");
string histname = "cThetaLabGate_"
- + to_string(minTheta+(i*gatesize)) + "-"
- + to_string(minTheta+((i+1)*gatesize));
+ + to_string((int) (minTheta+(i*gatesize))) + "-"
+ + to_string((int) (minTheta+((i+1)*gatesize)));
+ cout << "Loading " << histname << endl;
TList *list = (TList*)file->Get("GateThetaLabHistograms");
TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+// file->Close();
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize){
+ TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_ReadMe.root","READ");
+ string histname = "cPSpaceThetaLabGate_"
+ + to_string((int) (minTheta+(i*gatesize))) + "-"
+ + to_string((int) (minTheta+((i+1)*gatesize)));
+ cout << "Loading " << histname << endl;
+ TList *list = (TList*)file->Get("GatePhaseSpaceThetaLabHistograms");
+ TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+ file->Close();
return hist;
}
@@ -389,7 +476,7 @@ TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
}
////////////////////////////////////////////////////////////////////////////////
-double Chi2(TGraph* theory , TGraphErrors* exper){
+double Chi2(TGraph* theory, TGraphErrors* exper){
double Chi2 = 0 ;
double chi;
//for(int i = 1 ; i < exper->GetN() ; i++){
@@ -403,7 +490,6 @@ double Chi2(TGraph* theory , TGraphErrors* exper){
return Chi2;
}
-
////////////////////////////////////////////////////////////////////////////////
double ToMininize(const double* parameter){
static int f = 0 ;
diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C
index 9bfeae5f0fcf362f3fb75bd67c4c75aaf64623df..a7d5258fd3971e624b963a07a451ed37e3ea7d93 100755
--- a/Projects/e793s/macro/DrawPlots.C
+++ b/Projects/e793s/macro/DrawPlots.C
@@ -1,5 +1,4 @@
#include "GausFit.h"
-//#include "KnownPeakFitter_With47K.h"
#include "KnownPeakFitter.h"
#include "DrawPlots.h"
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index ebd6eaf71bcea2123cf557dbb47dd108b60db351..3e7ee21499cb88287dc4e92b08d8007981bc4543 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -1226,12 +1226,12 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
for (int i=0; i<numGates; i++){
double minTheta = startTheta + (i * gatesize);
- string title = to_string(minTheta)+" < ThetaLab < "+to_string(minTheta+gatesize);
+ string title = to_string((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize));
string gating = core
+ to_string(minTheta)
+ " && ThetaLab < "
+ to_string(minTheta+gatesize);
- string histname = "cThetaLabGate_" + to_string(minTheta) + "-" + to_string(minTheta+gatesize);
+ string histname = "cThetaLabGate_" + to_string((int) minTheta) + "-" + to_string((int) (minTheta+gatesize));
string draw = "Ex>>" + histname + "(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaLabGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000);
@@ -1239,6 +1239,7 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
TH1F* Ex_ThetaLabGate = (TH1F*) gDirectory->Get(histname.c_str());
Ex_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
Ex_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
+ Ex_ThetaLabGate->Sumw2();
Ex_ThetaLabGate->SetTitle(title.c_str());
list->Add(Ex_ThetaLabGate);
delete cEx_ThetaLabGate;
@@ -1248,3 +1249,39 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
list->Write("GateThetaLabHistograms",TObject::kSingleKey);
file->ls();
}
+
+void GatePhaseSpaceByThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize){
+ string core = "EventWeight*(Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > ";
+ string ytitle = "Counts / " + to_string(binsize) + " MeV";
+ double gatesize = (finishTheta-startTheta)/numGates;
+ TList* list = new TList();
+
+ TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_02Mar_47Kdp_PhaseSpace.root","READ");
+ TTree* PSTree = (TTree*) psfile->FindObjectAny("PhysicsTree");
+
+ for (int i=0; i<numGates; i++){
+ double minTheta = startTheta + (i * gatesize);
+ string title = to_string((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize));
+ string gating = core
+ + to_string(minTheta)
+ + " && ThetaLab < "
+ + to_string(minTheta+gatesize)
+ + ")";
+ string histname = "cPSpaceThetaLabGate_" + to_string((int) minTheta) + "-" + to_string((int) (minTheta+gatesize));
+ string draw = "Ex>>" + histname + "(" + to_string(10.0/binsize) + ",-1,9)";
+
+ TCanvas *cPSpace_ThetaLabGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000);
+ PSTree->Draw(draw.c_str(),gating.c_str(),"colz");
+ TH1F* PSpace_ThetaLabGate = (TH1F*) gDirectory->Get(histname.c_str());
+ PSpace_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
+ PSpace_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
+ PSpace_ThetaLabGate->Sumw2();
+ PSpace_ThetaLabGate->SetTitle(title.c_str());
+ list->Add(PSpace_ThetaLabGate);
+ delete cPSpace_ThetaLabGate;
+ }
+
+ TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms.root","RECREATE");
+ list->Write("GatePhaseSpaceThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index 3bacf35cd1125efa8167ade966514355f428d456..c70c4d6d06a68935caa978f520bb7eca77588790 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -60,6 +60,7 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
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;
@@ -73,7 +74,7 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
FlatBG->SetBinError(i,ErrBG);
}
hist->Add(FlatBG,-1);
- /**/
+ **/
//Build individual peak fit functions
string nameBase = "Peak ";
diff --git a/Projects/e793s/macro/ThreeBodyBreakup.h b/Projects/e793s/macro/ThreeBodyBreakup.h
index 5bdbe366e5519a29bc968ee0aecbc199d02dda8c..3df02dc9cc8c2d2b1211a0a78ac0c01cafa63e61 100644
--- a/Projects/e793s/macro/ThreeBodyBreakup.h
+++ b/Projects/e793s/macro/ThreeBodyBreakup.h
@@ -8,8 +8,12 @@ Double_t f_semi(Double_t *x, Double_t *par){
return f;
}
-
void ThreeBodyBreakup(){
+ TFile *fweight = new TFile("../../../Outputs/Simulation/EventWeight_25Feb.root");
+ TCanvas *cweight = (TCanvas*)fweight->Get("c1");
+ TH1F *hweight = (TH1F*)cweight->GetPrimitive("hEventWeight");
+// hweight->Draw();
+
TGenPhaseSpace PS1p;
//47K(d,p)47K+p+n
@@ -42,10 +46,11 @@ void ThreeBodyBreakup(){
double ex=R.ReconstructRelativistic(ELab,theta);
//cout << ELab << " " << ex << endl;
Ex1p->Fill(ex,weight);
- }
+ }
TCanvas* c_Ex1p = new TCanvas("c_Ex1p","c_Ex1p",1000,500);
Ex1p->Draw();
+ Ex1p->Multiply(hweight);
// TCanvas* c_semi = new TCanvas("c_semi","c_semi",1000,500);
TF1 *f1 = new TF1("f_semi",f_semi,-10,40,3);