From 50daad27a95f3bc2dc974cc3f36f41a82e482154 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Wed, 22 Jul 2020 09:37:06 +0200
Subject: [PATCH] * Strasse vertex reconstruction

---
 NPLib/TrackReconstruction/NPTrackingUtility.h |  5 +--
 Projects/Strasse/Analysis.cxx                 | 14 ++++++--
 Projects/Strasse/Analysis.h                   |  2 ++
 Projects/Strasse/Control.cxx                  | 32 +++++++++++++++----
 Projects/Strasse/strasse_optimized.detector   |  6 ++--
 5 files changed, 43 insertions(+), 16 deletions(-)

diff --git a/NPLib/TrackReconstruction/NPTrackingUtility.h b/NPLib/TrackReconstruction/NPTrackingUtility.h
index 113f4efea..cfc895251 100644
--- a/NPLib/TrackReconstruction/NPTrackingUtility.h
+++ b/NPLib/TrackReconstruction/NPTrackingUtility.h
@@ -25,11 +25,8 @@ namespace NPL{
   double MinimumDistance(const TVector3& v1,const TVector3& v2, const TVector3& w1, const TVector3& w2, TVector3& BestPosition){
   TVector3 v = v2-v1;
   TVector3 w = w2-w1;
-  // let be n perpendicular to both line
-  TVector3 n = v.Cross(w);
-
   // Finding best position
-  TVector3 e = v2-w2;
+  TVector3 e = v1-w1;
   double A = -(v.Mag2()*w.Mag2()-(v.Dot(w)*v.Dot(w)));
   double s = (-v.Mag2()*(w.Dot(e))+(v.Dot(e))*(w.Dot(v)))/A;
   double t = (w.Mag2()*(v.Dot(e))-(w.Dot(e)*w.Dot(v)))/A;
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index 2ad26b6c1..39d571973 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -26,6 +26,7 @@ using namespace std;
 #include"NPDetectorManager.h"
 #include"NPOptionManager.h"
 #include"NPFunction.h"
+#include"NPTrackingUtility.h"
 ////////////////////////////////////////////////////////////////////////////////
 Analysis::Analysis(){
 }
@@ -76,9 +77,11 @@ void Analysis::TreatEvent(){
       }
     
     // computing minimum distance of the two lines
-    TVector3 a;
-    TVector3 b;
-    
+    TVector3 Vertex;
+    Distance = NPL::MinimumDistance(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex);
+    VertexX=Vertex.X();
+    VertexY=Vertex.Y();
+    VertexZ=Vertex.Z();
     }
 }
 
@@ -91,6 +94,10 @@ void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
   RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
   RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+  RootOutput::getInstance()->GetTree()->Branch("VerteX",&VertexX,"VertexX/D");
+  RootOutput::getInstance()->GetTree()->Branch("VerteY",&VertexY,"VertexY/D");
+  RootOutput::getInstance()->GetTree()->Branch("VerteZ",&VertexZ,"VertexZ/D");
+  RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
   RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
   RootOutput::getInstance()->GetTree()->Branch("ReactionConditions","TReactionConditions",&RC);
 }
@@ -109,6 +116,7 @@ void Analysis::ReInitValue(){
   VertexX=-1000;
   VertexY=-1000;
   VertexZ=-1000;
+  Distance=-1000;
 }
 
 
diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h
index febf3d5e4..a676d48c8 100644
--- a/Projects/Strasse/Analysis.h
+++ b/Projects/Strasse/Analysis.h
@@ -57,6 +57,8 @@ class Analysis: public NPL::VAnalysis{
   double VertexX;
   double VertexY;
   double VertexZ;
+  double Distance;
+
 
   NPL::Reaction* myReaction;
   //	Energy loss table: the G4Table are generated by the simulation
diff --git a/Projects/Strasse/Control.cxx b/Projects/Strasse/Control.cxx
index 0d0d9840f..ed0922e1c 100644
--- a/Projects/Strasse/Control.cxx
+++ b/Projects/Strasse/Control.cxx
@@ -13,7 +13,7 @@ void Control(){
   PhysicsTree->Draw("InnerPosZ:fDetected_Position_Z",cond.c_str(),"col") ; 
   cInner->cd(4);
   PhysicsTree->Draw("InnerPosY:InnerPosX:InnerPosZ", cond.c_str(),"");            
-*/
+
   TCanvas* cOuter = new TCanvas("ControlOuter","ControlOuter",1000,1000);
   cOuter->Divide(2,2);
   cond = "OuterPosX!=-1000";
@@ -25,15 +25,35 @@ void Control(){
   PhysicsTree->Draw("OuterPosZ:fDetected_Position_Z[3]",cond.c_str(),"col") ; 
   cOuter->cd(4);
   PhysicsTree->Draw("OuterPosY:OuterPosX:OuterPosZ", cond.c_str(),"");            
+*/
+  TCanvas* cVertex = new TCanvas("ControlVertex","ControlVertex",1000,1000);
+  cVertex->Divide(2,2);
+  cond = "VertexX!=-1000";
+  cVertex->cd(1);
+  PhysicsTree->Draw("VertexX:fRC_Vertex_Position_X",cond.c_str(),"col") ; 
+  TLine* lx = new TLine(-15,-15,15,15);
+  lx->Draw();
+  cVertex->cd(2);
+  PhysicsTree->Draw("VertexY:fRC_Vertex_Position_Y",cond.c_str(),"col") ; 
+  TLine* ly = new TLine(-15,-15,15,15);
+  ly->Draw();
+  cVertex->cd(3);
+  PhysicsTree->Draw("VertexZ:fRC_Vertex_Position_Z",cond.c_str(),"col") ; 
+  TLine* lz = new TLine(-80,-80,80,80);
+  lz->Draw();
+  cVertex->cd(4);
+
+  PhysicsTree->Draw("Distance", cond.c_str(),"");            
+  //PhysicsTree->Draw("VertexY:VertexX:VertexZ", cond.c_str(),"");            
 
 
+ 
+/*
   TCanvas* c2 = new TCanvas("Control 2", "Control2",500,500,2000,1000);  
   c2->Divide(2,1);
   c2->cd(1);
-  PhysicsTree->Draw("OuterPosY:OuterPosX",cond.c_str());            
-  PhysicsTree->Draw("InnerPosY:InnerPosX",cond.c_str(),"same");            
-
+  PhysicsTree->Draw("VertexY:VertexX:VertexZ",cond.c_str());            
   c2->cd(2);
-  PhysicsTree->Draw("OuterPosY:OuterPosZ",cond.c_str());            
-  PhysicsTree->Draw("InnerPosY:InnerPosZ",cond.c_str(),"same");            
+
+*/
 }
diff --git a/Projects/Strasse/strasse_optimized.detector b/Projects/Strasse/strasse_optimized.detector
index 986bac241..201a72f39 100644
--- a/Projects/Strasse/strasse_optimized.detector
+++ b/Projects/Strasse/strasse_optimized.detector
@@ -98,7 +98,7 @@ Strasse Outer
 Strasse Chamber
   Z= -30 mm
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Catana Dummy
- Z= 300 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Catana Dummy
+% Z= 300 mm
 
-- 
GitLab