From 4e9e4ac9cc376bff0d506e07bad9ac811f499c6c Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Fri, 9 Oct 2020 18:24:05 +0200
Subject: [PATCH] *Progress on Samurai analysis
---
NPLib/Core/NPCalibrationManager.cxx | 22 ++++
NPLib/Core/NPCalibrationManager.h | 1 +
.../Detectors/Samurai/TSamuraiFDC2Physics.cxx | 105 ++++++++++++++----
NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h | 6 +
4 files changed, 114 insertions(+), 20 deletions(-)
diff --git a/NPLib/Core/NPCalibrationManager.cxx b/NPLib/Core/NPCalibrationManager.cxx
index 45864d6e0..7f25d9777 100644
--- a/NPLib/Core/NPCalibrationManager.cxx
+++ b/NPLib/Core/NPCalibrationManager.cxx
@@ -356,6 +356,28 @@ bool CalibrationManager::ApplyThreshold(const std::string& ParameterPath, const
return false;
}
+/////////////////////////////////////////////////////////////////////////////////////////////
+double CalibrationManager::ApplySigmoid(const std::string& ParameterPath, const double& RawValue) const{
+ std::map< std::string , std::vector<double> >::const_iterator it ;
+ static std::map< std::string , std::vector<double> >::const_iterator ite = fCalibrationCoeff.end();
+ // Find the good parameter in the Map
+ // Using Find method of stl is the fastest way
+ it = fCalibrationCoeff.find(ParameterPath) ;
+ // If the find methods return the end iterator it's mean the parameter was not found
+ if(it == ite ){
+ return RawValue ;
+ }
+
+ std::vector<double> Coeff = it->second ;
+ // Check that the number of coeff is ok
+ if(Coeff.size()!=3){
+ return RawValue ;
+ }
+
+ return (Coeff[0]/(exp(Coeff[1]*(Coeff[2]-RawValue))+1));
+
+
+ }
/////////////////////////////////////////////////////////////////////////////////////////////
double CalibrationManager::GetPedestal(const std::string& ParameterPath) const{
return GetValue(ParameterPath,0);
diff --git a/NPLib/Core/NPCalibrationManager.h b/NPLib/Core/NPCalibrationManager.h
index b1981418f..cc44bcd8d 100644
--- a/NPLib/Core/NPCalibrationManager.h
+++ b/NPLib/Core/NPCalibrationManager.h
@@ -63,6 +63,7 @@ class CalibrationManager{
double ApplyResistivePositionCalibrationDebug (const std::string& ParameterPath , const double& RawValue) const ;
bool ApplyThreshold (const std::string& ParameterPath, const double& RawValue) const ;
+ double ApplySigmoid(const std::string& ParameterPath, const double& RawValue) const ;
double GetPedestal (const std::string& ParameterPath) const ;
double GetValue (const std::string& ParameterPath,const unsigned int& order) const ;
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
index a9200971b..3bfbe8da1 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
@@ -36,6 +36,10 @@
#include "NPSystemOfUnits.h"
// ROOT
#include "TChain.h"
+#include "TVector2.h"
+#include "Math/Minimizer.h"
+#include "Math/Functor.h"
+#include "Math/Factory.h"
using namespace NPUNITS;
///////////////////////////////////////////////////////////////////////////
@@ -46,7 +50,7 @@ ClassImp(TSamuraiFDC2Physics)
m_PreTreatedData = new TSamuraiFDC2Data ;
m_EventPhysics = this ;
//m_Spectra = NULL;
- ToTThreshold = 0;
+ ToTThreshold = 100;
}
///////////////////////////////////////////////////////////////////////////
@@ -58,16 +62,30 @@ void TSamuraiFDC2Physics::BuildSimplePhysicalEvent(){
void TSamuraiFDC2Physics::BuildPhysicalEvent(){
PreTreat();
RemoveNoise();
-
- // one map per detector
- map<unsigned int, vector<double> > X ;
- map<unsigned int, vector<double> > Z ;
- map<unsigned int, vector<double> > R ;
+ // Map[detector & plane angle, vector of spatial information]
+ map<std::pair<unsigned int,double>, vector<double> > X ;
+ map<std::pair<unsigned int,double>, vector<double> > Z ;
+ map<std::pair<unsigned int,double>, vector<double> > R ;
unsigned int size = Detector.size();
for(unsigned int i = 0 ; i < size ; i++){
-
-
+ if(Matched[i]){
+ int det = Detector[i];
+ int layer = Layer[i];
+ int wire = Wire[i];
+ SamuraiDCIndex idx(det,layer,wire);
+ std::pair<unsigned int, double> p(det,Wire_T[idx]);
+ X[p].push_back(Wire_X[idx]);
+ Z[p].push_back(Wire_Z[idx]);
+ R[p].push_back(DriftLength[i]);
+ }
+ }
+
+ // Reconstruct the vector for each of the plane of each of the detector
+ double dirX,dirZ,refX;
+ for(auto it = X.begin();it!=X.end();++it){
+ if(it->first.second==0)
+ Track2D(X[it->first],Z[it->first],R[it->first],dirX,dirZ,PosX);
}
return;
@@ -75,12 +93,49 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::Track2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& dirX, double& dirZ,double& refX ){
+ fitX=X;
+ fitZ=Z;
+ fitR=R;
+ // assume all X,Z,R of same size
+ unsigned int size = X.size();
+ // Define the starting point of the fit: a straight line passing through the
+ // the two first and last wire
+ // z = ax+b -> x=(z-b)/a
+ double a = fitZ[size-1]-fitZ[0]/(fitX[size-1]-fitX[0]);
+ double b = fitZ[0]-a*fitX[0];
+ double parameter[2]={a,b};
+ static ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ static ROOT::Math::Functor f(this,&TSamuraiFDC2Physics::SumD,2);
+ min->SetFunction(f);
+ min->SetVariable(0,"a",parameter[0], 1);
+ min->SetVariable(1,"b",parameter[1], 1);
- }
+ min->Minimize();
+ const double *xs = min->X();
+ refX=(-xs[1])/xs[0];
+}
+////////////////////////////////////////////////////////////////////////////////
+double TSamuraiFDC2Physics::SumD(const double* parameter ){
+ unsigned int size = fitX.size();
+ // Compute the sum P of the distance between the circle and the track
+ double P = 0;
+ for(unsigned int i = 0 ; i < size ; i++){
+ double x=(fitZ[i]-parameter[1])/parameter[0];
+ P+= sqrt(abs( (x-fitX[i])*(x-fitX[i])-fitR[i]*fitR[i]));
+ }
+ return P;
+ }
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::PreTreat(){
ClearPreTreatedData();
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+ static string channel;
+ // one map per detector
+ map<unsigned int, vector<double> > X ;
+ map<unsigned int, vector<double> > Z ;
+ map<unsigned int, vector<double> > R ;
+
unsigned int size = m_EventData->Mult();
for(unsigned int i = 0 ; i < size ; i++){
// EDGE=1 is the leading edge, IE, the real time.
@@ -114,7 +169,8 @@ void TSamuraiFDC2Physics::PreTreat(){
Wire.push_back(wire);
Time.push_back(time);
ToT.push_back(etime-time);
- SamuraiDCIndex idx(det,layer,wire);
+ channel="SamuraiFDC2/L" + NPL::itoa(layer);
+ DriftLength.push_back(Cal->ApplySigmoid(channel,etime));
}
}
@@ -151,18 +207,20 @@ void TSamuraiFDC2Physics::RemoveNoise(){
if(det==adet && blayer && abs(wire-awire)<=1){
match=true;
break;
- }
}
-
- if(match)
- Matched.push_back(true);
- else
- Matched.push_back(false);
}
+
+ if(match)
+ Matched.push_back(true);
+ else
+ Matched.push_back(false);
+ }
return;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::Clear(){
+ PosX=-10000;
+ DriftLength.clear();
Detector.clear();
Layer.clear();
Wire.clear();
@@ -219,10 +277,10 @@ void TSamuraiFDC2Physics::AddDC(string name, NPL::XmlParser& xml){
cout << "ERROR: Unknown layer orientation for Samurai FDC2"<< endl;
exit(1);
}
- SamuraiDCIndex idx(det,layer,wire);
- Wire_X[idx]=X;
- Wire_Z[idx]=Z;
- Wire_T[idx]=T;
+ SamuraiDCIndex idx(det,layer,wire);
+ Wire_X[idx]=X;
+ Wire_Z[idx]=Z;
+ Wire_T[idx]=T;
}
}
}
@@ -265,6 +323,13 @@ void TSamuraiFDC2Physics::WriteSpectra(){
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::AddParameterToCalibrationManager(){
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ // each layer
+ for( int l = 0 ; l < 14 ; ++l){
+ Cal->AddParameter("SamuraiFDC2", "L"+ NPL::itoa(l),"FDC2_L"+ NPL::itoa(l));
+ }
+
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
index 7a36606aa..46484cb74 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
@@ -87,6 +87,7 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
public:
// Provide Physical Multiplicity
+ vector<double> DriftLength;
vector<int> Detector;
vector<int> Layer;
vector<int> Wire;
@@ -97,6 +98,7 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
vector<TVector3> ParticleDirection;
vector<TVector3> MiddlePosition;
+ double PosX;
public:
// Projected position at given Z plan
TVector3 ProjectedPosition(double Z);
@@ -113,6 +115,10 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
// Construct the 2D track and ref position at Z=0 based on X,Z and Radius provided
void Track2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& dirX, double& dirZ,double& refX );
+ double SumD(const double* parameter );
+ vector<double> fitX;
+ vector<double> fitZ;
+ vector<double> fitR;
public: // Innherited from VDetector Class
--
GitLab