diff --git a/NPLib/Detectors/Catana/TCatanaPhysics.cxx b/NPLib/Detectors/Catana/TCatanaPhysics.cxx
index edb4a8f5681c9c4b77fca050a55cacad477dfefe..f3618f7ee564cf9e3d8101ade2011e28962c3446 100644
--- a/NPLib/Detectors/Catana/TCatanaPhysics.cxx
+++ b/NPLib/Detectors/Catana/TCatanaPhysics.cxx
@@ -37,6 +37,7 @@ using namespace std;
#include "NPOptionManager.h"
#include "NPSystemOfUnits.h"
using namespace NPUNITS;
+#include "NPTrackingUtility.h"
// ROOT
#include "TChain.h"
@@ -46,14 +47,14 @@ ClassImp(TCatanaPhysics)
///////////////////////////////////////////////////////////////////////////
TCatanaPhysics::TCatanaPhysics()
- : m_EventData(new TCatanaData),
- m_PreTreatedData(new TCatanaData),
- m_EventPhysics(this),
- m_Spectra(0),
- m_E_RAW_Threshold(0), // adc channels
- m_E_Threshold(0), // MeV
- m_NumberOfDetectors(0) {
-}
+ : m_EventData(new TCatanaData),
+ m_PreTreatedData(new TCatanaData),
+ m_EventPhysics(this),
+ m_Spectra(0),
+ m_E_RAW_Threshold(0), // adc channels
+ m_E_Threshold(0), // MeV
+ m_NumberOfDetectors(0) {
+ }
///////////////////////////////////////////////////////////////////////////
@@ -65,13 +66,17 @@ void TCatanaPhysics::AddDetector(double X, double Y, double Z, double Theta, dou
m_Phi[ID]=Phi;
m_Type[ID]=Type;
}
-
+
///////////////////////////////////////////////////////////////////////////
void TCatanaPhysics::BuildSimplePhysicalEvent() {
BuildPhysicalEvent();
}
+///////////////////////////////////////////////////////////////////////////
+TVector3 TCatanaPhysics::GetPositionOfInteraction(int& i){
+ return m_Position[DetectorNumber[i]];
+}
///////////////////////////////////////////////////////////////////////////
void TCatanaPhysics::ReadCSV(string path){
std::ifstream csv(path);
@@ -86,11 +91,11 @@ void TCatanaPhysics::ReadCSV(string path){
// ignore first line
getline(csv,buffer);
while(csv >> ID >> buffer >> type >> buffer >> layer >> buffer >> X >> buffer >> Y >> buffer >> Z >> buffer >> Theta >> buffer >> Phi){
- if(type<6)
+ if(type<6)
AddDetector(X,Y,Z,Theta*deg,Phi*deg,ID,type);
- else{
- // ignore other type for which I don't have the geometry
- }
+ else{
+ // ignore other type for which I don't have the geometry
+ }
}
return;
@@ -114,6 +119,25 @@ void TCatanaPhysics::BuildPhysicalEvent() {
}
}
+///////////////////////////////////////////////////////////////////////////
+unsigned int TCatanaPhysics::FindClosestHitToLine(const TVector3& v1, const TVector3& v2,double& d){
+
+ d = 1e32;
+ unsigned result = 0;
+ unsigned int size = DetectorNumber.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ double current_d = NPL::MinimumDistancePointLine(v1,v2,m_Position[DetectorNumber[i]]) ;
+ if(current_d < d){
+ d=current_d;
+ result=i;
+ }
+ }
+
+ if(d==1e32)
+ d=-1000;
+
+ return result;
+}
///////////////////////////////////////////////////////////////////////////
void TCatanaPhysics::PreTreat() {
// This method typically applies thresholds and calibrations
@@ -243,7 +267,7 @@ void TCatanaPhysics::ReadConfiguration(NPL::InputParser parser){
exit(1);
}
}
-
+
// Type 1
blocks = parser.GetAllBlocksWithTokenAndValue("Catana","Detector");
if(NPOptionManager::getInstance()->GetVerboseLevel())
@@ -370,14 +394,14 @@ NPL::VDetector* TCatanaPhysics::Construct() {
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy_Catana{
- public:
- proxy_Catana(){
- NPL::DetectorFactory::getInstance()->AddToken("Catana","Catana");
- NPL::DetectorFactory::getInstance()->AddDetector("Catana",TCatanaPhysics::Construct);
- }
-};
+ class proxy_Catana{
+ public:
+ proxy_Catana(){
+ NPL::DetectorFactory::getInstance()->AddToken("Catana","Catana");
+ NPL::DetectorFactory::getInstance()->AddDetector("Catana",TCatanaPhysics::Construct);
+ }
+ };
-proxy_Catana p_Catana;
+ proxy_Catana p_Catana;
}
diff --git a/NPLib/Detectors/Catana/TCatanaPhysics.h b/NPLib/Detectors/Catana/TCatanaPhysics.h
index 9ad23779bd19029065560e329755fbb4caa32b7d..f9fc98d14c71bc3b39a9c9502031cf4d4c9742fd 100644
--- a/NPLib/Detectors/Catana/TCatanaPhysics.h
+++ b/NPLib/Detectors/Catana/TCatanaPhysics.h
@@ -69,7 +69,7 @@ class TCatanaPhysics : public TObject, public NPL::VDetector {
/// A usefull method to bundle all operation to add a detector
void AddDetector(double X, double Y, double Z, double Theta, double Phi, int ID, int Type);
void ReadCSV(string path);
-
+
// Position method and variable
public:
map<int,TVector3> m_Position;//!
@@ -78,6 +78,8 @@ class TCatanaPhysics : public TObject, public NPL::VDetector {
map<int,int> m_Type;//!
TVector3 m_Ref;//!
TVector3 GetPositionOfInteraction(int& i);//!
+ // Return index of the closest hit to line defined by v1 and v2
+ unsigned int FindClosestHitToLine(const TVector3& v1, const TVector3& v2, double& d);
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
diff --git a/NPLib/TrackReconstruction/NPTrackingUtility.h b/NPLib/TrackReconstruction/NPTrackingUtility.h
index 4eb1538666230303e9842d370f898e0fa8bbb7bd..ee7c30502399c3bbae7ac174869277243d923457 100644
--- a/NPLib/TrackReconstruction/NPTrackingUtility.h
+++ b/NPLib/TrackReconstruction/NPTrackingUtility.h
@@ -18,11 +18,12 @@
*****************************************************************************/
namespace NPL{
+ //////////////////////////////////////////////////////////////////////////////
// return the minimum distance between v and w defined respectively by points
// v1,v2 and w1 w1
// Also compute the best crossing position BestPosition, i.e. average position
// at the minimum distance.
- double MinimumDistance(const TVector3& v1,const TVector3& v2, const TVector3& w1, const TVector3& w2, TVector3& BestPosition, TVector3& delta){
+ double MinimumDistanceTwoLines(const TVector3& v1,const TVector3& v2, const TVector3& w1, const TVector3& w2, TVector3& BestPosition, TVector3& delta){
TVector3 v = v2-v1;
TVector3 w = w2-w1;
// Finding best position
@@ -36,6 +37,17 @@ namespace NPL{
delta = (v1+t*v-w1-s*w);
return d;
}
+
+ //////////////////////////////////////////////////////////////////////////////
+ // return the minimum distance between the line defines by v1,v2 and the point
+ // in space x
+ // demo is here: https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+ double MinimumDistancePointLine(const TVector3& v1, const TVector3& v2, const TVector3& x){
+ TVector3 w1 = x-v1;
+ TVector3 w2 = x-v2;
+ TVector3 w = w1.Cross(w2);
+ return w.Mag()/(v2-v1).Mag();
+ }
}
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index c58c88c9bb2493cc64564a82cc048bd6c1b04912..d22b19e089253998e9e38f22ba52b1e56a0a5447 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -45,6 +45,7 @@ void Analysis::Init(){
InitInputBranch();
Strasse = (TStrassePhysics*) m_DetectorManager -> GetDetector("Strasse");
+ Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana");
// reaction properties
myQFS = new NPL::QFS();
myQFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
@@ -57,7 +58,12 @@ void Analysis::Init(){
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
// EnergyLoss Tables
string BeamName = NPL::ChangeNameToG4Standard(myBeam->GetName());
- BeamTarget = NPL::EnergyLoss(BeamName+"_"+TargetMaterial+".G4table","G4Table",10000);
+ BeamTarget = NPL::EnergyLoss(BeamName+"_"+TargetMaterial+".G4table","G4Table",100);
+ protonTarget = NPL::EnergyLoss("proton_"+TargetMaterial+".G4table","G4Table",100);
+ protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100);
+ protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100);
+
+
}
////////////////////////////////////////////////////////////////////////////////
@@ -87,15 +93,24 @@ void Analysis::TreatEvent(){
// computing minimum distance of the two lines
TVector3 Vertex;
TVector3 delta;
- Distance = NPL::MinimumDistance(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex,delta);
+ Distance = NPL::MinimumDistanceTwoLines(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex,delta);
VertexX=Vertex.X();
VertexY=Vertex.Y();
VertexZ=Vertex.Z();
deltaX=delta.X();
deltaY=delta.Y();
deltaZ=delta.Z();
+
+ // Look for associated Catana event
+ double d1,d2;
+ unsigned int i1,i2;
+ i1 = Catana->FindClosestHitToLine(InnerPos1,OuterPos1,d1);
+ i2 = Catana->FindClosestHitToLine(InnerPos2,OuterPos2,d2);
+ if(i1!=i2){
+ E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]);
+ E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]);
+ }
}
-
//double thickness_before = 0;
//double EA_vertex = BeamTarget.Slow(InitialBeamEnergy,thickness_before,0);
@@ -110,6 +125,31 @@ void Analysis::TreatEvent(){
//Ex = TMath::Sqrt( EB*EB - PB.Mag2() ) - myQFS->GetNucleusB()->Mass();
}
+////////////////////////////////////////////////////////////////////////////////
+double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir,const double& Ecatana){
+ TVector3 Normal = TVector3(0,0,1);
+ Normal.SetPhi(dir.Phi());
+ double Theta = dir.Angle(Normal);
+ // Catana Al housing
+ double E = protonAl.EvaluateInitialEnergy(Ecatana,0.5*mm,Theta);
+ // Strasse Chamber
+ E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta);
+ // Outer Barrel
+ E = protonSi.EvaluateInitialEnergy(E,400*micrometer,Theta);
+ // Inner Barrel
+ E = protonSi.EvaluateInitialEnergy(E,200*micrometer,Theta);
+ // LH2 target
+ static TVector3 x1;
+ x1= x0+dir;
+ TVector3 T1(0,30,0);
+ TVector3 T2(0,30,1);
+ T1.SetPhi(dir.Phi());
+ T2.SetPhi(dir.Phi());
+ TVector3 Vertex,delta;
+ double d = NPL::MinimumDistanceTwoLines(x0,x1,T1,T2,Vertex,delta);
+ E = protonTarget.EvaluateInitialEnergy(E,d,Theta);
+ }
+
////////////////////////////////////////////////////////////////////////////////
TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj){
@@ -127,7 +167,8 @@ void Analysis::End(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
- RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("E1",&E1,"E1/D");
+ RootOutput::getInstance()->GetTree()->Branch("E2",&E2,"E2/D");
RootOutput::getInstance()->GetTree()->Branch("Theta12",&Theta12,"Theta12/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("VertexX",&VertexX,"VertexX/D");
@@ -153,7 +194,8 @@ void Analysis::InitInputBranch(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
Ex = -1000 ;
- ELab = -1000;
+ E1= -1000;
+ E2 = -1000;
Theta12 = -1000;
ThetaCM = -1000;
VertexX=-1000;
diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h
index 8092e2804ac7dddc34d16af9f902e64a7556b5b1..08494fc00f4f3de403a4d087a2de1b8ea8005531 100644
--- a/Projects/Strasse/Analysis.h
+++ b/Projects/Strasse/Analysis.h
@@ -27,6 +27,7 @@
#include"RootOutput.h"
#include"RootInput.h"
#include "TStrassePhysics.h"
+#include "TCatanaPhysics.h"
#include "TInitialConditions.h"
#include "TInteractionCoordinates.h"
#include "TReactionConditions.h"
@@ -52,7 +53,8 @@ class Analysis: public NPL::VAnalysis{
private:
double Ex;
- double ELab;
+ double E1;
+ double E2;
double Theta12;
double ThetaCM;
double VertexX;
@@ -74,6 +76,11 @@ class Analysis: public NPL::VAnalysis{
NPL::QFS* myQFS;
// Energy loss table: the G4Table are generated by the simulation
EnergyLoss BeamTarget;
+ EnergyLoss protonTarget;
+ EnergyLoss protonAl;
+ EnergyLoss protonSi;
+ double ReconstructProtonEnergy(const TVector3& x0,const TVector3& dir,const double& Ecatana);
+
TVector3 BeamImpact;
double TargetThickness ;
@@ -82,6 +89,7 @@ class Analysis: public NPL::VAnalysis{
// intermediate variable
TRandom3 Rand ;
TStrassePhysics* Strasse;
+ TCatanaPhysics* Catana;
TInitialConditions* IC ;
TInteractionCoordinates* DC;
TReactionConditions* RC;