Commit 964aaae2 authored by Adrien Matta's avatar Adrien Matta
Browse files

* Progress on SamuraiFieldMap:

        - infer Brho based on fdc0/2 position and direction
parent fb0186b7
......@@ -30,18 +30,63 @@ using namespace NPUNITS;
using namespace std;
ClassImp(SamuraiFieldMap);
////////////////////////////////////////////////////////////////////////////////
double SamuraiFieldMap::FindBrho(TVector3& p_fdc0,TVector3& d_fdc0,TVector3& p_fdc2,TVector3& d_fdc2){
if(!m_BrhoScan)
BrhoScan(2.5,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;
}
////////////////////////////////////////////////////////////////////////////////
TGraph* SamuraiFieldMap::BrhoScan(double min, double max,double step){
if(m_BrhoScan)
delete m_BrhoScan;
m_BrhoScan=new TGraph;
unsigned int size = (max-min)/step;
m_BrhoScan->Set(size);
unsigned int i=0;
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);
}
return m_BrhoScan;
}
////////////////////////////////////////////////////////////////////////////////
TVector3 SamuraiFieldMap::PropagateToFDC2(TVector3 pos, TVector3 dir, double angle,double R){
TVector3 SamuraiFieldMap::PropagateToFDC2(TVector3 pos, TVector3 dir){
// go to FDC2 frame reference
pos.RotateY(-angle);
dir.RotateY(-angle);
pos.RotateY(-m_fdc2angle);
dir.RotateY(-m_fdc2angle);
double deltaZ=R-pos.Z();
double deltaZ=m_fdc2R-pos.Z();
dir*=deltaZ/dir.Z();
pos+=dir;
cout << pos.Z() << " " << R << endl;
pos.SetX(pos.X());
pos.RotateY(angle);
pos.RotateY(m_fdc2angle);
return pos;
}
......@@ -112,7 +157,7 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double rmax, double Brho, TVe
count++;
}
imp=imp.Unit();
pos = PropagateToFDC2(pos, imp, (59.930-90.0)*deg, 3686.77 + 880.745/2.);
pos = PropagateToFDC2(pos, imp);
pos.RotateY(-m_angle);
track.push_back(pos);
......
......@@ -26,12 +26,13 @@
#include<vector>
#include<map>
#include"TObject.h"
#include"TGraph.h"
#include"TVector3.h"
#include "NPParticle.h"
class SamuraiFieldMap{
public:
SamuraiFieldMap(){};
SamuraiFieldMap(){m_BrhoScan=NULL;};
SamuraiFieldMap(std::string file);
~SamuraiFieldMap(){};
......@@ -69,7 +70,20 @@ class SamuraiFieldMap{
// 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(TVector3 pos, TVector3 dir,double angle,double R);
private:
double m_fdc2angle;
double m_fdc2R;
public:
void SetFDC2Angle(double angle){m_fdc2angle=angle;};
void SetFDC2R(double R){m_fdc2R=R;};
TVector3 PropagateToFDC2(TVector3 pos, TVector3 dir);
public:
TGraph* BrhoScan(double min,double max,double step);
double FindBrho(TVector3& p_fdc0,TVector3& d_fdc0,TVector3& p_fdc2,TVector3& d_fdc2);
private:
TGraph* m_BrhoScan;
//
ClassDef(SamuraiFieldMap,1);
};
......
......@@ -6,7 +6,6 @@ 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());
}
......@@ -18,6 +17,10 @@ void testB(){
auto c = new TCanvas("trajectory","trajectory",1000,1000);
double angle = 30*deg;
field.LoadMap(angle,"field_map/180702-2,40T-3000.table.bin",10);
field.SetFDC2Angle((59.930-90.0)*deg);
field.SetFDC2R(3686.77 + 880.745/2.);
auto scan = field.BrhoScan(2.5,10,0.1);
unsigned int size = 1000;
vector<float> pos = {0,0,0};
auto h = new TH2F("h","h", size,-4000,4000,size,-4000,4000);
......@@ -35,12 +38,13 @@ void testB(){
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);
TVector3 p(1,1,-3500); TVector3 d(0.01,-0.01,1); double b = 5.481923;
DrawT(p,d,b);
std::vector< TVector3 > track = field.Propagate(3000,b,p,d);
cout << field.FindBrho(p,d,track.back(),d)<< endl;;
double rFDC2 = 3686.77 + 880.745/2.;
double phiFDC2 = (59.930-90.0-angle/deg)*deg;
TVector3 C(-252 ,0, rFDC2 );
......@@ -64,4 +68,11 @@ void testB(){
mag->SetLineWidth(4);
mag->SetFillStyle(0);
mag->Draw();
new TCanvas();
scan->Draw("ac");
cout << scan->Eval(700-252) << endl;
}
void testB2(){
SamuraiFieldMap field;
double angle = 30*deg;
field.LoadMap(angle,"field_map/180702-2,40T-3000.table.bin",10);
field.SetFDC2Angle((59.930-90.0)*deg);
field.SetFDC2R(3686.77 + 880.745/2.);
auto scan = field.BrhoScan(2.5,10,0.1);
auto h=new TH1D("dB","dB",1000,-2e-3,2e-3);
double r = gRandom->Uniform(-1,1);
for(unsigned int i = 0 ; i < 1000 ; i++){
TVector3 p(gRandom->Uniform(-5,5),gRandom->Uniform(-5,5),-3500); TVector3 d(gRandom->Uniform(-0.05,0.05),gRandom->Uniform(-0.05,0.05),1);
double b = gRandom->Uniform(3,7);
std::vector< TVector3 > track = field.Propagate(3000,b,p,d);
h->Fill(field.FindBrho(p,d,track.back(),d)-b);
}
h->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