From 13949b0fd5bef59ef99c3fe430f1c8800e450a80 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Fri, 24 Jul 2020 17:36:14 +0200
Subject: [PATCH] * adding delta information in tracking reconstruction

---
 NPLib/TrackReconstruction/NPTrackingUtility.h |  6 ++---
 Projects/Strasse/Analysis.cxx                 | 22 ++++++++++++++-----
 Projects/Strasse/Analysis.h                   |  3 +++
 Projects/Strasse/Control.cxx                  | 12 +++++++++-
 4 files changed, 34 insertions(+), 9 deletions(-)

diff --git a/NPLib/TrackReconstruction/NPTrackingUtility.h b/NPLib/TrackReconstruction/NPTrackingUtility.h
index cfc895251..4eb153866 100644
--- a/NPLib/TrackReconstruction/NPTrackingUtility.h
+++ b/NPLib/TrackReconstruction/NPTrackingUtility.h
@@ -22,7 +22,7 @@ namespace NPL{
   // 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){
+  double MinimumDistance(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
@@ -31,9 +31,9 @@ namespace NPL{
   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;
   double d = sqrt((e+v*t-w*s).Mag2());
-  
+ 
   BestPosition = 0.5*(v1+t*v+w1+s*w);
-  
+  delta = (v1+t*v-w1-s*w);
   return d;
   }
 }
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index 58ea82dff..229b6e1e6 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -57,7 +57,7 @@ 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",10000);
 } 
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -85,10 +85,14 @@ void Analysis::TreatEvent(){
     
     // computing minimum distance of the two lines
     TVector3 Vertex;
-    Distance = NPL::MinimumDistance(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex);
+    TVector3 delta;
+    Distance = NPL::MinimumDistance(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex,delta);
     VertexX=Vertex.X();
     VertexY=Vertex.Y();
     VertexZ=Vertex.Z();
+    deltaX=delta.X();
+    deltaY=delta.Y();
+    deltaZ=delta.Z();
   }
 
     //double thickness_before = 0;
@@ -125,9 +129,13 @@ 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("VertexX",&VertexX,"VertexX/D");
+  RootOutput::getInstance()->GetTree()->Branch("VertexY",&VertexY,"VertexY/D");
+  RootOutput::getInstance()->GetTree()->Branch("VertexZ",&VertexZ,"VertexZ/D");
+  RootOutput::getInstance()->GetTree()->Branch("deltaX",&deltaX,"deltaX/D");
+  RootOutput::getInstance()->GetTree()->Branch("deltaY",&deltaY,"deltaY/D");
+  RootOutput::getInstance()->GetTree()->Branch("deltaZ",&deltaZ,"deltaZ/D");
+
   RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
   RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
   RootOutput::getInstance()->GetTree()->Branch("ReactionConditions","TReactionConditions",&RC);
@@ -147,6 +155,10 @@ void Analysis::ReInitValue(){
   VertexX=-1000;
   VertexY=-1000;
   VertexZ=-1000;
+  deltaX=-1000;
+  deltaY=-1000;
+  deltaZ=-1000;
+
   Distance=-1000;
 }
 
diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h
index 39579ba16..c1e0c3f38 100644
--- a/Projects/Strasse/Analysis.h
+++ b/Projects/Strasse/Analysis.h
@@ -58,6 +58,9 @@ class Analysis: public NPL::VAnalysis{
     double VertexX;
     double VertexY;
     double VertexZ;
+    double deltaX;
+    double deltaY;
+    double deltaZ;
     double Distance;
     TLorentzVector LV_A;
     TLorentzVector LV_T;
diff --git a/Projects/Strasse/Control.cxx b/Projects/Strasse/Control.cxx
index 97c177d30..2627945eb 100644
--- a/Projects/Strasse/Control.cxx
+++ b/Projects/Strasse/Control.cxx
@@ -40,7 +40,17 @@ void Control(){
   PhysicsTree->Draw("Distance>>hd(500,0,80)", cond.c_str(),"");            
   //PhysicsTree->Draw("VertexY:VertexX:VertexZ", cond.c_str(),"");            
 
-
+  TCanvas* cdelta= new TCanvas("ControlDelta","ControlDelta",1000,1000);
+  cdelta->Divide(2,2);
+  cond = "deltaX!=-1000";
+  cdelta->cd(1);
+  PhysicsTree->Draw("deltaX>>dx(500,-2,2)",cond.c_str(),"col") ; 
+  cdelta->cd(2);
+  PhysicsTree->Draw("deltaY>>dy(500,-2,2)",cond.c_str(),"col") ; 
+  cdelta->cd(3);
+  PhysicsTree->Draw("deltaZ>>dz(5000,-20,20)",cond.c_str(),"col") ; 
+  cdelta->cd(4);
+  PhysicsTree->Draw("Distance>>hd(500,0,80)", cond.c_str(),"");           
  
 /*
   TCanvas* c2 = new TCanvas("Control 2", "Control2",500,500,2000,1000);  
-- 
GitLab