Skip to content

GitLab

Commit 1fc0a534 authored by Adrien Matta's avatar Adrien Matta ☠
Browse files

* Progress on S034 FDC alignement

parent 94dc0140
Pipeline #121931 passed with stages
in 7 minutes and 1 second
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1424 additions and 728 deletions
NPLib/Detectors/Minos/TMinosPhysics.cxx
View file @ 1fc0a534
......@@ -88,11 +88,12 @@ void TMinosPhysics::BuildPhysicalEvent() {
void TMinosPhysics::PreTreat() {
// apply thresholds and calibration
static unsigned int sizePad,sizeQ ;
sizePad = m_EventData->GetPadMult();
static unsigned short PadNumber;
static double Q,T;
static auto cal= CalibrationManager::getInstance();
static string cal_v, cal_o;
sizePad = m_EventData->GetPadMult();
if(sizePad>20){
for(unsigned int i = 0 ; i < sizePad ; i++){
vector<unsigned short>* Charge = m_EventData->GetChargePtr(i);
......
NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
View file @ 1fc0a534
......@@ -50,10 +50,10 @@ double SamuraiFieldMap::Delta(const double* parameter){
// Move the fdc2 pos from lab frame to fdc2 frame
pos.back().RotateY(-m_fdc2angle+m_angle);
double d = (pos.back().X()-m_FitPosFDC2.X())*(pos.back().X()-m_FitPosFDC2.X());
//diff = pos.back()-m_FitPosFDC2;
return d;
//return diff.Mag2();
//double d = (pos.back().X()-m_FitPosFDC2.X())*(pos.back().X()-m_FitPosFDC2.X());
// return d;
diff = pos.back()-m_FitPosFDC2;
return diff.Mag2();
}
////////////////////////////////////////////////////////////////////////////////
......
NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
View file @ 1fc0a534
......@@ -53,6 +53,19 @@ ClassImp(TSamuraiBDCPhysics)
PowerThreshold=5;
}
///////////////////////////////////////////////////////////////////////////
TVector3 TSamuraiBDCPhysics::GetPos(unsigned int det){
TVector3 res(-10000,-10000,-10000);
unsigned int size = PosX.size();
for(unsigned int i = 0 ; i < size ; i++){
if(Detector[i]==det){
res = TVector3(PosX[i],PosY[i],PosZ[i]);
}
}
return res;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::BuildSimplePhysicalEvent(){
BuildPhysicalEvent();
......@@ -325,6 +338,10 @@ void TSamuraiBDCPhysics::ReadConfiguration(NPL::InputParser parser){
m_invertY[det] = invertY;
m_invertD[det] = invertD;
}
else{
cout << " --- ERROR : BDC block wrongly formatted" << endl;
exit(1);
}
}
#if __cplusplus > 199711L && NPMULTITHREADING
......
NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
View file @ 1fc0a534
......@@ -97,6 +97,13 @@ class TSamuraiBDCPhysics : public TObject, public NPL::VDetector{
std::vector<TVector3> Dir;
std::vector<int> PileUp;
private: // offset and inversion
std::map<unsigned int, TVector3> m_offset;//!
std::map<unsigned int, bool> m_invertX;//!
std::map<unsigned int, bool> m_invertY;//!
std::map<unsigned int, bool> m_invertD;//!
public:
// Projected position at given Z plan
TVector3 ProjectedPosition(int Detector, double Z);
......@@ -182,25 +189,20 @@ class TSamuraiBDCPhysics : public TObject, public NPL::VDetector{
TSamuraiBDCData* GetRawData() const {return m_EventData;}
TSamuraiBDCData* GetPreTreatedData() const {return m_PreTreatedData;}
double GetPosX(unsigned int det) {return PosX[det];}
double GetPosY(unsigned int det) {return PosY[det];}
double GetThetaX(unsigned int det){return ThetaX[det];}
double GetPhiY(unsigned int det) {return PhiY[det];}
double GetDevX(unsigned int det) {return devX[det];}
double GetDevY(unsigned int det) {return devY[det];}
int GetPileUp(unsigned int det){return PileUp[det];}
TVector3 GetPos(unsigned int det);
double GetPosX(unsigned int i) {return PosX[i];}
double GetPosY(unsigned int i) {return PosY[i];}
double GetThetaX(unsigned int i){return ThetaX[i];}
double GetPhiY(unsigned int i) {return PhiY[i];}
double GetDevX(unsigned int i) {return devX[i];}
double GetDevY(unsigned int i) {return devY[i];}
int GetPileUp(unsigned int i){return PileUp[i];}
private: // Root Input and Output tree classes
TSamuraiBDCData* m_EventData;//!
TSamuraiBDCData* m_PreTreatedData;//!
TSamuraiBDCPhysics* m_EventPhysics;//!
private: // offset and inversion
std::map<unsigned int, TVector3> m_offset;//!
std::map<unsigned int, bool> m_invertX;//!
std::map<unsigned int, bool> m_invertY;//!
std::map<unsigned int, bool> m_invertD;//!
private: // Spectra Class
// TSamuraiBDCSpectra* m_Spectra; // !
......
Projects/S034/Analysis.cxx
View file @ 1fc0a534
......@@ -37,18 +37,19 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
Minos= (TMinosPhysics*) m_DetectorManager->GetDetector("Minos");
//Nebula = (TNebulaPhysics*) m_DetectorManager->GetDetector("NEBULA");
BDC = (TSamuraiBDCPhysics*) m_DetectorManager->GetDetector("SAMURAIBDC");
FDC0 = (TSamuraiFDC0Physics*) m_DetectorManager->GetDetector("SAMURAIFDC0");
FDC2 = (TSamuraiFDC2Physics*) m_DetectorManager->GetDetector("SAMURAIFDC2");
Hodo = (TSamuraiHodoscopePhysics*) m_DetectorManager->GetDetector("SAMURAIHOD");
m_field.LoadMap(30*deg,"field_map/180702-2,40T-3000.table.bin",10);
m_field.SetFDC2Angle((59.930-90.0)*deg);
m_field.SetFDC2R(FDC2->GetOffset().Z());
InitOutputBranch();
InitInputBranch();
Minos= (TMinosPhysics*) m_DetectorManager->GetDetector("Minos");
//Nebula = (TNebulaPhysics*) m_DetectorManager->GetDetector("NEBULA");
BDC = (TSamuraiBDCPhysics*) m_DetectorManager->GetDetector("SAMURAIBDC");
FDC0 = (TSamuraiFDC0Physics*) m_DetectorManager->GetDetector("SAMURAIFDC0");
FDC2 = (TSamuraiFDC2Physics*) m_DetectorManager->GetDetector("SAMURAIFDC2");
Hodo = (TSamuraiHodoscopePhysics*) m_DetectorManager->GetDetector("SAMURAIHOD");
m_field.LoadMap(30*deg,"field_map/180702-2,40T-3000.table.bin",10);
m_field.SetFDC2Angle((59.930-90.0)*deg);
m_field.SetFDC2R(FDC2->GetOffset().Z());
InitOutputBranch();
InitInputBranch();
// for fdc/bdc alignement
file.open("Calibration/Pos/fdc.txt");
}
////////////////////////////////////////////////////////////////////////////////
......@@ -56,28 +57,60 @@ void Analysis::TreatEvent(){
Clear();
//cout << Trigger << " " ;
Trigger=Trigger&0x00ff;
//cout << Trigger << endl;
//cout << Trigger << endl;
// Compute Brho
//
if( FDC2->PosX>-1500 && FDC2->PosX<1000
&& FDC2->PosY>-500 && FDC2->PosY<500
&& FDC0->PosX>-80 && FDC0->PosX<80
&& FDC0->PosY>-80 && FDC0->PosY<80
&& FDC0->Dir.Z()>0.8) { // if both are correctly build
// Compute ThetaX and PhiY using Minos vertex and FDC0 XY
double FDC0_ThetaX = FDC0->ThetaX;
double FDC0_PhiY = FDC0->PhiY;
if(Minos->Z_Vertex>0){
FDC0_ThetaX = atan((FDC0->PosX-Minos->X_Vertex)/(1283.7-Minos->Z_Vertex));
FDC0_PhiY = atan((FDC0->PosY-Minos->Y_Vertex)/(1283.7-Minos->Z_Vertex));
}
//double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX/*-252.55*/, FDC2->ThetaX};
&& FDC0->PosY>-80 && FDC0->PosY<80 // both FDC ok
&& (Minos->Tracks_P0.size()==1)|| (Minos->Tracks_P0.size()==2)) { // p,pn or p,2p
// Compute ThetaX and PhiY using Minos vertex and FDC0 X
double FDC0_ThetaX = FDC0->ThetaX;
double FDC0_PhiY = FDC0->PhiY;
TVector3 FDC0_Dir = FDC0->Dir;
// Check if both BDC are reconstructed
//
TVector3 BDC1=BDC->GetPos(1);
TVector3 BDC2=BDC->GetPos(2);
//double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX};
if( BDC1.Z()!=-10000 && BDC2.Z()!=-10000){
TVector3 Vertex;
TVector3 P1 = Minos->Tracks_P0[0]+Minos->Tracks_Dir[0];
TVector3 delta;
MinimumDistanceTwoLines(BDC1,BDC2,
Minos->Tracks_P0[0], P1,
Vertex, delta) ;
FDC0_Dir= FDC0->GetPos()-Vertex;
FDC0_Dir=FDC0_Dir.Unit();
TVector3 BDCDir=BDC2-BDC1;
BDCDir=BDCDir.Unit();
BDCDir*=(Vertex.Z()-BDC2.Z())/BDCDir.Z();
//cout << "BDC" << endl;
BDCX=(BDC2+BDCDir).X();
BDCY=(BDC2+BDCDir).Y();
// Z relative to Minos entrance
Z=Vertex.Z()+4657.39;
//double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX};
// Brho=r_fit(brho_param);
Brho=m_field.FindBrho(FDC0->GetPos(),FDC0->Dir,FDC2->GetPos(),TVector3(0,0,1));
if(FDC0->GetPos().X()>-10000 && FDC0_Dir.Z()>0.8){
FDC0_ThetaX = atan((FDC0->PosX-Vertex.X())/(1283.7-Z));
FDC0_PhiY = atan((FDC0->PosY-Vertex.Y())/(1283.7-Z));
double brho_param[6]={FDC0->PosX-1.77, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX+252.55, FDC2->ThetaX};
BrhoP=r_fit(brho_param);
Brho=m_field.FindBrho(FDC0->GetPos(),FDC0_Dir,FDC2->GetPos(),TVector3(0,0,1));
// cout << Brho-BrhoP << endl;
}
/* // Calib//////////////////////////////////////////////////////////////////
static int count=0;
if(FDC2->PosX-252.55>0&&FDC0->GetPos().X()>-10000 && FDC0->Dir.Z()>0.8){
file << FDC0->GetPos().X() <<" " << FDC0->GetPos().Y() << " " << FDC0->GetPos().Z() <<" " << FDC0_Dir.X() <<" " << FDC0_Dir.Y() << " " << FDC0_Dir.Z()<< " " ;
file << FDC2->GetPos().X() <<" " << FDC2->GetPos().Y() << " " << FDC2->GetPos().Z() <<" " << FDC2->Dir.X() <<" " << FDC2->Dir.Y() << " " << FDC2->Dir.Z()<< endl;
count ++;
}
if(count>10000)
exit(1);*/
}
}
}
......@@ -88,18 +121,26 @@ void Analysis::End(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Clear(){
Brho=-1000;
BDCX=-1000;
BDCY=-1000;
BrhoP=-1000;
Beta_f=-1000;
Z=-1000;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Brho",&Brho,"Brho/D");
RootOutput::getInstance()->GetTree()->Branch("BrhoP",&BrhoP,"BrhoP/D");
RootOutput::getInstance()->GetTree()->Branch("BDCX",&BDCX,"BDCX/D");
RootOutput::getInstance()->GetTree()->Branch("BDCY",&BDCY,"BDCY/D");
RootOutput::getInstance()->GetTree()->Branch("Z",&Z,"Z/D");
RootOutput::getInstance()->GetTree()->Branch("Beta_f",&Beta_f,"Beta_f/D");
RootOutput::getInstance()->GetTree()->Branch("Trigger",&Trigger,"Trigger/I");
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchAddress("Trigger",&Trigger);
RootInput::getInstance()->GetChain()->SetBranchAddress("Trigger",&Trigger);
}
// -*- mode: c++ -*-
......@@ -213,7 +254,7 @@ static double e_gCoefficient[] = {
-0.188488,
-0.675742,
-0.181284
};
};
// Assignment to error coefficients vector.
static double e_gCoefficientRMS[] = {
......@@ -293,7 +334,7 @@ static double e_gCoefficientRMS[] = {
0.547759,
2.82987,
2.05714
};
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
......@@ -393,16 +434,16 @@ double e_fit(double *x) {
double v = 1 + 2. / (e_gXMax[j] - e_gXMin[j]) * (x[j] - e_gXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
......@@ -525,7 +566,7 @@ static double p_gCoefficient[] = {
-0.145767,
-1.13422,
0.865558
};
};
// Assignment to error coefficients vector.
static double p_gCoefficientRMS[] = {
......@@ -605,7 +646,7 @@ static double p_gCoefficientRMS[] = {
3.21076e-11,
8.11735e-11,
7.22006e-11
};
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
......@@ -705,16 +746,16 @@ double p_fit(double *x) {
double v = 1 + 2. / (p_gXMax[j] - p_gXMin[j]) * (x[j] - p_gXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
......@@ -837,7 +878,7 @@ static double r_gCoefficient[] = {
-0.0019449,
-0.0151334,
0.0115488
};
};
// Assignment to error coefficients vector.
static double r_gCoefficientRMS[] = {
......@@ -917,7 +958,7 @@ static double r_gCoefficientRMS[] = {
3.21076e-11,
8.11735e-11,
7.22006e-11
};
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
......@@ -1017,16 +1058,16 @@ double r_fit(double *x) {
double v = 1 + 2. / (r_gXMax[j] - r_gXMin[j]) * (x[j] - r_gXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
......@@ -1150,7 +1191,7 @@ static double l_gCoefficient[] = {
-25.8219,
2.18449,
-0.917242
};
};
// Assignment to error coefficients vector.
static double l_gCoefficientRMS[] = {
......@@ -1231,7 +1272,7 @@ static double l_gCoefficientRMS[] = {
31.7885,
20.2812,
22.29
};
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
......@@ -1332,16 +1373,16 @@ double l_fit(double *x) {
double v = 1 + 2. / (l_gXMax[j] - l_gXMin[j]) * (x[j] - l_gXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
......@@ -1465,7 +1506,7 @@ static double t_gCoefficient[] = {
-0.0377476,
-0.150836,
0.0571845
};
};
// Assignment to error coefficients vector.
static double t_gCoefficientRMS[] = {
......@@ -1546,7 +1587,7 @@ static double t_gCoefficientRMS[] = {
0.930273,
2.8245,
2.14083
};
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
......@@ -1647,16 +1688,16 @@ double t_fit(double *x) {
double v = 1 + 2. / (t_gXMax[j] - t_gXMin[j]) * (x[j] - t_gXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p3 = 2 * v * p2 - p1;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
......@@ -1668,10 +1709,8 @@ double t_fit(double *x) {
}
// EOF for simtree_tof.C
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
// Construct Method to be pass to the AnalysisFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
return (NPL::VAnalysis*) new Analysis();
......@@ -1681,13 +1720,13 @@ NPL::VAnalysis* Analysis::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy{
public:
proxy(){
NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
}
};
class proxy_analysis{
public:
proxy_analysis(){
NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
}
};
proxy p;
proxy_analysis p_analysis;
}
Projects/S034/Analysis.h
View file @ 1fc0a534
......@@ -29,7 +29,7 @@
#include"TSamuraiFDC2Physics.h"
#include"TSamuraiHodoscopePhysics.h"
#include"SamuraiFieldMap.h"
#include<fstream>
class Analysis: public NPL::VAnalysis{
public:
Analysis();
......@@ -50,14 +50,16 @@ class Analysis: public NPL::VAnalysis{
TSamuraiFDC2Physics* FDC2;
TSamuraiHodoscopePhysics* Hodo;
SamuraiFieldMap m_field ;
ofstream file;
private: // output variable
double Brho;
double Brho,BrhoP,BDCX,BDCY,Z;
double Beta_f;
int Trigger;
public:
void Clear();
void InitOutputBranch();
void InitInputBranch();
};
// use for broh calculation
......
Projects/S034/Calibration/Pos/FDC.cxx 0 → 100644
View file @ 1fc0a534
vector<TVector3*> pos0;
vector<TVector3*> dir0;
vector<TVector3*> pos2;
vector<TVector3*> dir2;
SamuraiFieldMap field;
////////////////////////////////////////////////////////////////////////////////
double devB(const double* parameter){
// Return the standard deviation in Brho
unsigned int size = pos0.size();
TVector3 o0(parameter[0],parameter[1],parameter[2]);
TVector3 o2(parameter[3],parameter[4],parameter[5]);
TVector3 p0,p2;
field.SetFDC2R(parameter[5]);
double Brho;
static auto h = new TH1D("h","h", 1000,0,10);
h->Reset();
for(unsigned int i = 0 ; i < 1000 ; i++){
p0=*(pos0[i])+o0;
p2=*(pos2[i])+o2;
if(dir0[i]->Z()>0.8)
h->Fill(field.FindBrho(p0,*dir0[i],p2,*dir2[i]));
}
o0.Print();
o2.Print();
cout << h->GetStdDev() << endl;
return h->GetStdDev();
}
////////////////////////////////////////////////////////////////////////////////
void LoadFile(){
ifstream file("Calibration/Pos/fdc.txt");
if(!file.is_open()){
cout << "fail to load file" << endl;
exit(1);
}
else {
cout << "Success opening file" << endl;
}
double x0,y0,z0;
double dx0,dy0,dz0;
double x2,y2,z2;
double dx2,dy2,dz2;
while(file >> x0 >> y0 >> z0 >> dx0 >> dy0 >> dz0 >> x2 >> y2 >> z2 >> dx2 >> dy2 >> dz2){
auto p0 = new TVector3(x0,y0,z0);
auto d0 = new TVector3(dx0,dy0,dz0);
auto p2 = new TVector3(x2,y2,z2);
auto d2 = new TVector3(dx2,dy2,dz2);
pos0.push_back(p0);
dir0.push_back(d0);
pos2.push_back(p2);
dir2.push_back(d2);
}
cout << " Val " << pos0.size() << " Value Loaded" << endl;
file.close();
}
////////////////////////////////////////////////////////////////////////////////
void FDC(){
// Data
LoadFile();
// Field map
field.LoadMap(30*deg,"field_map/180702-2,40T-3000.table.bin",10);
field.SetFDC2Angle((59.930-90.0)*deg);
double parameter[6]={0,0,-3456.52,-252.55,0,4123.47};
cout << devB(parameter) << endl;
// Minimizer
auto min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
auto func=ROOT::Math::Functor(&devB,6);
min->SetFunction(func);
min->SetPrintLevel(0);