From deed2a2cbca8700cca62da09ea8836e5e7b6b6e3 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Mon, 20 Jul 2020 18:24:48 +0200
Subject: [PATCH] * working progress on strasse analysis

---
 NPLib/Detectors/Strasse/TStrassePhysics.cxx | 34 ++++++++----------
 Projects/Strasse/Analysis.cxx               | 40 ++++++++++++++++++---
 Projects/Strasse/Analysis.h                 | 26 +++++---------
 3 files changed, 58 insertions(+), 42 deletions(-)

diff --git a/NPLib/Detectors/Strasse/TStrassePhysics.cxx b/NPLib/Detectors/Strasse/TStrassePhysics.cxx
index a0274a76c..a4543ebf8 100644
--- a/NPLib/Detectors/Strasse/TStrassePhysics.cxx
+++ b/NPLib/Detectors/Strasse/TStrassePhysics.cxx
@@ -358,22 +358,14 @@ void TStrassePhysics::BuildPhysicalEvent() {
       vector<TVector2> inner = MatchInner();
       vector<TVector2> outer = MatchOuter();
   
-      EventMultiplicity = inner.size();
       
       for(unsigned int i=0; i<inner.size(); i++){
         int N = m_PreTreatedData->GetInner_TE_DetectorNbr(inner[i].X());
-        int T = m_PreTreatedData->GetInner_TE_StripNbr(inner[i].X());
-        int L = m_PreTreatedData->GetInner_LE_StripNbr(inner[i].Y());
+        int innerT = m_PreTreatedData->GetInner_TE_StripNbr(inner[i].X());
+        int innerL = m_PreTreatedData->GetInner_LE_StripNbr(inner[i].Y());
   
         double TE = m_PreTreatedData->GetInner_TE_Energy(inner[i].X());
-        DetectorNumber.push_back(N);
-        InnerStripT.push_back(T);
-        InnerStripL.push_back(L);
-        DE.push_back(TE);
-        InnerPosX.push_back(GetInnerPositionOfInteraction(i).x());
-        InnerPosY.push_back(GetInnerPositionOfInteraction(i).y());
-        InnerPosZ.push_back(GetInnerPositionOfInteraction(i).z());
-        // look for outer  
+               // look for outer  
         double outerE = 0;
         int outerT=0;
         int outerL=0;
@@ -384,22 +376,24 @@ void TStrassePhysics::BuildPhysicalEvent() {
             outerL = m_PreTreatedData->GetOuter_LE_StripNbr(outer[j].Y());
             }
         }
+
         if(outerE){
+          EventMultiplicity++;
+          DetectorNumber.push_back(N);
+          InnerStripT.push_back(innerT);
+          InnerStripL.push_back(innerL);
+          DE.push_back(TE);
+          InnerPosX.push_back(GetInnerPositionOfInteraction(i).x());
+          InnerPosY.push_back(GetInnerPositionOfInteraction(i).y());
+          InnerPosZ.push_back(GetInnerPositionOfInteraction(i).z());
+
           OuterStripT.push_back(outerT);
           OuterStripL.push_back(outerL);
           E.push_back(outerE);
           OuterPosX.push_back(GetOuterPositionOfInteraction(i).x());
           OuterPosY.push_back(GetOuterPositionOfInteraction(i).y());
           OuterPosZ.push_back(GetOuterPositionOfInteraction(i).z());
-          }
-        else{
-          OuterStripT.push_back(-1000);
-          OuterStripL.push_back(-1000);
-          E.push_back(-1000);
-          OuterPosX.push_back(-1000);
-          OuterPosY.push_back(-1000);
-          OuterPosZ.push_back(-1000);
-          }
+        }
       }
 
     }
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index ac3ada148..9dcea923c 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -35,7 +35,11 @@ Analysis::~Analysis(){
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::Init(){
-  IC= new TInteractionCoordinates;
+  IC= new TInitialConditions;
+  DC= new TInteractionCoordinates;
+  RC= new TReactionConditions;
+  
+
   InitOutputBranch();
   InitInputBranch();
   
@@ -51,6 +55,30 @@ void Analysis::Init(){
 void Analysis::TreatEvent(){
   // Reinitiate calculated variable
   ReInitValue();
+  unsigned int size = Strasse->GetEventMultiplicity();
+  if(size==2){ // 2 proton detected
+    // Proton 1
+    TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0);
+    TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0);
+    TVector3 Proton1 = OuterPos1-InnerPos1;
+    // Proton 2
+    TVector3 InnerPos2 = Strasse->GetInnerPositionOfInteraction(1);
+    TVector3 OuterPos2 = Strasse->GetOuterPositionOfInteraction(1);
+    TVector3 Proton2 = OuterPos2-InnerPos2;
+
+    double deltaPhi = abs(Proton1.Phi()/deg-Proton2.Phi()/deg);
+    double sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg
+    
+    // reject event that make no physical sense
+    if(deltaPhi<170 && sumTheta<80){
+      return
+      }
+    
+    // computing minimum distance of the two lines
+    TVector3 a;
+    TVector3 b;
+    
+    }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -62,13 +90,14 @@ 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("InteractionCoordinates","TInteractionCoordinates",&IC);
+  RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
+  RootOutput::getInstance()->GetTree()->Branch("ReactionConditions","TReactionConditions",&RC);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitInputBranch(){
-
-    RootInput:: getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&IC);
+    RootInput:: getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&DC);
+    RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&RC);
 }
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::ReInitValue(){
@@ -76,6 +105,9 @@ void Analysis::ReInitValue(){
   ELab = -1000;
   ThetaLab = -1000;
   ThetaCM = -1000;
+  VertexX=-1000;
+  VertexY=-1000;
+  VertexZ=-1000;
 }
 
 
diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h
index ff841d675..febf3d5e4 100644
--- a/Projects/Strasse/Analysis.h
+++ b/Projects/Strasse/Analysis.h
@@ -29,6 +29,7 @@
 #include "TStrassePhysics.h"
 #include "TInitialConditions.h"
 #include "TInteractionCoordinates.h"
+#include "TReactionConditions.h"
 #include <TRandom3.h>
 #include <TVector3.h>
 #include <TMath.h>
@@ -53,8 +54,11 @@ class Analysis: public NPL::VAnalysis{
   double ELab;
   double ThetaLab;
   double ThetaCM;
+  double VertexX;
+  double VertexY;
+  double VertexZ;
+
   NPL::Reaction* myReaction;
-TInitialConditions* myInit ;
   //	Energy loss table: the G4Table are generated by the simulation
   EnergyLoss LightCD2;
   EnergyLoss LightAl;
@@ -67,23 +71,9 @@ TInitialConditions* myInit ;
   double OriginalBeamEnergy ; // AMEV
                                                            // intermediate variable
   TRandom3 Rand ;
-  int DetectorNumber  ;
-  double ThetaNormalTarget;
-  double ThetaM2Surface ;
-  double Si_E_M2 ;
-  double CsI_E_M2  ;
-  double Energy ;
-  double E_M2 ;
-  
-  double ThetaSharcSurface ;
-  double X_Sharc ;
-  double Y_Sharc ;
-  double Z_Sharc  ;
-  double Si_E_Sharc ;
-  double E_Sharc ;
-  double Si_X_Sharc ;
-  double Si_Y_Sharc ;
   TStrassePhysics* Strasse;
-  TInteractionCoordinates* IC;
+  TInitialConditions* IC ;
+  TInteractionCoordinates* DC;
+  TReactionConditions* RC;
 };
 #endif
-- 
GitLab