Magnetic Field Adaptability for Different Geometries

TrackPerf/EfficiencyHists.hxx

@@ -5,6 +5,10 @@
 #include <TVector2.h>
 #include <TMath.h>
 
+// DD4hep
+#include <DD4hep/Detector.h>
+#include <DD4hep/DD4hepUnits.h>
+
 namespace EVENT {
 class Track;
 class MCParticle;
@@ -22,12 +26,10 @@ class EfficiencyHists {
 
   // Fill histograms with a single track
   void fillMC(const EVENT::MCParticle* track, bool passed);
-  void fillTrack(const EVENT::Track* track, bool passed);
+  void fillTrack(const EVENT::Track* track, bool passed,
+                 dd4hep::Detector* lcdd);
 
  private:
-  //! magnetic field to use for curvature -> pT conversion
-  float _Bz = 3.57;
-
   //! Efficiency plots
   TEfficiency* h_effpt;
   TEfficiency* h_efftheta;

TrackPerf/ResoHists.hxx

@@ -4,6 +4,12 @@
 #include <TH2.h>
 #include <TVector2.h>
 
+// DD4hep
+#include <DD4hep/Detector.h>
+#include <DD4hep/DD4hepUnits.h>
+
+
+
 namespace EVENT {
 class Track;
 class MCParticle;
@@ -20,12 +26,11 @@ class ResoHists {
   ResoHists();
 
   // Fill histograms with a single track
-  void fill(const EVENT::Track* track, const EVENT::MCParticle* particle);
+  void fill(const EVENT::Track* track,
+            const EVENT::MCParticle* particle,
+            dd4hep::Detector* lcdd);
 
  private:
-  //! magnetic field to use for curvature -> pT conversion
-  float _Bz = 3.57;
-
   //! Reconstructed track pT
   TH2* h_track_truth_pt;
   TH1* h_reso_pt_rel;

TrackPerf/TrackHists.hxx

@@ -3,6 +3,10 @@
 #include <TH1.h>
 #include <TH2.h>
 
+// DD4hep
+#include <DD4hep/Detector.h>
+#include <DD4hep/DD4hepUnits.h>
+
 namespace EVENT {
 class Track;
 }
@@ -18,12 +22,10 @@ class TrackHists {
   TrackHists();
 
   // Fill histograms with a single track
-  void fill(const EVENT::Track* track);
+  void fill(const EVENT::Track* track,
+            dd4hep::Detector* lcdd);
 
  private:
-  //! magnetic field to use for curvature -> pT conversion
-  float _Bz = 3.57;
-
   //! Reconstructed track pT
   TH1* h_pt;
   TH1* h_lambda;

TrackPerf/TrackPerfHistProc.hxx

@@ -4,6 +4,9 @@
 
 #include <marlin/Processor.h>
 
+// DD4hep
+#include <DD4hep/Detector.h>
+
 namespace TrackPerf {
 class TrackHists;
 class TruthHists;
@@ -49,6 +52,10 @@ class TrackPerfHistProc : public marlin::Processor {
   //! Track to MC truth match collection
   std::string _trkMatchColName{};
 
+  //! Magnetic field to use for curvature -> pT conversion
+  dd4hep::Detector* _lcdd;
+
+
   //! Determination of good vs bad match
   float _matchProb = 0.5;
 

src/EfficiencyHists.cxx

@@ -17,8 +17,16 @@ EfficiencyHists::EfficiencyHists(bool effi) {
   }
 }
 
-void EfficiencyHists::fillTrack(const EVENT::Track* track, bool passed) {
-  float pt = fabs(0.3 * _Bz / track->getOmega() / 1000);
+void EfficiencyHists::fillTrack(const EVENT::Track* track, bool passed,
+				dd4hep::Detector* lcdd) {
+  //TODO: This assumes uniform magnetic field
+  const double position[3] = {0, 0, 0};           // position to calculate magnetic field (here, the origin)
+  double magneticFieldVector[3] = {0, 0, 0};      // initialise object to hold magnetic field
+  lcdd->field().magneticField(
+        position, magneticFieldVector); // get the magnetic field vector from DD4hep
+  float Bz = magneticFieldVector[2] / dd4hep::tesla;
+
+  float pt = fabs(0.3 * Bz / track->getOmega() / 1000);
   float theta = TMath::Pi() - std::atan(track->getTanLambda());
 
   h_effpt->Fill(passed, pt);

src/ResoHists.cxx

@@ -23,8 +23,16 @@ ResoHists::ResoHists() {
 }
 
 void ResoHists::fill(const EVENT::Track* track,
-                     const EVENT::MCParticle* particle) {
-  float track_pt = fabs(0.3 * _Bz / track->getOmega() / 1000);
+                     const EVENT::MCParticle* particle, 
+		     dd4hep::Detector* lcdd) {
+  //TODO: This assumes uniform magnetic field
+  const double position[3] = {0, 0, 0};           // position to calculate magnetic field (here, the origin)
+  double magneticFieldVector[3] = {0, 0, 0};      // initialise object to hold magnetic field
+  lcdd->field().magneticField(
+        position, magneticFieldVector); // get the magnetic field vector from DD4hep
+  float Bz = magneticFieldVector[2] / dd4hep::tesla;
+
+  float track_pt = fabs(0.3 * Bz / track->getOmega() / 1000);
   float track_lambda = std::atan(track->getTanLambda());
 
   const double* mom = particle->getMomentum();

src/TrackHists.cxx

@@ -32,8 +32,16 @@ TrackHists::TrackHists() {
                20, -0.5, 19.5);
 }
 
-void TrackHists::fill(const EVENT::Track* track) {
-  float pt = fabs(0.3 * _Bz / track->getOmega() / 1000);
+void TrackHists::fill(const EVENT::Track* track, 
+		      dd4hep::Detector* lcdd) {
+  //TODO: This assumes uniform magnetic field
+  const double position[3] = {0, 0, 0};           // position to calculate magnetic field (here, the origin)
+  double magneticFieldVector[3] = {0, 0, 0};      // initialise object to hold magnetic field
+  lcdd->field().magneticField(
+        position, magneticFieldVector); // get the magnetic field vector from DD4hep
+  float Bz = magneticFieldVector[2] / dd4hep::tesla;
+
+  float pt = fabs(0.3 * Bz / track->getOmega() / 1000);
   h_pt->Fill(pt);
 
   float lambda = std::atan(track->getTanLambda());

src/TrackPerfHistProc.cxx

@@ -6,6 +6,10 @@
 #include <EVENT/MCParticle.h>
 #include <EVENT/Track.h>
 
+// DD4hep
+#include <DD4hep/DD4hepUnits.h>
+#include <DD4hep/Detector.h>
+
 #include <AIDA/ITree.h>
 #include <marlin/AIDAProcessor.h>
 #include <set>
@@ -76,6 +80,8 @@ void TrackPerfHistProc::init() {
   tree->cd("../efficiency");
   _effiPlots = std::make_shared<TrackPerf::EfficiencyHists>(true);
   _fakePlots = std::make_shared<TrackPerf::EfficiencyHists>(false);
+
+  _lcdd = &dd4hep::Detector::getInstance();
 }
 
 void TrackPerfHistProc::processRunHeader(LCRunHeader* /*run*/) {}
@@ -135,7 +141,7 @@ void TrackPerfHistProc::processEvent(LCEvent* evt) {
         static_cast<const EVENT::Track*>(trkCol->getElementAt(i));
 
     trkSet.insert(trk);
-    _allTracks->fill(trk);
+    _allTracks->fill(trk, _lcdd);
   }
   h_number_of_tracks->Fill(trkSet.size());
 
@@ -160,11 +166,11 @@ void TrackPerfHistProc::processEvent(LCEvent* evt) {
 
     if (rel->getWeight() > _matchProb) {
       if (trkSet.find(trk) != trkSet.end()) {
-        _realTracks->fill(trk);
+        _realTracks->fill(trk, _lcdd);
         _realTruths->fill(mcp);
 	_effiPlots->fillMC(mcp, true);
-        _realReso->fill(trk, mcp);
-	_fakePlots->fillTrack(trk, false);
+        _realReso->fill(trk, mcp, _lcdd);
+	_fakePlots->fillTrack(trk, false, _lcdd);
 
         mcpSet.erase(mcp);
         trkSet.erase(trk);
@@ -179,8 +185,8 @@ void TrackPerfHistProc::processEvent(LCEvent* evt) {
     _effiPlots->fillMC(mcp, false);
   }
   for (const EVENT::Track* trk : trkSet) {
-    _fakeTracks->fill(trk);
-    _fakePlots->fillTrack(trk, true);
+    _fakeTracks->fill(trk, _lcdd);
+    _fakePlots->fillTrack(trk, true, _lcdd);
   }
   h_number_of_fakes->Fill(trkSet.size());
 }