Skip to content
Snippets Groups Projects
Commit 017cd118 authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

* Progress on s034 analysis 

        - BDC alignment on MINOS/SAMURAI
parent 1fc0a534
No related branches found
No related tags found
No related merge requests found
Pipeline #122270 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
+ 29
− 26
View file @ 017cd118
......@@ -58,7 +58,6 @@ double SamuraiFieldMap::Delta(const double* parameter){
////////////////////////////////////////////////////////////////////////////////
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;
......@@ -123,7 +122,8 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV
N.SetBrho(Brho);
// track result
std::vector< TVector3 > track;
static std::vector< TVector3 > track;
track.clear();
// starting point of the track
if(store){
......@@ -132,11 +132,13 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV
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>m_Rmax && count<limit){
......@@ -156,17 +158,18 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV
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 = 1*nanosecond;
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);
......@@ -193,26 +196,26 @@ std::vector< TVector3 > SamuraiFieldMap::Propagate(double Brho, TVector3 pos, TV
////////////////////////////////////////////////////////////////////////////////
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<double> 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));
......
Projects/S034/Analysis.cxx
+ 26
− 17
View file @ 017cd118
......@@ -49,7 +49,7 @@ void Analysis::Init(){
InitOutputBranch();
InitInputBranch();
// for fdc/bdc alignement
file.open("Calibration/Pos/fdc.txt");
//file.open("Calibration/Pos/bdc.txt");
}
////////////////////////////////////////////////////////////////////////////////
......@@ -64,13 +64,9 @@ void Analysis::TreatEvent(){
&& FDC2->PosY>-500 && FDC2->PosY<500
&& FDC0->PosX>-80 && FDC0->PosX<80
&& FDC0->PosY>-80 && FDC0->PosY<80 // both FDC ok
&& (Minos->Tracks_P0.size()==1)|| (Minos->Tracks_P0.size()==2)) { // p,pn or p,2p
&& (Minos->Tracks_P0.size()>1)) { // 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);
......@@ -81,7 +77,7 @@ void Analysis::TreatEvent(){
MinimumDistanceTwoLines(BDC1,BDC2,
Minos->Tracks_P0[0], P1,
Vertex, delta) ;
FDC0_Dir= FDC0->GetPos()-Vertex;
TVector3 FDC0_Dir= FDC0->GetPos()-Vertex;
FDC0_Dir=FDC0_Dir.Unit();
TVector3 BDCDir=BDC2-BDC1;
BDCDir=BDCDir.Unit();
......@@ -93,23 +89,36 @@ void Analysis::TreatEvent(){
Z=Vertex.Z()+4657.39;
//double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX};
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};
if(FDC0_Dir.Z()>0.6){
double FDC0_ThetaX = atan((FDC0->PosX-Vertex.X())/(1254.39-Z));
double FDC0_PhiY = atan((FDC0->PosY-Vertex.Y())/(1254.39-Z));
double brho_param[6]={FDC0->PosX, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX+252.416, 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;
// Calib//////////////////////////////////////////////////////////////////
/* static int count=0;
if(Minos->Delta_Vertex < 5 && FDC2->PosX-252.55>0&&FDC0->GetPos().X()>-10000 && FDC0->Dir.Z()>0.9 && Minos->Z_Vertex>0&& sqrt(Minos->X_Vertex*Minos->X_Vertex+Minos->Y_Vertex*Minos->Y_Vertex)<15){
file << FDC0->GetPos().X() <<" " << FDC0->GetPos().Y() << " " << FDC0->GetPos().Z() <<" " ;
file << Minos->X_Vertex <<" " << Minos->Y_Vertex << " " << Minos->Z_Vertex << " " ;
file << FDC2->GetPos().X() <<" " << FDC2->GetPos().Y() << " " << FDC2->GetPos().Z() <<" " << FDC2->Dir.X() <<" " << FDC2->Dir.Y() << " " << FDC2->Dir.Z()<< endl;
count ++;
}
if(count>1000)
exit(1);
*/
/* static int count=0;
if(Minos->Delta_Vertex < 5 && sqrt(Minos->X_Vertex*Minos->X_Vertex+Minos->Y_Vertex*Minos->Y_Vertex)<15 && Minos->Z_Vertex>-4650){
file << BDC1.X() <<" " << BDC1.Y()<< " " << BDC1.Z() <<" " ;
file << BDC2.X() <<" " << BDC2.Y()<< " " << BDC2.Z() <<" " ;
file << Minos->X_Vertex <<" " << Minos->Y_Vertex << " " << Minos->Z_Vertex << endl ;
count ++;
}
if(count>10000)
exit(1);*/
exit(1);
*/
}
}
}
......
Projects/S034/Calibration/Pos/FDC.cxx
+ 59
− 35
View file @ 017cd118
vector<TVector3*> pos0;
vector<TVector3*> dir0;
vector<TVector3*> posM;
vector<TVector3*> pos2;
vector<TVector3*> dir2;
SamuraiFieldMap field;
......@@ -10,23 +10,28 @@ 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;
TVector3 oM(parameter[1],parameter[2],parameter[3]);
TVector3 o0(parameter[4],parameter[5],parameter[6]);
TVector3 o2(parameter[7],parameter[8],parameter[9]);
TVector3 pM,p0,d0,p2;
oM.Print();
o0.Print();
o2.Print();
field.SetFDC2R(parameter[6]);
static auto h = new TH1D("h","h", 1000,0,10);
h->Reset();
for(unsigned int i = 0 ; i < 1000 ; i++){
for(unsigned int i = 0 ; i < size ; i++){
pM=*(posM[i]);
pM.RotateZ(parameter[0]*M_PI/180.);
pM+=oM;
p0=*(pos0[i])+o0;
p2=*(pos2[i])+o2;
if(dir0[i]->Z()>0.8)
h->Fill(field.FindBrho(p0,*dir0[i],p2,*dir2[i]));
d0=(p0-pM).Unit();
// p2.Print();
if(d0.Z()>0.9)
h->Fill(field.FindBrho(p0,d0,p2,*dir2[i]));
}
o0.Print();
o2.Print();
cout << h->GetStdDev() << endl;
cout << h->GetStdDev() << " " << parameter[0] << endl;
return h->GetStdDev();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -41,18 +46,18 @@ void LoadFile(){
cout << "Success opening file" << endl;
}
double xm,ym,zm;
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){
while(file >> x0 >> y0 >> z0 >> xm >> ym >> zm >> x2 >> y2 >> z2 >> dx2 >> dy2 >> dz2){
auto p0 = new TVector3(x0,y0,z0);
auto d0 = new TVector3(dx0,dy0,dz0);
auto pM = new TVector3(xm,ym,zm);
auto p2 = new TVector3(x2,y2,z2);
auto d2 = new TVector3(dx2,dy2,dz2);
pos0.push_back(p0);
dir0.push_back(d0);
posM.push_back(pM);
pos2.push_back(p2);
dir2.push_back(d2);
}
......@@ -69,28 +74,47 @@ void FDC(){
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);
auto func=ROOT::Math::Functor(&devB,10);
min->SetFunction(func);
min->SetPrintLevel(0);
min->SetPrecision(1e-3);
min->SetLimitedVariable(0,"X0",0,1,-10,10);
min->SetLimitedVariable(1,"Y0",0,1,-10,10);
min->SetLimitedVariable(2,"Z0",-3372.07,1,-3500,-3400);
min->SetLimitedVariable(3,"X2",-252.55,1,-260,-240);
min->SetLimitedVariable(4,"Y2",0,1,-10,10);
min->SetLimitedVariable(5,"Z2",4123.47,1,4000,4200);
min->SetPrecision(1e-6);
double parameter[10]={40.6,0,0,-4657.39, 0, 0,-3372.07,-252.55,0,4123.47};
devB(parameter);
// min->SetLimitedVariable(0,"AM",parameter[0],1,-90,90);
min->SetFixedVariable(0,"AM",parameter[0]);
min->SetLimitedVariable(1,"XM",parameter[1],1,-10,10);
min->SetLimitedVariable(2,"YM",parameter[2],1,-10,10);
min->SetLimitedVariable(3,"ZM",parameter[3],1,-4650,-4660);
min->SetLimitedVariable(4,"X0",parameter[4],1,-10,10);
min->SetLimitedVariable(5,"Y0",parameter[5],1,-10,10);
min->SetLimitedVariable(6,"Z0",parameter[6],1,-3370,-3375);
min->SetLimitedVariable(7,"X2",parameter[7],1,-260,-240);
min->SetLimitedVariable(8,"Y2",parameter[8],1,-10,10);
min->SetLimitedVariable(9,"Z2",parameter[9],1,4120,4125);
/*min->SetFixedVariable(1,"XM",parameter[1]);
min->SetFixedVariable(2,"YM",parameter[2]);
min->SetFixedVariable(3,"ZM",parameter[3]);
min->SetFixedVariable(4,"X0",parameter[4]);
min->SetFixedVariable(5,"Y0",parameter[5]);
min->SetFixedVariable(6,"Z0",parameter[6]);
min->SetFixedVariable(7,"X2",parameter[7]);
min->SetFixedVariable(8,"Y2",parameter[8]);
min->SetFixedVariable(9,"Z2",parameter[9]);
*/
min->Minimize();
const double* x = min->X();
cout << "X0 =" << x[0]<<endl;
cout << "Y0 =" << x[1]<<endl;
cout << "Z0 =" << x[2]<<endl;
cout << "X2 =" << x[3]<<endl;
cout << "Y2 =" << x[4]<<endl;
cout << "Z2 =" << x[5]<<endl;
cout << "AM =" << x[0]<<endl;
cout << "XM =" << x[1]<<endl;
cout << "YM =" << x[2]<<endl;
cout << "ZM =" << x[3]<<endl;
cout << "X0 =" << x[4]<<endl;
cout << "Y0 =" << x[5]<<endl;
cout << "Z0 =" << x[6]<<endl;
cout << "X2 =" << x[7]<<endl;
cout << "Y2 =" << x[8]<<endl;
cout << "Z2 =" << x[9]<<endl;
cout << "Minimum: " << devB(x) << endl;
}
Projects/S034/macro/rigz.cxx
+ 19
− 12
View file @ 017cd118
......@@ -4,9 +4,9 @@ void rigz(){
auto fz = new TFile("root/zaihong/run0582_RIG20210424_6He.root");
auto tz = (TTree*) fz->FindObjectAny("rig");
auto fl = new TFile("root/analysis/Results582.root");
auto tl = (TTree*) fl->FindObjectAny("ResultTree");
tl->AddFriend(tz);
auto fl = new TFile("root/analysis/Result582.root");
auto tl = (TTree*) fl->FindObjectAny("PhysicsTree");
tl->AddFriend(tz);
double FDC0_X,FDC0_Y,FDC2_X,FDC2_ThetaX,beta;
int FragID;
tz->SetBranchAddress("FDC0_X",&FDC0_X);
......@@ -24,11 +24,11 @@ void rigz(){
auto b = new TH1D("rig","rig",1000,2,8);
unsigned int entries = tz->GetEntries();
for(unsigned int i = 0 ; i < entries ; i++){
for(unsigned int i = 0 ; i < entries; i++){
tz->GetEntry(i);
double brho_param[6]={FDC0_X/*+1.77*/, FDC0_Y, 0, 0, FDC2_X/*-252.55*/, FDC2_ThetaX};
double Brho=r_fit(brho_param);
if(Brho>2&& Brho<8 && FragID>0 && FragID<27){
if(FragID>0 && FragID<27){
h->Fill(Brho);
// compute Brho based on beta and FragID
double rig ;
......@@ -48,18 +48,25 @@ void rigz(){
He6.SetBeta(beta);
rig = He6.GetBrho();
}
b->Fill(rig);
if(rig>2 && rig<8)
b->Fill(rig);
}
}
h->Scale(1./h->Integral());
h->Draw(); h->SetLineColor(kBlack);
b->Draw("same");
// h->Scale(1./h->Integral());
// h->Draw(); h->SetLineColor(kBlack);
b->Draw("");
b->Scale(1./b->Integral());
b->SetLineColor(kOrange+7);b->SetLineWidth(4);
cout << tl->Draw("Brho>>g","Brho>2 && Brho<8 &&FragID>0 && FragID<27","same") << endl;
b->SetLineColor(kOrange+7);b->SetLineWidth(2);
cout << "ratio Z: " << b->Integral(b->FindBin(3),b->FindBin(4.3))/ b->Integral(b->FindBin(4.5),b->FindBin(6.5));
cout << tl->Draw("Brho>>g(1000,2,8)","Brho>2 && Brho<8 &&FragID>0 && FragID<27","same") << endl;
cout << tl->Draw("BrhoP>>hp(1000,2,8)","Brho>2 && Brho<8 &&FragID>0 && FragID<27","same") << endl;
auto g = (TH1*) gDirectory->FindObjectAny("g");
g->SetLineColor(kAzure+7);g->SetLineWidth(4);
auto hp = (TH1*) gDirectory->FindObjectAny("hp");
g->SetLineColor(kAzure+7);g->SetLineWidth(2);
g->Scale(1./g->Integral());
cout << "ratio N: " << g->Integral(g->FindBin(3),g->FindBin(4.3))/ g->Integral(g->FindBin(4.5),g->FindBin(6.5));
hp->Scale(1./hp->Integral());
hp->SetLineColor(kBlack);b->SetLineWidth(2);
auto l = new TLine(3.62,0,3.62,h->GetMaximum());l->Draw();
l = new TLine(5.53,0,5.53,h->GetMaximum());l->Draw();
l = new TLine(5.48,0,5.48,h->GetMaximum());l->Draw();
......
Projects/S034/s034.detector
+ 7
− 11
View file @ 017cd118
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIBDC 1
XML= db/SAMURAIBDC1.xml
% Offset= -0.69 +1.1 -6876.34 mm
Offset= 0 0 -6876.34 mm
Offset= 0.959479 -0.17903 -6876.34 mm
InvertX= 0
InvertY= 0
InvertD= 1
......@@ -11,15 +10,16 @@ SAMURAIBDC 1
SAMURAIBDC 2
XML= db/SAMURAIBDC2.xml
% Offset= +0.25 +1.6 -5876.7 mm
Offset= 0 0 -5876.7 mm
% Offset= 0 0 -5876.7 mm
Offset= -0.0686464 0.294288 -5876.7 mm
InvertX= 0
InvertY= 0
InvertD= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minos
Position= 0 0 -4657.39 mm
ZRotation= 35 deg
Position= 0.345399 1.02061 -4650 mm
ZRotation= 40.6 deg
TargetLength= 151.72 mm
TargetMaterial= LH2
CellMaterial= Mylar
......@@ -33,10 +33,7 @@ Minos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIFDC0
XML= db/SAMURAIFDC0_20200109.xml
%Offset= 3.005 1.864 0 mm
% Offset= 1.77 0 -3372.07 mm
Offset= 4.52872 -9.97899 -3403.03 mm
Offset= 0 0 0 mm
Offset= -0.00666226 0.102191 -3370.01 mm
InvertX= 1
InvertY= 0
InvertD= 1
......@@ -44,8 +41,7 @@ SAMURAIFDC0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIFDC2
XML= db/SAMURAIFDC2.xml
Offset= -255.13 9.86883 4168.46 mm
% Offset= -252.55 0 4123.47 mm
Offset= -252.416 -0.228477 4122.57 mm
InvertX= 0
InvertY= 1
InvertD= 1
......
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