From 737addac37e6f994c66dbd193ac852c35d23a343 Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Thu, 3 Jun 2021 16:02:09 +0100
Subject: [PATCH] * Revamped beam spot minimiser - Now using ROOT::Math
minimizer. - Requires ROOT's minuit2 build flag to be on.
---
Projects/e793s/macro/BeamSpot/EventReader.C | 5 +-
.../e793s/macro/BeamSpot/MinimizeBeamSpot.cxx | 132 +++++++++++++++++
.../e793s/macro/BeamSpot/MinimizeBeamSpot.h | 135 ++++++++++++++++++
3 files changed, 269 insertions(+), 3 deletions(-)
create mode 100644 Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
create mode 100644 Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
diff --git a/Projects/e793s/macro/BeamSpot/EventReader.C b/Projects/e793s/macro/BeamSpot/EventReader.C
index ff1ec0d0f..e88744fa2 100755
--- a/Projects/e793s/macro/BeamSpot/EventReader.C
+++ b/Projects/e793s/macro/BeamSpot/EventReader.C
@@ -20,8 +20,7 @@ using namespace std;
void EventReader(){
// Read ROOT file and pull the tree
- //auto DataFile = new TFile("../../../../Outputs/Analysis/47K_RawEnergyBranch_Run63_May11.root", "READ");
- auto DataFile = new TFile("../../../../Outputs/Analysis/47K_RawEnergyBranch_Full_May11.root", "READ");
+ auto DataFile = new TFile("../../../../Outputs/Analysis/47K_Full_02June.root", "READ");
auto PhysTree = (TTree*) DataFile->FindObjectAny("PhysicsTree");
// Initilise the Mugast branch
@@ -56,7 +55,7 @@ void EventReader(){
// Open output file
ofstream outfile;
- outfile.open("./XYZE_gammaGated_Full.txt");
+ outfile.open("./XYZE_Full_02June.txt");
// Loop on entries
for(unsigned int i=0; i<numEntries; i++){
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
new file mode 100644
index 000000000..a89711353
--- /dev/null
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
@@ -0,0 +1,132 @@
+#include "MinimizeBeamSpot.h"
+
+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]);
+ unsigned int size = pos.size();
+
+ double dE,Theta;
+ TVector3 dir;
+
+ //Initilize histogram
+ h->Reset();
+ h1->Reset();
+ h2->Reset();
+ h3->Reset();
+ h4->Reset();
+ h5->Reset();
+ h7->Reset();
+
+ //Loop over events
+ for(unsigned int i = 0 ; i < size ; i++){
+ //Particle path vector
+ dir=*(pos[i])-offset;
+
+ //Detected energy, and angle of particle leaving target
+ double Theta= dir.Angle(TVector3(0,0,1));
+ double Energy = energy[i];
+
+ //NOTE!!! Not calucualting energy loss in Al???
+ //Energy loss in target
+ Energy=CD2.EvaluateInitialEnergy(
+ Energy, //energy after leaving target
+ 0.5*parameter[4]*micrometer, //pass through half target
+ Theta); //angle leaving target
+
+ //Final value of Ex
+ double Ex = reaction.ReconstructRelativistic(Energy,Theta);
+
+ //Fill histogram with ERROR in Ex!
+ //h->Fill(Ex-refE);
+ 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] << " @" << i << endl;
+ return 1; // Exit code
+ }
+
+ }
+ //End loop over events
+
+ //Write vals to screen
+ cout << "Mean: " << h->GetMean()
+ << "\t StdDev: " << h->GetStdDev()
+ << "\t Thickness??: " << parameter[4]
+ << endl;
+
+ //Draw histogram(s)
+ h->Draw();
+ if(flagDraw){ InitiliseCanvas(); }
+
+ //Adapt the metric as needed
+ return sqrt( pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
+}
+////////////////////////////////////////////////////////////////////////////////
+void MinimizeBeamSpot(){
+
+ // Read data in
+ LoadFile();
+
+ // Start with beam (0,0,0) and 4.7um 0.5mg/c2 target
+ double parameter[4] = {0.0, 0.0, 0.0, 4.7};
+ devE(parameter);
+
+ // Function with 4 parameter XYZ and Target thickness
+ auto func = ROOT::Math::Functor(&devE,4);
+
+ // 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.7-3.0, 4.7+3.0);
+
+ // Don't draw iterations of minimizer
+ flagDraw = 0;
+
+ // Shrink it, babeyyy
+ minim->Minimize();
+
+ // Draw minimal value
+ flagDraw = 1;
+
+ // Pull values from minimizer
+ const double* x = minim->X();
+ 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 << "Minimum: " << devE(x) << endl;
+ cout << "========================================" << endl;
+}
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
new file mode 100644
index 000000000..932bc8ead
--- /dev/null
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
@@ -0,0 +1,135 @@
+#include "Math/Minimizer.h"
+#include "Math/Factory.h"
+#include "Math/Functor.h"
+#include "TRandom2.h"
+#include "TError.h"
+#include <iostream>
+
+double refE = 0.143; // the energy of the selected states
+vector<TVector3*> pos;
+vector<double> energy;
+vector<int> detnum;
+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","h", 80,-1.,1.);
+static auto h1 = new TH1D("h1","h1", 40,-1.,1.);
+static auto h2 = new TH1D("h2","h2", 40,-1.,1.);
+static auto h3 = new TH1D("h3","h3", 40,-1.,1.);
+static auto h4 = new TH1D("h4","h4", 40,-1.,1.);
+static auto h5 = new TH1D("h5","h5", 40,-1.,1.);
+static auto h7 = new TH1D("h7","h7", 40,-1.,1.);
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadFile(){
+ // Open XYZE gamma-gated file
+ ifstream file("XYZE_Full_02June.txt");
+ if(!file.is_open()){
+ cout << "fail to load file" << endl;
+ exit(1);
+ }
+ else {
+ cout << "Success opening file" << endl;
+ }
+
+ // Read in
+ int mg;
+ double x,y,z,e;
+ while(file >> mg >> x >> y >> z >> e ){
+ auto p = new TVector3(x,y,z);
+ detnum.push_back(mg);
+ pos.push_back(p);
+ energy.push_back(e);
+ }
+ file.close();
+}
+////////////////////////////////////////////////////////////////////////////////
+void InitiliseCanvas(){
+ TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
+ gStyle->SetOptStat(0);
+ canv->Divide(2,1,0.005,0.005,0);
+ canv->cd(1)->SetLeftMargin(0.15);
+ canv->cd(1)->SetBottomMargin(0.15);
+ gPad->SetTickx();
+ gPad->SetTicky();
+ canv->cd(2)->SetLeftMargin(0.15);
+ canv->cd(2)->SetBottomMargin(0.15);
+ gPad->SetTickx();
+ gPad->SetTicky();
+
+ 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();
+ h1->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG2 -----
+ h2->SetStats(0);
+ h2->SetLineColor(kOrange);
+ h2->SetFillStyle(3244);
+ h2->SetFillColor(kOrange);
+ h2->Draw("same");
+ h2->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- MG3 -----
+ h3->SetStats(0);
+ h3->SetLineColor(kGreen);
+ h3->SetFillStyle(3344);
+ h3->SetFillColor(kGreen);
+ h3->Draw("same");
+ h3->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // ----- 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
+
+ // ----- 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
+
+ // ----- MG7 -----
+ h7->SetStats(0);
+ h7->SetLineColor(kViolet);
+ h7->SetFillStyle(3644);
+ h7->SetFillColor(kViolet);
+ h7->Draw("same");
+ h7->Fit("gaus","WQ"); //add N to stop it drawing
+
+ // 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");
+ legend->AddEntry(h3,"MUGAST 3","f");
+ legend->AddEntry(h4,"MUGAST 4","f");
+ legend->AddEntry(h5,"MUGAST 5","f");
+ legend->AddEntry(h7,"MUGAST 7","f");
+ legend->Draw();
+
+ // ----- ALL -----
+ canv->cd(2);
+ h->SetStats(0);
+ h->GetXaxis()->SetTitle("Ex [MeV]");
+ h->GetYaxis()->SetTitle("Counts");
+ h->Draw();
+ h->Fit("gaus", "WQ");
+ gPad->Update();
+}
+////////////////////////////////////////////////////////////////////////////////
+
--
GitLab