diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index 0f8ab929b7c59268a856a923ae3266a4e93e9592..646433aa6522c3e35dba24f969132881d5715e9c 100755
--- a/Projects/e793s/Analysis.cxx
+++ b/Projects/e793s/Analysis.cxx
@@ -32,6 +32,7 @@ using namespace std;
#include "NPAnalysisFactory.h"
#include "NPDetectorManager.h"
#include "NPOptionManager.h"
+#include "macro/DefineColours.h"
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
@@ -56,6 +57,7 @@ void Analysis::Init() {
} else {
cout << " == == == == PHASE SPACE == == == ==" << endl;
isSim=false;
+ //isSim=true;
isPhaseSpace=true;
}
@@ -76,8 +78,8 @@ void Analysis::Init() {
if (i==6){int j=7;}
string name1 = base1 + to_string(j);
string name2 = base2 + to_string(j);
- ThetaCM_detected_MGX[i] = new TH1F(name1.c_str(),name1.c_str(),180,0,180);
- ThetaLab_detected_MGX[i] = new TH1F(name2.c_str(),name2.c_str(),180,0,180);
+ ThetaCM_detected_MGX[i] = new TH1F(name1.c_str(),name1.c_str(),NumThetaAngleBins,0,180); //900 bins for 0.2 angular bin width
+ ThetaLab_detected_MGX[i] = new TH1F(name2.c_str(),name2.c_str(),NumThetaAngleBins,0,180);
}
// Initilize MMX histograms
@@ -87,8 +89,8 @@ void Analysis::Init() {
int j=i+1;
string name1 = base3 + to_string(j);
string name2 = base4 + to_string(j);
- ThetaCM_detected_MMX[i] = new TH1F(name1.c_str(),name1.c_str(),180,0,180);
- ThetaLab_detected_MMX[i] = new TH1F(name2.c_str(),name2.c_str(),180,0,180);
+ ThetaCM_detected_MMX[i] = new TH1F(name1.c_str(),name1.c_str(),NumThetaAngleBins,0,180);
+ ThetaLab_detected_MMX[i] = new TH1F(name2.c_str(),name2.c_str(),NumThetaAngleBins,0,180);
}
}
@@ -134,7 +136,13 @@ void Analysis::Init() {
FinalBeamEnergy = BeamTargetELoss.Slow(OriginalBeamEnergy, 0.5*TargetThickness, 0);
reaction.SetBeamEnergy(FinalBeamEnergy);
- cout << "Beam energy at mid-target: " << FinalBeamEnergy << endl;
+ cout << "\033[91m Running for reaction "
+ << reaction.GetParticle1()->GetName() << "("
+ << reaction.GetParticle2()->GetName() << ","
+ << reaction.GetParticle3()->GetName() << ")"
+ << reaction.GetParticle4()->GetName() << endl;
+
+ cout << "\033[36m Beam energy at mid-target: " << FinalBeamEnergy << "\033[37m"<< endl;
// initialize various parameters
Rand = TRandom3();
@@ -266,9 +274,19 @@ void Analysis::TreatEvent(){
// if(!isSim){
// Evaluate energy using the thickness
- elab_tmp = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface);
+ elab_tmp = LightAl.EvaluateInitialEnergy(
+ Energy,
+ 0.4*micrometer,
+ ThetaM2Surface);
+ ELoss_Al.push_back(Energy-elab_tmp);
+ double elab_tmp2 = elab_tmp;
// Target Correction
- elab_tmp = LightTarget.EvaluateInitialEnergy(elab_tmp, 0.5*TargetThickness, ThetaNormalTarget);
+ elab_tmp = LightTarget.EvaluateInitialEnergy(
+ elab_tmp,
+ 0.5*TargetThickness,
+ ThetaNormalTarget);
+ ELoss_Target.push_back(elab_tmp2-elab_tmp);
+ ELoss.push_back(Energy-elab_tmp);
// } else {elab_tmp = Energy;}
ELab.push_back(elab_tmp);
@@ -346,10 +364,14 @@ void Analysis::TreatEvent(){
Energy, //particle energy after Al
0.4*micrometer, //thickness of Al
ThetaMGSurface); //angle of impingement
+ ELoss_Al.push_back(Energy-elab_tmp);
+ double elab_tmp2 = elab_tmp;
elab_tmp = LightTarget.EvaluateInitialEnergy(
elab_tmp, //particle energy after leaving target
TargetThickness*0.5, //distance passed through target
ThetaNormalTarget); //angle of exit from target
+ ELoss_Target.push_back(elab_tmp2-elab_tmp);
+ ELoss.push_back(Energy-elab_tmp);
} else { //TESTING DIFFERENT ENERGY LOSSES IN SIMULATION
elab_tmp = Energy; //so I can add and remove sections
//elab_tmp = LightSi.EvaluateInitialEnergy(
@@ -368,6 +390,14 @@ void Analysis::TreatEvent(){
ELab.push_back(elab_tmp);
+ //cout << "===============" << endl;
+ //cout << "RawE:\t" << RawEnergy.back() << endl;
+ //cout << "ELAl:\t" << ELoss_Al.back() << endl;
+ //cout << "ELCD:\t" << ELoss_Target.back() << endl;
+ //cout << "ELTt:\t" << ELoss.back() << endl;
+ //cout << "ELab:\t" << ELab.back() << endl;
+
+
// Part 3 : Excitation Energy Calculation
//if(!isSim){ //TESTING!!!!
Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp));
@@ -512,17 +542,45 @@ void Analysis::TreatEvent(){
////////////////////////////////////////////////////////////////////////////////
void FillSolidAngles(TH1F* hSA, TH1F* hDet, TH1F* hEmm){
+ if(!filledCline){
+ cout << RED << "FILLING CLINE VECTOR! SHOULD ONLY OCCUR ONCE!" << RESET << endl;
+ for(int t=0; t<NumThetaAngleBins; t++){
+ double angleMin = (t)*(180./NumThetaAngleBins);
+ double angleMax = (t+1)*(180./NumThetaAngleBins);
+ cout << " Angle " << angleMin
+ << " to " << angleMax
+ << endl;
+ clineVal.push_back(2.0*M_PI*(cos(angleMin*degtorad) - cos(angleMax*degtorad)));
+ clineX.push_back(angleMin+((angleMax-angleMin)/2.0));
+ filledCline=true;
+ }
+ }
+
for (int i = 0; i < hSA->GetNbinsX(); ++i){
double val = hDet->GetBinContent(i) / hEmm->GetBinContent(i);
double valerr = val * sqrt(
pow(hDet->GetBinError(i) / hDet->GetBinContent(i), 2) +
pow(hEmm->GetBinError(i) / hEmm->GetBinContent(i), 2) );
if (isnan(val)) { val = 0; valerr = 0; }
+ val *= clineVal.at(i);
+ valerr *= clineVal.at(i);
+
hSA->SetBinContent(i, val);
hSA->SetBinError(i, valerr);
}
}
+
+//void DivideByCline(TH1F* histo){
+//
+// for(int b=0; b<NumThetaAngleBins; b++){
+// int store = histo->GetBinContent(b);
+// histo->SetBinContent(b,(int) store/clineVal.at(b));
+/// }
+//}
+
+
+
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
cout << endl << "\e[1;32m GATCONF Statistics " << endl ;
@@ -532,7 +590,7 @@ void Analysis::End(){
cout << endl ;
if(isSim && !isPhaseSpace){
-
+
//TObjArray HistList(0);
TList *HistList = new TList();
@@ -556,14 +614,19 @@ void Analysis::End(){
Efficiency_CM_MM->Divide(ThetaCM_emmitted);
Efficiency_Lab_MM->Divide(ThetaLab_emmitted);
- double dt_MM = 180./Efficiency_Lab_MM->GetNbinsX();
- cout << "Angular infinitesimal (MM) = " << dt_MM << "deg " << endl;
- auto Cline_MM = new TF1("Cline_MM",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt_MM,M_PI),0,180);
+ //double dt_MM = 180./Efficiency_Lab_MM->GetNbinsX();
+ //cout << "Angular infinitesimal (MM) = " << dt_MM << "deg " << endl;
+ //auto Cline_MM = new TF1("Cline_MM",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt_MM,M_PI),0,180);
- SolidAngle_CM_MM->Divide(ThetaCM_emmitted);
- SolidAngle_CM_MM->Divide(Cline_MM,1);
- SolidAngle_Lab_MM->Divide(ThetaLab_emmitted);
- SolidAngle_Lab_MM->Divide(Cline_MM,1);
+ FillSolidAngles(SolidAngle_CM_MM, ThetaCM_detected_MM, ThetaCM_emmitted);
+ FillSolidAngles(SolidAngle_Lab_MM, ThetaLab_detected_MM, ThetaLab_emmitted);
+
+ //SolidAngle_CM_MM->Divide(ThetaCM_emmitted);
+ //SolidAngle_CM_MM->Divide(Cline_MM,1);
+ //SolidAngle_CM_MM->Divide(Cline);
+ //SolidAngle_Lab_MM->Divide(ThetaLab_emmitted);
+ //SolidAngle_Lab_MM->Divide(Cline_MM,1);
+ //SolidAngle_Lab_MM->Divide(Cline);
HistList->Add(ThetaCM_emmitted);
HistList->Add(ThetaLab_emmitted);
@@ -573,7 +636,8 @@ void Analysis::End(){
HistList->Add(Efficiency_Lab_MM);
HistList->Add(SolidAngle_CM_MM);
HistList->Add(SolidAngle_Lab_MM);
- HistList->Add(Cline_MM);
+ //HistList->Add(Cline_MM);
+ //HistList->Add(Cline);
// MUGAST
@@ -596,16 +660,17 @@ void Analysis::End(){
Efficiency_CM_MG->Divide(ThetaCM_emmitted);
Efficiency_Lab_MG->Divide(ThetaLab_emmitted);
- double dt_MG = 180./Efficiency_Lab_MG->GetNbinsX();
- cout << "Angular infinitesimal (MG) = " << dt_MG << "deg " << endl;
- auto Cline_MG = new TF1("Cline_MG",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt_MG,M_PI),0,180);
+ //double dt_MG = 180./Efficiency_Lab_MG->GetNbinsX();
+ //cout << "Angular infinitesimal (MG) = " << dt_MG << "deg " << endl;
+ //auto Cline_MG = new TF1("Cline_MG",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt_MG,M_PI),0,180);
/* Testing method for better errors in SolidAngle histograms */
FillSolidAngles(SolidAngle_CM_MG, ThetaCM_detected_MG, ThetaCM_emmitted);
- SolidAngle_CM_MG->Divide(Cline_MG,1);
-
+ //SolidAngle_CM_MG->Divide(Cline_MG,1);
+ //SolidAngle_CM_MG->Divide(Cline);
FillSolidAngles(SolidAngle_Lab_MG, ThetaLab_detected_MG, ThetaLab_emmitted);
- SolidAngle_Lab_MG->Divide(Cline_MG,1);
+ //SolidAngle_Lab_MG->Divide(Cline_MG,1);
+ //SolidAngle_Lab_MG->Divide(Cline);
HistList->Add(ThetaCM_detected_MG);
HistList->Add(ThetaLab_detected_MG);
@@ -613,7 +678,7 @@ void Analysis::End(){
HistList->Add(Efficiency_Lab_MG);
HistList->Add(SolidAngle_CM_MG);
HistList->Add(SolidAngle_Lab_MG);
- HistList->Add(Cline_MG);
+ //HistList->Add(Cline_MG);
// MUGAST INDIVIDUAL
TH1F *SolidAngle_CM_MGX[6];
@@ -626,7 +691,8 @@ void Analysis::End(){
SolidAngle_CM_MGX[i]->SetName(name.c_str());
SolidAngle_CM_MGX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_CM_MGX[i], ThetaCM_detected_MGX[i], ThetaCM_emmitted);
- SolidAngle_CM_MGX[i]->Divide(Cline_MG,1);
+ //SolidAngle_CM_MGX[i]->Divide(Cline_MG,1);
+ //SolidAngle_CM_MGX[i]->Divide(Cline);
}
TH1F *SolidAngle_Lab_MGX[6];
@@ -639,7 +705,8 @@ void Analysis::End(){
SolidAngle_Lab_MGX[i]->SetName(name.c_str());
SolidAngle_Lab_MGX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_Lab_MGX[i], ThetaLab_detected_MGX[i], ThetaLab_emmitted);
- SolidAngle_Lab_MGX[i]->Divide(Cline_MG,1);
+ //SolidAngle_Lab_MGX[i]->Divide(Cline_MG,1);
+ //SolidAngle_Lab_MGX[i]->Divide(Cline);
}
// MUST2 INDIVIDUAL
@@ -652,7 +719,8 @@ void Analysis::End(){
SolidAngle_CM_MMX[i]->SetName(name.c_str());
SolidAngle_CM_MMX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_CM_MMX[i], ThetaCM_detected_MMX[i], ThetaCM_emmitted);
- SolidAngle_CM_MMX[i]->Divide(Cline_MM,1);
+ //SolidAngle_CM_MMX[i]->Divide(Cline_MM,1);
+ //SolidAngle_CM_MMX[i]->Divide(Cline);
}
TH1F *SolidAngle_Lab_MMX[6];
@@ -664,7 +732,8 @@ void Analysis::End(){
SolidAngle_Lab_MMX[i]->SetName(name.c_str());
SolidAngle_Lab_MMX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_Lab_MMX[i], ThetaLab_detected_MMX[i], ThetaLab_emmitted);
- SolidAngle_Lab_MMX[i]->Divide(Cline_MM,1);
+ //SolidAngle_Lab_MMX[i]->Divide(Cline_MM,1);
+ //SolidAngle_Lab_MMX[i]->Divide(Cline);
}
for(int i=0; i<6; i++){HistList->Add(ThetaCM_detected_MGX[i]);}
@@ -677,6 +746,14 @@ void Analysis::End(){
for(int i=0; i<5; i++){HistList->Add(SolidAngle_CM_MMX[i]);}
for(int i=0; i<5; i++){HistList->Add(SolidAngle_Lab_MMX[i]);}
+ auto clineValGraph = new TGraph(NumThetaAngleBins);
+ clineValGraph->SetName("clineValGraph");
+ clineValGraph->SetTitle("clineValGraph");
+ for(int b=0; b<NumThetaAngleBins; b++){
+ clineValGraph->SetPoint(b,clineX.at(b),clineVal.at(b));
+ }
+ HistList->Add(clineValGraph);
+
auto HistoFile = new TFile("~/Programs/nptool/Projects/e793s/SolidAngle_HistFile_New.root","RECREATE");
HistList->Write();
HistoFile->Close();
@@ -706,6 +783,9 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
RootOutput::getInstance()->GetTree()->Branch("Ecm",&Ecm);
RootOutput::getInstance()->GetTree()->Branch("RawEnergy",&RawEnergy);
+ RootOutput::getInstance()->GetTree()->Branch("ELoss_Al",&ELoss_Al);
+ RootOutput::getInstance()->GetTree()->Branch("ELoss_Target",&ELoss_Target);
+ RootOutput::getInstance()->GetTree()->Branch("ELoss",&ELoss);
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab);
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM);
@@ -967,6 +1047,9 @@ void Analysis::ReInitValue(){
EAgata = -1000;
ELab.clear();
RawEnergy.clear();
+ ELoss_Al.clear();
+ ELoss_Target.clear();
+ ELoss.clear();
ThetaLab.clear();
PhiLab.clear();
ThetaCM.clear();
@@ -991,6 +1074,8 @@ void Analysis::ReInitValue(){
thetalab_tmp = 0;
philab_tmp = 0;
+ filledCline=false;
+
MG_T=-1000;
MG_E=-1000;
MG_X=-1000;
diff --git a/Projects/e793s/Analysis.h b/Projects/e793s/Analysis.h
index c6894b6d4614d26b709a8c2076316cf69d9c75c2..3a4ad387239dcee1c7c97e2b288ebb1b8712894f 100755
--- a/Projects/e793s/Analysis.h
+++ b/Projects/e793s/Analysis.h
@@ -42,13 +42,20 @@
#include <TMath.h>
#include <bitset>
-auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",180,0,180);
-auto ThetaCM_detected_MM = new TH1F("ThetaCM_detected_MM","ThetaCM_detected_MM",180,0,180);
-auto ThetaCM_detected_MG = new TH1F("ThetaCM_detected_MG","ThetaCM_detected_MG",180,0,180);
+int NumThetaAngleBins = 1800;// 180 = 1 deg, 360 = 0.5 deg, 900 = 0.2 deg, 1800 = 0.1 deg
-auto ThetaLab_emmitted = new TH1F("ThetaLab_emmitted","ThetaLab_emmitted",180,0,180);
-auto ThetaLab_detected_MM = new TH1F("ThetaLab_detected_MM","ThetaLab_detected_MM",180,0,180);
-auto ThetaLab_detected_MG = new TH1F("ThetaLab_detected_MG","ThetaLab_detected_MG",180,0,180);
+auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",NumThetaAngleBins,0,180);
+auto ThetaCM_detected_MM = new TH1F("ThetaCM_detected_MM","ThetaCM_detected_MM",NumThetaAngleBins,0,180);
+auto ThetaCM_detected_MG = new TH1F("ThetaCM_detected_MG","ThetaCM_detected_MG",NumThetaAngleBins,0,180);
+
+auto ThetaLab_emmitted = new TH1F("ThetaLab_emmitted","ThetaLab_emmitted",NumThetaAngleBins,0,180);
+auto ThetaLab_detected_MM = new TH1F("ThetaLab_detected_MM","ThetaLab_detected_MM",NumThetaAngleBins,0,180);
+auto ThetaLab_detected_MG = new TH1F("ThetaLab_detected_MG","ThetaLab_detected_MG",NumThetaAngleBins,0,180);
+//auto HistCline_MM = new TH1F("HistCline_MM","HistClint_MM",NumThetaAngleBins,0,180);
+//auto HistCline_MG = new TH1F("HistCline_MG","HistClint_MG",NumThetaAngleBins,0,180);
+double degtorad = M_PI/180.;
+vector<double> clineVal, clineX;
+bool filledCline;
TH1F *ThetaCM_detected_MGX[6];
TH1F *ThetaCM_detected_MMX[5];
@@ -94,6 +101,9 @@ class Analysis: public NPL::VAnalysis{
std::vector<double> Ex;
std::vector<double> Ecm;
std::vector<double> RawEnergy;
+ std::vector<double> ELoss_Al;
+ std::vector<double> ELoss_Target;
+ std::vector<double> ELoss;
std::vector<double> ThetaLab;
std::vector<double> PhiLab;
std::vector<double> ThetaCM;
diff --git a/Projects/e793s/exp.sh b/Projects/e793s/exp.sh
index 263ee445adae3059accc1316e8a6ab05e4fed0ca..4d78921e4c30fc4245971f171cf26ac3d907b17c 100755
--- a/Projects/e793s/exp.sh
+++ b/Projects/e793s/exp.sh
@@ -2,5 +2,26 @@
cd ~/Programs/nptool/Projects/e793s;
cmake ./;
make -j6;
-#npanalysis --definition Exp -R RunToTreat_PartII.txt -E Reaction/47Kdp_08Nov.reaction -D Detector/mugast_08Nov.detector -C Calibration.txt -O $1;
-npanalysis --definition Exp -R RunToTreat_PartI.txt -E Reaction/47Kdp_08Nov.reaction -D Detector/mugast_18May22_TestingThickerTarget.detector -C Calibration.txt -O $1;
+
+#====================================================
+#rfile='Reaction/47Kdp_08Nov.reaction'
+#rfile='Reaction/47Kdd_08Nov.reaction'
+#rfile='Reaction/47Kpp_08Nov.reaction'
+#rfile='Reaction/47Kdt_08Nov.reaction'
+#----------------------------------------------------
+#rfile='Reaction/47Kdp_11Jul22.reaction'
+#rfile='Reaction/47Kdt_11Jul22.reaction'
+rfile='Reaction/47Kdd_11Jul22.reaction'
+#rfile='Reaction/47Kpp_11Jul22.reaction'
+#====================================================
+#dfile='Detector/mugast_08Nov.detector'
+#----------------------------------------------------
+dfile='Detector/mugast_11Jul22.detector'
+#====================================================
+
+
+ptI='_PartI'
+ptII='_PartII'
+
+npanalysis --definition Exp -R RunToTreat_PartI.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptI;
+npanalysis --definition Exp -R RunToTreat_PartII.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptII;
diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
index 4608b0578473ec44320cdb6134988495c90e3e9a..d35f3ebe395fcefd3a76537f03d1de09283ab42f 100644
--- a/Projects/e793s/macro/CS2.h
+++ b/Projects/e793s/macro/CS2.h
@@ -16,7 +16,7 @@ int indexE;
double globalS, globalSerr;
/* Output volume toggle */
-bool loud = 0;
+bool loud = 1;
/* Scale method toggle */
bool scaleTogether = 1;
@@ -25,22 +25,95 @@ bool scaleTogether = 1;
string orbitalname;
string orbital;
+
+////////////////////////////////////////////////////////////////////////////////
+void canclone(TCanvas* major, int padNum, string name){
+ string minorName = "./CS2_Figures/" + name + ".root";
+ cout << minorName << endl;
+ TFile* minorFile = new TFile(minorName.c_str(),"READ");
+ TCanvas* minor = (TCanvas*) minorFile->FindObjectAny("c_peakFits");
+
+ major->cd(padNum);
+ TPad *pad = (TPad*)gPad;
+ minor->cd();
+ TObject *obj, *clone;
+ pad->Range(minor->GetX1(),minor->GetY1(),minor->GetX2(),minor->GetY2());
+ pad->SetTickx(minor->GetTickx());
+ pad->SetTicky(minor->GetTicky());
+ pad->SetGridx(minor->GetGridx());
+ pad->SetGridy(minor->GetGridy());
+ pad->SetLogx(minor->GetLogx());
+ pad->SetLogy(minor->GetLogy());
+ pad->SetLogz(minor->GetLogz());
+ pad->SetBorderSize(minor->GetBorderSize());
+ pad->SetBorderMode(minor->GetBorderMode());
+ minor->TAttLine::Copy((TAttLine&)*pad);
+ minor->TAttFill::Copy((TAttFill&)*pad);
+ minor->TAttPad::Copy((TAttPad&)*pad);
+
+ TIter next(minor->GetListOfPrimitives());
+ gROOT->SetSelectedPad(pad);
+ while ((obj=next())) {
+ clone = obj->Clone();
+ obj->Draw("SAME");
+ pad->GetListOfPrimitives()->Add(clone,obj->GetDrawOption());
+ }
+ pad->Modified();
+ pad->Update();
+ major->cd(padNum);
+ pad->Draw();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void CS_Diagnosis(){
+
+ auto majorCanv = new TCanvas("CompareCanv","CompareCanv",1500,1500);
+ majorCanv->Divide(3,3);
+ canclone(majorCanv, 1, "c_peakFits_110_112");
+ canclone(majorCanv, 2, "c_peakFits_115_118");
+ canclone(majorCanv, 3, "c_peakFits_120_122");
+ canclone(majorCanv, 4, "c_peakFits_125_128");
+ canclone(majorCanv, 5, "c_peakFits_130_132");
+ canclone(majorCanv, 6, "c_peakFits_135_138");
+ canclone(majorCanv, 7, "c_peakFits_140_142");
+ canclone(majorCanv, 8, "c_peakFits_145_148");
+ canclone(majorCanv, 9, "c_peakFits_150_152");
+
+}
+
////////////////////////////////////////////////////////////////////////////////
void CS(){
/* Overload function */
- cout << " Inputs:\n Experimental...\n\t - Energy of state\n\t - Spin of state" << endl;
- cout << " Theory...\n\t - Orbital l\n\t - Orbital j\n\t - Orbital n\n\n" << endl;
-
- cout << " 0f5/2 | ----- | --?-- | l=3, j=2.5, n=0" << endl;
- cout << " | | |" << endl;
- cout << " 1p1/2 | ----- | --?-- | l=1, j=0.5, n=1" << endl;
- cout << " 1p3/2 | ----- | --?-- | l=1, j=1.5, n=1" << endl;
- cout << " 0f7/2 | ----- | ===== |" << endl;
- cout << " | | |" << endl;
- cout << " 0d3/2 | ----- | ===== |" << endl;
- cout << " 1s1/2 | --x-- | ===== |" << endl;
- cout << " | | |" << endl;
- cout << " | p | n |" << endl;
+/*
+ cout << " Inputs:\n Experimental...\n\t - Energy of state\n\t - Spin of state" << endl;
+ cout << " Theory...\n\t - Orbital l\n\t - Orbital j\n\t - Orbital n\n\n" << endl;
+
+ cout << " 0f5/2 | ----- | --?-- | l=3, j=2.5, n=0" << endl;
+ cout << " | | |" << endl;
+ cout << " 1p1/2 | ----- | --?-- | l=1, j=0.5, n=1" << endl;
+ cout << " 1p3/2 | ----- | --?-- | l=1, j=1.5, n=1" << endl;
+ cout << " 0f7/2 | ----- | ===== |" << endl;
+ cout << " | | |" << endl;
+ cout << " 0d3/2 | ----- | ===== |" << endl;
+ cout << " 1s1/2 | --x-- | ===== |" << endl;
+ cout << " | | |" << endl;
+ cout << " | p | n |" << endl;
+*/
+
+ cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
+ cout << "---- 0.143, p3/2 = CS(0.143, 2, 1, 1.5) "<< endl;
+ cout << "---- 0.279, p3/2 = CS(0.279, 2, 1, 1.5) "<< endl;
+ cout << "---- 0.728, f7/2 = CS(0.728, 3, 3, 3.5) "<< endl;
+ cout << "---- 0.968, p1/2 = CS(0.968, 0, 1, 0.5) "<< endl;
+ cout << "---- 1.410, p3/2 = CS(1.410, 1, 1, 1.5) "<< endl;
+ cout << "---- 1.981, p3/2 = CS(1.981, 1, 1, 0.5) "<< endl;
+ cout << "---- 2.410, p3/2 = CS(2.410, 0, 1, 0.5) "<< endl;
+ cout << "---- 3.2 , f7/2 = CS(3.2 , 3, 3, 3.5) "<< endl;
+ cout << "---- 3.6 , f5/2 = CS(3.6 , 3, 3, 2.5) "<< endl;
+ cout << "---- 3.8 , f5/2 = CS(3.8 , 3, 3, 2.5) "<< endl;
+ cout << "---- 4.1 , f5/2 = CS(4.1 , 3, 3, 2.5) "<< endl;
+ cout << "---- 4.4 , f5/2 = CS(4.4 , 3, 3, 2.5) "<< endl;
+
}
@@ -49,8 +122,11 @@ void CS(double Energy, double Spin, double spdf, double angmom){
// p3/2 -> 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.8 +- 1.3
-
+ //double ElasticNorm = 5.8, ElasticNormErr = 1.3; // DeuteronNorm in elastics, 5.8 +- 1.3
+ //double ElasticNorm = 2.70, ElasticNormErr = 0.46; // For thicker
+ //double ElasticNorm = 2.8, ElasticNormErr = 0.7; // DeuteronNorm in elastics, 2.8 +- 0.7
+ double ElasticNorm = 2.363, ElasticNormErr = 0.000; // DeuteronNorm in elastics, reconstructed from determined thickness
+
orbitalname.clear();
orbital.clear();
if(spdf==1){
@@ -61,7 +137,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
orbitalname="p_{3/2}";
orbital="p32";
} else {
- orbitalname="?????";
+ orbitalname;
orbital="???";
}
} else if (spdf==3){
@@ -82,7 +158,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){
/* Reduce by factor of 10,000 */
ElasticNorm /= 10000.;
+ //ElasticNorm /= 1000.;
ElasticNormErr /= 10000.;
+ //ElasticNormErr /= 1000.;
double nodes;
if(spdf==1){
@@ -122,12 +200,18 @@ void CS(double Energy, double Spin, double spdf, double angmom){
}
/* Solid Angle (from simulation) */
- //auto file = new TFile("../SolidAngle_HistFile_15Feb_47Kdp.root");
- auto file = new TFile("../SolidAngle_HistFile_17May22_47Kdp_0143.root");
- //auto file = new TFile("../SolidAngle_HistFile_06May_WideStripMatching_LargeRun.root");
+ /* 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_HistFile_New.root");
+ /* ADD OPTION TO CHANGE SOLID ANGLE FILE DEPENDING ON PEAK!!!!*/
TH1F* SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_Lab_MG");
TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
- SolidAngle->Draw();
+ SolidAngle->Draw("HIST");
+ SolidAngle->GetXaxis()->SetRangeUser(100.,160.);
/* (canvas deleted after Area/SA calculation) */
/* Area of experimental peaks */
@@ -155,6 +239,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){
int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001);
int binmax = SolidAngle->FindBin(areaArray[i][4]-0.0001);
+//cout << "binmin " << binmin << " to binmax " << binmax << endl;
+
anglecentres.push_back(((areaArray[i][4]-areaArray[i][3])/2.)+areaArray[i][3]);
anglewidth.push_back(areaArray[i][4]-areaArray[i][3]);
@@ -202,7 +288,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
<< endl;
}
}
- delete c_SolidAngle;
+ //delete c_SolidAngle;
/* Graph of Area/Solid Angle*/
TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
@@ -349,6 +435,13 @@ void CS(double Energy, double Spin, double spdf, double angmom){
c_Chi2Min->SaveAs(savestring1.c_str());
c_Chi2Min->SaveAs(savestring2.c_str());
+ cout << YELLOW
+ << " ----- C2S = (2J+1)S = (2*" << (int)Spin << "+1)S = ("
+ << (int)((2.*Spin)+1.) << ")S = "
+ << ((2.*Spin)+1.) * globalS << " ----- "
+ << RESET << endl;
+
+
//delete c_AoSA;
//delete c_dSdO;
}
@@ -358,19 +451,23 @@ vector<vector<double>> GetExpDiffCross(double Energy){
cout << "========================================================" << endl;
vector<vector<double>> AllPeaks_OneGate;
vector<vector<double>> OnePeak_AllGates;
- int numbins = 10;
- double x[numbins], y[numbins];
+ /****CHANGE ANGLE GATING****/
+ int numAngleBins = 20;
+ double widthAngleBins = 2.5;
+ double firstAngle = 105.;
+ /***************************/
+ double x[numAngleBins], y[numAngleBins];
//TList* list = new TList();
/* Determine scaling factor for PhaseSpace */
TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
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.);
+ TH1F* baseEx = PullThetaLabHist(0,firstAngle,widthAngleBins);
+ TH1F* basePS = PullPhaseSpaceHist(0,firstAngle,widthAngleBins);
+ for(int i=1; i<numAngleBins;i++){
+ TH1F* addEx = PullThetaLabHist(i,firstAngle,widthAngleBins); baseEx->Add(addEx,1.);
+ TH1F* addPS = PullPhaseSpaceHist(i,firstAngle,widthAngleBins); basePS->Add(addPS,1.);
}
/* Subtract flat background equal to smallest bin in range */
@@ -385,7 +482,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Begin scaling within range, track changes */
basePS->Scale(0.1);
trackScale = 0.1;
- int numbinsScale = baseEx->GetNbinsX();
+ int numAngleBinsScale = baseEx->GetNbinsX();
int nbinlow = basePS->FindBin(4.); int nbinhigh = basePS->FindBin(8.0);
for(int b=nbinlow; b<nbinhigh; b++){
if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
@@ -405,15 +502,14 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
// TEMPORARY!!! REMOVE LAST THREE BINS ON HIGH ENERGY STATES!!!
- if(means[indexE] > 3.0){numbins-=3;}
+ if(means[indexE] > 3.0){numAngleBins-=3;}
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- for(int i=0; i<numbins;i++){
- double bin = 5.;
- double min = 105. + (i*bin);
- double max = min + bin;
+ for(int i=0; i<numAngleBins;i++){
+ double min = firstAngle + (i*widthAngleBins);
+ double max = min + widthAngleBins;
cout << "===================================" << endl;
cout << "min: " << min << " max: " << max << endl;
@@ -424,8 +520,8 @@ 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.);
+ TH1F* gate = PullThetaLabHist(i,firstAngle,widthAngleBins);
+ TH1F* pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins);
/* Scale the Phase Space at this angle... */
/* ... for all angles together */
@@ -438,8 +534,8 @@ vector<vector<double>> GetExpDiffCross(double Energy){
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++){
+ int numAngleBins = gate->GetNbinsX();
+ for(int b=0; b<numAngleBins; b++){
if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
" > " << gate->GetBinContent(b);
}
@@ -471,14 +567,22 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* 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);
+
+ //AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0); cout << "!!!!!!! NO SLIDING SHIFT!!!!!" << endl;
+ AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0); cout << "!!!!!!! WITH A VARIABLE SLIDING SHIFT!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
+ //double slideshift = 0.00218*(((max-min)/2.)+min) - 0.29645; AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, slideshift); cout << "!!!!!!! WITH A FIXED SLIDING SHIFT!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
/* Write PS-subtracted spectrum to list */
//list->Add(gate);
//list->Add(c_peakFits);
+ gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
string filename = "./CS2_Figures/";
filename += tmp2 + ".root";
c_peakFits->SaveAs(filename.c_str());
+ //auto tempfile = new TFile(filename.c_str(),"UPDATE"); //Reopen newly made file
+ //auto templist = new TList();
+ //templist->Add();
+
/* Check correct OneGate vector is selected */
cout << "area of " << means[indexE] << " = "
@@ -505,7 +609,14 @@ vector<vector<double>> GetExpDiffCross(double Energy){
////////////////////////////////////////////////////////////////////////////////
TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
//TFile* file = new TFile("GateThetaLabHistograms.root","READ");
- TFile* file = new TFile("GateThetaLabHistograms_ReadMe.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_20May22.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_0p05BinWidth.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_ReadMe.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_26May22_v2.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_26May22_v2_bin0p02.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_2p5degAngles_26May22v2.root","READ");
+ //TFile* file = new TFile("GateThetaLabHistograms_11Apr22_20angles.root","READ");
+ TFile* file = new TFile("GateThetaLabHistograms_11Jul22.root","READ");
string histname = "cThetaLabGate_"
+ to_string((int) (minTheta+(i*gatesize))) + "-"
+ to_string((int) (minTheta+((i+1)*gatesize)));
@@ -518,7 +629,9 @@ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
////////////////////////////////////////////////////////////////////////////////
TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize){
- TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_ReadMe.root","READ");
+ //TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_ReadMe.root","READ");
+ //TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_2p5degAngles_26May22v2.root","READ");
+ TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_11Jul22.root","READ");
string histname = "cPSpaceThetaLabGate_"
+ to_string((int) (minTheta+(i*gatesize))) + "-"
+ to_string((int) (minTheta+((i+1)*gatesize)));
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 8106fb61de1bf7c081152caa123552e5a6a6daf0..559e798365138c3479ee7cadc7e112d6c8229331 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -2,7 +2,6 @@
#include "NPReaction.h"
#include <string>
#include <sstream>
-
using namespace std;
TChain* chain=NULL ;
@@ -10,6 +9,7 @@ char cond[1000];
NPL::Reaction Cadp("47Ca(d,p)48Ca@355");
NPL::Reaction Scdp("47Sc(d,p)48Sc@355");
+NPL::Reaction K46dp("46K(d,p)47K@355");
NPL::Reaction Kdp("47K(d,p)48K@355");
NPL::Reaction Kdt("47K(d,t)46K@355");
@@ -22,8 +22,8 @@ NPL::Reaction Tidd("47Ti(d,d)47Ti@355");
NPL::Reaction Ti12C12C("47Ti(12C,12C)47Ti@355");
void KnownLines_Ex(bool isVertical, double rangemin, double rangemax, Style_t lType, Color_t lColour);
-void AddGammaLinesMG(TH1F* hist, double particle, double ymax);
-void AddPlacedGammasMG(TH1F* hist, double ymax);
+void AddGammaLines(TH1F* hist, double particle, double ymax);
+void AddPlacedGammas(TH1F* hist, double ymax);
double tCentre;
double tRange;
@@ -48,23 +48,83 @@ TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventLi
void LoadChain47Kdp(){
vector<string> files;
+// files.push_back("../../../Outputs/Analysis/OriginalValues_ptI.root");
+// files.push_back("../../../Outputs/Analysis/OriginalValues_ptII.root");
+// files.push_back("../../../Outputs/Analysis/OriginalValues_ptIII.root");
+
//files.push_back("../../../Outputs/Analysis/47Kdp_08Nov_PartI.root");
//files.push_back("../../../Outputs/Analysis/47Kdp_08Nov_PartII.root");
- //files.push_back("../../../Outputs/Analysis/47Kdp_08Apr_PartI.root");
- //files.push_back("../../../Outputs/Analysis/47Kdp_08Apr_PartII.root");
+ /* With thresholds, strip matching, and bad strips out */
+ //files.push_back("../../../Outputs/Analysis/47Kdp_11Apr22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_11Apr22_PartII.root");
+
+ /* As above, with thicker target! */
+ //files.push_back("../../../Outputs/Analysis/47Kdp_26May22_v2_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_26May22_v2_PartII.root");
+
+ /* Fix Z=0 */
+ //files.push_back("../../../Outputs/Analysis/47Kdp_27May22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_27May22_PartII.root");
+
+ /* Not fixed, but Z near 0 */
+ //files.push_back("../../../Outputs/Analysis/47Kdp_01Jul22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_01Jul22_PartII.root");
+
+ /* New target thickness analysis */
+ files.push_back("../../../Outputs/Analysis/47Kdp_11Jul22_PartI.root");
+ files.push_back("../../../Outputs/Analysis/47Kdp_11Jul22_PartII.root");
- /* With thresholds, strip mathcing, and bad strips out */
- files.push_back("../../../Outputs/Analysis/47Kdp_11Apr22_PartI.root");
- files.push_back("../../../Outputs/Analysis/47Kdp_11Apr22_PartII.root");
chain = Chain("PhysicsTree",files,true);
}
void LoadChain47Kdt(){
vector<string> files;
- files.push_back("../../../Outputs/Analysis/47Kdt_13May22_PartI.root");
- files.push_back("../../../Outputs/Analysis/47Kdt_13May22_PartII.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_13May22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_13May22_PartII.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_25May22_TestingDiffTimeCalib_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_25May22_TestingDiffTimeCalib_PartII.root");
+
+ /* Offset MM1 timing by -5 */
+ //files.push_back("../../../Outputs/Analysis/47Kdt_22Jun22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_22Jun22_PartII.root");
+
+ /* As above, with 01Jul22 XYZT */
+ //files.push_back("../../../Outputs/Analysis/47Kdt_22Jun22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_22Jun22_PartII.root");
+
+ /* As above, with 11Jul22 XYZT */
+ //files.push_back("../../../Outputs/Analysis/47Kdt_11Jul22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_11Jul22_PartII.root");
+
+ files.push_back("../../../Outputs/Analysis/47Kdt_14Jul22_PartI.root");
+ files.push_back("../../../Outputs/Analysis/47Kdt_14Jul22_PartII.root");
+
+
+ //files.push_back("../../../Outputs/Analysis/TESTING_PartII.root");
+
+ chain = Chain("PhysicsTree",files,true);
+}
+
+void LoadChain47Kdd(){
+ vector<string> files;
+ //files.push_back("../../../Outputs/Analysis/47Kdd_08Nov_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdd_08Nov_PartII.root");
+
+ files.push_back("../../../Outputs/Analysis/47Kdd_11Jul22_PartI.root");
+ files.push_back("../../../Outputs/Analysis/47Kdd_11Jul22_PartII.root");
+
+ chain = Chain("PhysicsTree",files,true);
+}
+
+void LoadChain47Kpp(){
+ vector<string> files;
+ //files.push_back("../../../Outputs/Analysis/47Kpp_08Nov_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kpp_08Nov_PartII.root");
+
+ files.push_back("../../../Outputs/Analysis/47Kpp_11Jul22_PartI.root");
+ files.push_back("../../../Outputs/Analysis/47Kpp_11Jul22_PartII.root");
chain = Chain("PhysicsTree",files,true);
}
@@ -100,6 +160,16 @@ void DrawParticleStates(TCanvas* canvas){
gs->SetLineColor(kGreen);
gs->SetLineStyle(7);
gs->Draw();
+
+ TLine** lines = new TLine*[numPeaks];
+ for(int j = 0; j < numPeaks; j++){
+ lines[j] = new TLine(means[j],0.0,means[j],max);
+ lines[j]->SetLineStyle(8);
+ lines[j]->SetLineColor(kGray);
+ lines[j]->Draw();
+ }
+
+/*
TLine *l0143 = new TLine(0.143, 0.0, 0.143, max);
l0143->SetLineStyle(kDashed);
l0143->Draw();
@@ -139,6 +209,8 @@ void DrawParticleStates(TCanvas* canvas){
TLine *l4510 = new TLine(4.51, 0.0, 4.51, max);
l4510->SetLineStyle(kDotted);
l4510->Draw("same");
+
+*/
}
void plot_kine(NPL::Reaction r, double Ex,Color_t c,int w, int s){
@@ -171,10 +243,7 @@ void AddTiStates(double E){
string timegate; /* defined by choice of dp or dt */
string det_gate; /* defined by choice of dp or dt */
-//string mg_gate = "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8";
-//string mm_gate = "MUST2.TelescopeNumber>0";
-//string mm_lt5_gate = "MUST2.TelescopeNumber<5";
-//string mm_eq5_gate = "MUST2.TelescopeNumber==5";
+string reactionName; /* defined by choice of dp or dt */
string exclBmDcy = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003";
void Draw_1DGamma(){
@@ -207,6 +276,70 @@ void Draw_1DGamma_DetGate(){
Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
}
+void Draw_1DGamma_DetGate_x10x100(){
+ string gate = timegate
+ + " && " + det_gate
+ + " && " + exclBmDcy;
+
+ TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
+ gStyle->SetOptStat(0);
+ chain->Draw("AddBack_EDC>>Eg(1250,0,5)",gate.c_str(),"");
+ TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
+ Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ Eg->GetYaxis()->SetTitle("Counts / 0.004 MeV");
+
+ for(int b = 126; b<501; b++){
+ Eg->SetBinContent(b,(Eg->GetBinContent(b) * 10.));
+ }
+ for(int b = 501; b<1250; b++){
+ Eg->SetBinContent(b,(Eg->GetBinContent(b) * 100.));
+ }
+
+ Eg->Draw();
+
+ TLine *x10 = new TLine(0.5, 0.0, 0.5, 3000.);
+ x10->SetLineColor(kBlack); x10->SetLineStyle(7);
+ x10->Draw("same");
+ TLine *x100 = new TLine(2.0, 0.0, 2.0, 3000.);
+ x100->SetLineColor(kBlack); x100->SetLineStyle(7);
+ x100->Draw("same");
+
+}
+
+void Draw_1DGamma_DetGate_SplitCanv(){
+ string gate = timegate
+ + " && " + det_gate
+ + " && " + exclBmDcy;
+
+ TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
+ gStyle->SetOptStat(0);
+ cEg->Divide(2,1,0.005,0.005,0);
+ cEg->cd(1);
+ gStyle->SetPadLeftMargin(0.0);
+ gStyle->SetPadRightMargin(0.0);
+ gPad->SetTickx();
+ gPad->SetTicky();
+ cEg->cd(2);
+ gStyle->SetPadLeftMargin(0.0);
+ gStyle->SetPadRightMargin(0.0);
+ gPad->SetTickx();
+ gPad->SetTicky();
+
+ cEg->cd(1);
+ chain->Draw("AddBack_EDC>>Eg(5000,0,5)",gate.c_str(),"");
+ TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
+ Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
+ Eg->GetXaxis()->SetRangeUser(0.,2.);
+
+ cEg->cd(2);
+ chain->Draw("AddBack_EDC>>Eg2(1250,0,5)",gate.c_str(),"");
+ TH1F* Eg2 = (TH1F*) gDirectory->Get("Eg2");
+ Eg2->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ Eg2->GetYaxis()->SetTitle("Counts / 0.004 MeV");
+ Eg2->GetXaxis()->SetRangeUser(2.,4.5);
+}
+
void Load_1DGamma_MG(){
TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
TH1F *hEgMG = new TH1F("hEg","Loaded 1D Gamma Spectrum, MG gated",600,-15,15);
@@ -261,14 +394,35 @@ void Draw_2DParticleGamma(){
+ " && Ex@.size()==1";
TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000);
- chain->Draw("AddBack_EDC:Ex>>ExEg(600,-15,15,2500,0,5)", gate.c_str(), "colz");
+ gStyle->SetOptStat(0);
+ //chain->Draw("AddBack_EDC:Ex>>ExEg(600,-15,15,2500,0,5)", gate.c_str(), "colz");
+ chain->Draw("Ex:AddBack_EDC>>ExEg(5000,0,5,3000,-15,15)", gate.c_str(), "");
TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg");
- ExEg->SetTitle("Ex-Egamma");
- ExEg->GetXaxis()->SetTitle("Ex [MeV]");
- ExEg->GetYaxis()->SetTitle("Eg [MeV]");
- TLine *XeqY = new TLine(0,0,9,9);
- XeqY->SetLineColor(kRed);
- XeqY->Draw();
+ ExEg->SetTitle("");
+ //ExEg->GetXaxis()->SetTitle("Ex [MeV]");
+ ExEg->GetYaxis()->SetTitle("Ex [MeV]");
+ //ExEg->GetYaxis()->SetTitle("Eg [MeV]");
+ ExEg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ ExEg->GetYaxis()->SetRangeUser(-1.0,7.0);
+ ExEg->Draw();
+ TLine *XeqY = new TLine(0,0,5,5);
+ XeqY->SetLineColor(kRed);
+ XeqY->Draw("same");
+ TLine *Sn = new TLine(0,4.644,4.644,4.644);
+ Sn->SetLineColor(kRed);Sn->SetLineStyle(2);
+ Sn->Draw("same");
+ TLatex *TSn = new TLatex(.5,.5,"S_{n}");
+ TSn->SetTextColor(kRed);
+ TSn->SetTextSize(0.05);
+ TSn->SetX(2.50);
+ TSn->SetY(4.90);
+ TSn->Draw("same");
+ TLatex *Texeg = new TLatex(.5,.5,"Ex = E_{#gamma}");
+ Texeg->SetTextColor(kRed);
+ Texeg->SetTextSize(0.05);
+ Texeg->SetX(2.35);
+ Texeg->SetY(1.50);
+ Texeg->Draw("same");
}
void Load_2DGammaGamma(){
@@ -281,15 +435,20 @@ void Load_2DGammaGamma(){
}
void Draw_2DGammaGamma_ExcludeBeam(){
+ string gate = timegate
+ + " && " + exclBmDcy;
+
TCanvas *cEgEg = new TCanvas("cEgEg","cEgEg",1000,1000);
- chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(999,0.005,5,999,0.005,5)","abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003","colz");
+ chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(2500,0,5,2500,0,5)",gate.c_str(),"colz");
TH1F* EgEg = (TH1F*) gDirectory->Get("EgEg");
- chain->Draw("AddBack_EDC2:AddBack_EDC>>EgEg2(999,0.005,5,999,0.005,5)","abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003","colz");
+ chain->Draw("AddBack_EDC2:AddBack_EDC>>EgEg2(2500,0,5,2500,0,5)",gate.c_str(),"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->GetXaxis()->SetTitle("Eg [Counts / 0.002 MeV]");
+ EgEg->GetYaxis()->SetTitle("Eg [Counts / 0.002 MeV]");
+ EgEg->GetXaxis()->SetRangeUser(0.005,5.0);
+ EgEg->GetYaxis()->SetRangeUser(0.005,5.0);
EgEg->Draw("colz");
TLine *XeqY = new TLine(0,0,5,5);
XeqY->SetLineColor(kRed);
@@ -315,14 +474,16 @@ void gg(){
void Draw_2DGammaGamma_TimeGated(){
string gate = timegate;
TCanvas *cEgEg = new TCanvas("cEgEg","cEgEg",1000,1000);
- chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(999,0.005,5,999,0.005,5)",gate.c_str(),"colz");
+ chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(2500,0,5,2500,0,5)",gate.c_str(),"colz");
TH1F* EgEg = (TH1F*) gDirectory->Get("EgEg");
- chain->Draw("AddBack_EDC2:AddBack_EDC>>EgEg2(999,0.005,5,999,0.005,5)",gate.c_str(),"colz");
+ chain->Draw("AddBack_EDC2:AddBack_EDC>>EgEg2(2500,0,5,2500,0,5)",gate.c_str(),"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->GetXaxis()->SetTitle("Eg [Counts / 0.002 MeV]");
+ EgEg->GetYaxis()->SetTitle("Eg [Counts / 0.002 MeV]");
+ EgEg->GetXaxis()->SetRangeUser(0.005,5.0);
+ EgEg->GetYaxis()->SetRangeUser(0.005,5.0);
EgEg->Draw("colz");
TLine *XeqY = new TLine(0,0,5,5);
XeqY->SetLineColor(kRed);
@@ -425,8 +586,6 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg, double
DrawParticleStates(cEx_Gate);
}
-
-
void GateParticle_SeeGamma(double particle, double width){
gStyle->SetOptStat("nemMrRi");
@@ -452,7 +611,8 @@ void GateParticle_SeeGamma(double particle, double width){
EgGate->Draw();
limit->Draw();
- AddGammaLinesMG(EgGate, particle, cEg_Gate->GetUymax());
+
+ AddGammaLines(EgGate, particle, cEg_Gate->GetUymax());
}
@@ -1082,13 +1242,27 @@ void ExThetaLab(double gamma, double width){
void ELabThetaLab(){
TCanvas *cELabTLaab = new TCanvas("cELabTLab","cELabTLab",1000,1000);
gStyle->SetOptStat(0);
- chain->Draw("ELab:ThetaLab>>hKine(360,0,180,500,0,10)","abs(T_MUGAST_VAMOS-2700)<400","col");
+
+ string gate = timegate
+ + " && " + det_gate;
+ if(reactionName=="47K(d,t)"){
+ gate = gate + " && cutTritons";
+ }
+
+
+cout << "test" << endl;
+ chain->Draw("ELab:ThetaLab>>hKine(360,0,180,500,0,10)",gate.c_str(),"col");
+cout << "test2" << endl;
TH2F* hKine = (TH2F*) gDirectory->Get("hKine");
+cout << "test3" << endl;
hKine->SetTitle("");
hKine->GetXaxis()->SetTitle("#theta_{lab} [deg]");
hKine->GetYaxis()->SetTitle("E_{lab} [MeV]");
plot_kine(Kdt, 0.000, kBlack, 2, 1);
+ plot_kine(Kdd, 0.000, kRed, 2, 1);
+ plot_kine(Kpp, 0.000, kRed, 2, 1);
+
plot_kine(Kdp, 0.000, kBlack, 2, 1);
plot_kine(Kdp, 4.644, kBlack, 2, 1);
@@ -1111,6 +1285,7 @@ void ELabThetaLab(){
plot_kine(Cadp, 0.000, kRed, 2, 1);
plot_kine(Tidp, 0.000, kBlue, 2, 1);
plot_kine(Scdp, 0.000, kGreen, 2, 1);
+ plot_kine(K46dp, 0.000, kViolet, 2, 1);
//plot_kine(Tidp, 5.652, kBlack, 2, 6); //strongest populated state according to PDBarnes(1965)
}
@@ -1155,7 +1330,6 @@ void MM5_ExThetaLab(){
"colz");
}
-
void thickness(){
std::ifstream infile("thicknessTheory4.txt");
@@ -1508,12 +1682,8 @@ double expPyErr2[15]={0.5,
}
+/*
void ExThetaAnalysis(double gamma, double width, int version){
-// int version;
-// bool running = 1;
-
-// cout << "Constructing Ex:ThetaLab..." << endl;
-
string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ to_string(gamma) + ") < " + to_string(width);
@@ -1568,12 +1738,32 @@ void ExThetaAnalysis(double gamma, double width, int version){
}//else{running=0;}
// }
-
}
+*/
+void ExTheta_Analysis(double gamma, double width){
+ string gate = timegate
+ + " && " + det_gate
+ + " && Ex@.size()==1"
+ + " && abs(AddBack_EDC-"
+ + to_string(gamma) + ")<"
+ + to_string(width)
+ ;
+ string gateLow = gate + " && abs(ThetaLab-117.5)<12.5";
+ string gateHigh = gate + " && abs(ThetaLab-142.5)<12.5";
+ TCanvas *diagnosis_ExTheta = new TCanvas("diagnosis_ExTheta","diagnosis_ExTheta",1000,1000);
+ chain->Draw("Ex>>thetaHistLow(200,-1,9)", gateLow.c_str(),"");
+ TH2F* thetaHistLow = (TH2F*) gDirectory->Get("thetaHistLow");
+ thetaHistLow->SetLineColor(kBlue);
+ chain->Draw("Ex>>thetaHistHigh(200,-1,9)", gateHigh.c_str(),"");
+ TH2F* thetaHistHigh = (TH2F*) gDirectory->Get("thetaHistHigh");
+ thetaHistHigh->SetLineColor(kRed);
+ thetaHistLow->Draw();
+ thetaHistHigh->Draw("SAME");
+}
void ExMugast_ForPoster(){
@@ -1685,12 +1875,16 @@ void AGATA_efficiency(double Energy_keV){
}
void ElasticsGate(double EMin, double EMax){
- string gates = "abs(T_MUGAST_VAMOS-2700)<400 && MUST2.TelescopeNumber==5 && ELab > "
- + to_string(EMin)
- + " && ELab < "
- + to_string(EMax);
+ double width = EMax-EMin;
+ double centre = EMin+(0.5*width);
+// string gates = "abs(T_MUGAST_VAMOS-2700)<400 && MUST2.TelescopeNumber==5 && abs(ELab - "
+ string gate = timegate + "&&" + det_gate
+ + "&& abs(ELab - "
+ + to_string(centre)
+ + ")< "
+ + to_string(width);
- chain->Draw("ThetaLab>>hist(80,50,90)", gates.c_str(), "");
+ chain->Draw("ThetaLab>>hist(160,50,90)", gate.c_str(), "");
}
void GateThetaCM(double minTheta, double maxTheta, double binsize){
@@ -1766,7 +1960,10 @@ void GateThetaLab_AllOverlaid(){
}
void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize){
- string core = "abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > ";
+ string core = timegate
+ + " && " + det_gate
+ + " && Ex@.size()==1 && ThetaLab > ";
+// string core = "abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > ";
string ytitle = "Counts / " + to_string(binsize) + " MeV";
double gatesize = (finishTheta-startTheta)/numGates;
TList* list = new TList();
@@ -1797,13 +1994,135 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
file->ls();
}
+void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize, double gateGammaE, double gateGammaWdth){
+ string core = timegate
+ + " && " + det_gate
+ + " && abs(AddBack_EDC-" + to_string(gateGammaE) + ")<" + to_string(gateGammaWdth)
+ + " && Ex@.size()==1 && ThetaLab > ";
+// string core = "abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0 && abs(AddBack_EDC-" + to_string(gateGammaE) + ")<" + to_string(gateGammaWdth) + " && 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((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize));
+ string gating = core
+ + to_string(minTheta)
+ + " && ThetaLab < "
+ + to_string(minTheta+gatesize);
+ string histname = "cThetaLabGate_" + to_string((int) minTheta) + "-" + to_string((int) (minTheta+gatesize));
+ string draw = "Ex>>" + histname + "(" + to_string(30.0/binsize) + ",-15,15)";
+
+ 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->Sumw2();
+ Ex_ThetaLabGate->SetTitle(title.c_str());
+ list->Add(Ex_ThetaLabGate);
+ delete cEx_ThetaLabGate;
+ }
+
+ TFile* file = new TFile("GateThetaLabHistograms_GammaGate.root","RECREATE");
+ list->Write("GateThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
+
+void CompareThetaLabGatesFromFile(){
+ TFile* oldF = new TFile("GateThetaLabHistograms_11Apr22_20angles.root","READ");
+ TList* oldL = (TList*) oldF->FindObjectAny("GateThetaLabHistograms");
+
+ TFile* newF = new TFile("GateThetaLabHistograms_11Jul22.root","READ");
+ TList* newL = (TList*) newF->FindObjectAny("GateThetaLabHistograms");
+
+ double minTheta=105.;
+ double gatesize=2.5;
+ int numGates=20;
+
+ TCanvas* cOldNew = new TCanvas("cOldNew","cOldNew",2000,1000);
+ gStyle->SetOptStat(0);
+ cOldNew->Divide(3,2,0.01,0.01,0);
+ /*SetPad(xlow, ylow, xupp, yupp) */
+ cOldNew->cd(1)->SetPad(0., 0.5, 0.5, 1.0);
+ cOldNew->cd(2)->SetPad(0.5, 0.5, 0.75, 1.0 );
+ cOldNew->cd(3)->SetPad(0.75, 0.5, 1, 1.0 );
+ cOldNew->cd(4)->SetPad(0., 0., 0.5, 0.5);
+ cOldNew->cd(5)->SetPad(0.5, 0., 0.75, 0.5 );
+ cOldNew->cd(6)->SetPad(0.75, 0., 1, 0.5 );
+
+ for (int i=0; i<numGates; i++){
+ string hname = "cThetaLabGate_"
+ + to_string((int) (minTheta+(i*gatesize)))
+ + "-"
+ + to_string((int) (minTheta+((i+1.)*gatesize)));
+
+ cout << " Looking for " << hname << endl;
+
+ TH1F* oldH = (TH1F*) oldL->FindObject(hname.c_str());
+ TH1F* newH = (TH1F*) newL->FindObject(hname.c_str());
+
+ //would be faster to set some of these on just the last one but speed not issue here
+ cOldNew->cd(1);
+ oldH->SetLineColor(30+i);
+ oldH->Rebin(2);
+ oldH->GetYaxis()->SetTitle("Counts / 0.100000 MeV");
+ oldH->GetXaxis()->SetRangeUser(-1,7);
+ oldH->GetYaxis()->SetRangeUser(0,400);
+ oldH->SetTitle("08Nov22: 1.3008 um CD2");
+ oldH->Draw("HIST SAME");
+
+ cOldNew->cd(2);
+ TH1F* oldH1 = (TH1F*) oldH->Clone();
+ oldH1->SetName("oldH1");
+ oldH1->GetXaxis()->SetRangeUser(-1,1);
+ oldH1->SetTitle("");
+ oldH1->Draw("HIST SAME");
+
+ cOldNew->cd(3);
+ TH1F* oldH2 = (TH1F*) oldH->Clone();
+ oldH2->GetYaxis()->SetRangeUser(0,250);
+ oldH2->SetName("oldH2");
+ oldH2->GetXaxis()->SetRangeUser(1.5,2.5);
+ oldH2->SetTitle("");
+ oldH2->Draw("HIST SAME");
+
+ cOldNew->cd(4);
+ newH->GetXaxis()->SetRangeUser(-1,7);
+ newH->GetYaxis()->SetRangeUser(0,400);
+ newH->SetLineColor(30+i);
+ newH->SetTitle("11Jul22: 2.8798 um CD2");
+ newH->Draw("HIST SAME");
+
+ cOldNew->cd(5);
+ TH1F* newH1 = (TH1F*) newH->Clone();
+ newH1->SetName("newH1");
+ newH1->GetXaxis()->SetRangeUser(-1,1);
+ newH1->SetTitle("");
+ newH1->Draw("HIST SAME");
+
+ cOldNew->cd(6);
+ TH1F* newH2 = (TH1F*) newH->Clone();
+ newH2->GetYaxis()->SetRangeUser(0,250);
+ newH2->SetName("newH2");
+ newH2->GetXaxis()->SetRangeUser(1.5,2.5);
+ newH2->SetTitle("");
+ newH2->Draw("HIST SAME");
+
+ }
+
+
+}
+
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");
+ //TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_02Mar_47Kdp_PhaseSpace.root","READ");
+ TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_PhaseSpace_11Jul22.root","READ");
TTree* PSTree = (TTree*) psfile->FindObjectAny("PhysicsTree");
for (int i=0; i<numGates; i++){
@@ -1874,7 +2193,6 @@ void GammaSub_Actual_ExcludeBeamDecay(){
Eg->Draw();
}
-
/*
void gg(){
static vector<double> *AddBack_EDC;
@@ -2047,18 +2365,20 @@ cout << "THIS IS OLD!!!! UPDATE WITH THE BEAM EXCLUSION!!!" << endl;
}//for i
}
-
void gg(){
- cout << "LOADING FILES: 47Kdp_11Apr22_PartI & II" << endl;
+ //cout << "LOADING FILES: 47Kdp_11Apr22_PartI & II" << endl;
+ cout << "LOADING FILES: 47Kdp_11Jul22_PartI & II" << endl;
auto h=new TH2F("gg","gg",1000,0,10,1000,0,10);
- auto DataFile = new TFile("../../../Outputs/Analysis/47Kdp_11Apr22_PartI.root", "READ");
+ //auto DataFile = new TFile("../../../Outputs/Analysis/47Kdp_11Apr22_PartI.root", "READ");
+ auto DataFile = new TFile("../../../Outputs/Analysis/47Kdp_11Jul22_PartI.root", "READ");
auto chain = (TTree*) DataFile->FindObjectAny("PhysicsTree");
ggLoad(chain, h);
auto h2=new TH2F("gg","gg",1000,0,10,1000,0,10);
- auto DataFile2 = new TFile("../../../Outputs/Analysis/47Kdp_11Apr22_PartII.root", "READ");
+ //auto DataFile2 = new TFile("../../../Outputs/Analysis/47Kdp_11Apr22_PartII.root", "READ");
+ auto DataFile2 = new TFile("../../../Outputs/Analysis/47Kdp_11Jul22_PartII.root", "READ");
auto chain2 = (TTree*) DataFile->FindObjectAny("PhysicsTree");
ggLoad(chain2, h2);
@@ -2141,10 +2461,7 @@ void ggGater(TH2F* h, double E, double gate){
}
-
-
void Figure_Eg_MG(){
-
chain->Draw("AddBack_EDC>>Eg(5000,0,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8");
TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
@@ -2187,3 +2504,191 @@ void Figure_Eg_MG(){
}
+void Figure_GateGamma_SeeParticle(double gamma, double width, double bg, double widthbg, double gammaBinWidth, double particleBinWidth){
+ string gating = timegate + "&&" + det_gate + " && Ex@.size()==1 && abs(AddBack_EDC-"
+ + to_string(gamma)
+ + ")<"
+ + to_string(width);
+ string bggate = timegate + "&&" + det_gate + " && Ex@.size()==1 && abs(AddBack_EDC-"
+ + to_string(bg)
+ + ")<"
+ + to_string(widthbg);
+
+ string gammagate = timegate
+ + " && " + det_gate
+ + " && " + exclBmDcy;
+
+ double ratio = width/widthbg;
+
+ ////////////////////////////////////////////////
+
+ TCanvas *cFig_GGSP = new TCanvas("cFig_GGSP","cFig_GGSP",1000,1000);
+ gStyle->SetOptStat(0);
+ cFig_GGSP->Divide(1,2,0.005,0.005,0);
+ cFig_GGSP->cd(1);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.01);
+ gPad->SetTickx();
+ gPad->SetTicky();
+ cFig_GGSP->cd(2);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.01);
+ gPad->SetTickx();
+ gPad->SetTicky();
+
+ cFig_GGSP->cd(2);
+ string draw1 = "Ex>>hEx(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(draw1.c_str(),gating.c_str(),"");
+ TH1F* ExGate = (TH1F*) gDirectory->Get("hEx");
+ ExGate->GetXaxis()->SetTitle("Ex [MeV]");
+ string nameExYaxis = "Counts / " + to_string(particleBinWidth) + " MeV";
+ ExGate->GetYaxis()->SetTitle(nameExYaxis.c_str());
+ ExGate->SetLineColor(kGreen);
+ ExGate->SetFillColor(kGreen);
+ ExGate->SetFillStyle(3154);
+ ExGate->SetTitle("");
+ ExGate->GetXaxis()->SetTitleSize(0.05);
+ ExGate->GetXaxis()->SetLabelSize(0.05);
+ ExGate->GetYaxis()->SetTitleSize(0.05);
+ ExGate->GetYaxis()->SetLabelSize(0.05);
+ ExGate->GetYaxis()->SetTitleOffset(0.5);
+ ExGate->GetYaxis()->SetNdivisions(520);
+ ExGate->GetXaxis()->SetRangeUser(-1.,7.);
+
+ string draw2 = "Ex>>hExBG(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(draw2.c_str(),bggate.c_str(),"same hist");
+ TH1F* ExBG = (TH1F*) gDirectory->Get("hExBG");
+ ExBG->Scale(ratio);
+ ExBG->SetLineColor(kRed);
+ ExBG->SetFillColor(kRed);
+ ExBG->SetFillStyle(3345);
+ ExBG->SetTitle("");
+ ExBG->GetXaxis()->SetRangeUser(-1.,7.);
+ ExBG->Draw("same hist");
+
+ cFig_GGSP->Update();
+ double maxEx = cFig_GGSP->cd(2)->GetUymax();
+ TLine *Sn = new TLine(4.644, 0.0, 4.644, maxEx);
+ Sn->SetLineColor(kBlack); Sn->SetLineStyle(7);
+ Sn->Draw();
+ TLine *gs = new TLine(0.000, 0.0, 0.000, maxEx);
+ gs->SetLineColor(kBlack); gs->SetLineStyle(1);
+ gs->Draw();
+
+ cFig_GGSP->cd(1);
+ string drawg = "AddBack_EDC>>hEg(" + to_string(5./gammaBinWidth) + ",0,5)";
+ chain->Draw(drawg.c_str(),gammagate.c_str(),"");
+ TH1F* hEg = (TH1F*) gDirectory->Get("hEg");
+ hEg->SetLineColor(kBlack);
+ string nameEgYaxis = "Counts / " + to_string(gammaBinWidth) + " MeV";
+ hEg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ hEg->GetYaxis()->SetTitle(nameEgYaxis.c_str());
+ double zoomMin = min((double) floor((gamma-width)*10)/10., (double) floor((bg-widthbg)*10)/10.);
+ double zoomMax = max((double) ceil((gamma+width)*10)/10., (double) ceil((bg+widthbg)*10)/10.);
+ hEg->GetXaxis()->SetRangeUser(zoomMin,zoomMax);
+ hEg->SetTitle("");
+ hEg->GetXaxis()->SetTitleSize(0.05);
+ hEg->GetXaxis()->SetLabelSize(0.05);
+ hEg->GetYaxis()->SetTitleSize(0.05);
+ hEg->GetYaxis()->SetLabelSize(0.05);
+ hEg->GetYaxis()->SetTitleOffset(0.5);
+ hEg->Draw();
+
+ cFig_GGSP->Update();
+ double max = cFig_GGSP->cd(1)->GetUymax();
+ TLine* gateL = new TLine(gamma-width,0,gamma-width,max);
+ gateL->SetLineColor(kGreen); gateL->SetLineStyle(10);
+ TLine* gateH = new TLine(gamma+width,0,gamma+width,max);
+ gateH->SetLineColor(kGreen); gateH->SetLineStyle(10);
+ TLine* bgL = new TLine(bg-widthbg,0,bg-widthbg,max);
+ bgL->SetLineColor(kRed); bgL->SetLineStyle(10);
+ TLine* bgH = new TLine(bg+widthbg,0,bg+widthbg,max);
+ bgH->SetLineColor(kRed); bgH->SetLineStyle(10);
+ gateL->Draw("SAME");
+ gateH->Draw("SAME");
+ bgL->Draw("SAME");
+ bgH->Draw("SAME");
+
+}
+
+
+
+
+
+
+
+
+
+
+
+void Figure_TopGamma_BottomParticle(double gammaBinWidth, double particleBinWidth){
+ string gating = timegate + "&&" + det_gate + " && Ex@.size()==1";
+ string gammagate = timegate
+ + " && " + det_gate
+ + " && " + exclBmDcy;
+
+ ////////////////////////////////////////////////
+
+ TCanvas *cFig_GGSP = new TCanvas("cFig_GGSP","cFig_GGSP",1000,1000);
+ gStyle->SetOptStat(0);
+ cFig_GGSP->Divide(1,2,0.005,0.005,0);
+ cFig_GGSP->cd(1);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.01);
+ gPad->SetTickx();
+ gPad->SetTicky();
+ cFig_GGSP->cd(2);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.01);
+ gPad->SetTickx();
+ gPad->SetTicky();
+
+ cFig_GGSP->cd(2);
+ string draw1 = "Ex>>hEx(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(draw1.c_str(),gating.c_str(),"");
+ TH1F* hEx = (TH1F*) gDirectory->Get("hEx");
+ hEx->GetXaxis()->SetTitle("Ex [MeV]");
+ string nameExYaxis = "Counts / " + to_string(particleBinWidth) + " MeV";
+ hEx->GetYaxis()->SetTitle(nameExYaxis.c_str());
+ //ExGate->SetLineColor(kGreen);
+ //ExGate->SetFillColor(kGreen);
+ //ExGate->SetFillStyle(3154);
+ hEx->SetTitle("");
+ hEx->GetXaxis()->SetTitleSize(0.05);
+ hEx->GetXaxis()->SetLabelSize(0.05);
+ hEx->GetYaxis()->SetTitleSize(0.05);
+ hEx->GetYaxis()->SetLabelSize(0.05);
+ hEx->GetYaxis()->SetTitleOffset(0.5);
+ hEx->GetYaxis()->SetNdivisions(520);
+ hEx->GetXaxis()->SetRangeUser(-1.,7.);
+ hEx->Draw();
+ FitKnownPeaks(hEx);
+
+ cFig_GGSP->Update();
+ double maxEx = cFig_GGSP->cd(2)->GetUymax();
+ TLine *Sn = new TLine(4.644, 0.0, 4.644, maxEx);
+ Sn->SetLineColor(kBlack); Sn->SetLineStyle(7);
+ Sn->Draw();
+ TLine *gs = new TLine(0.000, 0.0, 0.000, maxEx);
+ gs->SetLineColor(kBlack); gs->SetLineStyle(1);
+ gs->Draw();
+
+ cFig_GGSP->cd(1);
+ string drawg = "AddBack_EDC>>hEg(" + to_string(5./gammaBinWidth) + ",0,5)";
+ chain->Draw(drawg.c_str(),gammagate.c_str(),"");
+ TH1F* hEg = (TH1F*) gDirectory->Get("hEg");
+ hEg->SetLineColor(kBlack);
+ string nameEgYaxis = "Counts / " + to_string(gammaBinWidth) + " MeV";
+ hEg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ hEg->GetYaxis()->SetTitle(nameEgYaxis.c_str());
+ hEg->SetTitle("");
+ hEg->GetXaxis()->SetTitleSize(0.05);
+ hEg->GetXaxis()->SetLabelSize(0.05);
+ hEg->GetYaxis()->SetTitleSize(0.05);
+ hEg->GetYaxis()->SetLabelSize(0.05);
+ hEg->GetYaxis()->SetTitleOffset(0.5);
+ hEg->Draw();
+
+ cFig_GGSP->Update();
+}
+
diff --git a/Projects/e793s/macro/GausFit.h b/Projects/e793s/macro/GausFit.h
index 252ef56f84fe5adaede0245d942d45b6f2ef17be..700b1c6d471754f82fd49425dc2b73540e23226c 100755
--- a/Projects/e793s/macro/GausFit.h
+++ b/Projects/e793s/macro/GausFit.h
@@ -13,13 +13,12 @@ Double_t pi = 3.14159265358979323846;
////////////////////////////////////////////////////////
////////////////////////////////////////////////////////
+/*
void DoubleGaus(TH1F* hist){
bool repeat=true, bgbool = true;
int repeatInt;
double minFit, maxFit, mean1, mean2;
- double binWidth = hist->GetXaxis()->GetBinWidth(3);
-
while (repeat){
cout << "====================================================================" << endl;
cout << " Input range to fit:" << endl;
@@ -34,6 +33,19 @@ void DoubleGaus(TH1F* hist){
cout << " Background, yes or no?" << endl;
cin >> bgbool;
+ DoubleGausNumbs(hist, minFit, maxFit, mean1, mean2, bgbool);
+
+}
+*/
+
+vector<double> DoubleGausNumbs(TH1F* hist, double minFit, double maxFit, double mean1, double mean2, bool bgbool){
+ bool repeat=true;
+ int repeatInt;
+ vector<double> areasOut;
+
+// while (repeat){ /* Comment out for mass runs */
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+
//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);
@@ -61,22 +73,26 @@ void DoubleGaus(TH1F* hist){
"Background");
f1->SetLineColor(kBlack);
+ /* OPTION */ //bg->FixParameter(0,0.0); f1->FixParameter(6,0.0);
+
+ /* OPTION */ //g1->FixParameter(0,0.0); f1->FixParameter(0,0.0);
+ /* OPTION */ g1->FixParameter(2,0.6); f1->FixParameter(2,0.6);
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->SetParameter(2, 0.14);
- g1->SetParLimits(2, 0.05, 0.30);
- //g1->SetParLimits(2, 0.13, 0.30);
+ //g1->SetParameter(2, 0.13);
+ //g1->SetParameter(2, 0.6);//FOR ELASTICS
+ // //g1->SetParLimits(2, 0.05, 0.20);
+ // g1->SetParLimits(2, 0.4, 1.0);//FOR ELASTICS
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);
+ //g2->SetParameter(2, 0.13);
+ g2->SetParameter(2, 0.6);//FOR ELASTICS
+ //g2->SetParLimits(2, 0.05, 0.20);
+ g2->SetParLimits(2, 0.4, 1.0);//FOR ELASTICS
hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit);
@@ -127,25 +143,33 @@ void DoubleGaus(TH1F* hist){
cout << fixed << setprecision(5);
- cout << "\033[91m Mean: \t" << finalPar[1]
+ areasOut.push_back(finalPar[0]/binWidth);
+ areasOut.push_back((finalPar[0]/binWidth) * (finalErr[0]/finalPar[0]));
+ areasOut.push_back(finalPar[3]/binWidth);
+ areasOut.push_back((finalPar[3]/binWidth) * (finalErr[3]/finalPar[3]));
+
+ cout << RED;
+ cout << " Mean: \t" << finalPar[1]
<< "\t +- " << finalErr[1]
<< endl;
- cout << "\033[91m Sigm: \t" << finalPar[2]
+ cout << " Sigm: \t" << finalPar[2]
<< "\t +- " << finalErr[2]
<< endl;
- cout << "\033[91m Area: \t" << finalPar[0]/binWidth
+ cout << " Area: \t" << finalPar[0]/binWidth
<< "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
<< endl;
- cout << "\033[36m Mean: \t" << finalPar[4]
+ cout << BLUE;
+ cout << " Mean: \t" << finalPar[4]
<< "\t +- " << finalErr[4]
<< endl;
- cout << "\033[36m Sigm: \t" << finalPar[5]
+ cout << " Sigm: \t" << finalPar[5]
<< "\t +- " << finalErr[5]
<< endl;
- cout << "\033[36m Area: \t" << finalPar[3]/binWidth
+ cout << " Area: \t" << finalPar[3]/binWidth
<< "\t +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3])
<< endl;
+ cout << RESET;
TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
@@ -159,14 +183,48 @@ void DoubleGaus(TH1F* hist){
gPad->Modified();
gPad->Update();
+// cout << "\33[37m Refit? " << endl;
+// cin >> repeatInt;
+// if(repeatInt!=1){ repeat=false; }
+// } /* Comment out for mass runs */
+
+ return areasOut;
+
+}
+
+
+
+void DoubleGaus(TH1F* hist){
+ bool repeat=true;
+ int repeatInt;
+ bool bgbool = true;
+ double minFit, maxFit, mean1, mean2;
+
+ TCanvas* canvGausFit = new TCanvas("canvGausFit","canvGausFit",1000,1000);
+ hist->Draw();
+
+ 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 << " Background, yes or no?" << endl;
+ cin >> bgbool;
+
+ vector<double> areas = DoubleGausNumbs(hist, minFit, maxFit, mean1, mean2, bgbool);
+
cout << "\33[37m Refit? " << endl;
cin >> repeatInt;
if(repeatInt!=1){ repeat=false; }
}
-
}
-
////////////////////////////////////////////////////////
////////////////////////////////////////////////////////
@@ -339,19 +397,19 @@ void TripleGaus(TH1F* hist){
*/
g1->SetParameter(0, 100);
- g1->SetParLimits(0, 0.0, 500.0);
+ g1->SetParLimits(0, 0.0, 1000.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->SetParLimits(0, 0.0, 1000.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->SetParLimits(0, 0.0, 1000.0);
g3->SetParameter(1, mean3);
g3->SetParLimits(1, mean3-1.0, mean3+1.0);
g3->SetParameter(2, 0.13);
@@ -668,4 +726,286 @@ void SingleGaus(TH1F* hist, bool isGamma){
+void DoubleGausForElasticsFitting(TH1F* hist){
+ bool repeat=true, bgbool = true;
+ int repeatInt;
+ double minFit, maxFit, mean1, mean2;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+ cout << "Bin Width: " << binWidth << endl;
+
+
+ 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 << " 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))",
+ 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.30);
+ //g1->SetParameter(2, 0.14);
+ g1->SetParLimits(2, 0.10, 0.80);
+ //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.30);
+ //g2->SetParameter(2, 0.14);
+ g2->SetParLimits(2, 0.10, 0.80);
+ //g2->SetParLimits(2, 0.13, 0.30);
+ bg->SetParameter(0,25.);
+ bg->SetParLimits(0,1.,75.);
+
+ 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);
+
+ /* 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();
+
+ 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 SingleGausForElasticsFitting(TH1F* hist){
+ bool repeat=true, bgbool = true;
+ int repeatInt;
+ double minFit, maxFit, mean1;
+
+ double binWidth = hist->GetXaxis()->GetBinWidth(3);
+ cout << "Bin Width: " << binWidth << endl;
+
+
+ while (repeat){
+ cout << "====================================================================" << endl;
+ cout << " Input range to fit:" << endl;
+ cout << " Min = ";
+ cin >> minFit;
+ cout << " Max = ";
+ cin >> maxFit;
+ cout << " Peak 1 = ";
+ cin >> mean1;
+ 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))",
+ minFit, maxFit);
+ g1->SetLineColor(kRed);
+ g1->SetLineStyle(2);
+ //TF1 *g2 = new TF1 ("m1", equation, minFit, maxFit);
+ 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("double_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, 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);
+ bg->SetParameter(0,25.);
+ bg->SetParLimits(0,1.,50.);
+
+ hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit);
+
+ Double_t par[7];
+ g1->GetParameters(&par[0]);
+ bg->GetParameters(&par[3]);
+ 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(3,0.);}
+
+ 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");
+
+
+/* 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;
+
+ 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(bg,"Background","l");
+ legend->Draw();
+
+ //cGate->Draw("SAME");
+ gPad->Modified();
+ gPad->Update();
+
+ cout << "\33[37m Refit? " << endl;
+ cin >> repeatInt;
+ if(repeatInt!=1){ repeat=false; }
+ }
+
+}
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index 9cfc1427d4e097163cbb289e082dbf0435229554..b38eeb6ef85174db1f2ae27717120445f5a92ea2 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -3,24 +3,25 @@
#include <cmath>
#include "stdlib.h"
-const int numPeaks = 16;//15;
+const int numPeaks = 17;
array<double,numPeaks> means = { 0.000,
0.143,
0.279,
0.728,
0.968,
1.410,
- 1.981,
+ 1.981,//1.952,//1.981,
2.412,
2.910,
- 3.246,
+ 3.253,
3.605,
- 3.797,//Split in two?
- 3.876,//Split in two?
- 4.055,//4.1,
- 4.38,
- 4.51//,
- //5.24
+ 3.795,//Split in two?
+ 3.870,//Split in two?
+ 4.045,//4.1,
+ 4.393,
+ //4.51//,
+ 5.15,
+ 5.82
};
array<double,27> knowngammas = { 0.143,
@@ -76,16 +77,18 @@ Double_t f_full(Double_t *x, Double_t *par) {
float xx = x[0];
double result, norm;
// Flat background
- result = par[0];
+// result = par[0];
+ result = 0;
// Add N peaks
for(int pk=0; pk<numPeaks; 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[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){
+vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){
double minFit=-1.0, maxFit=8.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
double sigma = 0.14;
@@ -124,6 +127,17 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
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);
@@ -137,23 +151,19 @@ 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->SetParLimits((i*3)+3,10.0,1e5);
}
- //full->SetParameter(0,30.);
//full->SetParLimits(0,0.,40.); /* FOR TOTAL SPECTRUM FITTING */
- full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
- //full->SetParLimits(0,0.,1.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
- //full->FixParameter(9,0.); //??
-
- // Specific limits
- // Set max 279 counts to 100
- //full->FixParameter(9,0.0);
- //full->SetParLimits(9,0.0,1e2); // Doesnt work???
+ //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->FixParameter(0,slideshift); /* FOR WHEN PAR[0] IS FIXED ENERGY CENTRIOD SLIDER */
+ full->SetParameter(52,0.39); /* SET 5.8MeV SIGMA */
+ full->SetParLimits(52,0.34,0.44); /* SET 5.8MeV SIGMA */
+
allPeaks[14]->SetLineColor(kOrange);
//Fit full function to histogram
@@ -211,5 +221,5 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
void FitKnownPeaks(TH1F* hist){
//Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
- vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist);
+ vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0);
}
diff --git a/Projects/e793s/macro/Plots_47Kdp.C b/Projects/e793s/macro/Plots_47Kdp.C
index 7508c206b093f945b37809d4cd1fc8770900b605..e1d13083041c6b3ae464a3082b15c53c679a7e4f 100644
--- a/Projects/e793s/macro/Plots_47Kdp.C
+++ b/Projects/e793s/macro/Plots_47Kdp.C
@@ -1,3 +1,4 @@
+#include "DefineColours.h"
#include "GausFit.h"
#include "KnownPeakFitter.h"
#include "DrawPlots.h"
@@ -7,7 +8,7 @@
#include "ThreeBodyBreakup_FitPhaseSpace.h"
-void AddGammaLinesMG(TH1F* hist, double particle, double ymax){
+void AddGammaLines(TH1F* hist, double particle, double ymax){
string base = "sub ";
for(int i=1; i<means.size();i++){
@@ -22,7 +23,7 @@ void AddGammaLinesMG(TH1F* hist, double particle, double ymax){
}
}
-void AddPlacedGammasMG(TH1F* hist, double ymax){
+void AddPlacedGammas(TH1F* hist, double ymax){
hist->Draw();
for(int i=0; i<knowngammas.size();i++){
TLine *line = new TLine(knowngammas.at(i), 0.0, knowngammas.at(i), ymax);
@@ -41,24 +42,13 @@ void Plots_47Kdp(){
tCentre = 2700; tRange = 400;
timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange);
det_gate = "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8";
-
+ reactionName = "47K(d,p)";
+
cout << "==============================================" << endl;
cout << "=============== (d,p) reaction ===============" << endl;
cout << "==============================================" << endl;
cout << ""<< endl;
- cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
- cout << "---- 0.143, p3/2 = CS(0.143, 2, 1, 1.5) "<< endl;
- cout << "---- 0.279, p3/2 = CS(0.279, 2, 1, 1.5) "<< endl;
- cout << "---- 0.728, f7/2 = CS(0.728, 3, 3, 3.5) "<< endl;
- cout << "---- 0.968, p1/2 = CS(0.968, 0, 1, 0.5) "<< endl;
- cout << "---- 1.410, p3/2 = CS(1.410, 2, 1, 1.5) "<< endl;
- cout << "---- 1.981, p3/2 = CS(1.981, 2, 1, 1.5) "<< endl;
- cout << "---- 2.410, p3/2 = CS(2.410, 0, 1, 0.5) "<< endl;
- cout << "---- 3.2 , f7/2 = CS(3.2 , 3, 3, 3.5) "<< endl;
- cout << "---- 3.6 , f5/2 = CS(3.6 , 3, 3, 2.5) "<< endl;
- cout << "---- 3.8 , f5/2 = CS(3.8 , 3, 3, 2.5) "<< endl;
- cout << "---- 4.1 , f5/2 = CS(4.1 , 3, 3, 2.5) "<< endl;
- cout << "---- 4.4 , f5/2 = CS(4.4 , 3, 3, 2.5) "<< endl;
+ CS();
cout << ""<< endl;
cout << "==============================================" << endl;
diff --git a/Projects/e793s/macro/Plots_47Kdt.C b/Projects/e793s/macro/Plots_47Kdt.C
index 0fb0215d1160e95bc38af80399a2ad9b88ce130a..872ca63e4dcc6bbddcf1074b48438f8c718bd489 100644
--- a/Projects/e793s/macro/Plots_47Kdt.C
+++ b/Projects/e793s/macro/Plots_47Kdt.C
@@ -1,24 +1,89 @@
+#include "DefineColours.h"
#include "GausFit.h"
-//#include "KnownPeakFitter.h"
+#include "KnownPeakFitter.h"
#include "DrawPlots.h"
//#include "CS2.h"
//#include "ThreeBodyBreakup.h"
//#include "ThreeBodyBreakup_FitPhaseSpace.h"
+
+void AddGammaLines(TH1F* hist, double particle, double ymax){
+// string base = "sub ";
+//
+// for(int i=1; i<means.size();i++){
+// string name = base + to_string(means.at(i));
+// TLine *line = new TLine(particle-means.at(i), 0.0, particle-means.at(i), ymax);
+// line->SetLineColor(kBlack); line->SetLineStyle(kDotted);
+// line->Draw();
+// TText *text = new TText((1.-(means.at(i)/particle))*particle,0.8*ymax,name.c_str());
+// text->SetTextAngle(90);
+// //text->SetTextSize(40);
+// text->Draw();
+// }
+}
+
+void AddPlacedGammas(TH1F* hist, double ymax){
+// hist->Draw();
+// for(int i=0; i<knowngammas.size();i++){
+// TLine *line = new TLine(knowngammas.at(i), 0.0, knowngammas.at(i), ymax);
+// line->SetLineColor(kBlack); line->SetLineStyle(kDotted);
+// line->Draw();
+// }
+}
+
/* MAIN FUNCTION */
+void MM_Timing_Comparison(){
+ TCanvas* c = new TCanvas("cMMT","cMMT",1000,1000);
+
+
+ chain->Draw("MUST2.Si_T>>h1(200,600,700)",
+ "abs(T_MUGAST_VAMOS-2750)<350 && MUST2.TelescopeNumber==1 && MUST2.CsI_E>0",
+ "same");
+ TH1F* h1 = (TH1F*) gDirectory->Get("h1");
+ chain->Draw("MUST2.Si_T>>h2(200,600,700)",
+ "abs(T_MUGAST_VAMOS-2750)<350 && MUST2.TelescopeNumber==2 && MUST2.CsI_E>0",
+ "same");
+ TH1F* h2 = (TH1F*) gDirectory->Get("h2");
+ chain->Draw("MUST2.Si_T>>h3(200,600,700)",
+ "abs(T_MUGAST_VAMOS-2750)<350 && MUST2.TelescopeNumber==3 && MUST2.CsI_E>0",
+ "same");
+ TH1F* h3 = (TH1F*) gDirectory->Get("h3");
+ chain->Draw("MUST2.Si_T>>h4(200,600,700)",
+ "abs(T_MUGAST_VAMOS-2750)<350 && MUST2.TelescopeNumber==4 && MUST2.CsI_E>0",
+ "same");
+ TH1F* h4 = (TH1F*) gDirectory->Get("h4");
+
+ h1->SetLineColor(kRed);
+ h2->SetLineColor(kGreen);
+ h3->SetLineColor(kBlue);
+ h4->SetLineColor(kViolet);
+
+// h1->Draw();
+// h2->Draw("same");
+// h3->Draw("same");
+// h4->Draw("same");
+}
+
void Plots_47Kdt(){
+ /* Load graphical cut */
+ TFile gcIn("GraphicalCut_22Jun22.root");
+ TCutG* cutTritons = (TCutG*) gcIn.FindObjectAny("cutTritons");
+
LoadChain47Kdt();
gStyle->SetOptStat("nemMrRi");
tCentre = 2750; tRange = 350;
timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange);
det_gate = "MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<5";
+ reactionName = "47K(d,t)";
cout << "==============================================" << endl;
cout << "=============== (d,t) reaction ===============" << endl;
cout << "==============================================" << endl;
+ cout << " - (d,t) selection: 'cutTritons' " << endl;
+
}