Skip to content

GitLab

Commit f7d038b2 authored by Morfouace's avatar Morfouace
Browse files

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

parents 9f394fa2 18d09002
Pipeline #122428 passed with stages
in 3 minutes and 27 seconds
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1854 additions and 877 deletions
NPLib/Detectors/Minos/TMinosPhysics.cxx
View file @ f7d038b2
......@@ -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 @ f7d038b2
......@@ -22,6 +22,7 @@
#include "SamuraiFieldMap.h"
#include "NPPhysicalConstants.h"
#include "Math/Factory.h"
using namespace NPUNITS;
#include <cmath>
......@@ -31,33 +32,48 @@ using namespace std;
ClassImp(SamuraiFieldMap);
SamuraiFieldMap::SamuraiFieldMap(){
m_BrhoScan=NULL;
m_min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
m_func=ROOT::Math::Functor(this,&SamuraiFieldMap::Delta,1);
m_min->SetFunction(m_func);
m_min->SetPrintLevel(0);
}
////////////////////////////////////////////////////////////////////////////////
double SamuraiFieldMap::FindBrho(TVector3& p_fdc0,TVector3& d_fdc0,TVector3& p_fdc2,TVector3& d_fdc2){
if(!m_BrhoScan)
BrhoScan(2.5,10,0.1);
double SamuraiFieldMap::Delta(const double* parameter){
static vector<TVector3>pos ;
static TVector3 diff;
pos =Propagate(parameter[0],m_FitPosFDC0,m_FitDirFDC0,false);
// 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());
// return d;
diff = pos.back()-m_FitPosFDC2;
return diff.Mag2();
}
////////////////////////////////////////////////////////////////////////////////
double SamuraiFieldMap::FindBrho(TVector3 p_fdc0,TVector3 d_fdc0,TVector3 p_fdc2,TVector3 d_fdc2){
m_FitPosFDC0=p_fdc0;
m_FitDirFDC0=d_fdc0;
m_FitPosFDC2=p_fdc2;
m_FitDirFDC2=d_fdc2;
if(!m_BrhoScan)
BrhoScan(1,10,0.1);
// do a first guess based on fdc2 pos
double b = m_BrhoScan->Eval(p_fdc2.X());
double b0 = b;
vector<TVector3> pos=Propagate(3000,b,p_fdc0,d_fdc0);
double step = 1;
double d = (pos.back()-p_fdc2).Mag();
double dd=d;
short sign =1;
unsigned int limit =1000;
unsigned count=0;
while(step>1e-6 && count<limit){
dd=d;
b+=sign*step;
pos=Propagate(3000,b,p_fdc0,d_fdc0);
d = (pos.back()-p_fdc2).Mag();
if(d>=dd){
step/=10;
sign=-sign;
}
count++;
}
return b-sign*0.5*step;
double b0[1] ={m_BrhoScan->Eval(p_fdc2.X())};
m_min->Clear();
m_min->SetPrecision(1e-6);
m_min->SetMaxFunctionCalls(1000);
m_min->SetLimitedVariable(0,"B",b0[0],0.1,1,10);
m_min->Minimize();
return m_min->X()[0];
}
////////////////////////////////////////////////////////////////////////////////
......@@ -68,11 +84,16 @@ TGraph* SamuraiFieldMap::BrhoScan(double min, double max,double step){
unsigned int size = (max-min)/step;
m_BrhoScan->Set(size);
unsigned int i=0;
TVector3 p(0,0,-3500);
TVector3 d(0,0,1);
p.RotateY(m_angle);
d.RotateY(m_angle);
for(double b = min ; b < max ; b+=step){
vector<TVector3> pos= Propagate(3000,b,TVector3(0,0,0),TVector3(0,0,1));
pos.back().RotateY(-m_fdc2angle);
m_BrhoScan->SetPoint(i++,pos.back()[0],b);
vector<TVector3> pos= Propagate(b,p,d,false);
pos.back().RotateY(-m_fdc2angle);
m_BrhoScan->SetPoint(i++,pos.back().X(),b);
}
m_BrhoScan->Sort();
return m_BrhoScan;
}
......@@ -91,7 +112,7 @@ TVector3 SamuraiFieldMap::PropagateToFDC2(TVector3 pos, TVector3 dir){
}
////////////////////////////////////////////////////////////////////////////////
std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVector3 pos, TVector3 dir){
std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TVector3 dir,bool store){
pos.RotateY(m_angle);
dir.RotateY(m_angle);
dir=dir.Unit();
......@@ -102,57 +123,66 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
N.SetBrho(Brho);
// track result
std::vector< TVector3 > track;
// starting point of the track
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
static std::vector< TVector3 > track;
track.clear();
// starting point of the track
if(store){
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
}
dir=dir.Unit();
double r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z());
static double r;
r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z());
// number of step taken
unsigned int count = 0;
static unsigned int count,limit;
count = 0;
// maximum number of state before giving up
unsigned int limit = 1000;
limit = 1000;
// First propagate to r_max with one line
while(r>rmax && count<limit){
pos+=(r-rmax)/cos(dir.Theta())*dir.Unit();
while(r>m_Rmax && count<limit){
pos+=(r-m_Rmax)/cos(dir.Theta())*dir.Unit();
r= 1.01*sqrt(pos.X()*pos.X()+pos.Z()*pos.Z());
}
if(r<=rmax){ // success
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
if(r<=m_Rmax){ // success
if(store){
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
}
}
else {// failure
//cout << "Fail" << endl;
return track;
}
TVector3 xk1,xk2,xk3,xk4; // position
TVector3 pk1,pk2,pk3,pk4; // impulsion
double K = N.GetEnergy(); // kinetic energy
double m = N.Mass(); // mc2
double P = sqrt(K*K+2*K*m)/c_light; // P
double px = P*dir.X();//px
double py = P*dir.Y();//py
double pz = P*dir.Z();//pz
TVector3 imp = P*dir;
double h = 0.1*nanosecond; // typical step length ~1mm at beta 0.6
while(r<=rmax && count < limit){
static TVector3 xk1,xk2,xk3,xk4; // position
static TVector3 pk1,pk2,pk3,pk4; // impulsion
static TVector3 imp;
static double K,m,P,px,py,pz;
K = N.GetEnergy(); // kinetic energy
m = N.Mass(); // mc2
P = sqrt(K*K+2*K*m)/c_light; // P
px = P*dir.X();//px
py = P*dir.Y();//py
pz = P*dir.Z();//pz
imp = P*dir;
static double h = 1*nanosecond;
while(r<=m_Rmax && count < limit){
func(N, pos , imp , xk1, pk1);
func(N, pos+xk1*(h/2.), imp+pk1*(h/2.) , xk2, pk2);
func(N, pos+xk2*(h/2.), imp+pk2*(h/2.) , xk3, pk3);
func(N, pos+xk3*h , imp+pk3*h , xk4, pk4);
pos +=(xk1+2*xk2+2*xk3+xk4)*(h/6.);
imp +=(pk1+2*pk2+2*pk3+pk4)*(h/6.);
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
if(store){
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
}
r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z());
count++;
}
......@@ -160,33 +190,33 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
pos = PropagateToFDC2(pos, imp);
pos.RotateY(-m_angle);
track.push_back(pos);
return track;
}
////////////////////////////////////////////////////////////////////////////////
void SamuraiFieldMap::func(NPL::Particle& N, TVector3 pos, TVector3 imp, TVector3& new_pos, TVector3& new_imp){
double px,py,pz;
static double px,py,pz,vx,vy,vz,Bx,By,Bz,q,P2,D,m2c4;
static vector<double> B;
px=imp.X();
py=imp.Y();
pz=imp.Z();
double P2,D,m2c4;
P2=imp.Mag2(); // P2
m2c4 = N.Mass()*N.Mass();
D=sqrt(m2c4+P2*c_squared); // sqrt(m2c4+P2c2)
double vx=px*c_squared/D;// pxc * c / D = pxc2/D
double vy=py*c_squared/D;
double vz=pz*c_squared/D;
vx=px*c_squared/D;// pxc * c / D = pxc2/D
vy=py*c_squared/D;
vz=pz*c_squared/D;
new_pos.SetX(vx);
new_pos.SetY(vy);
new_pos.SetZ(vz);
vector<float> B = InterpolateB(pos);
double Bx= B[0];
double By= B[1];
double Bz= B[2];
double q = N.GetZ()*eplus; // issue with the tesla/coulomb definition
B = InterpolateB(pos);
Bx= B[0];
By= B[1];
Bz= B[2];
q = N.GetZ()*eplus; // issue with the tesla/coulomb definition
new_imp.SetX(q*(vy*Bz-vz*By));// q*pyc2*Bz/D -q*pzc2*By/D
new_imp.SetY(q*(vz*Bx-vx*Bz));
new_imp.SetZ(q*(vx*By-vy*Bx));
......@@ -202,11 +232,11 @@ void SamuraiFieldMap::LoadMap(double angle,std::string file,unsigned int bin){
LoadAscii(file);
}
////////////////////////////////////////////////////////////////////////////////
std::vector<float> SamuraiFieldMap::GetB(std::vector<float>& pos){
static vector<float> nullv ={0,0,0};
std::vector<double> SamuraiFieldMap::GetB(std::vector<double>& pos){
static vector<double> nullv ={0,0,0};
// the map is only 1/4 of the detector so we apply symetrie:
double x,y,z ;
if(pos[0]<0)
pos[0] = -pos[0];
......@@ -222,11 +252,11 @@ std::vector<float> SamuraiFieldMap::GetB(std::vector<float>& pos){
}
////////////////////////////////////////////////////////////////////////////////
std::vector<float> SamuraiFieldMap::InterpolateB(const std::vector<float>& pos){
static vector<float> nullv ={0,0,0};
std::vector<double> SamuraiFieldMap::InterpolateB(const std::vector<double>& pos){
static vector<double> nullv ={0,0,0};
// the map is only 1/4 of the detector so we apply symetrie:
double x,y,z ;
if(pos[0]>0)
x = pos[0];
else
......@@ -239,37 +269,37 @@ std::vector<float> SamuraiFieldMap::InterpolateB(const std::vector<float>& pos){
else
z = -pos[2];
// out of bound
// out of bound
if(x<m_x_min || x>m_x_max)
return nullv;
if(y<m_y_min || y>m_y_max)
return nullv;
if(z<m_z_min || z>m_z_max)
return nullv;
float xm = (float)((int)x/m_bin*m_bin);
float ym = (float)((int)y/m_bin*m_bin);
float zm = (float)((int)z/m_bin*m_bin);
vector<float> p0={xm,ym,zm};
vector<float> p1={xm+m_bin,ym,zm};
vector<float> p2={xm,ym+m_bin,zm};
vector<float> p3={xm,ym,zm+m_bin};
vector<float> p4={xm-m_bin,ym,zm};
vector<float> p5={xm,ym-m_bin,zm};
vector<float> p6={xm,ym,zm-m_bin};
vector<map<vector<float>,vector<float>>::iterator> it=
double xm = (double)((int)x/m_bin*m_bin);
double ym = (double)((int)y/m_bin*m_bin);
double zm = (double)((int)z/m_bin*m_bin);
vector<double> p0={xm,ym,zm};
vector<double> p1={xm+m_bin,ym,zm};
vector<double> p2={xm,ym+m_bin,zm};
vector<double> p3={xm,ym,zm+m_bin};
vector<double> p4={xm-m_bin,ym,zm};
vector<double> p5={xm,ym-m_bin,zm};
vector<double> p6={xm,ym,zm-m_bin};
vector<map<vector<double>,vector<double>>::iterator> it=
{ m_field.lower_bound(p0),
m_field.lower_bound(p1),m_field.lower_bound(p2),m_field.lower_bound(p3),
m_field.lower_bound(p4),m_field.lower_bound(p5),m_field.lower_bound(p6)};
float Bx=0;
float By=0;
float Bz=0;
float totalW=0;
double Bx=0;
double By=0;
double Bz=0;
double totalW=0;
auto end=m_field.end();
unsigned int size = it.size();
for(unsigned int i = 0 ; i < size; i++){
......@@ -284,7 +314,7 @@ std::vector<float> SamuraiFieldMap::InterpolateB(const std::vector<float>& pos){
totalW+=1./(d*d);
}
}
vector<float> res = {Bx/totalW,By/totalW,Bz/totalW};
vector<double> res = {Bx/totalW,By/totalW,Bz/totalW};
return res;
}
......@@ -297,7 +327,7 @@ void SamuraiFieldMap::LoadAscii(std::string file){
}
cout << "//////// Loading Ascii Samurai field map " << file << endl;
float x,y,z,Bx,By,Bz;
double x,y,z,Bx,By,Bz;
m_x_max=m_y_max=m_z_max=-1e32;
m_x_min=m_y_min=m_z_min=1e32;
......@@ -312,11 +342,11 @@ void SamuraiFieldMap::LoadAscii(std::string file){
while(in >> x >> y >> z >> Bx >> By >> Bz){
if(++count%50000==0)
cout << "\r - Loading " << count << " values " << flush;
vector<float> p = {x,y,z};
vector<double> p = {x,y,z};
Bx*=tesla;
By*=tesla;
Bz*=tesla;
vector<float> B = {Bx,By,Bz};
vector<double> B = {Bx,By,Bz};
m_field[p]=B;
if(x<m_x_min)
m_x_min=x;
......@@ -332,8 +362,10 @@ void SamuraiFieldMap::LoadAscii(std::string file){
m_z_max=z;
}
m_Rmax=m_x_max;
cout << "\r - " << count << " values loaded" << endl;
cout << " - min(" << m_x_min <<";"<< m_y_min <<";" << m_z_min<< ") max(" << m_x_max <<";"<< m_y_max <<";" << m_z_max<< ")" << endl;
cout << " - Rmax = " << m_Rmax << endl;
in.close();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -345,7 +377,7 @@ void SamuraiFieldMap::LoadBinary(std::string file){
}
cout << "//////// Loading Binary Samurai field map " << file << endl;
float x,y,z,Bx,By,Bz;
double x,y,z,Bx,By,Bz;
m_x_max=m_y_max=m_z_max=-1e32;
m_x_min=m_y_min=m_z_min=1e32;
......@@ -361,12 +393,11 @@ void SamuraiFieldMap::LoadBinary(std::string file){
in.read((char*)&Bx,sizeof(Bx));
in.read((char*)&By,sizeof(By));
in.read((char*)&Bz,sizeof(Bz));
vector<float> p = {x,y,z};
vector<double> p = {x,y,z};
Bx*=tesla;
By*=tesla;
Bz*=tesla;
vector<float> B = {Bx,By,Bz};
vector<double> B = {Bx,By,Bz};
m_field[p]=B;
if(x<m_x_min)
m_x_min=x;
......@@ -381,7 +412,10 @@ void SamuraiFieldMap::LoadBinary(std::string file){
if(z>m_z_max)
m_z_max=z;
}
m_Rmax=m_x_max;
cout << "\r - " << count << " values loaded" << endl;
cout << " - min(" << m_x_min <<";"<< m_y_min <<";" << m_z_min<< ") max(" << m_x_max <<";"<< m_y_max <<";" << m_z_max<< ")" << endl;
cout << " - Rmax = " << m_Rmax << endl;
in.close();
}
NPLib/Detectors/Samurai/SamuraiFieldMap.h
View file @ f7d038b2
......@@ -28,11 +28,14 @@
#include"TObject.h"
#include"TGraph.h"
#include"TVector3.h"
#include "Math/Minimizer.h"
#include "Math/Functor.h"
#include "NPParticle.h"
class SamuraiFieldMap{
public:
SamuraiFieldMap(){m_BrhoScan=NULL;};
SamuraiFieldMap();
SamuraiFieldMap(std::string file);
~SamuraiFieldMap(){};
......@@ -45,31 +48,32 @@ class SamuraiFieldMap{
private:
// map[Pos]=B;
std::map<std::vector<float>,std::vector<float>> m_field;
float m_x_max,m_y_max,m_z_max,m_x_min,m_y_min,m_z_min;
std::map<std::vector<double>,std::vector<double>> m_field;
double m_x_max,m_y_max,m_z_max,m_x_min,m_y_min,m_z_min;
int m_bin;
double m_angle;
double m_Rmax ;
public: // getting the field at a point in space
// return B at an existing point
std::vector<float> GetB(std::vector<float>& pos);
inline std::vector<float> GetB(float x,float y ,float z){
std::vector<float> pos = {x,y,z};
std::vector<double> GetB(std::vector<double>& pos);
inline std::vector<double> GetB(double x,double y ,double z){
std::vector<double> pos = {x,y,z};
return GetB(pos);
};
// interpolate B witin volume (0 outside volume)
std::vector<float> InterpolateB(const std::vector<float>& pos);
std::vector<double> InterpolateB(const std::vector<double>& pos);
// interpolate B witin volume (0 outside volume)
inline std::vector<float> InterpolateB(const TVector3& pos){
std::vector<float> p={(float)pos.X(),(float)pos.Y(),(float)pos.Z()};
inline std::vector<double> InterpolateB(const TVector3& pos){
std::vector<double> p={(double)pos.X(),(double)pos.Y(),(double)pos.Z()};
return InterpolateB(p);
};
public: // Propagation of a particule in the field
// return a 3D track of the particle in the field
std::vector< TVector3 > Propagate(double rmax, double Brho, TVector3 pos, TVector3 dir);
void func(NPL::Particle& N, TVector3 pos, TVector3 dir, TVector3& new_pos, TVector3& new_dir);
std::vector< TVector3 > Propagate(double Brho, TVector3 pos, TVector3 dir,bool store=true);
void func(NPL::Particle& N, TVector3 pos, TVector3 imp, TVector3& new_pos, TVector3& new_dir);
private:
double m_fdc2angle;
double m_fdc2R;
......@@ -80,10 +84,17 @@ class SamuraiFieldMap{
public:
TGraph* BrhoScan(double min,double max,double step);
double FindBrho(TVector3& p_fdc0,TVector3& d_fdc0,TVector3& p_fdc2,TVector3& d_fdc2);
double FindBrho(TVector3 p_fdc0,TVector3 d_fdc0,TVector3 p_fdc2,TVector3 d_fdc2);
private:
TGraph* m_BrhoScan;
ROOT::Math::Minimizer* m_min;
ROOT::Math::Functor m_func;
double Delta(const double* parameter);
TVector3 m_FitPosFDC0,m_FitDirFDC0,m_FitPosFDC2,m_FitDirFDC2;
//
ClassDef(SamuraiFieldMap,1);
};
......
NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
View file @ f7d038b2
......@@ -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 @ f7d038b2
......@@ -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: