Skip to content
Snippets Groups Projects
Commit 2325c1f4 authored by Cyril Lenain's avatar Cyril Lenain :surfer_tone3:
Browse files

Merge branch 'NPTool.2.dev' of gitlab.in2p3.fr:np/nptool into NPTool.2.dev

parents 47d2b7b5 0048c522
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1335 additions and 196 deletions
Projects/e793s/Analysis.cxx
+ 112
− 27
View file @ 2325c1f4
...@@ -32,6 +32,7 @@ using namespace std; ...@@ -32,6 +32,7 @@ using namespace std;
#include "NPAnalysisFactory.h" #include "NPAnalysisFactory.h"
#include "NPDetectorManager.h" #include "NPDetectorManager.h"
#include "NPOptionManager.h" #include "NPOptionManager.h"
#include "macro/DefineColours.h"
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){ Analysis::Analysis(){
...@@ -56,6 +57,7 @@ void Analysis::Init() { ...@@ -56,6 +57,7 @@ void Analysis::Init() {
} else { } else {
cout << " == == == == PHASE SPACE == == == ==" << endl; cout << " == == == == PHASE SPACE == == == ==" << endl;
isSim=false; isSim=false;
//isSim=true;
isPhaseSpace=true; isPhaseSpace=true;
} }
...@@ -76,8 +78,8 @@ void Analysis::Init() { ...@@ -76,8 +78,8 @@ void Analysis::Init() {
if (i==6){int j=7;} if (i==6){int j=7;}
string name1 = base1 + to_string(j); string name1 = base1 + to_string(j);
string name2 = base2 + to_string(j); string name2 = base2 + to_string(j);
ThetaCM_detected_MGX[i] = new TH1F(name1.c_str(),name1.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(),180,0,180); ThetaLab_detected_MGX[i] = new TH1F(name2.c_str(),name2.c_str(),NumThetaAngleBins,0,180);
} }
// Initilize MMX histograms // Initilize MMX histograms
...@@ -87,8 +89,8 @@ void Analysis::Init() { ...@@ -87,8 +89,8 @@ void Analysis::Init() {
int j=i+1; int j=i+1;
string name1 = base3 + to_string(j); string name1 = base3 + to_string(j);
string name2 = base4 + to_string(j); string name2 = base4 + to_string(j);
ThetaCM_detected_MMX[i] = new TH1F(name1.c_str(),name1.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(),180,0,180); ThetaLab_detected_MMX[i] = new TH1F(name2.c_str(),name2.c_str(),NumThetaAngleBins,0,180);
} }
} }
...@@ -134,7 +136,13 @@ void Analysis::Init() { ...@@ -134,7 +136,13 @@ void Analysis::Init() {
FinalBeamEnergy = BeamTargetELoss.Slow(OriginalBeamEnergy, 0.5*TargetThickness, 0); FinalBeamEnergy = BeamTargetELoss.Slow(OriginalBeamEnergy, 0.5*TargetThickness, 0);
reaction.SetBeamEnergy(FinalBeamEnergy); 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 // initialize various parameters
Rand = TRandom3(); Rand = TRandom3();
...@@ -266,9 +274,19 @@ void Analysis::TreatEvent(){ ...@@ -266,9 +274,19 @@ void Analysis::TreatEvent(){
// if(!isSim){ // if(!isSim){
// Evaluate energy using the thickness // 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 // 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;} // } else {elab_tmp = Energy;}
ELab.push_back(elab_tmp); ELab.push_back(elab_tmp);
...@@ -346,10 +364,14 @@ void Analysis::TreatEvent(){ ...@@ -346,10 +364,14 @@ void Analysis::TreatEvent(){
Energy, //particle energy after Al Energy, //particle energy after Al
0.4*micrometer, //thickness of Al 0.4*micrometer, //thickness of Al
ThetaMGSurface); //angle of impingement ThetaMGSurface); //angle of impingement
ELoss_Al.push_back(Energy-elab_tmp);
double elab_tmp2 = elab_tmp;
elab_tmp = LightTarget.EvaluateInitialEnergy( elab_tmp = LightTarget.EvaluateInitialEnergy(
elab_tmp, //particle energy after leaving target elab_tmp, //particle energy after leaving target
TargetThickness*0.5, //distance passed through target TargetThickness*0.5, //distance passed through target
ThetaNormalTarget); //angle of exit from 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 } else { //TESTING DIFFERENT ENERGY LOSSES IN SIMULATION
elab_tmp = Energy; //so I can add and remove sections elab_tmp = Energy; //so I can add and remove sections
//elab_tmp = LightSi.EvaluateInitialEnergy( //elab_tmp = LightSi.EvaluateInitialEnergy(
...@@ -368,6 +390,14 @@ void Analysis::TreatEvent(){ ...@@ -368,6 +390,14 @@ void Analysis::TreatEvent(){
ELab.push_back(elab_tmp); 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 // Part 3 : Excitation Energy Calculation
//if(!isSim){ //TESTING!!!! //if(!isSim){ //TESTING!!!!
Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp)); Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp));
...@@ -512,17 +542,45 @@ void Analysis::TreatEvent(){ ...@@ -512,17 +542,45 @@ void Analysis::TreatEvent(){
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void FillSolidAngles(TH1F* hSA, TH1F* hDet, TH1F* hEmm){ 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){ for (int i = 0; i < hSA->GetNbinsX(); ++i){
double val = hDet->GetBinContent(i) / hEmm->GetBinContent(i); double val = hDet->GetBinContent(i) / hEmm->GetBinContent(i);
double valerr = val * sqrt( double valerr = val * sqrt(
pow(hDet->GetBinError(i) / hDet->GetBinContent(i), 2) + pow(hDet->GetBinError(i) / hDet->GetBinContent(i), 2) +
pow(hEmm->GetBinError(i) / hEmm->GetBinContent(i), 2) ); pow(hEmm->GetBinError(i) / hEmm->GetBinContent(i), 2) );
if (isnan(val)) { val = 0; valerr = 0; } if (isnan(val)) { val = 0; valerr = 0; }
val *= clineVal.at(i);
valerr *= clineVal.at(i);
hSA->SetBinContent(i, val); hSA->SetBinContent(i, val);
hSA->SetBinError(i, valerr); 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(){ void Analysis::End(){
cout << endl << "\e[1;32m GATCONF Statistics " << endl ; cout << endl << "\e[1;32m GATCONF Statistics " << endl ;
...@@ -532,7 +590,7 @@ void Analysis::End(){ ...@@ -532,7 +590,7 @@ void Analysis::End(){
cout << endl ; cout << endl ;
if(isSim && !isPhaseSpace){ if(isSim && !isPhaseSpace){
//TObjArray HistList(0); //TObjArray HistList(0);
TList *HistList = new TList(); TList *HistList = new TList();
...@@ -556,14 +614,19 @@ void Analysis::End(){ ...@@ -556,14 +614,19 @@ void Analysis::End(){
Efficiency_CM_MM->Divide(ThetaCM_emmitted); Efficiency_CM_MM->Divide(ThetaCM_emmitted);
Efficiency_Lab_MM->Divide(ThetaLab_emmitted); Efficiency_Lab_MM->Divide(ThetaLab_emmitted);
double dt_MM = 180./Efficiency_Lab_MM->GetNbinsX(); //double dt_MM = 180./Efficiency_Lab_MM->GetNbinsX();
cout << "Angular infinitesimal (MM) = " << dt_MM << "deg " << endl; //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); //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); FillSolidAngles(SolidAngle_CM_MM, ThetaCM_detected_MM, ThetaCM_emmitted);
SolidAngle_CM_MM->Divide(Cline_MM,1); FillSolidAngles(SolidAngle_Lab_MM, ThetaLab_detected_MM, ThetaLab_emmitted);
SolidAngle_Lab_MM->Divide(ThetaLab_emmitted);
SolidAngle_Lab_MM->Divide(Cline_MM,1); //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(ThetaCM_emmitted);
HistList->Add(ThetaLab_emmitted); HistList->Add(ThetaLab_emmitted);
...@@ -573,7 +636,8 @@ void Analysis::End(){ ...@@ -573,7 +636,8 @@ void Analysis::End(){
HistList->Add(Efficiency_Lab_MM); HistList->Add(Efficiency_Lab_MM);
HistList->Add(SolidAngle_CM_MM); HistList->Add(SolidAngle_CM_MM);
HistList->Add(SolidAngle_Lab_MM); HistList->Add(SolidAngle_Lab_MM);
HistList->Add(Cline_MM); //HistList->Add(Cline_MM);
//HistList->Add(Cline);
// MUGAST // MUGAST
...@@ -596,16 +660,17 @@ void Analysis::End(){ ...@@ -596,16 +660,17 @@ void Analysis::End(){
Efficiency_CM_MG->Divide(ThetaCM_emmitted); Efficiency_CM_MG->Divide(ThetaCM_emmitted);
Efficiency_Lab_MG->Divide(ThetaLab_emmitted); Efficiency_Lab_MG->Divide(ThetaLab_emmitted);
double dt_MG = 180./Efficiency_Lab_MG->GetNbinsX(); //double dt_MG = 180./Efficiency_Lab_MG->GetNbinsX();
cout << "Angular infinitesimal (MG) = " << dt_MG << "deg " << endl; //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); //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 */ /* Testing method for better errors in SolidAngle histograms */
FillSolidAngles(SolidAngle_CM_MG, ThetaCM_detected_MG, ThetaCM_emmitted); 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); 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(ThetaCM_detected_MG);
HistList->Add(ThetaLab_detected_MG); HistList->Add(ThetaLab_detected_MG);
...@@ -613,7 +678,7 @@ void Analysis::End(){ ...@@ -613,7 +678,7 @@ void Analysis::End(){
HistList->Add(Efficiency_Lab_MG); HistList->Add(Efficiency_Lab_MG);
HistList->Add(SolidAngle_CM_MG); HistList->Add(SolidAngle_CM_MG);
HistList->Add(SolidAngle_Lab_MG); HistList->Add(SolidAngle_Lab_MG);
HistList->Add(Cline_MG); //HistList->Add(Cline_MG);
// MUGAST INDIVIDUAL // MUGAST INDIVIDUAL
TH1F *SolidAngle_CM_MGX[6]; TH1F *SolidAngle_CM_MGX[6];
...@@ -626,7 +691,8 @@ void Analysis::End(){ ...@@ -626,7 +691,8 @@ void Analysis::End(){
SolidAngle_CM_MGX[i]->SetName(name.c_str()); SolidAngle_CM_MGX[i]->SetName(name.c_str());
SolidAngle_CM_MGX[i]->SetTitle(name.c_str()); SolidAngle_CM_MGX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_CM_MGX[i], ThetaCM_detected_MGX[i], ThetaCM_emmitted); 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]; TH1F *SolidAngle_Lab_MGX[6];
...@@ -639,7 +705,8 @@ void Analysis::End(){ ...@@ -639,7 +705,8 @@ void Analysis::End(){
SolidAngle_Lab_MGX[i]->SetName(name.c_str()); SolidAngle_Lab_MGX[i]->SetName(name.c_str());
SolidAngle_Lab_MGX[i]->SetTitle(name.c_str()); SolidAngle_Lab_MGX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_Lab_MGX[i], ThetaLab_detected_MGX[i], ThetaLab_emmitted); 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 // MUST2 INDIVIDUAL
...@@ -652,7 +719,8 @@ void Analysis::End(){ ...@@ -652,7 +719,8 @@ void Analysis::End(){
SolidAngle_CM_MMX[i]->SetName(name.c_str()); SolidAngle_CM_MMX[i]->SetName(name.c_str());
SolidAngle_CM_MMX[i]->SetTitle(name.c_str()); SolidAngle_CM_MMX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_CM_MMX[i], ThetaCM_detected_MMX[i], ThetaCM_emmitted); 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]; TH1F *SolidAngle_Lab_MMX[6];
...@@ -664,7 +732,8 @@ void Analysis::End(){ ...@@ -664,7 +732,8 @@ void Analysis::End(){
SolidAngle_Lab_MMX[i]->SetName(name.c_str()); SolidAngle_Lab_MMX[i]->SetName(name.c_str());
SolidAngle_Lab_MMX[i]->SetTitle(name.c_str()); SolidAngle_Lab_MMX[i]->SetTitle(name.c_str());
FillSolidAngles(SolidAngle_Lab_MMX[i], ThetaLab_detected_MMX[i], ThetaLab_emmitted); 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]);} for(int i=0; i<6; i++){HistList->Add(ThetaCM_detected_MGX[i]);}
...@@ -677,6 +746,14 @@ void Analysis::End(){ ...@@ -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_CM_MMX[i]);}
for(int i=0; i<5; i++){HistList->Add(SolidAngle_Lab_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"); auto HistoFile = new TFile("~/Programs/nptool/Projects/e793s/SolidAngle_HistFile_New.root","RECREATE");
HistList->Write(); HistList->Write();
HistoFile->Close(); HistoFile->Close();
...@@ -706,6 +783,9 @@ void Analysis::InitOutputBranch(){ ...@@ -706,6 +783,9 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab); RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
RootOutput::getInstance()->GetTree()->Branch("Ecm",&Ecm); RootOutput::getInstance()->GetTree()->Branch("Ecm",&Ecm);
RootOutput::getInstance()->GetTree()->Branch("RawEnergy",&RawEnergy); 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("ThetaLab",&ThetaLab);
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab); RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab);
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM); RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM);
...@@ -967,6 +1047,9 @@ void Analysis::ReInitValue(){ ...@@ -967,6 +1047,9 @@ void Analysis::ReInitValue(){
EAgata = -1000; EAgata = -1000;
ELab.clear(); ELab.clear();
RawEnergy.clear(); RawEnergy.clear();
ELoss_Al.clear();
ELoss_Target.clear();
ELoss.clear();
ThetaLab.clear(); ThetaLab.clear();
PhiLab.clear(); PhiLab.clear();
ThetaCM.clear(); ThetaCM.clear();
...@@ -991,6 +1074,8 @@ void Analysis::ReInitValue(){ ...@@ -991,6 +1074,8 @@ void Analysis::ReInitValue(){
thetalab_tmp = 0; thetalab_tmp = 0;
philab_tmp = 0; philab_tmp = 0;
filledCline=false;
MG_T=-1000; MG_T=-1000;
MG_E=-1000; MG_E=-1000;
MG_X=-1000; MG_X=-1000;
......
Projects/e793s/Analysis.h
+ 16
− 6
View file @ 2325c1f4
...@@ -42,13 +42,20 @@ ...@@ -42,13 +42,20 @@
#include <TMath.h> #include <TMath.h>
#include <bitset> #include <bitset>
auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",180,0,180); int NumThetaAngleBins = 1800;// 180 = 1 deg, 360 = 0.5 deg, 900 = 0.2 deg, 1800 = 0.1 deg
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);
auto ThetaLab_emmitted = new TH1F("ThetaLab_emmitted","ThetaLab_emmitted",180,0,180); auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",NumThetaAngleBins,0,180);
auto ThetaLab_detected_MM = new TH1F("ThetaLab_detected_MM","ThetaLab_detected_MM",180,0,180); auto ThetaCM_detected_MM = new TH1F("ThetaCM_detected_MM","ThetaCM_detected_MM",NumThetaAngleBins,0,180);
auto ThetaLab_detected_MG = new TH1F("ThetaLab_detected_MG","ThetaLab_detected_MG",180,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_MGX[6];
TH1F *ThetaCM_detected_MMX[5]; TH1F *ThetaCM_detected_MMX[5];
...@@ -94,6 +101,9 @@ class Analysis: public NPL::VAnalysis{ ...@@ -94,6 +101,9 @@ class Analysis: public NPL::VAnalysis{
std::vector<double> Ex; std::vector<double> Ex;
std::vector<double> Ecm; std::vector<double> Ecm;
std::vector<double> RawEnergy; 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> ThetaLab;
std::vector<double> PhiLab; std::vector<double> PhiLab;
std::vector<double> ThetaCM; std::vector<double> ThetaCM;
......
Projects/e793s/exp.sh
+ 23
− 2
View file @ 2325c1f4
...@@ -2,5 +2,26 @@ ...@@ -2,5 +2,26 @@
cd ~/Programs/nptool/Projects/e793s; cd ~/Programs/nptool/Projects/e793s;
cmake ./; cmake ./;
make -j6; 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;
Projects/e793s/macro/CS2.h
+ 155
− 42
View file @ 2325c1f4
...@@ -16,7 +16,7 @@ int indexE; ...@@ -16,7 +16,7 @@ int indexE;
double globalS, globalSerr; double globalS, globalSerr;
/* Output volume toggle */ /* Output volume toggle */
bool loud = 0; bool loud = 1;
/* Scale method toggle */ /* Scale method toggle */
bool scaleTogether = 1; bool scaleTogether = 1;
...@@ -25,22 +25,95 @@ bool scaleTogether = 1; ...@@ -25,22 +25,95 @@ bool scaleTogether = 1;
string orbitalname; string orbitalname;
string orbital; 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(){ void CS(){
/* Overload function */ /* 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 << " 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 << " 0f5/2 | ----- | --?-- | l=3, j=2.5, n=0" << endl;
cout << " 1p1/2 | ----- | --?-- | l=1, j=0.5, n=1" << endl; cout << " | | |" << endl;
cout << " 1p3/2 | ----- | --?-- | l=1, j=1.5, n=1" << endl; cout << " 1p1/2 | ----- | --?-- | l=1, j=0.5, n=1" << endl;
cout << " 0f7/2 | ----- | ===== |" << endl; cout << " 1p3/2 | ----- | --?-- | l=1, j=1.5, n=1" << endl;
cout << " | | |" << endl; cout << " 0f7/2 | ----- | ===== |" << endl;
cout << " 0d3/2 | ----- | ===== |" << endl; cout << " | | |" << endl;
cout << " 1s1/2 | --x-- | ===== |" << endl; cout << " 0d3/2 | ----- | ===== |" << endl;
cout << " | | |" << endl; cout << " 1s1/2 | --x-- | ===== |" << endl;
cout << " | p | n |" << 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){ ...@@ -49,8 +122,11 @@ void CS(double Energy, double Spin, double spdf, double angmom){
// p3/2 -> spdf = 1, angmom = 1.5 // p3/2 -> spdf = 1, angmom = 1.5
// J0 is incident spin, which is 47K g.s. therefore J0 = 1/2 // J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
double J0 = 0.5; 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(); orbitalname.clear();
orbital.clear(); orbital.clear();
if(spdf==1){ if(spdf==1){
...@@ -61,7 +137,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){ ...@@ -61,7 +137,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
orbitalname="p_{3/2}"; orbitalname="p_{3/2}";
orbital="p32"; orbital="p32";
} else { } else {
orbitalname="?????"; orbitalname;
orbital="???"; orbital="???";
} }
} else if (spdf==3){ } else if (spdf==3){
...@@ -82,7 +158,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){ ...@@ -82,7 +158,9 @@ void CS(double Energy, double Spin, double spdf, double angmom){
/* Reduce by factor of 10,000 */ /* Reduce by factor of 10,000 */
ElasticNorm /= 10000.; ElasticNorm /= 10000.;
//ElasticNorm /= 1000.;
ElasticNormErr /= 10000.; ElasticNormErr /= 10000.;
//ElasticNormErr /= 1000.;
double nodes; double nodes;
if(spdf==1){ if(spdf==1){
...@@ -122,12 +200,18 @@ void CS(double Energy, double Spin, double spdf, double angmom){ ...@@ -122,12 +200,18 @@ void CS(double Energy, double Spin, double spdf, double angmom){
} }
/* Solid Angle (from simulation) */ /* Solid Angle (from simulation) */
//auto file = new TFile("../SolidAngle_HistFile_15Feb_47Kdp.root"); /* ADD OPTION TO CHANGE SOLID ANGLE FILE DEPENDING ON PEAK!!!!*/
auto file = new TFile("../SolidAngle_HistFile_17May22_47Kdp_0143.root"); //auto file = new TFile("../SolidAngle_HistFile_47Kdp_26May22_v2_0143.root");
//auto file = new TFile("../SolidAngle_HistFile_06May_WideStripMatching_LargeRun.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"); TH1F* SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_Lab_MG");
TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000); 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) */ /* (canvas deleted after Area/SA calculation) */
/* Area of experimental peaks */ /* Area of experimental peaks */
...@@ -155,6 +239,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){ ...@@ -155,6 +239,8 @@ void CS(double Energy, double Spin, double spdf, double angmom){
int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001); int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001);
int binmax = SolidAngle->FindBin(areaArray[i][4]-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]); anglecentres.push_back(((areaArray[i][4]-areaArray[i][3])/2.)+areaArray[i][3]);
anglewidth.push_back(areaArray[i][4]-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){ ...@@ -202,7 +288,7 @@ void CS(double Energy, double Spin, double spdf, double angmom){
<< endl; << endl;
} }
} }
delete c_SolidAngle; //delete c_SolidAngle;
/* Graph of Area/Solid Angle*/ /* Graph of Area/Solid Angle*/
TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000); 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){ ...@@ -349,6 +435,13 @@ void CS(double Energy, double Spin, double spdf, double angmom){
c_Chi2Min->SaveAs(savestring1.c_str()); c_Chi2Min->SaveAs(savestring1.c_str());
c_Chi2Min->SaveAs(savestring2.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_AoSA;
//delete c_dSdO; //delete c_dSdO;
} }
...@@ -358,19 +451,23 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -358,19 +451,23 @@ vector<vector<double>> GetExpDiffCross(double Energy){
cout << "========================================================" << endl; cout << "========================================================" << endl;
vector<vector<double>> AllPeaks_OneGate; vector<vector<double>> AllPeaks_OneGate;
vector<vector<double>> OnePeak_AllGates; vector<vector<double>> OnePeak_AllGates;
int numbins = 10; /****CHANGE ANGLE GATING****/
double x[numbins], y[numbins]; int numAngleBins = 20;
double widthAngleBins = 2.5;
double firstAngle = 105.;
/***************************/
double x[numAngleBins], y[numAngleBins];
//TList* list = new TList(); //TList* list = new TList();
/* Determine scaling factor for PhaseSpace */ /* Determine scaling factor for PhaseSpace */
TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000); TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
double trackScale = 0.0; double trackScale = 0.0;
if(scaleTogether){ if(scaleTogether){
TH1F* baseEx = PullThetaLabHist(0,105.,5.); TH1F* baseEx = PullThetaLabHist(0,firstAngle,widthAngleBins);
TH1F* basePS = PullPhaseSpaceHist(0,105.,5.); TH1F* basePS = PullPhaseSpaceHist(0,firstAngle,widthAngleBins);
for(int i=1; i<numbins;i++){ for(int i=1; i<numAngleBins;i++){
TH1F* addEx = PullThetaLabHist(i,105.,5.); baseEx->Add(addEx,1.); TH1F* addEx = PullThetaLabHist(i,firstAngle,widthAngleBins); baseEx->Add(addEx,1.);
TH1F* addPS = PullPhaseSpaceHist(i,105.,5.); basePS->Add(addPS,1.); TH1F* addPS = PullPhaseSpaceHist(i,firstAngle,widthAngleBins); basePS->Add(addPS,1.);
} }
/* Subtract flat background equal to smallest bin in range */ /* Subtract flat background equal to smallest bin in range */
...@@ -385,7 +482,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -385,7 +482,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Begin scaling within range, track changes */ /* Begin scaling within range, track changes */
basePS->Scale(0.1); basePS->Scale(0.1);
trackScale = 0.1; trackScale = 0.1;
int numbinsScale = baseEx->GetNbinsX(); int numAngleBinsScale = baseEx->GetNbinsX();
int nbinlow = basePS->FindBin(4.); int nbinhigh = basePS->FindBin(8.0); int nbinlow = basePS->FindBin(4.); int nbinhigh = basePS->FindBin(8.0);
for(int b=nbinlow; b<nbinhigh; b++){ for(int b=nbinlow; b<nbinhigh; b++){
if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){ if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
...@@ -405,15 +502,14 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -405,15 +502,14 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
// TEMPORARY!!! REMOVE LAST THREE BINS ON HIGH ENERGY STATES!!! // 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++){ for(int i=0; i<numAngleBins;i++){
double bin = 5.; double min = firstAngle + (i*widthAngleBins);
double min = 105. + (i*bin); double max = min + widthAngleBins;
double max = min + bin;
cout << "===================================" << endl; cout << "===================================" << endl;
cout << "min: " << min << " max: " << max << endl; cout << "min: " << min << " max: " << max << endl;
...@@ -424,8 +520,8 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -424,8 +520,8 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */ /* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
/* To change angle gates, run GateThetaLab_MultiWrite() */ /* To change angle gates, run GateThetaLab_MultiWrite() */
TH1F* gate = PullThetaLabHist(i,105.,5.); TH1F* gate = PullThetaLabHist(i,firstAngle,widthAngleBins);
TH1F* pspace = PullPhaseSpaceHist(i,105.,5.); TH1F* pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins);
/* Scale the Phase Space at this angle... */ /* Scale the Phase Space at this angle... */
/* ... for all angles together */ /* ... for all angles together */
...@@ -438,8 +534,8 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -438,8 +534,8 @@ vector<vector<double>> GetExpDiffCross(double Energy){
if(pspace->Integral() > 50.){ // Non-garbage histogram if(pspace->Integral() > 50.){ // Non-garbage histogram
pspace->Scale(0.01); pspace->Scale(0.01);
trackScale=0.01; trackScale=0.01;
int numbins = gate->GetNbinsX(); int numAngleBins = gate->GetNbinsX();
for(int b=0; b<numbins; b++){ for(int b=0; b<numAngleBins; b++){
if(loud){cout << " FROM " << pspace->GetBinContent(b) << if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
" > " << gate->GetBinContent(b); " > " << gate->GetBinContent(b);
} }
...@@ -471,14 +567,22 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -471,14 +567,22 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Retrieve array containing all fits, for one angle gate. * /* Retrieve array containing all fits, for one angle gate. *
* Specific peak of interest selected from the vector by * * Specific peak of interest selected from the vector by *
* global variable indexE */ * 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 */ /* Write PS-subtracted spectrum to list */
//list->Add(gate); //list->Add(gate);
//list->Add(c_peakFits); //list->Add(c_peakFits);
gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
string filename = "./CS2_Figures/"; string filename = "./CS2_Figures/";
filename += tmp2 + ".root"; filename += tmp2 + ".root";
c_peakFits->SaveAs(filename.c_str()); 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 */ /* Check correct OneGate vector is selected */
cout << "area of " << means[indexE] << " = " cout << "area of " << means[indexE] << " = "
...@@ -505,7 +609,14 @@ vector<vector<double>> GetExpDiffCross(double Energy){ ...@@ -505,7 +609,14 @@ vector<vector<double>> GetExpDiffCross(double Energy){
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
//TFile* file = new TFile("GateThetaLabHistograms.root","READ"); //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_" string histname = "cThetaLabGate_"
+ to_string((int) (minTheta+(i*gatesize))) + "-" + to_string((int) (minTheta+(i*gatesize))) + "-"
+ to_string((int) (minTheta+((i+1)*gatesize))); + to_string((int) (minTheta+((i+1)*gatesize)));
...@@ -518,7 +629,9 @@ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){ ...@@ -518,7 +629,9 @@ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
TH1F* PullPhaseSpaceHist(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_" string histname = "cPSpaceThetaLabGate_"
+ to_string((int) (minTheta+(i*gatesize))) + "-" + to_string((int) (minTheta+(i*gatesize))) + "-"
+ to_string((int) (minTheta+((i+1)*gatesize))); + to_string((int) (minTheta+((i+1)*gatesize)));
......
Projects/e793s/macro/DrawPlots.h
+ 560
− 55
View file @ 2325c1f4
This diff is collapsed.
Projects/e793s/macro/GausFit.h
+ 361
− 21
View file @ 2325c1f4
...@@ -13,13 +13,12 @@ Double_t pi = 3.14159265358979323846; ...@@ -13,13 +13,12 @@ Double_t pi = 3.14159265358979323846;
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
/*
void DoubleGaus(TH1F* hist){ void DoubleGaus(TH1F* hist){
bool repeat=true, bgbool = true; bool repeat=true, bgbool = true;
int repeatInt; int repeatInt;
double minFit, maxFit, mean1, mean2; double minFit, maxFit, mean1, mean2;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
while (repeat){ while (repeat){
cout << "====================================================================" << endl; cout << "====================================================================" << endl;
cout << " Input range to fit:" << endl; cout << " Input range to fit:" << endl;
...@@ -34,6 +33,19 @@ void DoubleGaus(TH1F* hist){ ...@@ -34,6 +33,19 @@ void DoubleGaus(TH1F* hist){
cout << " Background, yes or no?" << endl; cout << " Background, yes or no?" << endl;
cin >> bgbool; 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", equation, minFit, maxFit);
TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))", TF1 *g1 = new TF1 ("m1", "([0]/([2]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[2],2))",
minFit, maxFit); minFit, maxFit);
...@@ -61,22 +73,26 @@ void DoubleGaus(TH1F* hist){ ...@@ -61,22 +73,26 @@ void DoubleGaus(TH1F* hist){
"Background"); "Background");
f1->SetLineColor(kBlack); 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->SetParameter(0, 100);
g1->SetParLimits(0, 0.0, 500.0); g1->SetParLimits(0, 0.0, 500.0);
g1->SetParameter(1, mean1); g1->SetParameter(1, mean1);
g1->SetParLimits(1, mean1-0.5, mean1+0.5); g1->SetParLimits(1, mean1-0.5, mean1+0.5);
g1->SetParameter(2, 0.13); //g1->SetParameter(2, 0.13);
//g1->SetParameter(2, 0.14); //g1->SetParameter(2, 0.6);//FOR ELASTICS
g1->SetParLimits(2, 0.05, 0.30); // //g1->SetParLimits(2, 0.05, 0.20);
//g1->SetParLimits(2, 0.13, 0.30); // g1->SetParLimits(2, 0.4, 1.0);//FOR ELASTICS
g2->SetParameter(0, 100); g2->SetParameter(0, 100);
g2->SetParLimits(0, 0.0, 500.0); g2->SetParLimits(0, 0.0, 500.0);
g2->SetParameter(1, mean2); g2->SetParameter(1, mean2);
g2->SetParLimits(1, mean2-1.0, mean2+1.0); g2->SetParLimits(1, mean2-1.0, mean2+1.0);
g2->SetParameter(2, 0.13); //g2->SetParameter(2, 0.13);
//g2->SetParameter(2, 0.14); g2->SetParameter(2, 0.6);//FOR ELASTICS
g2->SetParLimits(2, 0.05, 0.30); //g2->SetParLimits(2, 0.05, 0.20);
//g2->SetParLimits(2, 0.13, 0.30); g2->SetParLimits(2, 0.4, 1.0);//FOR ELASTICS
hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit); hist->Fit(g1, "WWR", "", minFit, mean1+5);//maxFit);
...@@ -127,25 +143,33 @@ void DoubleGaus(TH1F* hist){ ...@@ -127,25 +143,33 @@ void DoubleGaus(TH1F* hist){
cout << fixed << setprecision(5); 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] << "\t +- " << finalErr[1]
<< endl; << endl;
cout << "\033[91m Sigm: \t" << finalPar[2] cout << " Sigm: \t" << finalPar[2]
<< "\t +- " << finalErr[2] << "\t +- " << finalErr[2]
<< endl; << endl;
cout << "\033[91m Area: \t" << finalPar[0]/binWidth cout << " Area: \t" << finalPar[0]/binWidth
<< "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0]) << "\t +- " << (finalPar[0]/binWidth) * (finalErr[0]/finalPar[0])
<< endl; << endl;
cout << "\033[36m Mean: \t" << finalPar[4] cout << BLUE;
cout << " Mean: \t" << finalPar[4]
<< "\t +- " << finalErr[4] << "\t +- " << finalErr[4]
<< endl; << endl;
cout << "\033[36m Sigm: \t" << finalPar[5] cout << " Sigm: \t" << finalPar[5]
<< "\t +- " << finalErr[5] << "\t +- " << finalErr[5]
<< endl; << endl;
cout << "\033[36m Area: \t" << finalPar[3]/binWidth cout << " Area: \t" << finalPar[3]/binWidth
<< "\t +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3]) << "\t +- " << (finalPar[3]/binWidth) * (finalErr[3]/finalPar[3])
<< endl; << endl;
cout << RESET;
TLegend *legend=new TLegend(0.7,0.7,0.9,0.9); TLegend *legend=new TLegend(0.7,0.7,0.9,0.9);
...@@ -159,14 +183,48 @@ void DoubleGaus(TH1F* hist){ ...@@ -159,14 +183,48 @@ void DoubleGaus(TH1F* hist){
gPad->Modified(); gPad->Modified();
gPad->Update(); 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; cout << "\33[37m Refit? " << endl;
cin >> repeatInt; cin >> repeatInt;
if(repeatInt!=1){ repeat=false; } if(repeatInt!=1){ repeat=false; }
} }
} }
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
...@@ -339,19 +397,19 @@ void TripleGaus(TH1F* hist){ ...@@ -339,19 +397,19 @@ void TripleGaus(TH1F* hist){
*/ */
g1->SetParameter(0, 100); g1->SetParameter(0, 100);
g1->SetParLimits(0, 0.0, 500.0); g1->SetParLimits(0, 0.0, 1000.0);
g1->SetParameter(1, mean1); g1->SetParameter(1, mean1);
g1->SetParLimits(1, mean1-1.0, mean1+1.0); g1->SetParLimits(1, mean1-1.0, mean1+1.0);
g1->SetParameter(2, 0.13); g1->SetParameter(2, 0.13);
g1->SetParLimits(2, 0.05, 0.30); g1->SetParLimits(2, 0.05, 0.30);
g2->SetParameter(0, 100); g2->SetParameter(0, 100);
g2->SetParLimits(0, 0.0, 500.0); g2->SetParLimits(0, 0.0, 1000.0);
g2->SetParameter(1, mean2); g2->SetParameter(1, mean2);
g2->SetParLimits(1, mean2-1.0, mean2+1.0); g2->SetParLimits(1, mean2-1.0, mean2+1.0);
g2->SetParameter(2, 0.13); g2->SetParameter(2, 0.13);
g2->SetParLimits(2, 0.05, 0.30); g2->SetParLimits(2, 0.05, 0.30);
g3->SetParameter(0, 100); g3->SetParameter(0, 100);
g3->SetParLimits(0, 0.0, 500.0); g3->SetParLimits(0, 0.0, 1000.0);
g3->SetParameter(1, mean3); g3->SetParameter(1, mean3);
g3->SetParLimits(1, mean3-1.0, mean3+1.0); g3->SetParLimits(1, mean3-1.0, mean3+1.0);
g3->SetParameter(2, 0.13); g3->SetParameter(2, 0.13);
...@@ -668,4 +726,286 @@ void SingleGaus(TH1F* hist, bool isGamma){ ...@@ -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; }
}
}
Projects/e793s/macro/KnownPeakFitter.h
+ 36
− 26
View file @ 2325c1f4
...@@ -3,24 +3,25 @@ ...@@ -3,24 +3,25 @@
#include <cmath> #include <cmath>
#include "stdlib.h" #include "stdlib.h"
const int numPeaks = 16;//15; const int numPeaks = 17;
array<double,numPeaks> means = { 0.000, array<double,numPeaks> means = { 0.000,
0.143, 0.143,
0.279, 0.279,
0.728, 0.728,
0.968, 0.968,
1.410, 1.410,
1.981, 1.981,//1.952,//1.981,
2.412, 2.412,
2.910, 2.910,
3.246, 3.253,
3.605, 3.605,
3.797,//Split in two? 3.795,//Split in two?
3.876,//Split in two? 3.870,//Split in two?
4.055,//4.1, 4.045,//4.1,
4.38, 4.393,
4.51//, //4.51//,
//5.24 5.15,
5.82
}; };
array<double,27> knowngammas = { 0.143, array<double,27> knowngammas = { 0.143,
...@@ -76,16 +77,18 @@ Double_t f_full(Double_t *x, Double_t *par) { ...@@ -76,16 +77,18 @@ Double_t f_full(Double_t *x, Double_t *par) {
float xx = x[0]; float xx = x[0];
double result, norm; double result, norm;
// Flat background // Flat background
result = par[0]; // result = par[0];
result = 0;
// Add N peaks // Add N peaks
for(int pk=0; pk<numPeaks; pk++){ for(int pk=0; pk<numPeaks; pk++){
result += (par[3+(pk*3)]/(par[1+(pk*3)]*sqrt(2*pi))) 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; 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 minFit=-1.0, maxFit=8.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3); double binWidth = hist->GetXaxis()->GetBinWidth(3);
double sigma = 0.14; double sigma = 0.14;
...@@ -124,6 +127,17 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){ ...@@ -124,6 +127,17 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth"); 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 //Build background function
TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit); TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
bg->SetLineColor(kGreen); bg->SetLineColor(kGreen);
...@@ -137,23 +151,19 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){ ...@@ -137,23 +151,19 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
for(int i=0; i<numPeaks; i++) { for(int i=0; i<numPeaks; i++) {
full->FixParameter((i*3)+1,sigma); full->FixParameter((i*3)+1,sigma);
full->FixParameter((i*3)+2,means.at(i)); 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->SetParameter((i*3)+3,1e1);
full->SetParLimits((i*3)+3,0.0,1e5); 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.,40.); /* FOR TOTAL SPECTRUM FITTING */
full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */ //full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
//full->SetParLimits(0,0.,1.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */ //full->FixParameter(0,0.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
//full->FixParameter(9,0.); //?? 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 */
// Specific limits //full->FixParameter(0,slideshift); /* FOR WHEN PAR[0] IS FIXED ENERGY CENTRIOD SLIDER */
// Set max 279 counts to 100 full->SetParameter(52,0.39); /* SET 5.8MeV SIGMA */
//full->FixParameter(9,0.0); full->SetParLimits(52,0.34,0.44); /* SET 5.8MeV SIGMA */
//full->SetParLimits(9,0.0,1e2); // Doesnt work???
allPeaks[14]->SetLineColor(kOrange); allPeaks[14]->SetLineColor(kOrange);
//Fit full function to histogram //Fit full function to histogram
...@@ -211,5 +221,5 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){ ...@@ -211,5 +221,5 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
void FitKnownPeaks(TH1F* hist){ void FitKnownPeaks(TH1F* hist){
//Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen //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);
} }
Projects/e793s/macro/Plots_47Kdp.C
+ 6
− 16
View file @ 2325c1f4
#include "DefineColours.h"
#include "GausFit.h" #include "GausFit.h"
#include "KnownPeakFitter.h" #include "KnownPeakFitter.h"
#include "DrawPlots.h" #include "DrawPlots.h"
...@@ -7,7 +8,7 @@ ...@@ -7,7 +8,7 @@
#include "ThreeBodyBreakup_FitPhaseSpace.h" #include "ThreeBodyBreakup_FitPhaseSpace.h"
void AddGammaLinesMG(TH1F* hist, double particle, double ymax){ void AddGammaLines(TH1F* hist, double particle, double ymax){
string base = "sub "; string base = "sub ";
for(int i=1; i<means.size();i++){ for(int i=1; i<means.size();i++){
...@@ -22,7 +23,7 @@ void AddGammaLinesMG(TH1F* hist, double particle, double ymax){ ...@@ -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(); hist->Draw();
for(int i=0; i<knowngammas.size();i++){ for(int i=0; i<knowngammas.size();i++){
TLine *line = new TLine(knowngammas.at(i), 0.0, knowngammas.at(i), ymax); TLine *line = new TLine(knowngammas.at(i), 0.0, knowngammas.at(i), ymax);
...@@ -41,24 +42,13 @@ void Plots_47Kdp(){ ...@@ -41,24 +42,13 @@ void Plots_47Kdp(){
tCentre = 2700; tRange = 400; tCentre = 2700; tRange = 400;
timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange); timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange);
det_gate = "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8"; det_gate = "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8";
reactionName = "47K(d,p)";
cout << "==============================================" << endl; cout << "==============================================" << endl;
cout << "=============== (d,p) reaction ===============" << endl; cout << "=============== (d,p) reaction ===============" << endl;
cout << "==============================================" << endl; cout << "==============================================" << endl;
cout << ""<< endl; cout << ""<< endl;
cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl; CS();
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;
cout << ""<< endl; cout << ""<< endl;
cout << "==============================================" << endl; cout << "==============================================" << endl;
......
Projects/e793s/macro/Plots_47Kdt.C
+ 66
− 1
View file @ 2325c1f4
#include "DefineColours.h"
#include "GausFit.h" #include "GausFit.h"
//#include "KnownPeakFitter.h" #include "KnownPeakFitter.h"
#include "DrawPlots.h" #include "DrawPlots.h"
//#include "CS2.h" //#include "CS2.h"
//#include "ThreeBodyBreakup.h" //#include "ThreeBodyBreakup.h"
//#include "ThreeBodyBreakup_FitPhaseSpace.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 */ /* 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(){ void Plots_47Kdt(){
/* Load graphical cut */
TFile gcIn("GraphicalCut_22Jun22.root");
TCutG* cutTritons = (TCutG*) gcIn.FindObjectAny("cutTritons");
LoadChain47Kdt(); LoadChain47Kdt();
gStyle->SetOptStat("nemMrRi"); gStyle->SetOptStat("nemMrRi");
tCentre = 2750; tRange = 350; tCentre = 2750; tRange = 350;
timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange); timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange);
det_gate = "MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<5"; det_gate = "MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<5";
reactionName = "47K(d,t)";
cout << "==============================================" << endl; cout << "==============================================" << endl;
cout << "=============== (d,t) reaction ===============" << endl; cout << "=============== (d,t) reaction ===============" << endl;
cout << "==============================================" << endl; cout << "==============================================" << endl;
cout << " - (d,t) selection: 'cutTritons' " << endl;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment