diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx
index 3453e61c045196982a3931d0f4e600a75095ef17..1397989adf15ca0298567fe5286717383b559e18 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.cxx
+++ b/NPLib/Detectors/Minos/TMinosPhysics.cxx
@@ -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);
diff --git a/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx b/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
index 9968dce6bffe5fb62a6f1200083a0a000b86fca5..6c8e29329247a3720daa4ec980b2f44cb2b31302 100644
--- a/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
+++ b/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
@@ -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();
}
diff --git a/NPLib/Detectors/Samurai/SamuraiFieldMap.h b/NPLib/Detectors/Samurai/SamuraiFieldMap.h
index 09e5e1b7093485fc4767f61d4b9253815a557937..31ec83a54498afb5806ff10c0433d6527e7bedec 100644
--- a/NPLib/Detectors/Samurai/SamuraiFieldMap.h
+++ b/NPLib/Detectors/Samurai/SamuraiFieldMap.h
@@ -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);
};
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
index 9017c919d68a4ae77948733498a935f7f4f2f195..93871fe64c681818377404fa69dc3207e13c8781 100644
--- a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
@@ -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
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
index 37219d89cf01762b1bbeacb8831c184f9b7e1a1e..9fae5422e4c3c4b3692f92a3dd499524c6764095 100644
--- a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
@@ -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: // 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;//!
-
private: // Spectra Class
// TSamuraiBDCSpectra* m_Spectra; // !
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
index b1538925e436aa0884dbf3863d915501a989166b..f26eeb65502ac4f037b114863d05837e57a745c6 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
@@ -469,6 +469,9 @@ void TSamuraiFDC0Physics::InitializeRootInputRaw(){
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::InitializeRootInputPhysics(){
+ TChain* inputChain = RootInput::getInstance()->GetChain() ;
+ inputChain->SetBranchStatus( "SamuraiFDC0" , true );
+ inputChain->SetBranchAddress( "SamuraiFDC0" , &m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
index ae43f580d56268cc9f617f23c2c923ecbfb27913..bb17e6e5214afd20147c60fc376ca83535f9bc86 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
@@ -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; // !
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
index 7f2b5a7fdc99ad3c85a2d9b684461f5010c01ca5..9cd4e25e963ed4c71419e714cfe6f08f7dc58583 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
@@ -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);
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
index 99c03a43b478000246707d5702446cfa61cadbc0..3e370bccd2757b5534ceb71200ef8344a60e2f20 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
@@ -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; // !
diff --git a/Projects/S034/Analysis.cxx b/Projects/S034/Analysis.cxx
index d09bff10f22df38d8a32ee620464d26a0eba6066..544f76824ebf0db7a597f97ded69ad4c056fbacf 100644
--- a/Projects/S034/Analysis.cxx
+++ b/Projects/S034/Analysis.cxx
@@ -37,16 +37,19 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
- Minos= (TMinosPhysics*) m_DetectorManager->GetDetector("Minos");
- //Nebula = (TNebulaPhysics*) m_DetectorManager->GetDetector("NEBULA");
- BDC = (TSamuraiBDCPhysics*) m_DetectorManager->GetDetector("SAMURAIBDC");
- 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);
-
- InitOutputBranch();
- InitInputBranch();
+ Minos= (TMinosPhysics*) m_DetectorManager->GetDetector("Minos");
+ //Nebula = (TNebulaPhysics*) m_DetectorManager->GetDetector("NEBULA");
+ BDC = (TSamuraiBDCPhysics*) m_DetectorManager->GetDetector("SAMURAIBDC");
+ FDC0 = (TSamuraiFDC0Physics*) m_DetectorManager->GetDetector("SAMURAIFDC0");
+ FDC2 = (TSamuraiFDC2Physics*) m_DetectorManager->GetDetector("SAMURAIFDC2");
+ Hodo = (TSamuraiHodoscopePhysics*) m_DetectorManager->GetDetector("SAMURAIHOD");
+ 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();
+ // for fdc/bdc alignement
+ //file.open("Calibration/Pos/bdc.txt");
}
////////////////////////////////////////////////////////////////////////////////
@@ -54,21 +57,69 @@ void Analysis::TreatEvent(){
Clear();
//cout << Trigger << " " ;
Trigger=Trigger&0x00ff;
-//cout << Trigger << endl;
+ //cout << Trigger << endl;
// Compute Brho
- if(FDC2->PosX>-10000 && FDC0->PosX>-10000 ){ // 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;
+ //
+ 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 // both FDC ok
+ && (Minos->Tracks_P0.size()>1)) { // p,pn or p,2p
+ // Compute ThetaX and PhiY using Minos vertex and FDC0 X
+ // Check if both BDC are reconstructed
+ TVector3 BDC1=BDC->GetPos(1);
+ TVector3 BDC2=BDC->GetPos(2);
+
+ if( BDC1.Z()!=-10000 && BDC2.Z()!=-10000){
+ TVector3 Vertex;
+ TVector3 P1 = Minos->Tracks_P0[0]+Minos->Tracks_Dir[0];
+ TVector3 delta;
+ MinimumDistanceTwoLines(BDC1,BDC2,
+ Minos->Tracks_P0[0], P1,
+ Vertex, delta) ;
+ TVector3 FDC0_Dir= FDC0->GetPos()-Vertex;
+ FDC0_Dir=FDC0_Dir.Unit();
+ TVector3 BDCDir=BDC2-BDC1;
+ BDCDir=BDCDir.Unit();
+ BDCDir*=(Vertex.Z()-BDC2.Z())/BDCDir.Z();
+ //cout << "BDC" << endl;
+ BDCX=(BDC2+BDCDir).X();
+ BDCY=(BDC2+BDCDir).Y();
+ // Z relative to Minos entrance
+ Z=Vertex.Z()+4657.39;
+ //double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX};
- if(Minos->Z_Vertex>0){
- FDC0_ThetaX = atan((FDC0->PosX-Minos->X_Vertex)/(1283.7-Minos->Z_Vertex));
- FDC0_PhiY = atan((FDC0->PosY-Minos->Y_Vertex)/(1283.7-Minos->Z_Vertex));
- }
- //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(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);
+ */
- double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, 0, 0, FDC2->PosX/*-252.55*/, FDC2->ThetaX};
- Brho=r_fit(brho_param);
+ }
}
}
@@ -79,18 +130,26 @@ void Analysis::End(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Clear(){
Brho=-1000;
+ BDCX=-1000;
+ BDCY=-1000;
+ BrhoP=-1000;
Beta_f=-1000;
+ Z=-1000;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Brho",&Brho,"Brho/D");
+ RootOutput::getInstance()->GetTree()->Branch("BrhoP",&BrhoP,"BrhoP/D");
+ RootOutput::getInstance()->GetTree()->Branch("BDCX",&BDCX,"BDCX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BDCY",&BDCY,"BDCY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Z",&Z,"Z/D");
RootOutput::getInstance()->GetTree()->Branch("Beta_f",&Beta_f,"Beta_f/D");
RootOutput::getInstance()->GetTree()->Branch("Trigger",&Trigger,"Trigger/I");
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
- RootInput::getInstance()->GetChain()->SetBranchAddress("Trigger",&Trigger);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Trigger",&Trigger);
}
// -*- mode: c++ -*-
@@ -204,7 +263,7 @@ static double e_gCoefficient[] = {
-0.188488,
-0.675742,
-0.181284
- };
+};
// Assignment to error coefficients vector.
static double e_gCoefficientRMS[] = {
@@ -284,7 +343,7 @@ static double e_gCoefficientRMS[] = {
0.547759,
2.82987,
2.05714
- };
+};
// Assignment to powers vector.
// The powers are stored row-wise, that is
@@ -384,16 +443,16 @@ double e_fit(double *x) {
double v = 1 + 2. / (e_gXMax[j] - e_gXMin[j]) * (x[j] - e_gXMax[j]);
// what is the power to use!
switch(power) {
- case 1: r = 1; break;
- case 2: r = v; break;
- default:
- p2 = v;
- for (k = 3; k <= power; k++) {
- p3 = p2 * v;
- p3 = 2 * v * p2 - p1;
- p1 = p2; p2 = p3;
- }
- r = p3;
+ case 1: r = 1; break;
+ case 2: r = v; break;
+ default:
+ p2 = v;
+ for (k = 3; k <= power; k++) {
+ p3 = p2 * v;
+ p3 = 2 * v * p2 - p1;
+ p1 = p2; p2 = p3;
+ }
+ r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
@@ -516,7 +575,7 @@ static double p_gCoefficient[] = {
-0.145767,
-1.13422,
0.865558
- };
+};
// Assignment to error coefficients vector.
static double p_gCoefficientRMS[] = {
@@ -596,7 +655,7 @@ static double p_gCoefficientRMS[] = {
3.21076e-11,
8.11735e-11,
7.22006e-11
- };
+};
// Assignment to powers vector.
// The powers are stored row-wise, that is
@@ -696,16 +755,16 @@ double p_fit(double *x) {
double v = 1 + 2. / (p_gXMax[j] - p_gXMin[j]) * (x[j] - p_gXMax[j]);
// what is the power to use!
switch(power) {
- case 1: r = 1; break;
- case 2: r = v; break;
- default:
- p2 = v;
- for (k = 3; k <= power; k++) {
- p3 = p2 * v;
- p3 = 2 * v * p2 - p1;
- p1 = p2; p2 = p3;
- }
- r = p3;
+ case 1: r = 1; break;
+ case 2: r = v; break;
+ default:
+ p2 = v;
+ for (k = 3; k <= power; k++) {
+ p3 = p2 * v;
+ p3 = 2 * v * p2 - p1;
+ p1 = p2; p2 = p3;
+ }
+ r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
@@ -828,7 +887,7 @@ static double r_gCoefficient[] = {
-0.0019449,
-0.0151334,
0.0115488
- };
+};
// Assignment to error coefficients vector.
static double r_gCoefficientRMS[] = {
@@ -908,7 +967,7 @@ static double r_gCoefficientRMS[] = {
3.21076e-11,
8.11735e-11,
7.22006e-11
- };
+};
// Assignment to powers vector.
// The powers are stored row-wise, that is
@@ -1008,16 +1067,16 @@ double r_fit(double *x) {
double v = 1 + 2. / (r_gXMax[j] - r_gXMin[j]) * (x[j] - r_gXMax[j]);
// what is the power to use!
switch(power) {
- case 1: r = 1; break;
- case 2: r = v; break;
- default:
- p2 = v;
- for (k = 3; k <= power; k++) {
- p3 = p2 * v;
- p3 = 2 * v * p2 - p1;
- p1 = p2; p2 = p3;
- }
- r = p3;
+ case 1: r = 1; break;
+ case 2: r = v; break;
+ default:
+ p2 = v;
+ for (k = 3; k <= power; k++) {
+ p3 = p2 * v;
+ p3 = 2 * v * p2 - p1;
+ p1 = p2; p2 = p3;
+ }
+ r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
@@ -1141,7 +1200,7 @@ static double l_gCoefficient[] = {
-25.8219,
2.18449,
-0.917242
- };
+};
// Assignment to error coefficients vector.
static double l_gCoefficientRMS[] = {
@@ -1222,7 +1281,7 @@ static double l_gCoefficientRMS[] = {
31.7885,
20.2812,
22.29
- };
+};
// Assignment to powers vector.
// The powers are stored row-wise, that is
@@ -1323,16 +1382,16 @@ double l_fit(double *x) {
double v = 1 + 2. / (l_gXMax[j] - l_gXMin[j]) * (x[j] - l_gXMax[j]);
// what is the power to use!
switch(power) {
- case 1: r = 1; break;
- case 2: r = v; break;
- default:
- p2 = v;
- for (k = 3; k <= power; k++) {
- p3 = p2 * v;
- p3 = 2 * v * p2 - p1;
- p1 = p2; p2 = p3;
- }
- r = p3;
+ case 1: r = 1; break;
+ case 2: r = v; break;
+ default:
+ p2 = v;
+ for (k = 3; k <= power; k++) {
+ p3 = p2 * v;
+ p3 = 2 * v * p2 - p1;
+ p1 = p2; p2 = p3;
+ }
+ r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
@@ -1456,7 +1515,7 @@ static double t_gCoefficient[] = {
-0.0377476,
-0.150836,
0.0571845
- };
+};
// Assignment to error coefficients vector.
static double t_gCoefficientRMS[] = {
@@ -1537,7 +1596,7 @@ static double t_gCoefficientRMS[] = {
0.930273,
2.8245,
2.14083
- };
+};
// Assignment to powers vector.
// The powers are stored row-wise, that is
@@ -1638,16 +1697,16 @@ double t_fit(double *x) {
double v = 1 + 2. / (t_gXMax[j] - t_gXMin[j]) * (x[j] - t_gXMax[j]);
// what is the power to use!
switch(power) {
- case 1: r = 1; break;
- case 2: r = v; break;
- default:
- p2 = v;
- for (k = 3; k <= power; k++) {
- p3 = p2 * v;
- p3 = 2 * v * p2 - p1;
- p1 = p2; p2 = p3;
- }
- r = p3;
+ case 1: r = 1; break;
+ case 2: r = v; break;
+ default:
+ p2 = v;
+ for (k = 3; k <= power; k++) {
+ p3 = p2 * v;
+ p3 = 2 * v * p2 - p1;
+ p1 = p2; p2 = p3;
+ }
+ r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
@@ -1659,10 +1718,8 @@ double t_fit(double *x) {
}
// EOF for simtree_tof.C
-
-
////////////////////////////////////////////////////////////////////////////////
-// Construct Method to be pass to the DetectorFactory //
+// Construct Method to be pass to the AnalysisFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
return (NPL::VAnalysis*) new Analysis();
@@ -1672,13 +1729,13 @@ NPL::VAnalysis* Analysis::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy{
- public:
- proxy(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
-};
+ class proxy_analysis{
+ public:
+ proxy_analysis(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
-proxy p;
+ proxy_analysis p_analysis;
}
diff --git a/Projects/S034/Analysis.h b/Projects/S034/Analysis.h
index e5ace49dd1dce98fb0a3ffd6c1386eb528669aa2..7dfeb6e35196e4257c880fd6e5e278d3b75355af 100644
--- a/Projects/S034/Analysis.h
+++ b/Projects/S034/Analysis.h
@@ -29,7 +29,7 @@
#include"TSamuraiFDC2Physics.h"
#include"TSamuraiHodoscopePhysics.h"
#include"SamuraiFieldMap.h"
-
+#include<fstream>
class Analysis: public NPL::VAnalysis{
public:
Analysis();
@@ -50,14 +50,16 @@ class Analysis: public NPL::VAnalysis{
TSamuraiFDC2Physics* FDC2;
TSamuraiHodoscopePhysics* Hodo;
SamuraiFieldMap m_field ;
+ ofstream file;
private: // output variable
- double Brho;
+ double Brho,BrhoP,BDCX,BDCY,Z;
double Beta_f;
int Trigger;
public:
void Clear();
void InitOutputBranch();
void InitInputBranch();
+
};
// use for broh calculation
diff --git a/Projects/S034/Calibration/Pos/BDC.cxx b/Projects/S034/Calibration/Pos/BDC.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bebc79d18e64ce1973b658c2d2f742624f7577c7
--- /dev/null
+++ b/Projects/S034/Calibration/Pos/BDC.cxx
@@ -0,0 +1,98 @@
+
+
+vector<TVector3*> pos1;
+vector<TVector3*> pos2;
+vector<TVector3*> posM;
+
+////////////////////////////////////////////////////////////////////////////////
+double devR(const double* parameter){
+ // Return the standard deviation in Brho
+ unsigned int size = pos1.size();
+
+ TVector3 o1(parameter[0],parameter[1],parameter[2]);
+ TVector3 o2(parameter[3],parameter[4],parameter[5]);
+
+ double dR;
+ TVector3 p1,p2,pB,pM;
+ static auto h = new TH1D("h","h", 1000,-0,100);
+ h->Reset();
+ for(unsigned int i = 0 ; i < size ; i++){
+ pM=*(posM[i]);
+ p1=*(pos1[i])+o1;
+ p2=*(pos2[i])+o2;
+ TVector3 BDCDir=p2-p1;
+ BDCDir=BDCDir.Unit();
+ BDCDir*=(pM.Z()-p2.Z())/BDCDir.Z();
+ pB=p2+BDCDir;
+ dR=(pB-pM).Mag();
+ //cout << (pB-pM).Mag2() << " " << dR << endl;
+ double R = sqrt(pB.X()*pB.X()+pB.Y()*pB.Y());
+ if(R<15)
+ h->Fill(dR);
+ }
+ cout << h->GetMean() << " " << h->GetStdDev() << " " << parameter[0] << endl;
+ h->Draw();
+ gPad->Update();
+ return (h->GetMean()+h->GetStdDev() );
+}
+////////////////////////////////////////////////////////////////////////////////
+void LoadFile(){
+
+ ifstream file("Calibration/Pos/bdc.txt");
+ if(!file.is_open()){
+ cout << "fail to load file" << endl;
+ exit(1);
+ }
+ else {
+ cout << "Success opening file" << endl;
+ }
+
+ double xm,ym,zm;
+ double x1,y1,z1;
+ double x2,y2,z2;
+
+ while(file >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> xm >> ym >> zm ){
+ auto p1 = new TVector3(x1,y1,z1);
+ auto p2 = new TVector3(x2,y2,z2);
+ auto pM = new TVector3(xm,ym,zm);
+ pos1.push_back(p1);
+ pos2.push_back(p2);
+ posM.push_back(pM);
+ }
+ file.close();
+}
+////////////////////////////////////////////////////////////////////////////////
+void BDC(){
+
+ // Data
+ LoadFile();
+
+ // Minimizer
+ auto min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ auto func=ROOT::Math::Functor(&devR,6);
+ min->SetFunction(func);
+ min->SetPrintLevel(0);
+ min->SetPrecision(1e-10);
+ double parameter[6]={0,0,-6876.34,0,0,-5876.7};
+ devR(parameter);
+
+ // min->SetLimitedVariable(0,"AM",parameter[0],1,-90,90);
+ min->SetLimitedVariable(0,"X1",parameter[0],1,-10,10);
+ min->SetLimitedVariable(1,"Y1",parameter[1],1,-10,10);
+// min->SetLimitedVariable(2,"Z1",parameter[2],1,parameter[2]-100,parameter[2]+100);
+ min->SetFixedVariable(2,"Z1",parameter[2]);
+ min->SetLimitedVariable(3,"X2",parameter[3],1,-10,10);
+ min->SetLimitedVariable(4,"Y2",parameter[4],1,-10,10);
+ //min->SetLimitedVariable(5,"Z2",parameter[5],1,parameter[5]-100,parameter[5]+100);
+ min->SetFixedVariable(5,"Z2",parameter[5]);
+ min->Minimize();
+
+ const double* x = min->X();
+ cout << "X1 =" << x[0]<<endl;
+ cout << "Y1 =" << x[1]<<endl;
+ cout << "Z1 =" << x[2]<<endl;
+ cout << "X2 =" << x[3]<<endl;
+ cout << "Y2 =" << x[4]<<endl;
+ cout << "Z2 =" << x[5]<<endl;
+ cout << "Minimum: " << devR(x) << endl;
+}
diff --git a/Projects/S034/Calibration/Pos/FDC.cxx b/Projects/S034/Calibration/Pos/FDC.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0bf889b1c58fb8242ddf302a35683f3aa2170bc4
--- /dev/null
+++ b/Projects/S034/Calibration/Pos/FDC.cxx
@@ -0,0 +1,120 @@
+
+
+vector<TVector3*> pos0;
+vector<TVector3*> posM;
+vector<TVector3*> pos2;
+vector<TVector3*> dir2;
+SamuraiFieldMap field;
+
+////////////////////////////////////////////////////////////////////////////////
+double devB(const double* parameter){
+ // Return the standard deviation in Brho
+ unsigned int size = pos0.size();
+ 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 < size ; i++){
+ pM=*(posM[i]);
+ pM.RotateZ(parameter[0]*M_PI/180.);
+ pM+=oM;
+ p0=*(pos0[i])+o0;
+ p2=*(pos2[i])+o2;
+ d0=(p0-pM).Unit();
+// p2.Print();
+ if(d0.Z()>0.9)
+ h->Fill(field.FindBrho(p0,d0,p2,*dir2[i]));
+ }
+ cout << h->GetStdDev() << " " << parameter[0] << endl;
+ return h->GetStdDev();
+}
+////////////////////////////////////////////////////////////////////////////////
+void LoadFile(){
+
+ ifstream file("Calibration/Pos/fdc.txt");
+ if(!file.is_open()){
+ cout << "fail to load file" << endl;
+ exit(1);
+ }
+ else {
+ cout << "Success opening file" << endl;
+ }
+
+ double xm,ym,zm;
+ double x0,y0,z0;
+ double x2,y2,z2;
+ double 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 pM = new TVector3(xm,ym,zm);
+ auto p2 = new TVector3(x2,y2,z2);
+ auto d2 = new TVector3(dx2,dy2,dz2);
+ pos0.push_back(p0);
+ posM.push_back(pM);
+ pos2.push_back(p2);
+ dir2.push_back(d2);
+ }
+ cout << " Val " << pos0.size() << " Value Loaded" << endl;
+ file.close();
+}
+////////////////////////////////////////////////////////////////////////////////
+void FDC(){
+
+ // Data
+ LoadFile();
+
+ // Field map
+ field.LoadMap(30*deg,"field_map/180702-2,40T-3000.table.bin",10);
+ field.SetFDC2Angle((59.930-90.0)*deg);
+
+ // Minimizer
+ auto min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ auto func=ROOT::Math::Functor(&devB,10);
+ min->SetFunction(func);
+ min->SetPrintLevel(0);
+ 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 << "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;
+}
diff --git a/Projects/S034/Calibration/Pos/fdc.txt b/Projects/S034/Calibration/Pos/fdc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Projects/S034/field_map/ascii2bin.cxx b/Projects/S034/field_map/ascii2bin.cxx
index ef2e7ea617bb37e4ca3eba1a43ee6088355af8f9..5b3153ea23430dd77c6ec2168f83ea1f3c4afe35 100644
--- a/Projects/S034/field_map/ascii2bin.cxx
+++ b/Projects/S034/field_map/ascii2bin.cxx
@@ -5,9 +5,9 @@ void ascii2bin(){
cout << "Error: failed to load samurai field map " << file << endl;
exit(1);
}
-
+
cout << "//////// Loading Samurai Field Map " << file << endl;
- float x,y,z,Bx,By,Bz;
+ double x,y,z,Bx,By,Bz;
std::string name = file+".bin";
ofstream out(name.c_str(),std::ofstream::binary);
@@ -19,19 +19,27 @@ void ascii2bin(){
getline(in,buffer);
}
- while(in >> x >> y >> z >> Bx >> By >> Bz){
+ while(!in.eof()){
if(++count%10000==0)
- cout << "\r - Loading " << count << " values " << flush;
- out.write((char*)&x,sizeof(x));
- out.write((char*)&y,sizeof(y));
- out.write((char*)&z,sizeof(z));
- out.write((char*)&Bx,sizeof(Bx));
- out.write((char*)&By,sizeof(By));
- out.write((char*)&Bz,sizeof(Bz));
- //out << x << y << z << Bx << By << Bz;
+ cout << "\r - Loading " << count << " values " << flush;
+
+ if(in >> x >> y >> z >> Bx >> By >> Bz){
+ out.write((char*)&x,sizeof(x));
+ out.write((char*)&y,sizeof(y));
+ out.write((char*)&z,sizeof(z));
+ out.write((char*)&Bx,sizeof(Bx));
+ out.write((char*)&By,sizeof(By));
+ out.write((char*)&Bz,sizeof(Bz));
+ }
+ else{
+ cout << x <<" " << y << " "<< z << " "<< Bx <<" "<< By << " " << Bz << endl;;
+ if(!in.eof())
+ in.clear();
+ }
+ // cout << x << endl;
}
- cout << "\r - " << count << " values loaded" << endl;
+ cout << "\r - " << count << " values loaded" << endl;
out.close();
in.close();
diff --git a/Projects/S034/macro/betaz.cxx b/Projects/S034/macro/betaz.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b6abb85a180d231f4e852a15661ae7b86bcbf9ac
--- /dev/null
+++ b/Projects/S034/macro/betaz.cxx
@@ -0,0 +1,49 @@
+void betaz(){
+
+ auto fz = new TFile("root/zaihong/run0582_RIG20210424_6He.root");
+ auto tz = (TTree*) fz->FindObjectAny("rig");
+ auto fl = new TFile("root/analysis/Result582.root");
+ auto tl = (TTree*) fl->FindObjectAny("PhysicsTree");
+ tl->AddFriend(tz);
+ double Brho,betaZ,beta;
+ int FragID;
+ tl->SetBranchAddress("Brho",&Brho);
+ tl->SetBranchAddress("fBeta",&betaZ);
+ tl->SetBranchAddress("FragID",&FragID);
+ NPL::Particle H2("2H");
+ NPL::Particle H3("3H");
+ NPL::Particle He4("4He");
+ NPL::Particle He6("6He");
+
+ auto z = new TH1D("betaZ","betaZ",100,0.45,0.6);
+ auto b = new TH1D("betaF","betaZ",100,0.45,0.6);
+ auto d = new TH1D("betaF","betaZ",1000,-0.01,0.01);
+ unsigned int entries = tl->GetEntries();
+
+ for(unsigned int i = 0 ; i < entries; i++){
+ tl->GetEntry(i);
+ if(FragID>0 && FragID<27){
+ // compute Brho based on beta and FragID
+ double rig ;
+ if(FragID==26){
+ He6.SetBrho(Brho);
+ beta= He6.GetBeta();
+ b->Fill(beta);
+ z->Fill(betaZ);
+ d->Fill(beta-betaZ);
+ }
+ }
+ }
+ // h->Scale(1./h->Integral());
+ // h->Draw(); h->SetLineColor(kBlack);
+ b->Draw("");
+ b->Scale(1./b->Integral());
+ b->SetLineColor(kAzure+7);b->SetLineWidth(2);
+ z->Scale(1./z->Integral());
+ z->SetLineColor(kOrange+7);z->SetLineWidth(2);
+ z->Draw("same");
+ new TCanvas();
+ d->Draw();
+ d->SetLineColor(kAzure+7);b->SetLineWidth(2);
+}
+
diff --git a/Projects/S034/macro/rigz.cxx b/Projects/S034/macro/rigz.cxx
index 22c4f402a3498c93454ffa4c71fc33b68f50da28..7da954ecf10aad232c0532b2dd77a3962a621bbc 100644
--- a/Projects/S034/macro/rigz.cxx
+++ b/Projects/S034/macro/rigz.cxx
@@ -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/test582.root");
+ 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);
@@ -20,73 +20,66 @@ void rigz(){
NPL::Particle He4("4He");
NPL::Particle He6("6He");
- auto h = new TH1D("brho","brho",1000,0,8);
- auto b = new TH1D("rig","rig",1000,0,8);
- auto b1 = new TH1D("rig1","rig1",1000,0,8);
- auto b2 = new TH1D("rig2","rig2",1000,0,8);
- auto b3 = new TH1D("rig3","rig3",1000,0,8);
- auto b4 = new TH1D("rig4","rig4",1000,0,8);
+ auto h = new TH1D("brho","brho",1000,2,8);
+ 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>0&&FragID>0 && FragID<27){
+ if(FragID>0 && FragID<27){
h->Fill(Brho);
// compute Brho based on beta and FragID
double rig ;
if(FragID==12){
H2.SetBeta(beta);
rig = H2.GetBrho();
- b1->Fill(rig);
}
else if(FragID==13){
H3.SetBeta(beta);
rig = H3.GetBrho();
- b2->Fill(rig);
}
else if(FragID==24){
He4.SetBeta(beta);
rig = He4.GetBrho();
- b3->Fill(rig);
}
else if(FragID==26){
He6.SetBeta(beta);
rig = He6.GetBrho();
- b4->Fill(rig);
}
- b->Fill(rig);
+ if(rig>2 && rig<8)
+ b->Fill(rig);
}
}
// h->Scale(1./h->Integral());
- h->Draw(); h->SetLineColor(kBlack);
- b->Draw("same");
- b->SetLineColor(kOrange+7);b->SetLineWidth(4);
- cout << tl->Draw("Brho>>g","Brho>0&&SamuraiFDC2.devX<10","same") << endl;
+// h->Draw(); h->SetLineColor(kBlack);
+ b->Draw("");
+ b->Scale(1./b->Integral());
+ 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);
-// g->Scale(1./g->Integral());
+ 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();
- new TCanvas();
- b->Draw();
- b1->Draw("same");
- b2->Draw("same");
- b3->Draw("same");
- b4->Draw("same");
}
-// EOF for ../Geant4/samurai/simtree_energy.C
// -*- mode: c++ -*-
//
-// File ../Geant4/samurai/simtree_momentum.C generated by TMultiDimFit::MakeRealCode
-// on Tue Jun 9 23:41:17 2020
+// File simtree_energy.C generated by TMultiDimFit::MakeRealCode
+// on Fri May 28 00:25:23 2021
// ROOT version 5.32/00
//
// This file contains the function
//
-// double p_fit(double *x);
+// double e_fit(double *x);
//
// For evaluating the parameterization obtained
// from TMultiDimFit and the point x
@@ -96,204 +89,573 @@ void rigz(){
//
// Static data variables
//
-static int p_gNVariables = 6;
-static int p_gNCoefficients = 57;
-static double p_gDMean = 524.342;
+static int e_gNVariables = 6;
+static int e_gNCoefficients = 76;
+static double e_gDMean = 69.5201;
// Assignment to mean vector.
-static double p_gXMean[] = {
- 0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+static double e_gXMean[] = {
+ 0.00503955, -0.0144156, -0.101042, -0.00374735, 0.252452, -0.512868 };
// Assignment to minimum vector.
-static double p_gXMin[] = {
- -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+static double e_gXMin[] = {
+ -64.3965, -52.0748, -0.0335722, -0.0401226, -1065.23, -299.807 };
// Assignment to maximum vector.
-static double p_gXMax[] = {
- 329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+static double e_gXMax[] = {
+ 52.3459, 55.8258, 0.0432727, 0.04243, 1147.99, 931.885 };
// Assignment to coefficients vector.
-static double p_gCoefficient[] = {
- -4150.11,
- -7349.66,
- -1376.07,
- 726.65,
- -604.755,
- 157.175,
- -7402.61,
- 137.906,
- 1266,
- 607.703,
- 9432.11,
- 2336.63,
- -12820,
- -981.145,
- -2219,
- -1146.1,
- 1069.02,
- -9607.95,
- -135.98,
- 236.78,
- -5477.97,
- 2213.85,
- 6282.03,
- 89.2823,
- -10.974,
- -104.594,
- -62.7333,
- -9.89892,
- -21.8523,
- -18.8223,
- 16.9863,
- 3546.43,
- -461.655,
- -17.9853,
- -2026.01,
- -134.173,
- 5565.29,
- 18164.2,
- 5239.56,
- -4192.18,
- 22.9582,
- -15479.1,
- -114.231,
- 2974.51,
- -10402.5,
- 2516.63,
- 25.4881,
- 303.465,
- 17.2719,
- 1202.72,
- 50.8832,
- -1.17495,
- 1284.76,
- 17.0545,
- -0.894684,
- -136.872,
- -134.371
+static double e_gCoefficient[] = {
+ -871.253,
+ -1919.06,
+ -788.479,
+ 2645.64,
+ 443.115,
+ 984.07,
+ 243.005,
+ 3552.25,
+ -157.556,
+ 2188.2,
+ 618.533,
+ 1176.32,
+ 1303.87,
+ -1056.05,
+ -772.684,
+ 6.86646,
+ 1343.46,
+ 3390.39,
+ 871.459,
+ 7388.91,
+ -599.788,
+ -4023.35,
+ -299.575,
+ 1670.13,
+ 0.620256,
+ 10.0888,
+ -0.557305,
+ -1.12306,
+ 9.51285,
+ 0.728805,
+ -0.620644,
+ -4.05522,
+ 0.633792,
+ -11.1877,
+ -4.56183,
+ -5.71034,
+ -4.45376,
+ 1.38591,
+ 1.18557,
+ -4.92592,
+ -5.16675,
+ -4.66713,
+ -1506.82,
+ -306.44,
+ 4269.09,
+ 2292.62,
+ 19.7398,
+ 2612.26,
+ 1682.7,
+ -1163.29,
+ -597.13,
+ -7838.7,
+ 14404.4,
+ 6607.08,
+ -532.085,
+ -317.801,
+ 681.575,
+ -546.576,
+ -775.207,
+ -3803.59,
+ -1160.61,
+ 0.0848953,
+ 4.22188,
+ -1.96972,
+ 0.794997,
+ -2260.46,
+ -1.21301,
+ -7409.46,
+ 269.548,
+ -0.00925447,
+ 2.77788,
+ 1.95742,
+ 0.453985,
+ -0.188488,
+ -0.675742,
+ -0.181284
};
// Assignment to error coefficients vector.
-static double p_gCoefficientRMS[] = {
- 5.0702e-12,
- 9.69107e-12,
- 3.89769e-11,
- 1.13183e-11,
- 1.53099e-11,
- 5.35431e-12,
- 8.94257e-12,
- 1.12034e-11,
- 2.16374e-11,
- 6.32208e-12,
- 6.87455e-11,
- 2.70028e-11,
- 1.70926e-11,
- 2.9263e-11,
- 1.4199e-11,
- 5.24502e-12,
- 6.53991e-12,
- 4.1161e-11,
- 1.19525e-11,
- 1.99626e-11,
- 3.47403e-11,
- 9.32059e-11,
- 1.19932e-10,
- 6.39059e-12,
- 2.49985e-11,
- 4.91273e-11,
- 1.92969e-11,
- 2.95574e-10,
- 1.39773e-11,
- 6.12734e-11,
- 3.16965e-11,
- 1.26652e-10,
- 1.15348e-11,
- 1.5818e-11,
- 1.00252e-11,
- 1.09154e-10,
- 2.29228e-10,
- 1.31396e-10,
- 5.16128e-11,
- 2.71397e-11,
- 2.1153e-10,
- 7.86728e-11,
- 9.25091e-12,
- 1.78175e-10,
- 6.64139e-11,
- 1.25002e-11,
- 8.61292e-11,
- 3.38293e-11,
- 1.16568e-11,
- 1.64392e-10,
- 1.42658e-11,
- 4.97717e-11,
- 1.20839e-11,
- 4.82003e-11,
- 1.9685e-11,
- 4.28747e-11,
- 1.11506e-11
+static double e_gCoefficientRMS[] = {
+ 66.6599,
+ 130.954,
+ 45.1508,
+ 146.993,
+ 49.7567,
+ 42.2773,
+ 8.59505,
+ 193.537,
+ 31.5244,
+ 92.7045,
+ 67.4399,
+ 97.1825,
+ 63.1346,
+ 63.8809,
+ 37.0527,
+ 0.0875267,
+ 99.2757,
+ 205.578,
+ 127.446,
+ 361.809,
+ 107.175,
+ 235.193,
+ 35.8634,
+ 105.166,
+ 2.78933,
+ 2.0882,
+ 3.14393,
+ 3.06267,
+ 8.25249,
+ 0.616472,
+ 1.47943,
+ 2.11076,
+ 0.930677,
+ 7.39957,
+ 2.26951,
+ 2.50688,
+ 0.95993,
+ 2.42862,
+ 0.393511,
+ 0.458493,
+ 2.0148,
+ 2.57894,
+ 72.1773,
+ 61.4191,
+ 180.529,
+ 189.341,
+ 7.78812,
+ 193.422,
+ 248.274,
+ 208.754,
+ 69.8817,
+ 458.225,
+ 704.917,
+ 400.536,
+ 30.129,
+ 31.4199,
+ 27.3927,
+ 40.2961,
+ 54.0649,
+ 207.698,
+ 79.9551,
+ 9.4176,
+ 1.89081,
+ 0.422792,
+ 1.61055,
+ 155.773,
+ 1.45066,
+ 404.667,
+ 36.5963,
+ 0.729026,
+ 2.55087,
+ 0.355723,
+ 1.83536,
+ 0.547759,
+ 2.82987,
+ 2.05714
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
-// p_ij = p_gPower[i * NVariables + j];
-static int p_gPower[] = {
+// p_ij = e_gPower[i * NVariables + j];
+static int e_gPower[] = {
1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 2,
2, 1, 1, 1, 1, 1,
1, 1, 1, 1, 2, 1,
- 1, 2, 1, 1, 1, 1,
1, 1, 1, 2, 1, 1,
1, 1, 2, 1, 1, 1,
1, 1, 1, 1, 1, 3,
1, 1, 1, 1, 2, 2,
1, 3, 1, 1, 1, 1,
2, 1, 2, 1, 1, 1,
+ 1, 1, 1, 2, 1, 2,
1, 2, 2, 1, 1, 1,
1, 1, 2, 1, 1, 2,
- 1, 2, 1, 1, 1, 2,
- 1, 1, 3, 1, 1, 1,
+ 2, 1, 1, 1, 1, 2,
3, 1, 1, 1, 1, 1,
1, 1, 1, 1, 3, 1,
2, 1, 1, 2, 1, 1,
1, 1, 1, 2, 2, 1,
+ 1, 2, 1, 2, 1, 1,
1, 1, 2, 1, 2, 1,
- 1, 2, 1, 1, 2, 1,
+ 1, 1, 2, 2, 1, 1,
2, 1, 1, 1, 2, 1,
- 2, 2, 1, 1, 1, 1,
1, 1, 1, 3, 1, 1,
1, 1, 1, 1, 2, 3,
2, 1, 1, 3, 1, 1,
1, 2, 1, 3, 1, 1,
- 2, 2, 1, 1, 2, 1,
+ 1, 1, 2, 3, 1, 1,
+ 1, 1, 3, 2, 1, 1,
+ 1, 2, 2, 2, 1, 1,
1, 3, 1, 1, 2, 1,
+ 1, 3, 1, 2, 1, 1,
1, 2, 2, 1, 2, 1,
1, 1, 3, 1, 2, 1,
2, 2, 1, 2, 1, 1,
+ 2, 1, 1, 2, 2, 1,
+ 1, 2, 1, 2, 2, 1,
+ 1, 1, 1, 3, 2, 1,
+ 3, 1, 1, 2, 1, 1,
+ 2, 1, 1, 1, 3, 1,
1, 1, 2, 1, 3, 1,
3, 2, 1, 1, 1, 1,
+ 1, 2, 3, 1, 1, 1,
3, 1, 1, 1, 1, 2,
- 2, 1, 3, 1, 1, 1,
- 2, 2, 1, 1, 1, 2,
+ 1, 3, 1, 1, 1, 2,
2, 1, 2, 1, 1, 2,
1, 2, 2, 1, 1, 2,
- 1, 1, 3, 1, 1, 2,
2, 2, 2, 1, 1, 1,
2, 1, 1, 2, 1, 2,
- 3, 1, 2, 1, 1, 1,
+ 1, 2, 1, 2, 1, 2,
+ 1, 1, 2, 2, 1, 2,
+ 1, 1, 1, 3, 1, 2,
2, 1, 1, 1, 2, 2,
- 1, 2, 1, 1, 2, 2,
- 1, 1, 1, 1, 3, 2,
+ 1, 1, 2, 1, 2, 2,
+ 1, 1, 1, 2, 2, 2,
2, 1, 1, 1, 1, 3,
1, 2, 1, 1, 1, 3,
+ 1, 1, 2, 1, 1, 3,
+ 1, 2, 1, 1, 1, 1,
+ 1, 2, 1, 1, 1, 2,
+ 1, 2, 1, 1, 2, 1,
+ 2, 2, 1, 1, 1, 1,
+ 2, 1, 2, 2, 1, 1,
+ 2, 2, 1, 1, 2, 1,
+ 1, 1, 1, 2, 3, 1,
+ 2, 1, 3, 1, 1, 1,
+ 2, 2, 1, 1, 1, 2,
+ 3, 1, 2, 1, 1, 1,
+ 1, 2, 1, 1, 2, 2,
+ 1, 1, 1, 2, 1, 3,
3, 1, 1, 1, 2, 1,
+ 1, 1, 2, 2, 2, 1,
+ 1, 2, 1, 1, 3, 1,
+ 2, 3, 1, 1, 1, 1,
+ 1, 1, 3, 1, 1, 1,
2, 1, 2, 1, 2, 1,
+ 1, 3, 2, 1, 1, 1
+};
+
+//
+// The function double e_fit(double *x)
+//
+double e_fit(double *x) {
+ double returnValue = e_gDMean;
+ int i = 0, j = 0, k = 0;
+ for (i = 0; i < e_gNCoefficients ; i++) {
+ // Evaluate the ith term in the expansion
+ double term = e_gCoefficient[i];
+ for (j = 0; j < e_gNVariables; j++) {
+ // Evaluate the polynomial in the jth variable.
+ int power = e_gPower[e_gNVariables * i + j];
+ double p1 = 1, p2 = 0, p3 = 0, r = 0;
+ double v = 1 + 2. / (e_gXMax[j] - e_gXMin[j]) * (x[j] - e_gXMax[j]);
+ // what is the power to use!
+ switch(power) {
+ case 1: r = 1; break;
+ case 2: r = v; break;
+ default:
+ p2 = v;
+ for (k = 3; k <= power; k++) {
+ p3 = p2 * v;
+ p3 = 2 * v * p2 - p1;
+ p1 = p2; p2 = p3;
+ }
+ r = p3;
+ }
+ // multiply this term by the poly in the jth var
+ term *= r;
+ }
+ // Add this term to the final result
+ returnValue += term;
+ }
+ return returnValue;
+}
+
+// EOF for simtree_energy.C
+// -*- mode: c++ -*-
+//
+// File simtree_momentum.C generated by TMultiDimFit::MakeRealCode
+// on Fri May 28 00:25:23 2021
+// ROOT version 5.32/00
+//
+// This file contains the function
+//
+// double p_fit(double *x);
+//
+// For evaluating the parameterization obtained
+// from TMultiDimFit and the point x
+//
+// See TMultiDimFit class documentation for more information
+//
+//
+// Static data variables
+//
+static int p_gNVariables = 6;
+static int p_gNCoefficients = 76;
+static double p_gDMean = 361.833;
+// Assignment to mean vector.
+static double p_gXMean[] = {
+ 0.00503955, -0.0144156, -0.101042, -0.00374735, 0.252452, -0.512868 };
+
+// Assignment to minimum vector.
+static double p_gXMin[] = {
+ -64.3965, -52.0748, -0.0335722, -0.0401226, -1065.23, -299.807 };
+
+// Assignment to maximum vector.
+static double p_gXMax[] = {
+ 52.3459, 55.8258, 0.0432727, 0.04243, 1147.99, 931.885 };
+
+// Assignment to coefficients vector.
+static double p_gCoefficient[] = {
+ -1679.5,
+ -3768.96,
+ -1773.35,
+ 5927.56,
+ -1052.24,
+ 786.483,
+ 2247.73,
+ 553.038,
+ 7920.2,
+ -326.613,
+ 1147.19,
+ 3010.2,
+ -1520.59,
+ -2398.64,
+ 12.0318,
+ 2542.51,
+ 7843.96,
+ 1874.12,
+ 17655.7,
+ -902.702,
+ -9667.32,
+ -2259.82,
+ -673.898,
+ 3744.62,
+ 23.4495,
+ -8.41591,
+ 23.1301,
+ 23.3445,
+ -0.850973,
+ -21.3182,
+ -5.51343,
+ -13.2501,
+ 0.0420299,
+ 4.8444,
+ 4.76409,
+ -11.8105,
+ -4406.58,
+ -637.169,
+ 8498.88,
+ 4695.17,
+ 4938.56,
+ 3623.63,
+ -1738.8,
+ -1348.34,
+ 34420.8,
+ 15289.1,
+ -1177.71,
+ -587.346,
+ 1529.61,
+ 450.733,
+ 4356.32,
+ 2406.82,
+ 0.828932,
+ -1494.72,
+ -7.21335,
+ -4.61121,
+ 2.36389,
+ -9.47311,
+ 9.7888,
+ -3.50166,
+ -12.2852,
+ -2912.85,
+ -0.277797,
+ 42.1753,
+ -1.26653,
+ -18825.4,
+ -8755.82,
+ 3.69156,
+ 8.16369,
+ -7.45288,
+ -0.273504,
+ -9.0052,
+ -17056.8,
+ -0.145767,
+ -1.13422,
+ 0.865558
+ };
+
+// Assignment to error coefficients vector.
+static double p_gCoefficientRMS[] = {
+ 1.69031e-11,
+ 3.29223e-11,
+ 6.40666e-11,
+ 2.89085e-11,
+ 6.58543e-11,
+ 8.69785e-11,
+ 7.0898e-11,
+ 3.57517e-11,
+ 5.62969e-11,
+ 1.89724e-11,
+ 1.69397e-10,
+ 1.38078e-10,
+ 1.28278e-10,
+ 1.24785e-10,
+ 2.44521e-11,
+ 3.47315e-10,
+ 1.44807e-10,
+ 2.73924e-10,
+ 1.18391e-10,
+ 3.51888e-10,
+ 1.21525e-10,
+ 2.55786e-10,
+ 1.82279e-11,
+ 6.11644e-11,
+ 7.41771e-11,
+ 9.78353e-11,
+ 1.16067e-09,
+ 1.31561e-09,
+ 3.16376e-11,
+ 1.03428e-09,
+ 6.01156e-10,
+ 3.13162e-11,
+ 9.7543e-11,
+ 9.29464e-11,
+ 8.80057e-11,
+ 1.16199e-10,
+ 4.98203e-10,
+ 3.6952e-11,
+ 4.19168e-10,
+ 5.29444e-10,
+ 6.7647e-10,
+ 5.33585e-10,
+ 6.85262e-10,
+ 3.55029e-11,
+ 2.30548e-10,
+ 2.81963e-10,
+ 1.35502e-10,
+ 1.39257e-10,
+ 1.49982e-10,
+ 1.83996e-10,
+ 2.15232e-10,
+ 2.71861e-10,
+ 1.88434e-11,
+ 1.92304e-11,
+ 7.69361e-11,
+ 1.02682e-10,
+ 3.52953e-10,
+ 6.53796e-10,
+ 5.36745e-10,
+ 1.27112e-10,
+ 7.41853e-11,
+ 3.7455e-11,
+ 7.64436e-11,
+ 9.64104e-10,
+ 8.7203e-11,
+ 2.36674e-10,
+ 1.34232e-10,
+ 6.85812e-11,
+ 5.67547e-10,
+ 4.46781e-10,
+ 3.2384e-11,
+ 7.34435e-11,
+ 2.61417e-10,
+ 3.21076e-11,
+ 8.11735e-11,
+ 7.22006e-11
+ };
+
+// Assignment to powers vector.
+// The powers are stored row-wise, that is
+// p_ij = p_gPower[i * NVariables + j];
+static int p_gPower[] = {
+ 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 2,
+ 2, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 2, 1,
+ 1, 2, 1, 1, 1, 1,
+ 1, 1, 1, 2, 1, 1,
+ 1, 1, 2, 1, 1, 1,
+ 1, 1, 1, 1, 1, 3,
+ 1, 1, 1, 1, 2, 2,
+ 1, 3, 1, 1, 1, 1,
+ 1, 1, 1, 2, 1, 2,
+ 1, 1, 2, 1, 1, 2,
+ 1, 2, 1, 1, 1, 2,
+ 2, 1, 1, 1, 1, 2,
+ 1, 1, 1, 1, 3, 1,
+ 2, 1, 1, 2, 1, 1,
+ 1, 1, 1, 2, 2, 1,
+ 1, 2, 1, 2, 1, 1,
+ 1, 1, 2, 1, 2, 1,
+ 1, 1, 2, 2, 1, 1,
+ 2, 1, 1, 1, 2, 1,
+ 2, 2, 1, 1, 1, 1,
+ 1, 1, 1, 3, 1, 1,
+ 1, 1, 1, 1, 2, 3,
+ 1, 2, 1, 3, 1, 1,
+ 1, 1, 3, 2, 1, 1,
+ 1, 2, 2, 2, 1, 1,
+ 2, 1, 2, 2, 1, 1,
+ 1, 3, 1, 1, 2, 1,
+ 2, 2, 1, 2, 1, 1,
+ 1, 2, 1, 2, 2, 1,
1, 1, 1, 3, 2, 1,
+ 3, 1, 1, 2, 1, 1,
2, 1, 1, 1, 3, 1,
- 1, 3, 1, 1, 1, 2,
- 2, 3, 1, 1, 1, 1,
1, 2, 1, 1, 3, 1,
+ 1, 1, 2, 1, 3, 1,
+ 2, 2, 1, 1, 1, 2,
+ 1, 3, 1, 1, 1, 2,
+ 2, 1, 2, 1, 1, 2,
+ 1, 2, 2, 1, 1, 2,
+ 2, 1, 1, 2, 1, 2,
+ 1, 2, 1, 2, 1, 2,
+ 1, 1, 2, 2, 1, 2,
+ 1, 1, 1, 3, 1, 2,
+ 1, 1, 2, 1, 2, 2,
+ 1, 1, 1, 2, 2, 2,
+ 2, 1, 1, 1, 1, 3,
+ 1, 2, 1, 1, 1, 3,
+ 1, 1, 2, 1, 1, 3,
+ 1, 1, 1, 2, 1, 3,
+ 2, 1, 2, 1, 1, 1,
+ 1, 2, 2, 1, 1, 1,
+ 1, 1, 3, 1, 1, 1,
+ 3, 1, 1, 1, 1, 1,
+ 1, 1, 2, 3, 1, 1,
+ 1, 3, 1, 2, 1, 1,
+ 2, 1, 2, 1, 2, 1,
+ 2, 1, 1, 2, 2, 1,
+ 1, 1, 2, 2, 2, 1,
+ 1, 1, 1, 2, 3, 1,
+ 3, 2, 1, 1, 1, 1,
+ 3, 1, 1, 1, 1, 2,
+ 2, 1, 3, 1, 1, 1,
+ 2, 2, 2, 1, 1, 1,
+ 3, 1, 2, 1, 1, 1,
+ 2, 1, 1, 1, 2, 2,
+ 1, 2, 1, 1, 2, 1,
+ 2, 1, 1, 3, 1, 1,
+ 2, 2, 1, 1, 2, 1,
+ 1, 2, 2, 1, 2, 1,
+ 1, 1, 3, 1, 2, 1,
1, 2, 3, 1, 1, 1,
- 1, 3, 2, 1, 1, 1
+ 1, 2, 1, 1, 2, 2,
+ 3, 1, 1, 1, 2, 1,
+ 1, 3, 2, 1, 1, 1,
+ 2, 3, 1, 1, 1, 1
};
//
@@ -332,11 +694,11 @@ double p_fit(double *x) {
return returnValue;
}
-// EOF for ../Geant4/samurai/simtree_momentum.C
+// EOF for simtree_momentum.C
// -*- mode: c++ -*-
//
-// File ../Geant4/samurai/simtree_rigidity.C generated by TMultiDimFit::MakeRealCode
-// on Tue Jun 9 23:41:17 2020
+// File simtree_rigidity.C generated by TMultiDimFit::MakeRealCode
+// on Fri May 28 00:25:23 2021
// ROOT version 5.32/00
//
// This file contains the function
@@ -352,140 +714,178 @@ double p_fit(double *x) {
// Static data variables
//
static int r_gNVariables = 6;
-static int r_gNCoefficients = 57;
-static double r_gDMean = 5.24705;
+static int r_gNCoefficients = 76;
+static double r_gDMean = 4.82778;
// Assignment to mean vector.
static double r_gXMean[] = {
- 0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+ 0.00503955, -0.0144156, -0.101042, -0.00374735, 0.252452, -0.512868 };
// Assignment to minimum vector.
static double r_gXMin[] = {
- -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+ -64.3965, -52.0748, -0.0335722, -0.0401226, -1065.23, -299.807 };
// Assignment to maximum vector.
static double r_gXMax[] = {
- 329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+ 52.3459, 55.8258, 0.0432727, 0.04243, 1147.99, 931.885 };
// Assignment to coefficients vector.
static double r_gCoefficient[] = {
- -41.5298,
- -73.5475,
- -13.7703,
- 7.27153,
- -6.05174,
- 1.57283,
- -74.0774,
- 1.38001,
- 12.6688,
- 6.08124,
- 94.3864,
- 23.3825,
- -128.289,
- -9.81824,
- -22.2054,
- -11.4689,
- 10.6976,
- -96.146,
- -1.36074,
- 2.36944,
- -54.8176,
- 22.1538,
- 62.8638,
- 0.893442,
- -0.109816,
- -1.04666,
- -0.627767,
- -0.0990577,
- -0.218674,
- -0.188353,
- 0.16998,
- 35.4888,
- -4.61974,
- -0.179978,
- -20.2741,
- -1.34266,
- 55.6915,
- 181.768,
- 52.4319,
- -41.9508,
- 0.229741,
- -154.898,
- -1.1431,
- 29.7657,
- -104.097,
- 25.1837,
- 0.255057,
- 3.03675,
- 0.172838,
- 12.0355,
- 0.509184,
- -0.0117577,
- 12.8565,
- 0.170664,
- -0.00895304,
- -1.36967,
- -1.34464
+ -22.4088,
+ -50.2876,
+ -23.6611,
+ 79.0889,
+ -14.0396,
+ 10.4937,
+ 29.9905,
+ 7.37895,
+ 105.676,
+ -4.35786,
+ 15.3065,
+ 40.1638,
+ -20.2886,
+ -32.0041,
+ 0.160535,
+ 33.9237,
+ 104.659,
+ 25.0055,
+ 235.572,
+ -12.0444,
+ -128.987,
+ -30.1518,
+ -8.99153,
+ 49.9628,
+ 0.312876,
+ -0.11229,
+ 0.308614,
+ 0.311475,
+ -0.0113542,
+ -0.28444,
+ -0.0735633,
+ -0.17679,
+ 0.000560787,
+ 0.0646367,
+ 0.0635652,
+ -0.157583,
+ -58.795,
+ -8.50147,
+ 113.397,
+ 62.6455,
+ 65.8931,
+ 48.3485,
+ -23.2001,
+ -17.9904,
+ 459.262,
+ 203.996,
+ -15.7137,
+ -7.8367,
+ 20.4089,
+ 6.01393,
+ 58.1245,
+ 32.1131,
+ 0.0110601,
+ -19.9434,
+ -0.0962446,
+ -0.0615253,
+ 0.0315403,
+ -0.126396,
+ 0.130608,
+ -0.0467211,
+ -0.163916,
+ -38.8649,
+ -0.00370653,
+ 0.562727,
+ -0.0168988,
+ -251.179,
+ -116.825,
+ 0.0492549,
+ 0.108925,
+ -0.0994405,
+ -0.00364924,
+ -0.120152,
+ -227.582,
+ -0.0019449,
+ -0.0151334,
+ 0.0115488
};
// Assignment to error coefficients vector.
static double r_gCoefficientRMS[] = {
- 5.0702e-12,
- 9.69107e-12,
- 3.89769e-11,
- 1.13183e-11,
- 1.53099e-11,
- 5.35431e-12,
- 8.94257e-12,
- 1.12034e-11,
- 2.16374e-11,
- 6.32208e-12,
- 6.87455e-11,
- 2.70028e-11,
- 1.70926e-11,
- 2.9263e-11,
- 1.4199e-11,
- 5.24502e-12,
- 6.53991e-12,
- 4.1161e-11,
- 1.19525e-11,
- 1.99626e-11,
- 3.47403e-11,
- 9.32059e-11,
- 1.19932e-10,
- 6.39059e-12,
- 2.49985e-11,
- 4.91273e-11,
- 1.92969e-11,
- 2.95574e-10,
- 1.39773e-11,
- 6.12734e-11,
- 3.16965e-11,
- 1.26652e-10,
- 1.15348e-11,
- 1.5818e-11,
- 1.00252e-11,
- 1.09154e-10,
- 2.29228e-10,
- 1.31396e-10,
- 5.16128e-11,
- 2.71397e-11,
- 2.1153e-10,
- 7.86728e-11,
- 9.25091e-12,
- 1.78175e-10,
- 6.64139e-11,
- 1.25002e-11,
- 8.61292e-11,
- 3.38293e-11,
- 1.16568e-11,
- 1.64392e-10,
- 1.42658e-11,
- 4.97717e-11,
- 1.20839e-11,
- 4.82003e-11,
- 1.9685e-11,
- 4.28747e-11,
- 1.11506e-11
+ 1.69031e-11,
+ 3.29223e-11,
+ 6.40666e-11,
+ 2.89085e-11,
+ 6.58543e-11,
+ 8.69785e-11,
+ 7.0898e-11,
+ 3.57517e-11,
+ 5.62969e-11,
+ 1.89724e-11,
+ 1.69397e-10,
+ 1.38078e-10,
+ 1.28278e-10,
+ 1.24785e-10,
+ 2.44521e-11,
+ 3.47315e-10,
+ 1.44807e-10,
+ 2.73924e-10,
+ 1.18391e-10,
+ 3.51888e-10,
+ 1.21525e-10,
+ 2.55786e-10,
+ 1.82279e-11,
+ 6.11644e-11,
+ 7.41771e-11,
+ 9.78353e-11,
+ 1.16067e-09,
+ 1.31561e-09,
+ 3.16376e-11,
+ 1.03428e-09,
+ 6.01156e-10,
+ 3.13162e-11,
+ 9.7543e-11,
+ 9.29464e-11,
+ 8.80057e-11,
+ 1.16199e-10,
+ 4.98203e-10,
+ 3.6952e-11,
+ 4.19168e-10,
+ 5.29444e-10,
+ 6.7647e-10,
+ 5.33585e-10,
+ 6.85262e-10,
+ 3.55029e-11,
+ 2.30548e-10,
+ 2.81963e-10,
+ 1.35502e-10,
+ 1.39257e-10,
+ 1.49982e-10,
+ 1.83996e-10,
+ 2.15232e-10,
+ 2.71861e-10,
+ 1.88434e-11,
+ 1.92304e-11,
+ 7.69361e-11,
+ 1.02682e-10,
+ 3.52953e-10,
+ 6.53796e-10,
+ 5.36745e-10,
+ 1.27112e-10,
+ 7.41853e-11,
+ 3.7455e-11,
+ 7.64436e-11,
+ 9.64104e-10,
+ 8.7203e-11,
+ 2.36674e-10,
+ 1.34232e-10,
+ 6.85812e-11,
+ 5.67547e-10,
+ 4.46781e-10,
+ 3.2384e-11,
+ 7.34435e-11,
+ 2.61417e-10,
+ 3.21076e-11,
+ 8.11735e-11,
+ 7.22006e-11
};
// Assignment to powers vector.
@@ -502,53 +902,72 @@ static int r_gPower[] = {
1, 1, 1, 1, 1, 3,
1, 1, 1, 1, 2, 2,
1, 3, 1, 1, 1, 1,
- 2, 1, 2, 1, 1, 1,
- 1, 2, 2, 1, 1, 1,
+ 1, 1, 1, 2, 1, 2,
1, 1, 2, 1, 1, 2,
1, 2, 1, 1, 1, 2,
- 1, 1, 3, 1, 1, 1,
- 3, 1, 1, 1, 1, 1,
+ 2, 1, 1, 1, 1, 2,
1, 1, 1, 1, 3, 1,
2, 1, 1, 2, 1, 1,
1, 1, 1, 2, 2, 1,
+ 1, 2, 1, 2, 1, 1,
1, 1, 2, 1, 2, 1,
- 1, 2, 1, 1, 2, 1,
+ 1, 1, 2, 2, 1, 1,
2, 1, 1, 1, 2, 1,
2, 2, 1, 1, 1, 1,
1, 1, 1, 3, 1, 1,
1, 1, 1, 1, 2, 3,
- 2, 1, 1, 3, 1, 1,
1, 2, 1, 3, 1, 1,
- 2, 2, 1, 1, 2, 1,
+ 1, 1, 3, 2, 1, 1,
+ 1, 2, 2, 2, 1, 1,
+ 2, 1, 2, 2, 1, 1,
1, 3, 1, 1, 2, 1,
- 1, 2, 2, 1, 2, 1,
- 1, 1, 3, 1, 2, 1,
2, 2, 1, 2, 1, 1,
+ 1, 2, 1, 2, 2, 1,
+ 1, 1, 1, 3, 2, 1,
+ 3, 1, 1, 2, 1, 1,
+ 2, 1, 1, 1, 3, 1,
+ 1, 2, 1, 1, 3, 1,
1, 1, 2, 1, 3, 1,
- 3, 2, 1, 1, 1, 1,
- 3, 1, 1, 1, 1, 2,
- 2, 1, 3, 1, 1, 1,
2, 2, 1, 1, 1, 2,
+ 1, 3, 1, 1, 1, 2,
2, 1, 2, 1, 1, 2,
1, 2, 2, 1, 1, 2,
- 1, 1, 3, 1, 1, 2,
- 2, 2, 2, 1, 1, 1,
2, 1, 1, 2, 1, 2,
- 3, 1, 2, 1, 1, 1,
- 2, 1, 1, 1, 2, 2,
- 1, 2, 1, 1, 2, 2,
- 1, 1, 1, 1, 3, 2,
+ 1, 2, 1, 2, 1, 2,
+ 1, 1, 2, 2, 1, 2,
+ 1, 1, 1, 3, 1, 2,
+ 1, 1, 2, 1, 2, 2,
+ 1, 1, 1, 2, 2, 2,
2, 1, 1, 1, 1, 3,
1, 2, 1, 1, 1, 3,
- 3, 1, 1, 1, 2, 1,
+ 1, 1, 2, 1, 1, 3,
+ 1, 1, 1, 2, 1, 3,
+ 2, 1, 2, 1, 1, 1,
+ 1, 2, 2, 1, 1, 1,
+ 1, 1, 3, 1, 1, 1,
+ 3, 1, 1, 1, 1, 1,
+ 1, 1, 2, 3, 1, 1,
+ 1, 3, 1, 2, 1, 1,
2, 1, 2, 1, 2, 1,
- 1, 1, 1, 3, 2, 1,
- 2, 1, 1, 1, 3, 1,
- 1, 3, 1, 1, 1, 2,
- 2, 3, 1, 1, 1, 1,
- 1, 2, 1, 1, 3, 1,
+ 2, 1, 1, 2, 2, 1,
+ 1, 1, 2, 2, 2, 1,
+ 1, 1, 1, 2, 3, 1,
+ 3, 2, 1, 1, 1, 1,
+ 3, 1, 1, 1, 1, 2,
+ 2, 1, 3, 1, 1, 1,
+ 2, 2, 2, 1, 1, 1,
+ 3, 1, 2, 1, 1, 1,
+ 2, 1, 1, 1, 2, 2,
+ 1, 2, 1, 1, 2, 1,
+ 2, 1, 1, 3, 1, 1,
+ 2, 2, 1, 1, 2, 1,
+ 1, 2, 2, 1, 2, 1,
+ 1, 1, 3, 1, 2, 1,
1, 2, 3, 1, 1, 1,
- 1, 3, 2, 1, 1, 1
+ 1, 2, 1, 1, 2, 2,
+ 3, 1, 1, 1, 2, 1,
+ 1, 3, 2, 1, 1, 1,
+ 2, 3, 1, 1, 1, 1
};
//
@@ -587,11 +1006,11 @@ double r_fit(double *x) {
return returnValue;
}
-// EOF for ../Geant4/samurai/simtree_rigidity.C
+// EOF for simtree_rigidity.C
// -*- mode: c++ -*-
//
-// File ../Geant4/samurai/simtree_length.C generated by TMultiDimFit::MakeRealCode
-// on Tue Jun 9 23:41:17 2020
+// File simtree_length.C generated by TMultiDimFit::MakeRealCode
+// on Fri May 28 00:25:23 2021
// ROOT version 5.32/00
//
// This file contains the function
@@ -607,138 +1026,180 @@ double r_fit(double *x) {
// Static data variables
//
static int l_gNVariables = 6;
-static int l_gNCoefficients = 56;
-static double l_gDMean = 9302.09;
+static int l_gNCoefficients = 77;
+static double l_gDMean = 9270.11;
// Assignment to mean vector.
static double l_gXMean[] = {
- 0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+ 0.00503955, -0.0144156, -0.101042, -0.00374735, 0.252452, -0.512868 };
// Assignment to minimum vector.
static double l_gXMin[] = {
- -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+ -64.3965, -52.0748, -0.0335722, -0.0401226, -1065.23, -299.807 };
// Assignment to maximum vector.
static double l_gXMax[] = {
- 329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+ 52.3459, 55.8258, 0.0432727, 0.04243, 1147.99, 931.885 };
// Assignment to coefficients vector.
static double l_gCoefficient[] = {
- -5063.74,
- -8210.4,
- 8531.82,
- -734.412,
- -360.971,
- -8616.61,
- 331.099,
- -854.689,
- 13989.9,
- -14949,
- -2334.33,
- 23725.3,
- -232.89,
- 863.204,
- -5135.67,
- -671.923,
- 193.744,
- -7892.44,
- 5470.48,
- 2106.04,
- 220.14,
- 392.627,
- -750.825,
- -13.1656,
- -15.4854,
- -15.0127,
- -692.948,
- -46.1235,
- -10.895,
- -2288.18,
- 13.9869,
- -7.15341,
- 0.581274,
- 42.6349,
- -8.22167,
- -116.588,
- 8175.36,
- 1080.47,
- 129.214,
- 26910.7,
- 7002.79,
- -4598.69,
- -14.7257,
- 2944.04,
- -467.226,
- -15008,
- 2105.01,
- 6.06541,
- 361.147,
- -1073.29,
- 3745.34,
- -284.423,
- 1666.15,
- -499.767,
- 8786.59,
- 119.522
+ -6590.16,
+ -13914.6,
+ -2835.19,
+ 11286.2,
+ -2318.17,
+ 1506.93,
+ 2806.41,
+ 1055.27,
+ 15246.9,
+ -582.244,
+ 14287.3,
+ 2024.12,
+ 5826.32,
+ 2883.06,
+ -3211.64,
+ -2565.14,
+ -3508.49,
+ -3945.73,
+ 81.3328,
+ 6976.65,
+ 15491.2,
+ 3340.02,
+ 32961.2,
+ -4257.58,
+ -18445.7,
+ -1204.58,
+ 7213.73,
+ 1.15164,
+ 56.7961,
+ -6.55251,
+ -3.1322,
+ 32.3364,
+ 17.5881,
+ 2.39586,
+ -17.4049,
+ 4.99685,
+ -44.3394,
+ -18.7891,
+ -32.238,
+ -26.3708,
+ 6.54775,
+ 7.12538,
+ -7694.27,
+ -1133.08,
+ 27850.4,
+ 11364,
+ -5002.96,
+ 13573.4,
+ 6437.24,
+ -8275.75,
+ -2420.82,
+ -5.63532,
+ -35904.8,
+ 64281.3,
+ 30190.8,
+ -1440.27,
+ 2620.77,
+ 1005.06,
+ -17864.9,
+ -5198.67,
+ 1.04599,
+ -1.61611,
+ -3.90405,
+ -15.7528,
+ -9.48741,
+ -25.0997,
+ -10129.3,
+ 2.95363,
+ 99.1311,
+ -34800.1,
+ -2335.5,
+ -8.51776,
+ 11.7194,
+ 9.32626,
+ -25.8219,
+ 2.18449,
+ -0.917242
};
// Assignment to error coefficients vector.
static double l_gCoefficientRMS[] = {
- 1092.36,
- 1874.66,
- 2316.91,
- 992.765,
- 293.062,
- 1327.01,
- 164.035,
- 1469.78,
- 1822.87,
- 2291.78,
- 587.934,
- 2842.24,
- 222.135,
- 300.266,
- 2772.37,
- 341.844,
- 284.914,
- 849.513,
- 1225.2,
- 5295.19,
- 203.834,
- 313.389,
- 883.734,
- 23.072,
- 38.4748,
- 15.8304,
- 233.869,
- 387.214,
- 119.356,
- 5419.4,
- 456.337,
- 6.95297,
- 5.94522,
- 60.7779,
- 168.339,
- 150.975,
- 1134.29,
- 329.863,
- 198.355,
- 3547.52,
- 1330.81,
- 1077.78,
- 656.951,
- 1360.04,
- 385.463,
- 1653.15,
- 541.211,
- 56.6718,
- 276.195,
- 707.914,
- 681.038,
- 733.694,
- 700.282,
- 173.163,
- 2230.31,
- 136.351
+ 647.441,
+ 1258.13,
+ 378.249,
+ 1228.53,
+ 338.526,
+ 439.055,
+ 444.451,
+ 68.299,
+ 1618.83,
+ 242.058,
+ 929.339,
+ 595.678,
+ 754.718,
+ 687.373,
+ 454.652,
+ 379.081,
+ 542.258,
+ 290.313,
+ 1.10705,
+ 744.15,
+ 1703.8,
+ 1010.19,
+ 2870.46,
+ 856.075,
+ 1864.98,
+ 296.741,
+ 874.288,
+ 33.2233,
+ 23.4443,
+ 37.7505,
+ 37.4312,
+ 100.661,
+ 24.1279,
+ 7.88562,
+ 26.9518,
+ 11.868,
+ 89.8789,
+ 28.9366,
+ 32.0384,
+ 12.1687,
+ 29.4468,
+ 4.91973,
+ 565.232,
+ 471.598,
+ 1808.9,
+ 1470.56,
+ 738.43,
+ 1449.61,
+ 1967.42,
+ 1666.8,
+ 577.982,
+ 17.9786,
+ 3633.48,
+ 5592.56,
+ 3319.49,
+ 257.515,
+ 239.116,
+ 312.322,
+ 1701.52,
+ 563.862,
+ 9.27568,
+ 114.437,
+ 17.2737,
+ 5.63337,
+ 5.32805,
+ 24.809,
+ 1098.33,
+ 25.1201,
+ 96.0282,
+ 3315.12,
+ 242.887,
+ 36.0287,
+ 32.4886,
+ 4.45102,
+ 31.7885,
+ 20.2812,
+ 22.29
};
// Assignment to powers vector.
@@ -749,58 +1210,79 @@ static int l_gPower[] = {
1, 1, 1, 1, 1, 2,
2, 1, 1, 1, 1, 1,
1, 1, 1, 1, 2, 1,
+ 1, 2, 1, 1, 1, 1,
1, 1, 1, 2, 1, 1,
1, 1, 2, 1, 1, 1,
1, 1, 1, 1, 1, 3,
1, 1, 1, 1, 2, 2,
+ 1, 3, 1, 1, 1, 1,
2, 1, 2, 1, 1, 1,
+ 1, 1, 1, 2, 1, 2,
+ 1, 2, 2, 1, 1, 1,
1, 1, 2, 1, 1, 2,
+ 1, 2, 1, 1, 1, 2,
1, 1, 3, 1, 1, 1,
2, 1, 1, 1, 1, 2,
3, 1, 1, 1, 1, 1,
1, 1, 1, 1, 3, 1,
2, 1, 1, 2, 1, 1,
1, 1, 1, 2, 2, 1,
+ 1, 2, 1, 2, 1, 1,
1, 1, 2, 1, 2, 1,
- 1, 2, 1, 1, 2, 1,
+ 1, 1, 2, 2, 1, 1,
2, 1, 1, 1, 2, 1,
- 2, 2, 1, 1, 1, 1,
1, 1, 1, 3, 1, 1,
1, 1, 1, 1, 2, 3,
2, 1, 1, 3, 1, 1,
- 3, 1, 1, 1, 2, 1,
+ 1, 2, 1, 3, 1, 1,
+ 1, 1, 2, 3, 1, 1,
+ 1, 1, 3, 2, 1, 1,
+ 1, 2, 2, 2, 1, 1,
2, 2, 1, 1, 2, 1,
1, 3, 1, 1, 2, 1,
- 1, 3, 1, 2, 1, 1,
1, 2, 2, 1, 2, 1,
1, 1, 3, 1, 2, 1,
2, 2, 1, 2, 1, 1,
+ 2, 1, 1, 2, 2, 1,
+ 1, 2, 1, 2, 2, 1,
1, 1, 1, 3, 2, 1,
+ 3, 1, 1, 2, 1, 1,
2, 1, 1, 1, 3, 1,
- 1, 2, 1, 1, 3, 1,
- 3, 2, 1, 1, 1, 1,
- 1, 2, 3, 1, 1, 1,
- 2, 1, 3, 1, 1, 1,
- 2, 2, 1, 1, 1, 2,
+ 3, 1, 1, 1, 1, 2,
1, 3, 1, 1, 1, 2,
- 1, 3, 2, 1, 1, 1,
2, 1, 2, 1, 1, 2,
1, 2, 2, 1, 1, 2,
1, 1, 3, 1, 1, 2,
- 2, 2, 2, 1, 1, 1,
+ 2, 1, 1, 2, 1, 2,
1, 2, 1, 2, 1, 2,
+ 1, 1, 2, 2, 1, 2,
+ 1, 1, 1, 3, 1, 2,
3, 1, 2, 1, 1, 1,
+ 2, 1, 1, 1, 2, 2,
+ 1, 1, 2, 1, 2, 2,
+ 1, 1, 1, 2, 2, 2,
+ 1, 2, 1, 1, 1, 3,
+ 1, 1, 2, 1, 1, 3,
+ 1, 1, 1, 2, 1, 3,
+ 1, 2, 1, 1, 2, 1,
+ 2, 2, 1, 1, 1, 1,
+ 3, 1, 1, 1, 2, 1,
+ 2, 1, 2, 2, 1, 1,
+ 1, 3, 1, 2, 1, 1,
+ 1, 1, 2, 1, 3, 1,
+ 1, 1, 1, 2, 3, 1,
+ 3, 2, 1, 1, 1, 1,
+ 2, 2, 1, 1, 1, 2,
+ 1, 3, 2, 1, 1, 1,
+ 2, 2, 2, 1, 1, 1,
1, 2, 1, 1, 2, 2,
- 1, 1, 1, 1, 3, 2,
- 2, 3, 1, 1, 1, 1,
2, 1, 1, 1, 1, 3,
- 1, 2, 1, 1, 1, 1,
- 1, 2, 2, 1, 1, 1,
- 1, 2, 1, 3, 1, 1,
2, 1, 2, 1, 2, 1,
- 1, 1, 2, 1, 3, 1,
- 2, 1, 1, 1, 2, 2,
- 1, 2, 1, 1, 1, 3
+ 1, 1, 2, 2, 2, 1,
+ 1, 2, 1, 1, 3, 1,
+ 1, 2, 3, 1, 1, 1,
+ 2, 1, 3, 1, 1, 1,
+ 2, 3, 1, 1, 1, 1
};
//
@@ -839,12 +1321,11 @@ double l_fit(double *x) {
return returnValue;
}
-
-// EOF for ../Geant4/samurai/simtree_length.C
+// EOF for simtree_length.C
// -*- mode: c++ -*-
//
-// File ../Geant4/samurai/simtree_tof.C generated by TMultiDimFit::MakeRealCode
-// on Tue Jun 9 23:41:17 2020
+// File simtree_tof.C generated by TMultiDimFit::MakeRealCode
+// on Fri May 28 00:25:23 2021
// ROOT version 5.32/00
//
// This file contains the function
@@ -860,138 +1341,180 @@ double l_fit(double *x) {
// Static data variables
//
static int t_gNVariables = 6;
-static int t_gNCoefficients = 56;
-static double t_gDMean = 63.5106;
+static int t_gNCoefficients = 77;
+static double t_gDMean = 88.0547;
// Assignment to mean vector.
static double t_gXMean[] = {
- 0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+ 0.00503955, -0.0144156, -0.101042, -0.00374735, 0.252452, -0.512868 };
// Assignment to minimum vector.
static double t_gXMin[] = {
- -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+ -64.3965, -52.0748, -0.0335722, -0.0401226, -1065.23, -299.807 };
// Assignment to maximum vector.
static double t_gXMax[] = {
- 329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+ 52.3459, 55.8258, 0.0432727, 0.04243, 1147.99, 931.885 };
// Assignment to coefficients vector.
static double t_gCoefficient[] = {
- 296.575,
- 509.773,
- -1298.15,
- -177.102,
- 38.3448,
- 495.622,
- -7.92753,
- -88.3852,
- -44.7681,
- -204.021,
- 916.644,
- 158.841,
- -954.907,
- 54.8928,
- -71.5493,
- -1160.79,
- -151.56,
- -6.48142,
- -20.8229,
- -95.551,
- 21.8388,
- 4.53363,
- -368.211,
- -11.1978,
- -0.122623,
- 0.899953,
- 0.664196,
- -4.26981,
- -2.76381,
- -353.964,
- -41.6085,
- -0.104552,
- 28.9895,
- 2.23953,
- 8.78264,
- 123.105,
- 3.83715,
- -338.603,
- -92.1278,
- -1248.05,
- -432.176,
- 299.722,
- 10.2068,
- 17.7278,
- -20.2405,
- -135.529,
- 644.109,
- -168.691,
- -0.0126702,
- -23.4315,
- -633.876,
- -171.644,
- -48.5535,
- 7.73714,
- -3.93291,
- 0.0801397
+ 282.419,
+ 600.69,
+ 199.214,
+ -610.381,
+ 16.8965,
+ -181.103,
+ -65.4683,
+ 27.7219,
+ -586.193,
+ -221.165,
+ 131.422,
+ 257.64,
+ 1.87553,
+ -199.501,
+ -2326.17,
+ 1262.92,
+ 88.6015,
+ -380.386,
+ -2.50251,
+ 0.192506,
+ -0.0682083,
+ -0.0605175,
+ -0.00933975,
+ 2.54256,
+ 2.13194,
+ 1.70986,
+ -0.734912,
+ 0.725959,
+ 0.566525,
+ 1.41423,
+ 262.324,
+ 54.1268,
+ -1143.9,
+ 255.544,
+ -5.4915,
+ -383.569,
+ 158.012,
+ 177.209,
+ 2458.18,
+ 2069.51,
+ -4536.21,
+ -1851.8,
+ 144.418,
+ 21.6852,
+ -147.382,
+ 49.0783,
+ -791.882,
+ 13.6185,
+ -269.147,
+ 75.7361,
+ 134.56,
+ -230.177,
+ -950.155,
+ 1062.38,
+ 81.3635,
+ -0.315919,
+ -1.87225,
+ -0.378093,
+ 0.88221,
+ -0.437349,
+ -0.563023,
+ 1.42746,
+ -525.101,
+ -447.633,
+ 0.297019,
+ 20.3115,
+ 195.852,
+ 0.424978,
+ -0.802113,
+ 0.778162,
+ 381.743,
+ -0.338946,
+ -0.109268,
+ -0.201413,
+ -0.0377476,
+ -0.150836,
+ 0.0571845
};
// Assignment to error coefficients vector.
static double t_gCoefficientRMS[] = {
- 90.0633,
- 153.552,
- 286.58,
- 57.4084,
- 10.6353,
- 98.8219,
- 12.2663,
- 110.128,
- 18.1831,
- 52.623,
- 177.01,
- 44.7459,
- 220.758,
- 29.2743,
- 22.5181,
- 280.17,
- 36.8213,
- 66.015,
- 23.1365,
- 90.5589,
- 10.4606,
- 23.5416,
- 99.1486,
- 50.5072,
- 1.90998,
- 3.21059,
- 1.31771,
- 31.4861,
- 9.39573,
- 67.4273,
- 33.8813,
- 0.497038,
- 12.9849,
- 4.9474,
- 13.3374,
- 41.0747,
- 11.7911,
- 86.8346,
- 27.0137,
- 280.477,
- 98.0117,
- 82.1288,
- 53.8616,
- 98.7791,
- 31.4828,
- 165.676,
- 125.204,
- 40.4042,
- 4.27406,
- 11.177,
- 143.484,
- 55.4336,
- 52.6846,
- 15.5019,
- 22.798,
- 0.581552
+ 51.6658,
+ 100.306,
+ 31.0942,
+ 100.473,
+ 36.5499,
+ 36.1352,
+ 5.48829,
+ 19.3199,
+ 74.8745,
+ 55.5697,
+ 31.1116,
+ 44.4306,
+ 0.110749,
+ 81.8808,
+ 231.073,
+ 150.427,
+ 24.4632,
+ 71.4537,
+ 2.1768,
+ 3.60448,
+ 0.792011,
+ 1.63035,
+ 3.6143,
+ 8.62908,
+ 3.21448,
+ 1.21207,
+ 0.490659,
+ 0.55933,
+ 0.532584,
+ 2.39502,
+ 45.341,
+ 37.6402,
+ 145.644,
+ 60.5957,
+ 9.27703,
+ 159.431,
+ 133.982,
+ 47.623,
+ 293.067,
+ 271.101,
+ 450.201,
+ 272.399,
+ 19.7108,
+ 21.145,
+ 19.4434,
+ 28.0278,
+ 132.355,
+ 49.4246,
+ 60.5858,
+ 37.5262,
+ 23.2998,
+ 58.9152,
+ 139.819,
+ 139.147,
+ 68.8468,
+ 3.6068,
+ 9.68357,
+ 10.9918,
+ 2.70767,
+ 1.19191,
+ 0.444182,
+ 3.0704,
+ 118.051,
+ 114.732,
+ 1.74076,
+ 25.8244,
+ 43.6475,
+ 3.15977,
+ 2.42206,
+ 2.89927,
+ 84.993,
+ 3.25617,
+ 1.97584,
+ 2.41757,
+ 0.930273,
+ 2.8245,
+ 2.14083
};
// Assignment to powers vector.
@@ -1005,55 +1528,76 @@ static int t_gPower[] = {
1, 1, 1, 2, 1, 1,
1, 1, 2, 1, 1, 1,
1, 1, 1, 1, 1, 3,
- 1, 1, 1, 1, 2, 2,
1, 3, 1, 1, 1, 1,
- 1, 2, 2, 1, 1, 1,
+ 2, 1, 2, 1, 1, 1,
1, 1, 2, 1, 1, 2,
1, 1, 3, 1, 1, 1,
2, 1, 1, 1, 1, 2,
- 3, 1, 1, 1, 1, 1,
1, 1, 1, 1, 3, 1,
- 2, 1, 1, 2, 1, 1,
- 1, 1, 1, 2, 2, 1,
1, 2, 1, 2, 1, 1,
1, 1, 2, 1, 2, 1,
2, 1, 1, 1, 2, 1,
1, 1, 1, 3, 1, 1,
1, 1, 1, 1, 2, 3,
- 2, 1, 1, 3, 1, 1,
1, 2, 1, 3, 1, 1,
- 3, 1, 1, 1, 2, 1,
- 2, 2, 1, 1, 2, 1,
+ 1, 1, 3, 2, 1, 1,
1, 3, 1, 1, 2, 1,
- 1, 2, 2, 1, 2, 1,
- 1, 1, 3, 1, 2, 1,
+ 1, 3, 1, 2, 1, 1,
+ 2, 1, 2, 1, 2, 1,
+ 2, 2, 1, 2, 1, 1,
1, 2, 1, 2, 2, 1,
1, 1, 1, 3, 2, 1,
- 1, 2, 1, 1, 3, 1,
+ 2, 1, 1, 1, 3, 1,
1, 1, 2, 1, 3, 1,
+ 1, 1, 1, 2, 3, 1,
3, 2, 1, 1, 1, 1,
- 1, 2, 3, 1, 1, 1,
3, 1, 1, 1, 1, 2,
- 2, 1, 3, 1, 1, 1,
- 2, 2, 1, 1, 1, 2,
1, 3, 1, 1, 1, 2,
2, 1, 2, 1, 1, 2,
- 1, 2, 2, 1, 1, 2,
1, 1, 3, 1, 1, 2,
2, 2, 2, 1, 1, 1,
1, 2, 1, 2, 1, 2,
- 3, 1, 2, 1, 1, 1,
+ 1, 1, 2, 2, 1, 2,
+ 1, 1, 1, 3, 1, 2,
2, 1, 1, 1, 2, 2,
1, 2, 1, 1, 2, 2,
- 1, 1, 1, 1, 3, 2,
- 2, 3, 1, 1, 1, 1,
+ 1, 1, 2, 1, 2, 2,
+ 1, 1, 1, 2, 2, 2,
+ 2, 1, 1, 1, 1, 3,
1, 2, 1, 1, 1, 3,
- 2, 1, 2, 1, 1, 1,
+ 1, 1, 2, 1, 1, 3,
+ 1, 2, 1, 1, 1, 1,
+ 1, 1, 1, 1, 2, 2,
+ 1, 1, 1, 2, 1, 2,
+ 1, 2, 2, 1, 1, 1,
+ 1, 2, 1, 1, 1, 2,
+ 3, 1, 1, 1, 1, 1,
+ 2, 1, 1, 2, 1, 1,
+ 1, 1, 1, 2, 2, 1,
+ 1, 2, 1, 1, 2, 1,
+ 1, 1, 2, 2, 1, 1,
+ 1, 1, 2, 3, 1, 1,
+ 1, 2, 2, 2, 1, 1,
+ 2, 1, 2, 2, 1, 1,
+ 1, 2, 2, 1, 2, 1,
+ 1, 1, 3, 1, 2, 1,
+ 1, 2, 1, 1, 3, 1,
+ 1, 2, 3, 1, 1, 1,
+ 1, 2, 2, 1, 1, 2,
+ 2, 1, 1, 2, 1, 2,
+ 3, 1, 2, 1, 1, 1,
+ 1, 1, 1, 2, 1, 3,
2, 2, 1, 1, 1, 1,
- 2, 1, 2, 1, 2, 1,
+ 2, 1, 1, 3, 1, 1,
+ 2, 2, 1, 1, 2, 1,
+ 2, 1, 1, 2, 2, 1,
+ 2, 2, 1, 1, 1, 2,
+ 1, 1, 2, 2, 2, 1,
+ 2, 1, 3, 1, 1, 1,
1, 3, 2, 1, 1, 1,
- 2, 1, 1, 1, 1, 3,
- 2, 1, 1, 1, 3, 1
+ 3, 1, 1, 1, 2, 1,
+ 3, 1, 1, 2, 1, 1,
+ 2, 3, 1, 1, 1, 1
};
//
@@ -1092,6 +1636,5 @@ double t_fit(double *x) {
return returnValue;
}
-// EOF for ../Geant4/samurai/simtree_tof.C
-
+// EOF for simtree_tof.C
diff --git a/Projects/S034/macro/testB.cxx b/Projects/S034/macro/testB.cxx
index 61f2f6df5a5884517bd60f3f00acac1f778e900c..78abab8dcc053d93e3086d101dbea9a4e53c8fdf 100644
--- a/Projects/S034/macro/testB.cxx
+++ b/Projects/S034/macro/testB.cxx
@@ -2,7 +2,7 @@
SamuraiFieldMap field;
void DrawT(TVector3 pos, TVector3 dir, double Brho){
-std::vector< TVector3 > track = field.Propagate(3000,Brho,pos,dir);
+std::vector< TVector3 > track = field.Propagate(Brho,pos,dir);
auto g = new TGraph();
unsigned int size = track.size();
g->Set(size);
@@ -40,9 +40,11 @@ void testB(){
DrawT(TVector3(0,0,-3500),TVector3(0,0,1),5.48);
DrawT(TVector3(0,0,-3500),TVector3(0,0,1),3.62);
+ DrawT(TVector3(),TVector3(0,0,1),3.62);
+
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);
+ std::vector< TVector3 > track = field.Propagate(b,p,d);
cout << field.FindBrho(p,d,track.back(),d)<< endl;;
double rFDC2 = 3686.77 + 880.745/2.;
diff --git a/Projects/S034/macro/testB2.cxx b/Projects/S034/macro/testB2.cxx
index 5265ce9e85e2e99bf35fface0345a48c20ac4ee0..ca9b10532ff4f624554365a175d43ace894b8f7a 100644
--- a/Projects/S034/macro/testB2.cxx
+++ b/Projects/S034/macro/testB2.cxx
@@ -2,17 +2,30 @@ void testB2(){
SamuraiFieldMap field;
double angle = 30*deg;
+ double fdc2angle = (59.930-90.0)*deg;
field.LoadMap(angle,"field_map/180702-2,40T-3000.table.bin",10);
- field.SetFDC2Angle((59.930-90.0)*deg);
+ field.SetFDC2Angle(fdc2angle);
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);
+ //auto scan = field.BrhoScan(2.5,10,0.1);
+// scan->Draw();
+ //auto h=new TH1D("dB","dB",1000,-2e-1,2e-1);
+ // new TCanvas();
+ auto h=new TH1D("dB","dB",1000,-0.1,0.1);
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);
+ for(unsigned int i = 0 ; i < 100 ; 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);
+ //
+ //TVector3 p(0,0,-3500);
+ //TVector3 d(0,0,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);
+ std::vector< TVector3 > track = field.Propagate(b,p,d,false);
+ // rotate from lab angle to FDC2 frame
+ track.back().RotateY(-fdc2angle+angle);
+ double br=field.FindBrho(p,d,track.back(),d);
+ double pc = 100*(br-b)/b;
+ cout <<" " << 100*(br-b)/b << " " << br << " " << b << endl;
+ h->Fill(pc);
}
h->Draw();
diff --git a/Projects/S034/s034.detector b/Projects/S034/s034.detector
index 6d99cb8d8220a1dd9233417b483fb2f34c2c22e0..f0c1a6b7d82a2e99446b87da9344154b19482b9e 100644
--- a/Projects/S034/s034.detector
+++ b/Projects/S034/s034.detector
@@ -1,24 +1,26 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIBDC 1
- XML= db/SAMURAIBDC1.xml
- Offset= 0 0 0 mm
- InvertX= 0
- InvertY= 0
- InvertD= 1
+ XML= db/SAMURAIBDC1.xml
+ Offset= 0.959479 -0.17903 -6876.34 mm
+ InvertX= 0
+ InvertY= 0
+ InvertD= 1
%%%%%%%%%%%%
SAMURAIBDC 2
- XML= db/SAMURAIBDC2.xml
- Offset= 0 0 0 mm
- InvertX= 0
- InvertY= 0
- IverteD= 1
+ XML= db/SAMURAIBDC2.xml
+% Offset= +0.25 +1.6 -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 0 mm
- ZRotation= 35 deg
- TargetLength= 152.76 mm
+ Position= 0.345399 1.02061 -4650 mm
+ ZRotation= 40.6 deg
+ TargetLength= 151.72 mm
TargetMaterial= LH2
CellMaterial= Mylar
TPCOnly= 0
@@ -31,7 +33,7 @@ Minos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIFDC0
XML= db/SAMURAIFDC0_20200109.xml
- Offset= 3.005 1.864 0 mm
+ Offset= -0.00666226 0.102191 -3370.01 mm
InvertX= 1
InvertY= 0
InvertD= 1
@@ -39,7 +41,7 @@ SAMURAIFDC0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIFDC2
XML= db/SAMURAIFDC2.xml
- Offset= 0 0 0 mm
+ Offset= -252.416 -0.228477 4122.57 mm
InvertX= 0
InvertY= 1
InvertD= 1
diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index 60fa8381bdb67c3c232c67f9e4a466fdef6f0fd4..302b3276c051d2626c9e7f431ba81adf0b635f97 100755
--- a/Projects/e793s/Analysis.cxx
+++ b/Projects/e793s/Analysis.cxx
@@ -138,7 +138,7 @@ void Analysis::TreatEvent(){
}
}
- TVector3 BeamDirection(XBeam,YBeam,1);
+ TVector3 BeamDirection(0.,0.,1.);
BeamImpact = TVector3(XBeam,YBeam,m_DetectorManager->GetTargetZ());
ParticleMult=M2->Si_E.size()+MG->DSSD_E.size();
diff --git a/Projects/e793s/macro/BeamSpot/BeamSpot.C b/Projects/e793s/macro/BeamSpot/BeamSpot.C
index 7d9eaa2cc525877cc534cbcc47cc4d391d9abb37..8bf8a71c3a0bc3740d2d779760a1ee6933b0e889 100755
--- a/Projects/e793s/macro/BeamSpot/BeamSpot.C
+++ b/Projects/e793s/macro/BeamSpot/BeamSpot.C
@@ -8,6 +8,8 @@
#include <fstream>
#include <vector>
#include <cmath>
+#include <string>
+#include <sstream>
//Root
#include <TVector3.h>
//NPTool
@@ -30,23 +32,23 @@ void BeamSpot(){
vector <double> Xd, Yd, Zd; //Vector of particle direction. Calculated as Xp-Xb, Yp-Yb...
ifstream MugastDataFile;
double ThetaNormalTarget;
-
+ TVector3 beamDir{ 0.0, 0.0, 1.0 };
gErrorIgnoreLevel = kWarning; // Suppress ".pdf created" lines
/*** ITERATIVE GRID CONTROLS ***/
/***** pos varied as offset ****/
- /**/ double xmin = +0.020; /**/
- /**/ double xmax = +0.080; /**/
- /**/ unsigned int xdiv = 12; /**/
+ /**/ double xmin = -8.000; /**/
+ /**/ double xmax = +2.000; /**/
+ /**/ unsigned int xdiv = 10; /**/
/**/ /**/
- /**/ double ymin = -0.050; /**/
- /**/ double ymax = +0.050; /**/
+ /**/ double ymin = -5.000; /**/
+ /**/ double ymax = +5.000; /**/
/**/ unsigned int ydiv = 10; /**/
/**/ /**/
- /**/ double zmin = -0.050; /**/
- /**/ double zmax = +0.050; /**/
- /**/ unsigned int zdiv = 2; /**/
+ /**/ double zmin = -5.000; /**/
+ /**/ double zmax = +5.000; /**/
+ /**/ unsigned int zdiv = 10; /**/
/**/ /**/
/***** thick varied as %ge *****/
/**/ unsigned int tmin = 7; /**/
@@ -56,19 +58,48 @@ void BeamSpot(){
/******* METRIC CONTROLS *******/
/**/ double peakE = 0.143; /**/
- /**/ double sigMultip = 0.1; /**/
+ /**/ double sigMultip = 0.05; /**/
/*******************************/
// File name controls
- const char* XYZE_file = "XYZE_gammaGated_Full_TestThetaNormalTarget.txt";
- const char* outputMetric = "output_Run63_metrics_ThetaNormal.txt";
- const char* outputHisto = "output_Run63_histograms_ThetaNormal.root";
+ string tag;
+ string addtag;
+ cout << "===================================================" << endl;
+ cout << " Enter the version tag " << endl;
+ cout << " (i.e. Run63_BeamDirectionCorrection)" << endl;
+ cout << " - - - - - - - - - - - - - - - - - - - - - - - - - " << endl;
+ getline(cin,tag);
+ cout << " - - - - - - - - - - - - - - - - - - - - - - - - - " << endl;
+ cout << " Additional output tagging? Coarse, fine? " << endl;
+ cout << " - - - - - - - - - - - - - - - - - - - - - - - - - " << endl;
+ getline(cin,addtag);
+
+ string XYZE_string = "XYZE_";
+ XYZE_string.append(tag);
+ XYZE_string.append(".txt");
+
+ string out_string = "output_";
+ out_string.append(tag);
+ out_string.append("_");
+ out_string.append(addtag);
+ string met_string = out_string;
+ met_string.append("_metrics.txt");
+ string hst_string = out_string;
+ hst_string.append("_histograms.root");
+
+ const char* XYZE_file = XYZE_string.c_str();
+ const char* outputMetric = met_string.c_str();
+ const char* outputHisto = hst_string.c_str();
// Calculate size of iteratve steps
double xstp = (xmax-xmin)/ ((double) xdiv);
double ystp = (ymax-ymin)/ ((double) ydiv);
double zstp = (zmax-zmin)/ ((double) zdiv);
- cout << "Xstp = " << xstp << " Ystp = " << ystp << " Zstp = " << zstp << endl;
+ cout << "===================================================" << endl;
+ cout << " X: " << xmin << " to " << xmax << " in steps of " << xstp << endl;
+ cout << " Y: " << ymin << " to " << ymax << " in steps of " << ystp << endl;
+ cout << " Z: " << zmin << " to " << zmax << " in steps of " << zstp << endl;
+ cout << "===================================================" << endl;
// Vectors of the normal for each detector. UPDATE WITH NEW POSITIONS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
TVector3 MugastNormal1{ -0.453915, +0.455463, -0.765842};
@@ -157,7 +188,6 @@ void BeamSpot(){
ThetaNormalTarget = tempTheta = ELab = Ex = 0.0;
-
switch(DetNum[i]){
case 1:
tempTheta = particleDir.Angle(MugastNormal1);
@@ -182,9 +212,9 @@ void BeamSpot(){
return; // Exit code
}
- // Change beam spot vector to inverse beam direction vector
- beamSpot.SetZ(-1.0);
- ThetaNormalTarget = particleDir.Angle(beamSpot);
+ // Change beamDir vector to inverse beam direction
+ beamDir.SetZ(-1.0);
+ ThetaNormalTarget = particleDir.Angle(beamDir);
//micrometer defined in NPSystemOfUnits.h
ELab = LightAl.EvaluateInitialEnergy(
@@ -196,9 +226,9 @@ void BeamSpot(){
0.5*TargetThickness, //pass through half target
ThetaNormalTarget); //angle leaving target
- // Change beam spot vector to beam direction vector
- beamSpot.SetZ(1.0);
- Ex = reaction.ReconstructRelativistic(ELab, particleDir.Angle(beamSpot));
+ // Change beamDir vector back to forward direction
+ beamDir.SetZ(1.0);
+ Ex = reaction.ReconstructRelativistic(ELab, particleDir.Angle(beamDir));
// Fill Ex histograms
tempHist->Fill(Ex);
diff --git a/Projects/e793s/macro/minimizer.cxx b/Projects/e793s/macro/minimizer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6fb73246063b27886d28852a0eb0d42d02da7843
--- /dev/null
+++ b/Projects/e793s/macro/minimizer.cxx
@@ -0,0 +1,92 @@
+
+double ref = 0.143; // the energy of the selected states
+vector<TVector3*> pos;
+vector<double> energy;
+NPL::EnergyLoss CD2("proton_CD2.G4table","G4Table",100);
+NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
+
+////////////////////////////////////////////////////////////////////////////////
+double devE(const double* parameter){
+
+ // My reaction at 8A MeV*47A=376 MeV
+ static NPL::Reaction reaction("47K(d,p)48K@376");
+ // Return the standard deviation in Brho
+ unsigned int size = pos1.size();
+
+ TVector3 offset(parameter[0],parameter[1],parameter[2]);
+
+ double dE,Theta;
+ TVector3 dir;
+ static auto h = new TH1D("h","h", 1000,-0,100);
+ h->Reset();
+ for(unsigned int i = 0 ; i < size ; i++){
+ dir=*(pos[i])-o1;
+ double Theta= dir.Angle(TVector3(0,0,1));
+ double Energy = energy[i];
+ // do some energy loss correction here:
+ // This need to take parameter[4]*0.5*micrometer as the target thickness
+ Energy=CD2.EvaluateInitialEnergy(Energy,0,5*parameter[4]*micrometer,Theta);
+
+ double Ex=reaction->ReconstructRelativistic(Energy,Theta);
+ //cout << (pB-pM).Mag2() << " " << dR << endl;
+ h->Fill(Ex-ref);
+ }
+ cout << h->GetMean() << " " << h->GetStdDev() << " " << parameter[4] << endl;
+ h->Draw();
+ gPad->Update();
+ // adapt the metric as needed
+ return (h->GetMean()+h->GetStdDev() );
+}
+////////////////////////////////////////////////////////////////////////////////
+void LoadFile(){
+
+ ifstream file("myFile.txt");
+ if(!file.is_open()){
+ cout << "fail to load file" << endl;
+ exit(1);
+ }
+ else {
+ cout << "Success opening file" << endl;
+ }
+
+ double x,y,z,e;
+
+ while(file >> x >> y >> z >> e ){
+ auto p1 = new TVector3(x,y,z);
+ double e;
+ pos.push_back(p1);
+ energy.push_back(e);
+ }
+ file.close();
+}
+////////////////////////////////////////////////////////////////////////////////
+void BDC(){
+
+ // Data
+ LoadFile();
+
+ // Minimizer
+ auto min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ // Function with 4 parameter XYZ and Target thickness
+ auto func=ROOT::Math::Functor(&devE,4);
+ min->SetFunction(func);
+ min->SetPrintLevel(0);
+ min->SetPrecision(1e-10);
+ // Start value: Center beam (0,0,0) and 4.7um target 0.5mg/cm2 target
+ double parameter[6]={0,0,0,4.7};
+ devR(parameter);
+
+ // min->SetLimitedVariable(0,"AM",parameter[0],1,-90,90);
+ min->SetLimitedVariable(0,"X",parameter[0],0.1,-10,10);
+ min->SetLimitedVariable(1,"Y",parameter[1],0.1,-10,10);
+ min->SetLimitedVariable(2,"Z",parameter[2],0.1,-5,5);
+ min->SetLimitedVariable(3,"T",parameter[3],0.1,parameter[2]-0.1*parameter[2],parameter[2]+0.1*parameter[2]);
+ min->Minimize();
+
+ const double* x = min->X();
+ cout << "X =" << x[0]<<endl;
+ cout << "Y =" << x[1]<<endl;
+ cout << "Z =" << x[2]<<endl;
+ cout << "Y =" << x[3]<<endl;
+ cout << "Minimum: " << devR(x) << endl;
+}