Commit 1b55825d authored by Charlie Paxman's avatar Charlie Paxman
Browse files

* Improving e793s beam spot minimizer

 -- NOTE: GetMean() and GetStdDev() return values inconsistent with the fitting protocol.
 -- Next update will implement TFitResultPtr to rectify this
parent f45dada2
Pipeline #123098 passed with stages
in 16 minutes and 29 seconds
......@@ -100,10 +100,47 @@ double devE(const double* parameter){
}
//End loop over events
//Select MG# being minimized
double mean = 0;
double stddev = 0;
switch(mgSelect){
case 0:
mean = h->GetMean();
stddev = h->GetStdDev();
break;
case 1:
mean = h1->GetMean();
stddev = h1->GetStdDev();
break;
case 2:
mean = h2->GetMean();
stddev = h2->GetStdDev();
break;
case 3:
mean = h3->GetMean();
stddev = h3->GetStdDev();
break;
case 4:
mean = h4->GetMean();
stddev = h4->GetStdDev();
break;
case 5:
mean = h5->GetMean();
stddev = h5->GetStdDev();
break;
case 6:
mean = h7->GetMean();
stddev = h7->GetStdDev();
break;
default:
cout << "ERROR:: Invalid MG# selection! -> " << mgSelect << endl;
return 1; // Exit code
}
//Write vals to screen
cout << "Mean: " << h->GetMean()
<< "\t StdDev: " << h->GetStdDev()
<< "\t Thickness: " << parameter[3] << " um"
cout << "Mean: " << mean
<< " StdDev: " << stddev
<< " Thickness: " << parameter[3] << " um"
<< endl;
//Draw histogram(s)
......@@ -118,7 +155,7 @@ double devE(const double* parameter){
*/
//Adapt the metric as needed
return sqrt( pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
return sqrt( pow(mean-refE,2) + pow(0.1*stddev,2) );
}
////////////////////////////////////////////////////////////////////////////////
void MinimizeBeamSpot(){
......@@ -126,6 +163,19 @@ void MinimizeBeamSpot(){
//Read data in
LoadFile();
//Output formatting
cout << fixed << showpoint << setprecision(6) << showpos;
//Read in
cout << "==================================================" << endl;
cout << "=--------- SELECT TELESCOPE TO MINIMIZE ---------=" << endl;
cout << "= Type MG# of telescope metric to use, or type 0 =" << endl;
cout << "= to use the sum of all MG's =" << endl;
cout << "==================================================" << endl;
cin >> mgSelect;
if(mgSelect==7){mgSelect=6;} // Correct the input for MG7
cout << "==================================================" << endl;
//Start with beam (0,0,0) and 4.76um 0.5mg/c2 target
double parameter[4] = {0.0, 0.0, 0.0, 4.76};
devE(parameter);
......@@ -159,11 +209,16 @@ void MinimizeBeamSpot(){
//Pull values from minimizer
const double* x = minim->X();
cout << "========================================" << endl;
cout << "==================================================" << endl;
cout << "=---------------- FINAL PEAK FITS ---------------=" << endl;
cout << "==================================================" << endl;
devE(x);
cout << "==================================================" << endl;
cout << "=------------ RESULTS OF MINIMIZATION -----------=" << endl;
cout << "==================================================" << endl;
cout << "\t\tX =" << x[0] << endl;
cout << "\t\tY =" << x[1] << endl;
cout << "\t\tZ =" << x[2] << endl;
cout << "\t\tT =" << x[3] << endl;
devE(x);
cout << "========================================" << endl;
cout << "==================================================" << endl;
}
......@@ -9,12 +9,13 @@ double refE = 0.143; // the energy of the selected states
vector<TVector3*> pos;
vector<double> energy;
vector<int> detnum;
unsigned int mgSelect;
NPL::EnergyLoss CD2("proton_CD2.G4table","G4Table",100);
NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
using namespace std;
bool flagDraw = 0;
static auto h = new TH1D("h","All MG#'s", 60,-1.,1.);
static auto h = new TH1D("h","All MG#'s", 80,-1.,1.);
static auto h1 = new TH1D("h1","Individual MG#'s", 40,-1.,1.);
static auto h2 = new TH1D("h2","h2", 40,-1.,1.);
static auto h3 = new TH1D("h3","h3", 40,-1.,1.);
......@@ -70,7 +71,9 @@ void InitiliseCanvas(){
h1->SetFillStyle(3244);
h1->SetFillColor(kRed);
h1->Draw();
h1->Fit("gaus","WQ"); //add N to stop it drawing
cout << "\n==================================================" << endl;
cout << "=---------------------- MG1 ---------------------=" << endl;
h1->Fit("gaus","W"); //N=stop drawing, Q=stop writing
// ----- MG2 -----
h2->SetStats(0);
......@@ -78,7 +81,9 @@ void InitiliseCanvas(){
h2->SetFillStyle(3244);
h2->SetFillColor(kOrange);
h2->Draw("same");
h2->Fit("gaus","WQ"); //add N to stop it drawing
cout << "\n==================================================" << endl;
cout << "=---------------------- MG2 ---------------------=" << endl;
h2->Fit("gaus","W");
// ----- MG3 -----
h3->SetStats(0);
......@@ -86,23 +91,29 @@ void InitiliseCanvas(){
h3->SetFillStyle(3344);
h3->SetFillColor(kGreen);
h3->Draw("same");
h3->Fit("gaus","WQ"); //add N to stop it drawing
cout << "\n==================================================" << endl;
cout << "=---------------------- MG3 ---------------------=" << endl;
h3->Fit("gaus","W");
// ----- MG4 -----
h4->SetStats(0);
h4->SetLineColor(kTeal);
h4->SetFillStyle(3444);
h4->SetFillColor(kTeal);
h4->Draw("same");
h4->Fit("gaus","WQ"); //add N to stop it drawing
cout << "\n==================================================" << endl;
cout << "=---------------------- MG4 ---------------------=" << endl;
h4->Fit("gaus","W");
// ----- MG5 -----
h5->SetStats(0);
h5->SetLineColor(kBlue);
h5->SetFillStyle(3544);
h5->SetFillColor(kBlue);
h5->Draw("same");
h5->Fit("gaus","WQ"); //add N to stop it drawing
cout << "\n==================================================" << endl;
cout << "=---------------------- MG5 ---------------------=" << endl;
h5->Fit("gaus","W");
// ----- MG7 -----
h7->SetStats(0);
......@@ -110,7 +121,9 @@ void InitiliseCanvas(){
h7->SetFillStyle(3644);
h7->SetFillColor(kViolet);
h7->Draw("same");
h7->Fit("gaus","WQ"); //add N to stop it drawing
cout << "\n==================================================" << endl;
cout << "=---------------------- MG7 ---------------------=" << endl;
h7->Fit("gaus","W");
// Format legend
auto legend = new TLegend(0.15,0.7,0.35,0.9);
......@@ -129,7 +142,9 @@ void InitiliseCanvas(){
h->GetXaxis()->SetTitle("Ex [MeV]");
h->GetYaxis()->SetTitle("Counts");
h->Draw();
h->Fit("gaus", "WQ");
cout << "\n==================================================" << endl;
cout << "=---------------------- SUM ---------------------=" << endl;
h->Fit("gaus", "W");
gPad->Update();
}
////////////////////////////////////////////////////////////////////////////////
......
Markdown is supported
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