Commit 31d48849 authored by Charlie Paxman's avatar Charlie Paxman
Browse files

* Progress on e793s beam spot minimizer

  - Added grid-based parameter initilization, to avoid getting stuck in local minima
  - This means looping, but with a max of ~15s for each iteration, it's not too terrible
parent 226ed843
Pipeline #123597 passed with stages
in 6 minutes and 21 seconds
#include "MinimizeBeamSpot.h"
double devE(const double* parameter){
////////////////////////////////////////////////////////////////////////////////
double devE(const double * parameter) {
//Beam energy: 7.7 [MeV/A] * 47 [A] = 361.9 [MeV]
static NPL::Reaction reaction("47K(d,p)48K@362");
//Beam spot offset
TVector3 offset(parameter[0],parameter[1],parameter[2]);
TVector3 offset(parameter[0], parameter[1], parameter[2]);
//Other variable initilizations
unsigned int size = pos.size();
TVector3 MugastNormal;
double dE,Theta;
double dE, Theta;
TVector3 dir;
//Initilize histogram
//canv->Clear();
//canv->ResetDrawn();
h->Reset();
h1->Reset();
h2->Reset();
h3->Reset();
h4->Reset();
h5->Reset();
h7->Reset();
h -> Reset();
h1 -> Reset();
h2 -> Reset();
h3 -> Reset();
h4 -> Reset();
h5 -> Reset();
h7 -> Reset();
//Initilize results array
// 7 => Sum in 0 and them MG's in 1-6
// 5 => Mean, MeanErr, StdDev, StdDevErr, Chi2/NDF
double FitResultMatrix[7][5];
//Loop over events
for(unsigned int i = 0 ; i < size ; i++){
for (unsigned int i = 0; i < size; i++) {
//Particle path vector
dir=*(pos[i])-offset;
dir = * (pos[i]) - offset;
//Define normal vector for the MG# of detection
DetectorSwitch(detnum[i], MugastNormal);
//Angle leaving target, angle entering MUGAST & energy deposited in MUGAST
double ThetaTarget = dir.Angle(TVector3(0,0,1));
double ThetaTarget = dir.Angle(TVector3(0, 0, 1));
double ThetaMugast = dir.Angle(MugastNormal);
double Energy = energy[i];
//Energy loss in Al
Energy=Al.EvaluateInitialEnergy(
Energy, //energy after Al
0.4*micrometer, //thickness of Al
ThetaMugast); //angle impinging on MUGAST
Energy = Al.EvaluateInitialEnergy(
Energy, //energy after Al
0.4 * micrometer, //thickness of Al
ThetaMugast //angle impinging on MUGAST
);
//Energy loss in target
Energy=CD2.EvaluateInitialEnergy(
Energy, //energy after leaving target
0.5*parameter[3]*micrometer, //pass through half target
ThetaTarget); //angle leaving target
Energy = CD2.EvaluateInitialEnergy(
Energy, //energy after leaving target
0.5 * parameter[3] * micrometer, //pass through half target
ThetaTarget //angle leaving target
);
//Final value of Ex
double Ex = reaction.ReconstructRelativistic(Energy,ThetaTarget);
double Ex = reaction.ReconstructRelativistic(Energy, ThetaTarget);
//Fill histograms with Ex
h->Fill(Ex);
h -> Fill(Ex);
DetectorSwitch(detnum[i], Ex);
}
//Initilise, Draw & Fit histograms
InitiliseCanvas(FitResultMatrix);
//Write vals to screen
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;
if (flagDraw) {
cout << "==================================================" << endl;
}
//Adapt the metric as needed
return sqrt( pow(FitResultMatrix[mgSelect][0]-refE,2) + pow(0.1*FitResultMatrix[mgSelect][2],2) );
double multiplier = 0.07; //0.1;
double metric = sqrt(pow(FitResultMatrix[mgSelect][0] - refE, 2)
+ pow(multiplier * FitResultMatrix[mgSelect][2], 2));
WriteToCout(FitResultMatrix[mgSelect], parameter[3], metric);
// double multiplier = 0.2;
// double metric = abs(FitResultMatrix[1][0]-refE) + multiplier*FitResultMatrix[1][2]
// + abs(FitResultMatrix[2][0]-refE) + multiplier*FitResultMatrix[2][2]
// + abs(FitResultMatrix[3][0]-refE) + multiplier*FitResultMatrix[3][2]
// + abs(FitResultMatrix[4][0]-refE) + multiplier*FitResultMatrix[4][2]
// + abs(FitResultMatrix[5][0]-refE) + multiplier*FitResultMatrix[5][2]
// + abs(FitResultMatrix[6][0]-refE) + multiplier*FitResultMatrix[6][2];
if (flagDraw) { WriteToFile(FitResultMatrix[mgSelect], parameter, metric); }
return metric;
}
////////////////////////////////////////////////////////////////////////////////
void MinimizeBeamSpot(){
////////////////////////////////////////////////////////////////////////////////
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
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);
//Function with 4 parameter XYZ and Target thickness
auto func = ROOT::Math::Functor(&devE,4);
//Initilise minimizer
auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
minim->SetPrintLevel(0);
minim->SetPrecision(1e-10);
//Set minimizer function
minim->SetFunction(func);
//Assign variable limits
minim->SetLimitedVariable(0,"X",parameter[0],0.01,-10,10);
minim->SetLimitedVariable(1,"Y",parameter[1],0.01,-10,10);
minim->SetLimitedVariable(2,"Z",parameter[2],0.01,-5,5);
minim->SetLimitedVariable(3,"T",parameter[3],0.01,4.76*0.5,4.76*1.5);
//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();
cout << "==================================================" << endl;
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] << " 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;
if (mgSelect >= 7) {
cout << "ERROR!! INVALID SELECTION" << endl;
return;
}
//Open output file
file.open("gridMinResults.txt", ios::out);
file << "minimX \tminimY \tminimZ \tminimT \tMetric \tMean \tMeanErr \tStdDev \tStdDvErr\tChi2/NDF" << endl;
file << setprecision(6) << fixed;
//Start timer
auto start = high_resolution_clock::now();
//TESTING: Grid method of avoiding local minima
for (int x = 0; x < 2; x++) {
for (int y = 0; y < 2; y++) {
for (int z = 0; z < 2; z++) {
for (int t = 0; t < 5; t++) {
//Start with beam (0,0,0) and 4.76um 0.5mg/c2 target
//double parameter[4] = {0.0, 0.0, 0.0, 4.76};
double parameter[4] = {
9.0 - (9.0 * x),
9.0 - (9.0 * y),
9.0 - (9.0 * z),
(1.5 - (0.25 * t)) * 4.75
};
//Don't draw iterations of minimizer
flagDraw = 0;
gROOT -> SetBatch(kTRUE);
//Initial pass through
devE(parameter);
//Function with 4 parameter XYZ and Target thickness
auto func = ROOT::Math::Functor( & devE, 4);
//Initilise minimizer
auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
minim -> SetMaxFunctionCalls(1000000); // used by Minuit and Minuit2
minim -> SetMaxIterations(1000000); // used by GSL
minim -> SetPrintLevel(3);
minim -> SetPrecision(1e-10);
//Set minimizer function
minim -> SetFunction(func);
//Assign variable limits
minim -> SetLimitedVariable(0, "X", parameter[0], 0.01, -10, 10);
minim -> SetLimitedVariable(1, "Y", parameter[1], 0.01, -10, 10);
minim -> SetLimitedVariable(2, "Z", parameter[2], 0.01, -5, 5);
minim -> SetLimitedVariable(3, "T", parameter[3], 0.01, 4.76 * 0.5, 4.76 * 1.5);
//Don't draw iterations of minimizer
flagDraw = 0;
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();
cout << "==================================================" << endl;
cout << "=---------------- FINAL PEAK FITS ---------------=" << endl;
cout << "==================================================" << endl;
devE(x);
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] << " 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;
}
}
}
}
//Stop timer
auto stop = high_resolution_clock::now();
//Close output file
file.close();
//Program runtime
auto duration = duration_cast < seconds > (stop - start);
cout << " Runtime = " << duration.count() << " s" << endl;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
#include "TRandom2.h"
#include "TError.h"
#include <iostream>
#include <chrono>
using namespace std;
using namespace std::chrono;
////////////////////////////////////////////////////////////////////////////////
/* Global */
//Various numbers and objects
double refE = 0.143; // the energy of the selected states
vector<TVector3*> pos;
vector<double> energy;
......@@ -12,17 +20,14 @@ vector<int> detnum;
unsigned int mgSelect = 10;
NPL::EnergyLoss CD2("proton_CD2.G4table","G4Table",100);
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);
//Output files
TFile *histfile = new TFile("./gridHistograms.root", "RECREATE");
ofstream file;
int writeCount = 0;
//Histograms
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.);
......@@ -34,7 +39,7 @@ static auto h7 = new TH1D("h7","h7", 40,-1.,1.);
////////////////////////////////////////////////////////////////////////////////
void LoadFile(){
// Open XYZE gamma-gated file
ifstream file("XYZE_Full_02June.txt");
ifstream file("XYZE_Full_09June_MG3.txt");
if(!file.is_open()){
cout << "fail to load file" << endl;
exit(1);
......@@ -94,7 +99,7 @@ void DetectorSwitch(int MG, TVector3& Normal ){
default:
cout << "ERROR:: Invalid DetNum " << MG << endl;
return 1; // Exit code
}
}
}
////////////////////////////////////////////////////////////////////////////////
//Overloaded function definiton; this is for filling individual Ex histograms
......@@ -121,7 +126,42 @@ void DetectorSwitch(int MG, double Ex){
default:
cout << "ERROR:: Invalid DetNum " << MG << endl;
return 1; // Exit code
}
}
}
////////////////////////////////////////////////////////////////////////////////
void WriteToCout(double* result, double thick, double metric){
cout
<< "Mean: "
<< result[0]
<< " +/- "
<< result[1]
<< " StdDev: "
<< result[2]
<< " +/- "
<< result[3]
<< " Thick: "
<< thick
<< " um"
<< " Chi2/NDF = "
<< result[4]
<< " Metric: "
<< metric
<< endl;
}
////////////////////////////////////////////////////////////////////////////////
void WriteToFile(double* result, const double* parameter, double metric){
file
<< parameter[0] << "\t"
<< parameter[1] << "\t"
<< parameter[2] << "\t"
<< parameter[3] << "\t"
<< metric << "\t"
<< result[0] << "\t"
<< result[1] << "\t"
<< result[2] << "\t"
<< result[3] << "\t"
<< result[4]
<< endl;
}
////////////////////////////////////////////////////////////////////////////////
void DrawOneHistogram(TH1D* hist, int mg, int colour, int fill, double *FitResultMatrixMG){
......@@ -178,7 +218,14 @@ void InitiliseCanvas(double FitResultMatrix[7][5]){
DrawOneHistogram(h4, 4, 840, 3444, FitResultMatrix[4]);
DrawOneHistogram(h5, 5, 600, 3544, FitResultMatrix[5]);
DrawOneHistogram(h7, 6, 880, 3644, FitResultMatrix[6]);
canv->Update();
TLine *line=new TLine(0.143,0.0,0.143,75.0);
line->SetLineColor(kBlack);
line->SetLineStyle(7);
line->Draw();
//Format legend
auto legend = new TLegend(0.15,0.7,0.35,0.9);
legend->AddEntry(h1,"MUGAST 1","f");
......@@ -200,7 +247,14 @@ void InitiliseCanvas(double FitResultMatrix[7][5]){
//Refresh
gPad->Update();
if(!flagDraw){delete canv;}
if(flagDraw){
//writeCount++;
histfile->cd();
canv->Write("Minimum#");
}
else {
delete canv;
}
}
////////////////////////////////////////////////////////////////////////////////
/***************************************************************************
* Charlie Paxman (cp00474@surrey.ac.uk) *
***************************************************************************/
//-------------------------------
//C++
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <iomanip>
//Root
//#include <TVector3.h>
//NPTool
//#include "NPEnergyLoss.h"
//#include "NPReaction.h"
//#include "NPSystemOfUnits.h"
using namespace std;
//-------------------------------
double Minimize(vector<double>, vector<vector<double>>, int, int);
//-------------------------------
void FindMinimum(){
ifstream metricsFile;
vector<double> T, X, Y, Z;
vector<vector<double>> mean, sigma, metric;
/* loop temps */
double tT, tX, tY, tZ;
double tmn, tsg, tmt;
vector<double> tmnV, tsgV, tmtV;
int lines = 0;
const char* metricsFilePath = "output_Full_02June_coarse_metrics.txt";
metricsFile.open(metricsFilePath, ios::in);
if(!metricsFile){
cout << "ERROR: File not opened." << endl;
return 1;
}
// Read in the output from BeamSpot.C
while(true){
// Reset temps
tT = tX = tY = tZ = tmn = tsg = tmt = 0.0;
tmnV.clear(); tsgV.clear(); tmtV.clear();
// Read from file
metricsFile >> tT >> tX >> tY >> tZ;
for(int i=0; i<7; i++){
tmn = tsg = tmt = 0.0;
metricsFile >> tmn >> tsg >> tmt;
tmnV.push_back(tmn);
tsgV.push_back(tsg);
tmtV.push_back(tmt);
}
// Push temps to vectors
T.push_back(tT);
X.push_back(tX);
Y.push_back(tY);
Z.push_back(tZ);
mean.push_back(tmnV);
sigma.push_back(tsgV);
metric.push_back(tmtV);
// Count lines
lines++;
// Break out
if(metricsFile.eof()){break;}
}
cout << "==================================================" << endl;
cout << "= Reading " << metricsFilePath << endl;
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;
unsigned int selectVar;
cin >> selectVar;
if(selectVar==7){selectVar=6;} // Correct the input for MG7
if(selectVar>=7){
cout << " FAIL! Invalid selection" << endl;
cout << " Exiting..." << endl;
return 1;
}
cout << fixed << showpoint;
cout << "> Thickness: " << setprecision(6)
<< Minimize(T, metric, lines, selectVar)
<< " mm"
<< endl;
cout << "> X: " << setprecision(6)
<< Minimize(X, metric, lines, selectVar)
<< " mm"
<< endl;
cout << "> Y: " << setprecision(6)
<< Minimize(Y, metric, lines, selectVar)
<< " mm"
<< endl;
cout << "> Z: " << setprecision(6)
<< Minimize(Z, metric, lines, selectVar)
<< " mm"
<< endl;
//TO DO:
// - Automatically open the histogram that corresponds to the values
// above, to save time and check for unreasonable results
}
/* * * * * * * * * * * * */
/* * * * * * * * * * * * */
/* * * * * * * * * * * * */
double Minimize(vector<double> V, vector<vector<double>> metric, int lines, int metSelect){
vector<double> minVector;
double minVal;
unsigned int minIndex;
//this is pretty inefficient, but there are no speed issues so oh well
for(int i=0; i<lines-1; i++){
minVector.push_back(metric[i][metSelect]);
}
minVal = *min_element(minVector.begin(), minVector.end());
minIndex = min_element(minVector.begin(), minVector.end()) - minVector.begin();
return V[minIndex];
}
/* * * * * * * * * * * * */
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