Commit db65046f authored by Charlie Paxman's avatar Charlie Paxman
Browse files

* Progress on e793s

 -- Now use FitResult mean & std dev for minimization, giving more accurate results
 -- Also improved the visual outputs a little
parent aab2e5c7
Pipeline #123186 passed with stages
in 4 minutes and 57 seconds
......@@ -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;
}
......@@ -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;}
}
////////////////////////////////////////////////////////////////////////////////
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