Commit 3025de0e authored by Adrien Matta's avatar Adrien Matta
Browse files

* Progress on SamuraiFieldMap

parent 9b1ab916
......@@ -19,6 +19,6 @@ add_custom_command(OUTPUT TSamuraiFDC0PhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR
add_library(NPSamurai SHARED SamuraiFieldMapDict.cxx SamuraiFieldMap.cxx TSamuraiHodoscopeData.cxx TSamuraiBDCData.cxx TSamuraiBDCPhysics.cxx TSamuraiBDCDataDict.cxx TSamuraiBDCPhysicsDict.cxx TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiHodoscopePhysics.cxx TSamuraiHodoscopePhysicsDict.cxx TSamuraiFDC2Physics.cxx TSamuraiHodoscopeDataDict.cxx TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC0Data.cxx TSamuraiFDC0DataDict.cxx TSamuraiFDC0Physics.cxx TSamuraiFDC0PhysicsDict.cxx)
target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore NPTrackReconstruction NPPhysics)
install(FILES SamuraiFieldMap.h TSamuraiBDCData.h TSamuraiBDCPhysics.h TSamuraiHodoscopeData.h TSamuraiHodoscopePhysics.h TSamuraiFDC2Data.h TSamuraiFDC2Physics.h TSamuraiFDC0Data.h TSamuraiFDC0Physics.h SamuraiDCIndex.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
......@@ -21,18 +21,115 @@
*****************************************************************************/
#include "SamuraiFieldMap.h"
#include "NPPhysicalConstants.h"
using namespace NPUNITS;
#include <fstream>
#include <iostream>
using namespace std;
ClassImp(SamuraiFieldMap);
////////////////////////////////////////////////////////////////////////////////
std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVector3 pos, TVector3 dir){
pos.RotateY(m_angle);
dir.RotateY(m_angle);
// Property of a particle with the correct Brho:
// We assume a 4He to compute v
// The choice of the particle is of no importance
static NPL::Particle N("4He");
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);
dir=dir.Unit();
double r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z());
// number of step taken
unsigned int count = 0;
// maximum number of state before giving up
unsigned int limit = 1000;
// First propagate to r_max with one line
while(r>rmax && count<limit){
pos+=(r-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);
}
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){
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);
r = sqrt(pos.X()*pos.X()+pos.Z()*pos.Z());
count++;
}
return track;
}
////////////////////////////////////////////////////////////////////////////////
SamuraiFieldMap::SamuraiFieldMap(std::string file){
LoadMap(file);
void SamuraiFieldMap::func(NPL::Particle& N, TVector3 pos, TVector3 imp, TVector3& new_pos, TVector3& new_imp){
double px,py,pz;
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;
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
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));
}
////////////////////////////////////////////////////////////////////////////////
void SamuraiFieldMap::LoadMap(std::string file){
void SamuraiFieldMap::LoadMap(double angle,std::string file,unsigned int bin){
m_bin=bin;
m_angle=angle;
if(file.find(".bin")!=std::string::npos)
LoadBinary(file);
else
......@@ -41,6 +138,15 @@ void SamuraiFieldMap::LoadMap(std::string file){
////////////////////////////////////////////////////////////////////////////////
std::vector<float> SamuraiFieldMap::GetB(std::vector<float>& pos){
static vector<float> 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];
if(pos[2]<0)
pos[2] = -pos[2];
auto it=m_field.find(pos);
if(it!=m_field.end()){
return it->second;
......@@ -50,21 +156,33 @@ std::vector<float> SamuraiFieldMap::GetB(std::vector<float>& pos){
}
////////////////////////////////////////////////////////////////////////////////
std::vector<float> SamuraiFieldMap::InterpolateB(std::vector<float>& pos){
std::vector<float> SamuraiFieldMap::InterpolateB(const std::vector<float>& pos){
static vector<float> nullv ={0,0,0};
double x = pos[0];
double y = pos[1];
double z = pos[2];
// 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
x = -pos[0];
y = pos[1];
if(pos[2]>0)
z = pos[2];
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;
m_bin = 10;
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);
......@@ -90,14 +208,14 @@ std::vector<float> SamuraiFieldMap::InterpolateB(std::vector<float>& pos){
unsigned int size = it.size();
for(unsigned int i = 0 ; i < size; i++){
if(it[i]!=end){
double d = 1e-6+sqrt( (pos[0]-it[i]->first[0])*(pos[0]-it[i]->first[0])+
(pos[1]-it[i]->first[1])*(pos[1]-it[i]->first[1])+
(pos[2]-it[i]->first[2])*(pos[2]-it[i]->first[2]));
Bx+=it[i]->second[0]/d;
By+=it[i]->second[1]/d;
Bz+=it[i]->second[2]/d;
totalW+=1./d;
double d = 1e-6+sqrt( (x-it[i]->first[0])*(x-it[i]->first[0])+
(y-it[i]->first[1])*(y-it[i]->first[1])+
(z-it[i]->first[2])*(z-it[i]->first[2]));
Bx+=it[i]->second[0]/(d*d);
By+=it[i]->second[1]/(d*d);
Bz+=it[i]->second[2]/(d*d);
totalW+=1./(d*d);
}
}
vector<float> res = {Bx/totalW,By/totalW,Bz/totalW};
......@@ -129,6 +247,9 @@ void SamuraiFieldMap::LoadAscii(std::string file){
if(++count%50000==0)
cout << "\r - Loading " << count << " values " << flush;
vector<float> p = {x,y,z};
Bx*=tesla;
By*=tesla;
Bz*=tesla;
vector<float> B = {Bx,By,Bz};
m_field[p]=B;
if(x<m_x_min)
......@@ -146,6 +267,7 @@ void SamuraiFieldMap::LoadAscii(std::string file){
}
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;
in.close();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -175,6 +297,9 @@ void SamuraiFieldMap::LoadBinary(std::string file){
in.read((char*)&Bz,sizeof(Bz));
vector<float> p = {x,y,z};
Bx*=tesla;
By*=tesla;
Bz*=tesla;
vector<float> B = {Bx,By,Bz};
m_field[p]=B;
if(x<m_x_min)
......
......@@ -22,11 +22,12 @@
* *
*****************************************************************************/
#include"TObject.h"
#include<string>
#include<vector>
#include<map>
#include"TObject.h"
#include"TVector3.h"
#include "NPParticle.h"
class SamuraiFieldMap{
public:
......@@ -35,7 +36,9 @@ class SamuraiFieldMap{
~SamuraiFieldMap(){};
public: // Map reading
void LoadMap(std::string file);
void LoadMap(double angle, std::string file, unsigned int bin);
private:
void LoadAscii(std::string file);
void LoadBinary(std::string file);
......@@ -44,16 +47,30 @@ class SamuraiFieldMap{
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;
int m_bin;
double m_angle;
public:
std::vector<float> GetB(std::vector<float>& pos);
std::vector<float> InterpolateB(std::vector<float>& pos);
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};
return GetB(pos);
};
// interpolate B witin volume (0 outside volume)
std::vector<float> InterpolateB(const std::vector<float>& 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()};
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);
//TVector3 PropagateToFDC2(double angle,double);
//
ClassDef(SamuraiFieldMap,1);
};
......
......@@ -118,6 +118,7 @@ namespace NPL {
double GetGamma() const {return fGamma;}
double GetVelocity() const {return fVelocity;}
TLorentzVector GetEnergyImpulsion() const {return fEnergyImpulsion;}
TVector3 GetImpulsion() const {return fEnergyImpulsion.Vect();}
double GetExcitationEnergy() const {return fExcitationEnergy;}
void SetName(const char* name) {fName = name;}
void SetZ(int charge) {fCharge = charge;}
......
This diff is collapsed.
SamuraiFieldMap field;
void DrawT(TVector3 pos, TVector3 dir, double Brho){
std::vector< TVector3 > track = field.Propagate(3000,Brho,pos,dir);
auto g = new TGraph();
unsigned int size = track.size();
g->Set(size);
cout << size << endl;
for(unsigned int i = 0 ; i < size ; i++){
g->SetPoint(i,-track[i].X(),track[i].Z());
}
g->SetLineWidth(6);
g->Draw("l");
}
////////////////////////////////////////////////////////////////////////////////
void testB(){
auto c = new TCanvas("trajectory","trajectory",1000,1000);
double angle = 0*deg;
field.LoadMap(angle,"field_map/180702-2,40T-3000.table.bin",10);
unsigned int size = 1000;
vector<float> pos = {0,0,0};
auto h = new TH2F("h","h", size,-4000,4000,size,-4000,4000);
for(int x = -3000 ; x < 3000 ; x+=10){
for(int z = -3000 ; z < 3000 ; z+=10){
pos[0]=x;pos[1]=0;pos[2]=z;
TVector3 p(x,0,z);
//float b = 1000*field.GetB(pos)[2];
float b = 1000*field.InterpolateB(p)[1];
p.RotateY(angle);
if(b!=0)
h->Fill(p.X(),p.Z(),b);
}
}
h->Draw("colz");
h->GetZaxis()->SetRangeUser(-1,3);
DrawT(TVector3(0,0,-3500),TVector3(0,0,1),5.48);
//DrawT(TVector3(0,0,-3500),TVector3(0.01,0.01,1),5);
//DrawT(TVector3(0,0,-3500),TVector3(-0.01,-0.01,1),5);
DrawT(TVector3(0,0,-3500),TVector3(0,0,1),3.62);
// DrawT(TVector3(0,0,-3500),TVector3(0.01,0.01,1),3.5);
// DrawT(TVector3(0,0,-3500),TVector3(-0.01,-0.01,1),3.5);
double rFDC2 = 3686.77 + 880.745/2.;
double phiFDC2 = (59.930-90.0-angle/deg)*deg;
TVector3 C(-252 ,0, rFDC2 );
TVector3 L(-252+1308 ,0, rFDC2 );
TVector3 R(-252-1308, 0, rFDC2 );
C.RotateY(phiFDC2);
L.RotateY(phiFDC2);
R.RotateY(phiFDC2);
auto FDC2 = new TLine(-L.X(),L.Z(),-R.X(),R.Z());
FDC2->SetLineWidth(4);
FDC2->SetLineColor(kRed);
FDC2->Draw();
auto CFDC2= new TMarker(-C.X(),C.Z(),20);
C.Print();
CFDC2->SetMarkerColor(kRed);
CFDC2->Draw();
auto mag = new TEllipse (0,0, 1000, 0);
mag->SetLineColor(kBlue);
mag->SetLineWidth(4);
mag->SetFillStyle(0);
mag->Draw();
}
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