From db65046fd0984e72f17bec09bfa7fa060ef59c46 Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Tue, 8 Jun 2021 10:21:24 +0100
Subject: [PATCH] * Progress on e793s -- Now use FitResult mean & std dev for
minimization, giving more accurate results -- Also improved the visual
outputs a little
---
.../e793s/macro/BeamSpot/MinimizeBeamSpot.cxx | 144 +++----------
.../e793s/macro/BeamSpot/MinimizeBeamSpot.h | 200 ++++++++++--------
2 files changed, 150 insertions(+), 194 deletions(-)
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
index d28d08da5..a551e0c5f 100755
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
@@ -14,6 +14,8 @@ double devE(const double* parameter){
TVector3 dir;
//Initilize histogram
+ //canv->Clear();
+ //canv->ResetDrawn();
h->Reset();
h1->Reset();
h2->Reset();
@@ -22,36 +24,15 @@ double devE(const double* parameter){
h5->Reset();
h7->Reset();
+ double FitResultMatrix[7][5];
+
//Loop over events
for(unsigned int i = 0 ; i < size ; i++){
//Particle path vector
dir=*(pos[i])-offset;
//Define normal vector for the MG# of detection
- //[[[[ UPDATE WITH NEW MG POSITIONS FROM SURVEY ]]]]
- switch(detnum[i]){
- case 1:
- MugastNormal.SetXYZ(-0.453915, +0.455463, -0.765842);
- break;
- case 2:
- MugastNormal.SetXYZ(-0.642828, +0.000000, -0.766010);
- break;
- case 3:
- MugastNormal.SetXYZ(-0.454594, -0.450670, -0.768271);
- break;
- case 4:
- MugastNormal.SetXYZ(-0.002437, -0.638751, -0.769409);
- break;
- case 5:
- MugastNormal.SetXYZ(+0.452429, -0.454575, -0.767248);
- break;
- case 7:
- MugastNormal.SetXYZ(+0.443072, +0.443265, -0.779232);
- break;
- default:
- cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << i << endl;
- return 1; // Exit code
- }
+ DetectorSwitch(detnum[i], MugastNormal);
//Angle leaving target, angle entering MUGAST & energy deposited in MUGAST
double ThetaTarget = dir.Angle(TVector3(0,0,1));
@@ -74,88 +55,26 @@ double devE(const double* parameter){
//Fill histograms with Ex
h->Fill(Ex);
- switch(detnum[i]){
- case 1:
- h1->Fill(Ex);
- break;
- case 2:
- h2->Fill(Ex);
- break;
- case 3:
- h3->Fill(Ex);
- break;
- case 4:
- h4->Fill(Ex);
- break;
- case 5:
- h5->Fill(Ex);
- break;
- case 7:
- h7->Fill(Ex);
- break;
- default:
- cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << i << endl;
- return 1; // Exit code
- }
- }
- //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
+ DetectorSwitch(detnum[i], Ex);
}
-
+
+ //Initilise, Draw & Fit histograms
+ InitiliseCanvas(FitResultMatrix);
+
//Write vals to screen
- cout << "Mean: " << mean
- << " StdDev: " << stddev
- << " Thickness: " << parameter[3] << " um"
+ if(flagDraw){cout << "==================================================" << endl;}
+ cout << "Mean: " << FitResultMatrix[mgSelect][0]
+ << " +/- "
+ << FitResultMatrix[mgSelect][1]
+ << " StdDev: " << FitResultMatrix[mgSelect][2]
+ << " +/- "
+ << FitResultMatrix[mgSelect][3]
+ << " Thick: " << parameter[3] << " um"
+ << " Fit Chi2/NDF = " << FitResultMatrix[mgSelect][4]
<< endl;
- //Draw histogram(s)
- h->Draw();
- if(flagDraw){ InitiliseCanvas(); }
-
- /*
- cout << pow(h->GetMean()-refE,2) << " + "
- << pow(0.1*h->GetStdDev(),2) << " = "
- << pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2)
- << endl;
- */
-
//Adapt the metric as needed
- return sqrt( pow(mean-refE,2) + pow(0.1*stddev,2) );
+ return sqrt( pow(FitResultMatrix[mgSelect][0]-refE,2) + pow(0.1*FitResultMatrix[mgSelect][2],2) );
}
////////////////////////////////////////////////////////////////////////////////
void MinimizeBeamSpot(){
@@ -172,8 +91,8 @@ void MinimizeBeamSpot(){
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
+ 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
@@ -183,9 +102,8 @@ void MinimizeBeamSpot(){
//Function with 4 parameter XYZ and Target thickness
auto func = ROOT::Math::Functor(&devE,4);
- //Minimizer
+ //Initilise minimizer
auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
-
minim->SetPrintLevel(0);
minim->SetPrecision(1e-10);
@@ -200,12 +118,15 @@ void MinimizeBeamSpot(){
//Don't draw iterations of minimizer
flagDraw = 0;
-
+ //canv->SetBatch(kTRUE);
+ gROOT->SetBatch(kTRUE);
+
//Shrink it, babeyyy
minim->Minimize();
//Draw minimal value
flagDraw = 1;
+ gROOT->SetBatch(kFALSE);
//Pull values from minimizer
const double* x = minim->X();
@@ -213,12 +134,15 @@ void MinimizeBeamSpot(){
cout << "=---------------- FINAL PEAK FITS ---------------=" << endl;
cout << "==================================================" << endl;
devE(x);
+// canv->DrawClone();
cout << "==================================================" << endl;
cout << "=------------ RESULTS OF MINIMIZATION -----------=" << endl;
+ if(mgSelect==6){mgSelect=7;} // Correct the input for MG7
+ cout << "=------------------- USING MG " << mgSelect << " -----------------=" << 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;
+ cout << "\t\tX = " << x[0] << " mm" << endl;
+ cout << "\t\tY = " << x[1] << " mm" << endl;
+ cout << "\t\tZ = " << x[2] << " mm" << endl;
+ cout << "\t\tT = " << x[3] << " um" << endl;
cout << "==================================================" << endl;
}
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
index 03eb5feee..ef2448e5e 100755
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
@@ -15,6 +15,14 @@ NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
using namespace std;
bool flagDraw = 0;
+
+//double FitResultMatrix[7][5];
+// 7 => Sum in 0 and them MG's in 1-6
+// 5 => Mean, MeanErr, StdDev, StdDevErr, Chi2/NDF
+
+
+//TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
+
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.);
@@ -47,23 +55,104 @@ void LoadFile(){
file.close();
}
////////////////////////////////////////////////////////////////////////////////
-void FillResultMatrix(double* matrix, TFitResultPtr fit){
+void FillMatrix(double* matrix, TFitResultPtr fit){
matrix[0] = fit->Parameter(1); //Mean
matrix[1] = fit->ParError(1);
matrix[2] = fit->Parameter(2); //StdDev
matrix[3] = fit->ParError(2);
matrix[4] = fit->Chi2()/fit->Ndf(); //Chi2/NDF
- cout << "\n Mean = " << matrix[0] << " +/- " << matrix[1] << endl;
- cout << " StdDev = " << matrix[2] << " +/- " << matrix[3] << endl;
- cout << " Chi2/NDF = " << matrix[4] << endl;
+ if(flagDraw){
+ cout << "\n Mean = " << matrix[0] << " +/- " << matrix[1] << endl;
+ cout << " StdDev = " << matrix[2] << " +/- " << matrix[3] << endl;
+ cout << " Chi2/NDF = " << matrix[4] << endl;
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+//[[[[ UPDATE WITH NEW MG POSITIONS FROM SURVEY ]]]]
+//Overloaded function definiton; this is for MUGAST Normal vectors
+void DetectorSwitch(int MG, TVector3& Normal ){
+ switch(MG){
+ case 1:
+ Normal.SetXYZ(-0.453915, +0.455463, -0.765842);
+ break;
+ case 2:
+ Normal.SetXYZ(-0.642828, +0.000000, -0.766010);
+ break;
+ case 3:
+ Normal.SetXYZ(-0.454594, -0.450670, -0.768271);
+ break;
+ case 4:
+ Normal.SetXYZ(-0.002437, -0.638751, -0.769409);
+ break;
+ case 5:
+ Normal.SetXYZ(+0.452429, -0.454575, -0.767248);
+ break;
+ case 7:
+ Normal.SetXYZ(+0.443072, +0.443265, -0.779232);
+ break;
+ default:
+ cout << "ERROR:: Invalid DetNum " << MG << endl;
+ return 1; // Exit code
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+//Overloaded function definiton; this is for filling individual Ex histograms
+void DetectorSwitch(int MG, double Ex){
+ switch(MG){
+ case 1:
+ h1->Fill(Ex);
+ break;
+ case 2:
+ h2->Fill(Ex);
+ break;
+ case 3:
+ h3->Fill(Ex);
+ break;
+ case 4:
+ h4->Fill(Ex);
+ break;
+ case 5:
+ h5->Fill(Ex);
+ break;
+ case 7:
+ h7->Fill(Ex);
+ break;
+ default:
+ cout << "ERROR:: Invalid DetNum " << MG << endl;
+ return 1; // Exit code
+ }
}
////////////////////////////////////////////////////////////////////////////////
-void InitiliseCanvas(){
- double FitResultMatrix[7][5];
- // 7 => Sum in 0 and them MG's in 1-6
- // 5 => Mean, MeanErr, StdDev, StdDevErr, Chi2/NDF
+void DrawOneHistogram(TH1D* hist, int mg, int colour, int fill, double *FitResultMatrixMG){
+ //Hist settings
+ hist->SetStats(0);
+ hist->SetLineColor(colour);
+ hist->SetFillStyle(fill);
+ hist->SetFillColor(colour);
+ hist->Draw("same");
+ if (flagDraw){
+ //Header
+ cout << noshowpos;
+ cout << "\n==================================================" << endl;
+ if (mg==6){
+ cout << "=---------------------- MG7 ---------------------=" << endl;
+ } else if (mg==0) {
+ cout << "=---------------------- SUM ---------------------=" << endl;
+ } else {
+ cout << "=---------------------- MG" << mg << " ---------------------=" << endl;
+ }
+ cout << showpos;
+ }
+
+ TFitResultPtr fit = hist->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
+ FillMatrix(FitResultMatrixMG,fit);
+}
+////////////////////////////////////////////////////////////////////////////////
+void InitiliseCanvas(double FitResultMatrix[7][5]){
+
+ //Canvas setup
TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
gStyle->SetOptStat(0);
canv->Divide(2,1,0.005,0.005,0);
@@ -75,79 +164,22 @@ void InitiliseCanvas(){
canv->cd(2)->SetBottomMargin(0.15);
gPad->SetTickx();
gPad->SetTicky();
-
+
+ //Histogram setup - Individual
canv->cd(1);
h1->SetMaximum(75.);
h1->GetXaxis()->SetTitle("Ex [MeV]");
h1->GetYaxis()->SetTitle("Counts");
-
- // ----- MG1 -----
- h1->SetStats(0);
- h1->SetLineColor(kRed);
- h1->SetFillStyle(3244);
- h1->SetFillColor(kRed);
- h1->Draw();
- cout << "\n==================================================" << endl;
- cout << "=---------------------- MG1 ---------------------=" << endl;
- TFitResultPtr f1 = h1->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[1],f1);
-
- // ----- MG2 -----
- h2->SetStats(0);
- h2->SetLineColor(kOrange);
- h2->SetFillStyle(3244);
- h2->SetFillColor(kOrange);
- h2->Draw("same");
- cout << "\n==================================================" << endl;
- cout << "=---------------------- MG2 ---------------------=" << endl;
- TFitResultPtr f2 = h2->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[2],f2);
-
- // ----- MG3 -----
- h3->SetStats(0);
- h3->SetLineColor(kGreen);
- h3->SetFillStyle(3344);
- h3->SetFillColor(kGreen);
- h3->Draw("same");
- cout << "\n==================================================" << endl;
- cout << "=---------------------- MG3 ---------------------=" << endl;
- TFitResultPtr f3 = h3->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[3],f3);
-
- // ----- MG4 -----
- h4->SetStats(0);
- h4->SetLineColor(kTeal);
- h4->SetFillStyle(3444);
- h4->SetFillColor(kTeal);
- h4->Draw("same");
- cout << "\n==================================================" << endl;
- cout << "=---------------------- MG4 ---------------------=" << endl;
- TFitResultPtr f4 = h4->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[4],f4);
- // ----- MG5 -----
- h5->SetStats(0);
- h5->SetLineColor(kBlue);
- h5->SetFillStyle(3544);
- h5->SetFillColor(kBlue);
- h5->Draw("same");
- cout << "\n==================================================" << endl;
- cout << "=---------------------- MG5 ---------------------=" << endl;
- TFitResultPtr f5 = h5->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[5],f5);
-
- // ----- MG7 -----
- h7->SetStats(0);
- h7->SetLineColor(kViolet);
- h7->SetFillStyle(3644);
- h7->SetFillColor(kViolet);
- h7->Draw("same");
- cout << "\n==================================================" << endl;
- cout << "=---------------------- MG7 ---------------------=" << endl;
- TFitResultPtr f7 = h7->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[6],f7);
-
- // Format legend
+ //Histogram draw - Individual
+ DrawOneHistogram(h1, 1, 632, 3244, FitResultMatrix[1]);
+ DrawOneHistogram(h2, 2, 800, 3244, FitResultMatrix[2]);
+ DrawOneHistogram(h3, 3, 416, 3344, FitResultMatrix[3]);
+ DrawOneHistogram(h4, 4, 840, 3444, FitResultMatrix[4]);
+ DrawOneHistogram(h5, 5, 600, 3544, FitResultMatrix[5]);
+ DrawOneHistogram(h7, 6, 880, 3644, FitResultMatrix[6]);
+
+ //Format legend
auto legend = new TLegend(0.15,0.7,0.35,0.9);
legend->AddEntry(h1,"MUGAST 1","f");
legend->AddEntry(h2,"MUGAST 2","f");
@@ -157,18 +189,18 @@ void InitiliseCanvas(){
legend->AddEntry(h7,"MUGAST 7","f");
legend->Draw();
- // ----- ALL -----
+ //Histogram setup - Sum
canv->cd(2);
- h->SetStats(0);
- h->SetLineColor(kBlack);
h->GetXaxis()->SetTitle("Ex [MeV]");
h->GetYaxis()->SetTitle("Counts");
- h->Draw();
- cout << "\n==================================================" << endl;
- cout << "=---------------------- SUM ---------------------=" << endl;
- TFitResultPtr f0 = h->Fit("gaus","WQS"); //N=stop drawing, Q=stop writing
- FillResultMatrix(FitResultMatrix[0],f0);
+
+ //Histogram draw - Sum
+ DrawOneHistogram(h, 0, 1, 0, FitResultMatrix[0]);
+
+ //Refresh
gPad->Update();
+
+ if(!flagDraw){delete canv;}
}
////////////////////////////////////////////////////////////////////////////////
--
GitLab