Commit ba1b7c9c authored by Adrien Matta's avatar Adrien Matta
Browse files

* Progress on S034 analysis

parent 43680dd2
......@@ -22,6 +22,7 @@
#include "SamuraiFieldMap.h"
#include "NPPhysicalConstants.h"
#include "Math/Factory.h"
using namespace NPUNITS;
#include <cmath>
......@@ -31,33 +32,56 @@ 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);
// cout << parameter[0] << " " << pos.back().X()-m_FitPosFDC2.X() << endl;
// m_FitPosFDC2.Print();
// pos.back().Print();
double d = (pos.back().X()-m_FitPosFDC2.X())*(pos.back().X()-m_FitPosFDC2.X());
//diff = pos.back()-m_FitPosFDC2;
//diff.Print();
return d;
//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-5);
m_min->SetMaxFunctionCalls(1000);
m_min->SetLimitedVariable(0,"B",b0[0],1,1,10);
m_min->Minimize();
//cout << Delta(m_min->X()) << endl;
// exit(1);
//cout << m_min->NCalls() << endl;
return m_min->X()[0];
}
////////////////////////////////////////////////////////////////////////////////
......@@ -68,10 +92,14 @@ 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);
}
return m_BrhoScan;
}
......@@ -91,7 +119,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();
......@@ -103,12 +131,13 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
// track result
std::vector< TVector3 > track;
// starting point of the track
pos.RotateY(-m_angle);
track.push_back(pos);
pos.RotateY(m_angle);
// 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());
// number of step taken
......@@ -117,15 +146,17 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
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();
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;
......@@ -143,16 +174,18 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
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){
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,7 +193,7 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
pos = PropagateToFDC2(pos, imp);
pos.RotateY(-m_angle);
track.push_back(pos);
return track;
}
......@@ -182,7 +215,7 @@ void SamuraiFieldMap::func(NPL::Particle& N, TVector3 pos, TVector3 imp, TVector
new_pos.SetX(vx);
new_pos.SetY(vy);
new_pos.SetZ(vz);
vector<float> B = InterpolateB(pos);
vector<double> B = InterpolateB(pos);
double Bx= B[0];
double By= B[1];
double Bz= B[2];
......@@ -202,11 +235,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 +255,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 +272,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 +317,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 +330,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 +345,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 +365,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 +380,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 +396,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 +415,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();
}
......@@ -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);
};
......
......@@ -469,6 +469,9 @@ void TSamuraiFDC0Physics::InitializeRootInputRaw(){
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::InitializeRootInputPhysics(){
TChain* inputChain = RootInput::getInstance()->GetChain() ;
inputChain->SetBranchStatus( "SamuraiFDC0" , true );
inputChain->SetBranchAddress( "SamuraiFDC0" , &m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
......
......@@ -82,6 +82,13 @@ class TSamuraiFDC0Physics : public TObject, public NPL::VDetector{
int MultMean;
int PileUp;
private: // offset and inversion
TVector3 m_offset;//!
bool m_invertX;//!
bool m_invertY;//!
bool m_invertD;//!
public:
// Projected position at given Z plan
TVector3 ProjectedPosition(double Z);
......@@ -168,6 +175,7 @@ class TSamuraiFDC0Physics : public TObject, public NPL::VDetector{
TSamuraiFDC0Data* GetRawData() const {return m_EventData;}
TSamuraiFDC0Data* GetPreTreatedData() const {return m_PreTreatedData;}
TVector3 GetPos(){return TVector3(PosX,PosY,m_offset.Z());}
double GetPosX(){return PosX;}
double GetPosY(){return PosY;}
double GetThetaX(){return ThetaX;}
......@@ -181,12 +189,6 @@ class TSamuraiFDC0Physics : public TObject, public NPL::VDetector{
TSamuraiFDC0Data* m_PreTreatedData;//!
TSamuraiFDC0Physics* m_EventPhysics;//!
private: // offset and inversion
TVector3 m_offset;//!
bool m_invertX;//!
bool m_invertY;//!
bool m_invertD;//!
private: // Spectra Class
// TSamuraiFDC0Spectra* m_Spectra; // !
......
......@@ -397,6 +397,9 @@ void TSamuraiFDC2Physics::ReadConfiguration(NPL::InputParser parser){
}
#endif
GetOffset().Print();
PosX=1;
cout << m_invertY << endl;
}
///////////////////////////////////////////////////////////////////////////
......@@ -494,6 +497,9 @@ void TSamuraiFDC2Physics::InitializeRootInputRaw(){
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::InitializeRootInputPhysics(){
TChain* inputChain = RootInput::getInstance()->GetChain() ;
inputChain->SetBranchStatus( "SamuraiFDC2" , true );
inputChain->SetBranchAddress( "SamuraiFDC2" , &m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
......
......@@ -79,12 +79,19 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
int Mult;
int MultMean;
int PileUp;
private: // offset and inversion
TVector3 m_offset;//!
bool m_invertX;//!
bool m_invertY;//!
bool m_invertD;//!
public:
// Projected position at given Z plan
TVector3 ProjectedPosition(double Z);
double ProjectedPositionX(double Z);
double ProjectedPositionY(double Z);
private: // Charateristic of the DC
void AddDC(std::string name, NPL::XmlParser&);//! take the XML file and fill in Wire_X and Layer_Angle
std::map<SamuraiDCIndex,double> Wire_X;//! X position of the wires
......@@ -157,6 +164,8 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
// Write Spectra to file
void WriteSpectra();
public: // Specific to SamuraiFDC2 Array
// Remove bad channel, calibrate the data and apply threshold
......@@ -165,6 +174,10 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
// Retrieve raw and pre-treated data
TSamuraiFDC2Data* GetRawData() const {return m_EventData;}
TVector3 GetPos(){return TVector3(PosX,PosY,m_offset.Z());}
TVector3 GetOffset(){return m_offset;}
bool GetInvertX(){return m_invertX;};
bool GetInvertY(){return m_invertY;};
double GetPosX(){return PosX;}
double GetPosY(){return PosY;}
double GetThetaX(){return ThetaX;}
......@@ -176,12 +189,6 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
TSamuraiFDC2Data* m_EventData;//!
TSamuraiFDC2Physics* m_EventPhysics;//!
private: // offset and inversion
TVector3 m_offset;//!
bool m_invertX;//!
bool m_invertY;//!
bool m_invertD;//!
private: // Spectra Class
// TSamuraiFDC2Spectra* m_Spectra; // !
......
......@@ -43,8 +43,10 @@ void Analysis::Init(){
FDC0 = (TSamuraiFDC0Physics*) m_DetectorManager->GetDetector("SAMURAIFDC0");
FDC2 = (TSamuraiFDC2Physics*) m_DetectorManager->GetDetector("SAMURAIFDC2");
Hodo = (TSamuraiHodoscopePhysics*) m_DetectorManager->GetDetector("SAMURAIHOD");
// m_field.LoadMap("field_map/180702-2,40T-3000.table.bin",10);
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();
}
......@@ -56,7 +58,12 @@ void Analysis::TreatEvent(){
Trigger=Trigger&0x00ff;
//cout << Trigger << endl;
// Compute Brho
if(FDC2->PosX>-10000 && FDC0->PosX>-10000 ){ // if both are correctly build
//
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;
......@@ -67,8 +74,10 @@ void Analysis::TreatEvent(){
}
//double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX/*-252.55*/, FDC2->ThetaX};
double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX};
Brho=r_fit(brho_param);
//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));
}