From ddc7530013abfdbc6191d747dcfb8b9f24dbeb6a Mon Sep 17 00:00:00 2001
From: Savannah Rose Shively <sshively@lxplus793.cern.ch>
Date: Wed, 15 Jun 2022 12:35:16 +0200
Subject: [PATCH] implement trees and track paramas

---
 .../python/SingleTrackExtrapolationConfig.py  |  6 +-
 .../src/SingleTrackExtrapolation.cxx          | 63 ++++++++++++++++---
 .../src/SingleTrackExtrapolation.h            | 22 +++++--
 3 files changed, 73 insertions(+), 18 deletions(-)

diff --git a/Tracking/Vertexing/SingleTrackExtrapolation/python/SingleTrackExtrapolationConfig.py b/Tracking/Vertexing/SingleTrackExtrapolation/python/SingleTrackExtrapolationConfig.py
index 6ab02b5f..7d313e63 100644
--- a/Tracking/Vertexing/SingleTrackExtrapolation/python/SingleTrackExtrapolationConfig.py
+++ b/Tracking/Vertexing/SingleTrackExtrapolation/python/SingleTrackExtrapolationConfig.py
@@ -33,8 +33,8 @@ if __name__ == "__main__":
 
    # Configure
     ConfigFlags.Input.Files = [f'/afs/cern.ch/work/s/sshively/private/MDC/samples/FaserMC-MDC_FS_Aee_100MeV_1Em5-110001-00000-00009-s0003-RDO-r0007-xAOD.root']
-    ConfigFlags.GeoModel.FaserVersion = "FASER-01"
-    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" 
+    ConfigFlags.GeoModel.FaserVersion = "FASERNU-03"
+    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" 
     ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" #is this necessary?
     ConfigFlags.GeoModel.Align.Dynamic = False
     ConfigFlags.lock()
@@ -53,4 +53,4 @@ if __name__ == "__main__":
     # replicaSvc.UseGeomSQLite = True
 
     # Execute and finish
-    sys.exit(int(acc.run(maxEvents=-1).isFailure()))
+    sys.exit(int(acc.run(maxEvents=10).isFailure()))
diff --git a/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.cxx b/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.cxx
index 568e209c..cc8dbcc7 100644
--- a/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.cxx
+++ b/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.cxx
@@ -19,16 +19,23 @@ SingleTrackExtrapolation::SingleTrackExtrapolation(const std::string &name, ISvc
 
 StatusCode SingleTrackExtrapolation::initialize() {
     ATH_CHECK( m_trackCollection.initialize() );
+    ATH_CHECK( m_mcEventKey.initialize() );
     ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
 
     m_tree = new TTree("trackCounts", "Track Counts");
 
-    m_tree->Branch("x", &m_x, "x/F");
-    m_tree->Branch("y", &m_y, "y/F");
-    m_tree->Branch("z", &m_z, "z/F");
-    m_tree->Branch("px", &m_px, "px/F");
-    m_tree->Branch("py", &m_py, "py/F");
-    m_tree->Branch("pz", &m_pz, "pz/F");
+    m_tree->Branch("event_counts", &m_event_counts, "event_counts/I");
+
+    m_tree->Branch("positiveMom",&m_positiveMom, "positiveMom/F");
+    m_tree->Branch("positiveDirX",&m_positiveDirX, "positiveDirX/F");
+    m_tree->Branch("positiveDirY",&m_positiveDirY, "positiveDirY/F");
+    m_tree->Branch("positiveDirZ",&m_positiveDirZ, "positiveDirZ/F");
+
+    m_tree->Branch("negativeMom",&m_negativeMom, "negativeMom/F");
+    m_tree->Branch("negativeDirX",&m_negativeDirX, "negativeDirX/F");
+    m_tree->Branch("negativeDirY",&m_negativeDirY, "negativeDirY/F");
+    m_tree->Branch("negativeDirZ",&m_negativeDirZ, "negativeDirZ/F");
+
     m_tree->Branch("eta", &m_eta, "eta/F");
     m_tree->Branch("theta", &m_theta, "theta/F");
     m_tree->Branch("pt", &m_pt, "pt/F");
@@ -52,12 +59,50 @@ StatusCode SingleTrackExtrapolation::execute(const EventContext& ctx) const {
     {
       track_count+=1;
       const Trk::TrackStateOnSurface* trackstate = (*track->trackStateOnSurfaces())[0];
-      const Trk::TrackParameters* parameters = trackstate->trackParameters();
-      //auto test = parameters-> #//how do i get charge and stuff
+      const Trk::TrackParameters* param = trackstate->trackParameters();
+
+      
+      auto position = param->position();
+      auto momentum = param->momentum();
+
+      double charge = param->charge();
 
-      //parameters->
+      if (charge>0)
+      {
+        m_positiveMom= sqrt(pow(momentum[0],2)+pow(momentum[1],2)+pow(momentum[2],2)); //momentum.normalize();
+        m_positiveDirX=position[0];
+        m_positiveDirY=position[1];
+        m_positiveDirZ=position[2];
+      }
+      else if (charge<0)
+      {
+        m_negativeMom= sqrt(pow(momentum[0],2)+pow(momentum[1],2)+pow(momentum[2],2));//momentum.normalize();
+        m_negativeDirX=position[0];
+        m_negativeDirY=position[1];
+        m_negativeDirZ=position[2];
+      }
+      else 
+      {
+        m_positiveMom=0;
+        m_positiveDirX=0;
+        m_positiveDirY=0;
+        m_positiveDirZ=0;
+
+        m_negativeMom=0;
+        m_negativeDirX=0;
+        m_negativeDirY=0;
+        m_negativeDirZ=0;
+      }
+      
+
+      m_event_counts=1;
+      m_tree->Fill();
 
     }
+
+    SG::ReadHandle<McEventCollection> h_mcEvents(m_mcEventKey);
+    if (h_mcEvents->size() == 0) return StatusCode::FAILURE;
+    //what now?? isnt this already iterating over events? wtf
 /*
 1. go through list of tracks
 2. save tracks in tree. importantly should include: charge of ptle, ptle ID, x,y,z, extrapolation into magnet
diff --git a/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.h b/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.h
index 5e9d812e..3c58e17f 100644
--- a/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.h
+++ b/Tracking/Vertexing/SingleTrackExtrapolation/src/SingleTrackExtrapolation.h
@@ -15,6 +15,7 @@
 #include "GaudiKernel/ServiceHandle.h"
 #include "GaudiKernel/ToolHandle.h"
 
+
 #include "TTree.h"
 #include "TFile.h"
 #include <TH1.h>
@@ -42,17 +43,26 @@ class SingleTrackExtrapolation : public AthReentrantAlgorithm, AthHistogramming
     const ServiceHandle<ITHistSvc>& histSvc() const;
 
     private:
+        
         const FaserSCT_ID* m_idHelper;
         SG::ReadHandleKey<TrackCollection> m_trackCollection {this, "TrackCollection", "CKFTrackCollection", "Input track collection name" };
+        SG::ReadHandleKey<McEventCollection> m_mcEventKey { this, "McEventCollection", "TruthEvent" };
 
         mutable TTree* m_tree;
 
-        mutable float m_x;
-        mutable float m_y;
-        mutable float m_z;
-        mutable float m_px;
-        mutable float m_py;
-        mutable float m_pz;
+        mutable int m_event_counts=0;
+
+        mutable float m_positiveMom;
+        mutable float m_positiveDirX;
+        mutable float m_positiveDirY;
+        mutable float m_positiveDirZ;
+
+        mutable float m_negativeMom;
+        mutable float m_negativeDirX;
+        mutable float m_negativeDirY;
+        mutable float m_negativeDirZ;
+
+
         mutable double m_chi2;
         mutable float m_eta;
         mutable float m_theta;
-- 
GitLab