Commit 5a8739f2 authored by flavigny's avatar flavigny
Browse files

Merge branch 'NPTool.2.dev' of https://gitlab.in2p3.fr/np/nptool into NPTool.2.dev

parents 6259db4a 7dcb3734
Pipeline #123247 passed with stages
in 3 minutes and 28 seconds
#include "NPXmlParser.h"
#include <stdio.h>
#include <iostream>
using namespace NPL;
using namespace NPL::XML;
////////////////////////////////////////////////////////////////////////////////
......@@ -113,6 +112,7 @@ void XmlParser::LoadNode(TXMLEngine* xml, XMLNodePointer_t node, Int_t level){
param=xml->GetNext(param);
}
std::string name = xml->GetNodeName(child);
b.SetName(name);
m_blocks[name].push_back(b);
child = xml->GetNext(child);
}
......
......@@ -4,6 +4,7 @@
#include<map>
#include<set>
#include<vector>
#include <iostream>
#include "TXMLEngine.h"
namespace NPL{
/////////////////////
......@@ -44,20 +45,25 @@ namespace NPL{
public:
Channel();
~Channel();
Channel(int device,int geo,int ch){
Channel(int device, int fpl, int detector, int geo,int ch){
m_device=device;
m_fpl=fpl;
m_detector=detector;
m_geo=geo;
m_ch=ch;
};
private:
int m_device;
int m_fpl;
int m_detector;
int m_geo;
int m_ch;
public:
int norme() const {return (m_device*1000000+m_geo*1000+m_ch);} ;
int norme() const {return (m_device*10000000000000+m_fpl*10000000000+m_detector*1000000+m_geo*1000+m_ch);} ;
void EraseGeoCh(){m_ch=-1;m_geo=-1;};// use for MRDC
void Print() const {std::cout << m_device << " " << m_fpl << " " << m_detector << " " << m_geo << " " << m_ch<< std::endl;}
public:
bool operator<(const Channel p2){
return this->norme()<p2.norme();
......@@ -82,7 +88,8 @@ namespace NPL{
double AsDouble(std::string name);
std::string AsString(std::string name);
void AddParameter(parameter p){ m_parameters.insert(p);};
std::string GetName(){return m_name;};
void SetName(std::string name) {m_name=name;};
private:
std::string m_name;
std::set<parameter> m_parameters;
......@@ -112,7 +119,7 @@ namespace NPL{
void SetBlock(const XML::Channel& c,XML::block* b){
m_Channels[c]=b;
}
std::map<XML::Channel,XML::block*> GetChannels(){return m_Channels;};
public:
std::vector<XML::block*> GetAllBlocksWithName(std::string name) ;
std::vector<std::string> GetAllBlocksName();
......
This diff is collapsed.
This diff is collapsed.
......@@ -7,4 +7,4 @@ db/SAMURAIFDC0_20200109.xml
db/SAMURAIHOD_s034_all40mV_s037_20170702.xml
db/SAMURAIBDC1.xml
db/SAMURAIBDC2.xml
db/SAMURAIPlastic.s034.2.xml
db/SAMURAIPlastic.s034.0.xml
......@@ -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,51 +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
}
DetectorSwitch(detnum[i], Ex);
}
//End loop over events
//Initilise, Draw & Fit histograms
InitiliseCanvas(FitResultMatrix);
//Write vals to screen
cout << "Mean: " << h->GetMean()
<< "\t StdDev: " << h->GetStdDev()
<< "\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)
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;
*/
//Adapt the metric as needed
return sqrt( pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
return sqrt( pow(FitResultMatrix[mgSelect][0]-refE,2) + pow(0.1*FitResultMatrix[mgSelect][2],2) );
}
////////////////////////////////////////////////////////////////////////////////
void MinimizeBeamSpot(){
......@@ -126,6 +82,19 @@ 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
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);
......@@ -133,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);
......@@ -150,20 +118,31 @@ 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();
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;
devE(x);
cout << "========================================" << endl;
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;
}
......@@ -9,12 +9,21 @@ 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);
using namespace std;
bool flagDraw = 0;
static auto h = new TH1D("h","All MG#'s", 60,-1.,1.);
//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.);
static auto h3 = new TH1D("h3","h3", 40,-1.,1.);
......@@ -46,7 +55,104 @@ void LoadFile(){
file.close();
}
////////////////////////////////////////////////////////////////////////////////
void InitiliseCanvas(){
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
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 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);
......@@ -58,61 +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();
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
//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");
......@@ -122,15 +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();
h->Fit("gaus", "WQ");
//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