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

* Small general updates to e793s

parent 7a741ffb
Pipeline #135433 passed with stages
in 8 minutes and 53 seconds
......@@ -45,11 +45,11 @@ void Analysis::Init() {
///////////////////////////////////////////////////////////////////////////////
// if(NPOptionManager::getInstance()->HasDefinition("simulation")){
cout << " == == == == SIMULATION == == == ==" << endl;
isSim=true;
// cout << " == == == == SIMULATION == == == ==" << endl;
// isSim=true;
// } else {
// cout << " == == == == EXPERIMENT == == == ==" << endl;
// isSim=false;
cout << " == == == == EXPERIMENT == == == ==" << endl;
isSim=false;
// }
agata_zShift=51*mm;
......@@ -58,6 +58,7 @@ void Analysis::Init() {
if(isSim){
Initial = new TInitialConditions();
ReactionConditions = new TReactionConditions();
RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
}
// initialize input and output branches
......@@ -138,8 +139,8 @@ void Analysis::Init() {
GATCONF_Counter[i] = 0 ;
}
ThetaCM_detected->Sumw2();
ThetaLab_detected->Sumw2();
// ThetaCM_detected->Sumw2();
// ThetaLab_detected->Sumw2();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -150,8 +151,11 @@ void Analysis::TreatEvent(){
if(isSim){
//ThetaCM_emmitted->Fill(Initial->GetThetaCM(0));
ThetaLab_emmitted->Fill(Initial->GetThetaLab_WorldFrame(0));
ThetaCM_emmitted->Fill(ReactionConditions->GetThetaCM());
//ThetaLab_emmitted->Fill(Initial->GetThetaLab_WorldFrame(0));
ThetaLab_emmitted->Fill(ReactionConditions->GetTheta(0));
/*
TVector3 labDir = Initial->GetParticleDirection(0);
TVector3 comDir = labDir;
......@@ -159,14 +163,14 @@ void Analysis::TreatEvent(){
comDir.SetX( labDir.X() / (2*comDir.Z()) );
comDir.SetY( labDir.Y() / (2*comDir.Z()) );
/*
cout << Initial->GetThetaLab_WorldFrame(0) << " "
<< labDir.Theta()/deg << " "
<< comDir.Theta()/deg << " "
<< endl;
*/
ThetaCM_emmitted->Fill(comDir.Theta()/deg);
*/
}
GATCONF_MASTER=ML->GetCalibratedValue("GATCONF_MASTER");
......@@ -177,16 +181,19 @@ void Analysis::TreatEvent(){
}
}
// if(isSim){
// OriginalELab = ReactionConditions->GetKineticEnergy(0);
// OriginalThetaLab = ReactionConditions->GetTheta(0);
// BeamEnergy = ReactionConditions->GetBeamEnergy();
// }
if(isSim){
OriginalELab = ReactionConditions->GetKineticEnergy(0);
OriginalThetaLab = ReactionConditions->GetTheta(0);
BeamEnergy = ReactionConditions->GetBeamEnergy();
}
TVector3 BeamDirection(0.,0.,1.);
BeamImpact = TVector3(XBeam,YBeam,m_DetectorManager->GetTargetZ());
ParticleMult=M2->Si_E.size()+MG->DSSD_E.size();
//// ParticleMult=M2->Si_E.size()+MG->DSSD_E.size();
ParticleMult=M2->Si_E.size();////+MG->DSSD_E.size();
//ParticleMult=M2->Si_E.size();
// FinalBeamEnergy=BeamCD2.Slow(OriginalBeamEnergy,0.5*TargetThickness,BeamDirection.Angle(TVector(0,0,1)));
//reaction.SetBeamEnergy(FinalBeamEnergy);
......@@ -203,11 +210,15 @@ void Analysis::TreatEvent(){
// MUST2
int TelescopeNumber = M2->TelescopeNumber[countMust2];
if(isSim){
if(TelescopeNumber==5){
//ThetaCM_detected->Fill(Initial->GetThetaCM(0));
ThetaLab_detected->Fill(Initial->GetThetaLab_WorldFrame(0));
ThetaCM_detected->Fill(ReactionConditions->GetThetaCM());
//ThetaLab_detected->Fill(Initial->GetThetaLab_WorldFrame(0));
ThetaLab_detected->Fill(ReactionConditions->GetTheta(0));
/*
TVector3 labDir = Initial->GetParticleDirection(0);
TVector3 comDir = labDir;
......@@ -215,16 +226,15 @@ void Analysis::TreatEvent(){
comDir.SetX( labDir.X() / (2*comDir.Z()) );
comDir.SetY( labDir.Y() / (2*comDir.Z()) );
/*
// Check the angles agree
cout << Initial->GetThetaLab_WorldFrame(0) << " "
<< labDir.Theta()/deg << " "
<< comDir.Theta()/deg << " "
<< endl;
*/
ThetaCM_detected->Fill(comDir.Theta()/deg);
*/
}
}
......@@ -330,7 +340,7 @@ void Analysis::TreatEvent(){
PhiLab.push_back(philab_tmp/deg);
}
/*
////////////////////////////////////////////////////////////////////////////
//////////////////////////////// LOOP on MUGAST ////////////////////////////
////////////////////////////////////////////////////////////////////////////
......@@ -342,6 +352,14 @@ void Analysis::TreatEvent(){
thetalab_tmp = 0;
philab_tmp = 0;
TVector3 HitDirection = MG->GetPositionOfInteraction(countMugast) - BeamImpact;
//TVector3 tempVector;
//if(MG->TelescopeNumber[0]==3){
// tempVector = {MG->GetPositionOfInteraction(countMugast).X(),
// MG->GetPositionOfInteraction(countMugast).Y(),
// MG->GetPositionOfInteraction(countMugast).Z()+1.0}; //add 1mm to MG3
// HitDirection = tempVector - BeamImpact;
// if(warning){cout << "!!! EDITING MG3 !!!" << endl; warning=false;}
//}
thetalab_tmp = HitDirection.Angle(BeamDirection);
philab_tmp = HitDirection.Phi();
......@@ -396,7 +414,8 @@ void Analysis::TreatEvent(){
}//end loop Mugast
*/
////////////////////////////////////////////////////////////////////////////
///////////////////////////////// LOOP on AGATA ////////////////////////////
......@@ -433,7 +452,6 @@ void Analysis::TreatEvent(){
}
*/
/********
// Agata add back is not always multiplicity 1 ?? NO, not necessarily!
for(int i=0; i<nbAdd; i++){
//if(nbAdd==1){
......@@ -478,9 +496,8 @@ void Analysis::TreatEvent(){
MWT[i] = (MW_T[i]-offset[j])/slope[j];
}
********/
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
......@@ -567,11 +584,11 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("Y",&Y);
RootOutput::getInstance()->GetTree()->Branch("Z",&Z);
RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
RootOutput::getInstance()->GetTree()->Branch("MG_T",MG_T);
RootOutput::getInstance()->GetTree()->Branch("MG_E",MG_E);
RootOutput::getInstance()->GetTree()->Branch("MG_X",MG_X);
RootOutput::getInstance()->GetTree()->Branch("MG_Y",MG_Y);
RootOutput::getInstance()->GetTree()->Branch("MG_D",MG_D);
//// RootOutput::getInstance()->GetTree()->Branch("MG_T",MG_T);
//// RootOutput::getInstance()->GetTree()->Branch("MG_E",MG_E);
//// RootOutput::getInstance()->GetTree()->Branch("MG_X",MG_X);
//// RootOutput::getInstance()->GetTree()->Branch("MG_Y",MG_Y);
//// RootOutput::getInstance()->GetTree()->Branch("MG_D",MG_D);
// Vamos
RootOutput::getInstance()->GetTree()->Branch("LTS",&LTS,"LTS/l");
......@@ -720,11 +737,14 @@ void Analysis::InitInputBranch(){
if(isSim){
//RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
//RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",
&Initial);
RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",
&Interaction);
//RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
//RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",
&ReactionConditions);
RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions",
&ReactionConditions);
}
}
......@@ -780,6 +800,7 @@ void Analysis::SetBranchStatus(){
if(isSim){
RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("InteractionCoordinates",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
}
......
......@@ -31,6 +31,7 @@
#include "TF1.h"
#include "TInitialConditions.h"
#include "TReactionConditions.h"
#include "TInteractionCoordinates.h"
#include "TMust2Physics.h"
#include "TMugastPhysics.h"
#include "TCATSPhysics.h"
......@@ -231,10 +232,14 @@ class Analysis: public NPL::VAnalysis{
//TCATSPhysics* CATS;
TModularLeafPhysics* ML;
bool warning=true;
//Simulation Stuff
bool isSim;
bool writetoscreen;
TInitialConditions* Initial;
TInteractionCoordinates* Interaction;
TReactionConditions* ReactionConditions;
// Beam object
......
......@@ -2,7 +2,7 @@
////////////////////////////////////////////////////////////////////////////////
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");
static NPL::Reaction reaction("47K(d,p)48K@355");
//Beam spot offset
TVector3 offset(parameter[0], parameter[1], parameter[2]);
......@@ -21,6 +21,18 @@ double devE(const double * parameter) {
h4 -> Reset();
h5 -> Reset();
h7 -> Reset();
hT -> Reset();
//Now that initial range is wide, crop to single peak
if(refE==0.143){
h->SetAxisRange(-1.0, +1.0, "X");
h1->SetAxisRange(-1.0, +1.0, "X");
h2->SetAxisRange(-1.0, +1.0, "X");
h3->SetAxisRange(-1.0, +1.0, "X");
h4->SetAxisRange(-1.0, +1.0, "X");
h5->SetAxisRange(-1.0, +1.0, "X");
h7->SetAxisRange(-1.0, +1.0, "X");
}
//Initilize results array
// 7 => Sum in 0 and them MG's in 1-6
......@@ -32,67 +44,114 @@ double devE(const double * parameter) {
//Particle path vector
dir = * (pos[i]) - offset;
//Set final beam energy (added 05July)
double FinalBeamEnergy = BeamTarget.Slow(
// reaction.GetBeamEnergy(),
354.75, //post-CATS beam energy
0.5*parameter[3]*micrometer, //thickness
0); //angle, probably
reaction.SetBeamEnergy(FinalBeamEnergy);
//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, 0.0, 1.0));
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 Al
0.4 * micrometer, //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);
//Fill histograms with Ex
h -> Fill(Ex);
//Fill histograms with
if(allButMG3){
if(detnum[i]!=3){
h -> Fill(Ex);
cout << "removed MG3 from sum!!" << endl;
}
} else {
h -> Fill(Ex);
}
DetectorSwitch(detnum[i], Ex);
hT -> Fill(ThetaTarget/deg,Ex);
}
//Initilise, Draw & Fit histograms
InitiliseCanvas(FitResultMatrix);
//Write vals to screen
if (flagDraw) {
cout << "==================================================" << endl;
}
//Adapt the metric as needed
double multiplier = 0.10; //0.07;
double metric = sqrt(pow(FitResultMatrix[mgSelect][0] - refE, 2)
+ pow(multiplier * FitResultMatrix[mgSelect][2], 2));
if (flagDraw) {cout << "==================================================" << endl;}
/*** Minimize by one peak ***/
/**/
double multiplier = 0.07;
double metric = abs(FitResultMatrix[mgSelect][0]-refE) + abs(multiplier*FitResultMatrix[mgSelect][2]);
/**/
/*** Minimize by all peaks ***/
/**
//double multiplier = 0.125;
double multiplier = 0.005;
double metric =
(1.0/6.0)*abs(FitResultMatrix[1][0]-refE)
+ (1.0/6.0)*abs(FitResultMatrix[2][0]-refE)
+ (1.0/6.0)*abs(FitResultMatrix[3][0]-refE)
+ (1.0/6.0)*abs(FitResultMatrix[4][0]-refE)
+ (1.0/6.0)*abs(FitResultMatrix[5][0]-refE)
+ (1.0/6.0)*abs(FitResultMatrix[6][0]-refE)
+ 1.0*abs(FitResultMatrix[0][0]-refE)
//+ multiplier*FitResultMatrix[1][2]
//+ multiplier*FitResultMatrix[2][2]
//+ multiplier*FitResultMatrix[3][2]
//+ multiplier*FitResultMatrix[4][2]
//+ multiplier*FitResultMatrix[5][2]
//+ multiplier*FitResultMatrix[6][2]
;
**/
/*** Minimize by variation in peaks ***/
/**
vector<double> r = {FitResultMatrix[1][0],
FitResultMatrix[2][0],
FitResultMatrix[3][0],
FitResultMatrix[4][0],
FitResultMatrix[5][0],
FitResultMatrix[6][0]};
double metric = (*max_element(r.begin(),r.end()) - *min_element(r.begin(),r.end()));
**/
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() {
////filename = "XYZE_ShiftMG3_GateOn1838.txt";
//filename = "XYZE_GateOn1838.txt";
//refE = 1.981; // the energy of the selected states
////filename = "XYZE_ShiftMG3_GateOn0143.txt";
filename = "XYZE_GateOn0143.txt";
refE = 0.143; // the energy of the selected states
//Read data in
LoadFile();
......@@ -102,16 +161,22 @@ void MinimizeBeamSpot() {
//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 << "= Type MG# of telescope metric to use, or: =" << endl;
cout << "= 0 = sum of all MG's =" << endl;
cout << "= 8 = sum of all MG's without MG3 =" << endl;
cout << "==================================================" << endl;
cin >> mgSelect;
if (mgSelect == 7) {
mgSelect = 6;
} // Correct the input for MG7
if (mgSelect == 8) {
allButMG3=true;
mgSelect = 0;
}
cout << "==================================================" << endl;
if (mgSelect >= 7) {
if (mgSelect >= 7 ) {
cout << "ERROR!! INVALID SELECTION" << endl;
return;
}
......@@ -125,18 +190,23 @@ void MinimizeBeamSpot() {
auto start = high_resolution_clock::now();
//TESTING: Grid method of avoiding local minima
// for (int x = 0; x < 3; x++) { //7; x++) {
// for (int y = 0; y < 3; y++){ //7; y++) {
for (int z = 0; z < 3; z++){ //7; z++) {
for (int t = 0; t < 9; t++) {
//for (int x = 0; x < 3; x++) { //7; x++) {
//for (int y = 0; y < 3; y++){ //7; y++) {
//for (int z = 0; z < 3; z++){ //7; z++) {
//for (int t = 0; t < 3; 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] = {
0.0, //9.0 - (x * 9.0), //3.0),
0.0, //9.0 - (y * 9.0), //3.0),
9.0 - (z * 9.0), //3.0),
(1.5 - (0.125 * t)) * 4.75
//-3.64, 0.35, -10., 10.
//-3.64, 0.35, +0.11, 3.02
//-3.8748, 0.0228, +1.016354, 1.650928
//-3.502156, 0.391660, 0.986920, 1.147161
//-4.509762, 0.179826, 1.386245, 1.398818
//-4.509762, 0.179826, 1.386245, 1.081334
0.0, 0.0, 0.0, 5.0
//-4.509762, 0.179826, 1.623007, 1.081334
//-3.9164, +0.0550, 1.0558, 1.7685
};
//Don't draw iterations of minimizer
......@@ -151,22 +221,23 @@ void MinimizeBeamSpot() {
//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 -> SetMaxFunctionCalls(100000000); // used by Minuit and Minuit2
minim -> SetMaxIterations(100000000); // used by GSL
minim -> SetPrintLevel(3);
minim -> SetPrecision(1e-10);
minim -> SetPrecision(1e-06);
//Set minimizer function
minim -> SetFunction(func);
//Assign variable limits
minim -> SetLimitedVariable(0, "X", parameter[0], 0.01, -0.5, 0.5);
minim -> SetLimitedVariable(1, "Y", parameter[1], 0.01, -0.5, 0.5);
//minim -> SetLimitedVariable(2, "Z", parameter[2], 0.01, -0.5, 0.5);
//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);
//minim -> SetLimitedVariable(0, "X", parameter[0], 0.10, -5.0, -3.5);
minim -> SetFixedVariable(0, "X", parameter[0]);
//minim -> SetLimitedVariable(1, "Y", parameter[1], 0.10, -0.5, +1.0);
minim -> SetFixedVariable(1, "Y", parameter[1]);
//minim -> SetLimitedVariable(2, "Z", parameter[2], 0.05, -0.0, +2.0);//-1.50, +1.50);
minim -> SetFixedVariable(2, "Z", parameter[2]);
//minim -> SetLimitedVariable(3, "T", parameter[3], 0.05, +0.5, +3.5);
minim -> SetFixedVariable(3, "T", parameter[3]);
//Don't draw iterations of minimizer
flagDraw = 0;
......@@ -195,10 +266,10 @@ void MinimizeBeamSpot() {
cout << "\t\tZ = " << x[2] << " mm" << endl;
cout << "\t\tT = " << x[3] << " um" << endl;
cout << "==================================================" << endl;
}
}
// }
// }
//}
//}
//}
//}
//Stop timer
auto stop = high_resolution_clock::now();
......
......@@ -13,14 +13,16 @@ 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;
vector<int> detnum;
unsigned int mgSelect = 10;
NPL::EnergyLoss CD2("proton_CD2.G4table","G4Table",100);
NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
NPL::EnergyLoss BeamTarget("K47_CD2.G4table","G4Table",100);
bool flagDraw = 0;
bool allButMG3 = 0;
//Output files
TFile *histfile = new TFile("./gridHistograms.root", "RECREATE");
......@@ -28,24 +30,25 @@ 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.);
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.);
string filename;
double refE; // the energy of the selected states
static auto h = new TH1D("h","All MG#'s", 40,-1.0,3.0);
static auto h1 = new TH1D("h1","Individual MG#'s", 40,-1.0,3.0);
static auto h2 = new TH1D("h2","h2", 40,-1.0,3.0);
static auto h3 = new TH1D("h3","h3", 40,-1.0,3.0);
static auto h4 = new TH1D("h4","h4", 40,-1.0,3.0);
static auto h5 = new TH1D("h5","h5", 40,-1.0,3.0);
static auto h7 = new TH1D("h7","h7", 40,-1.0,3.0);
static auto hT = new TH2F("hT","hT", 20,100.,160.,20,-0.1,0.4);
////////////////////////////////////////////////////////////////////////////////
void LoadFile(){
// Open XYZE gamma-gated file
ifstream file("XYZE_Full_09June_MG3.txt");
ifstream file(filename.c_str());
if(!file.is_open()){
cout << "fail to load file" << endl;
exit(1);
cout << "fail to load file " << filename << endl;
}
else {
cout << "Success opening file" << endl;
cout << "Success opening file " << filename << endl;
}
// Read in
......@@ -53,6 +56,10 @@ void LoadFile(){
double x,y,z,e;
while(file >> mg >> x >> y >> z >> e ){
auto p = new TVector3(x,y,z);
// if(mg==3){
// p->SetZ(z+1.0); //Edit MG3 position directly
// cout << "EDITING MG3 POSITION" << endl;
// }
detnum.push_back(mg);
pos.push_back(p);
energy.push_back(e);
......@@ -74,27 +81,26 @@ void FillMatrix(double* matrix, TFitResultPtr fit){
}
}
////////////////////////////////////////////////////////////////////////////////
//[[[[ 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);
Normal.SetXYZ(-0.454552, +0.454996, -0.765742);
break;
case 2:
Normal.SetXYZ(-0.642828, +0.000000, -0.766010);
Normal.SetXYZ(-0.641920, +0.002239, -0.766769);
break;
case 3:
Normal.SetXYZ(-0.454594, -0.450670, -0.768271);
Normal.SetXYZ(-0.455476, -0.452406, -0.766727);
break;
case 4:
Normal.SetXYZ(-0.002437, -0.638751, -0.769409);
Normal.SetXYZ(-0.003212, -0.641831, -0.766839);
break;
case 5:
Normal.SetXYZ(+0.452429, -0.454575, -0.767248);