* Progress on Stasse/Catana analysis with evaluation of Initial proton

energy

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;
 }
 

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

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();
+    }
 }
 
 

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;

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;