diff --git a/NPLib/Detectors/Minos/TMinosPhysics.h b/NPLib/Detectors/Minos/TMinosPhysics.h
index 8139e7536f65b548dd3249205db2722ea3c5459f..f9ce20ac5fb8f8261674c2cb3f9d6f89a8b55741 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.h
+++ b/NPLib/Detectors/Minos/TMinosPhysics.h
@@ -160,6 +160,9 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
double GetVertexX() {return X_Vertex;} //!
double GetVertexY() {return Y_Vertex;} //!
double GetVertexZ() {return Z_Vertex;} //!
+
+ double GetDeltaVertex() {return Delta_Vertex;} //!
+ double GetTheta12() {return Theta_12;} //!
int GetNbrOfTracks(){return Tracks_P0.size();}
diff --git a/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.h b/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.h
index 542198b365f50e80f2772585114fceccc41ba771..bacf5ceaff354383ef7500b2414577ec71a3eac7 100644
--- a/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.h
@@ -156,6 +156,11 @@ class TSamuraiHodoscopePhysics : public TObject, public NPL::VDetector {
// clear the pre-treated object
void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+
+
+ vector <double> GetCharge() {return Charge;}
+ vector <double> GetTime() {return Time;}
+ vector<int> GetID() {return ID;}
// give and external TSamuraiHodoscopeData object to TSamuraiHodoscopePhysics.
// needed for online analysis for example
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
old mode 100644
new mode 100755
index a89711353f530ab9e7352f448d580bba1890d714..26d119ac77efcb47c6020b7254afe83576c1d22d
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
@@ -6,8 +6,10 @@ double devE(const double* parameter){
//Beam spot offset
TVector3 offset(parameter[0],parameter[1],parameter[2]);
- unsigned int size = pos.size();
+ //Other variable initilizations
+ unsigned int size = pos.size();
+ TVector3 MugastNormal;
double dE,Theta;
TVector3 dir;
@@ -25,24 +27,53 @@ double devE(const double* parameter){
//Particle path vector
dir=*(pos[i])-offset;
- //Detected energy, and angle of particle leaving target
- double Theta= dir.Angle(TVector3(0,0,1));
+ //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
+ }
+
+ //Angle leaving target, angle entering MUGAST & energy deposited in MUGAST
+ double ThetaTarget = dir.Angle(TVector3(0,0,1));
+ double ThetaMugast = dir.Angle(MugastNormal);
double Energy = energy[i];
- //NOTE!!! Not calucualting energy loss in Al???
+ //Energy loss in Al
+ 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[4]*micrometer, //pass through half target
- Theta); //angle leaving target
+ 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,Theta);
+ double Ex = reaction.ReconstructRelativistic(Energy,ThetaTarget);
- //Fill histogram with ERROR in Ex!
- //h->Fill(Ex-refE);
- h->Fill(Ex);
-
+ //Fill histograms with Ex
+ h->Fill(Ex);
switch(detnum[i]){
case 1:
h1->Fill(Ex);
@@ -63,70 +94,76 @@ double devE(const double* parameter){
h7->Fill(Ex);
break;
default:
- cout << "ERROR! Invalid detnum: " << detnum[i] << " @" << i << endl;
+ cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << 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]
+ << "\t Thickness: " << parameter[3] << " um"
<< 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(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
}
////////////////////////////////////////////////////////////////////////////////
void MinimizeBeamSpot(){
- // Read data in
+ //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};
+ //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
+ //Function with 4 parameter XYZ and Target thickness
auto func = ROOT::Math::Functor(&devE,4);
- // Minimizer
+ //Minimizer
auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
minim->SetPrintLevel(0);
minim->SetPrecision(1e-10);
- // Set minimizer function
+ //Set minimizer function
minim->SetFunction(func);
- // Assign variable limits
+ //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);
+ minim->SetLimitedVariable(3,"T",parameter[3],0.01,4.76*0.5,4.76*1.5);
- // Don't draw iterations of minimizer
+ //Don't draw iterations of minimizer
flagDraw = 0;
- // Shrink it, babeyyy
+ //Shrink it, babeyyy
minim->Minimize();
- // Draw minimal value
+ //Draw minimal value
flagDraw = 1;
- // Pull values from minimizer
+ //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;
+ devE(x);
cout << "========================================" << endl;
}
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
old mode 100644
new mode 100755
index 932bc8eadf5916f3b0ec0e4c1f97ec1c867e5c62..5994406553c00aa43a6cce6ec12b5a62b08ea547
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
@@ -14,8 +14,8 @@ 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 h = new TH1D("h","All MG#'s", 60,-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.);
static auto h4 = new TH1D("h4","h4", 40,-1.,1.);
@@ -125,6 +125,7 @@ void InitiliseCanvas(){
// ----- ALL -----
canv->cd(2);
h->SetStats(0);
+ h->SetLineColor(kBlack);
h->GetXaxis()->SetTitle("Ex [MeV]");
h->GetYaxis()->SetTitle("Counts");
h->Draw();