Update to run with k4FWCore#173

key4hep · Feb 17, 2024 · 6b06055 · 6b06055
1 parent 9f4ba89
commit 6b06055
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 61 deletions.
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -23,7 +23,7 @@
 #include "edm4hep/SimTrackerHit.h"
 #include "edm4hep/TrackerHitPlaneCollection.h"
 
-#include "Gaudi/Accumulators/Histogram.h"
+#include "Gaudi/Accumulators/RootHistogram.h"
 #include "Gaudi/Histograming/Sink/Utils.h"
 
 #include "DD4hep/BitFieldCoder.h"
@@ -57,27 +57,17 @@ DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     throw std::runtime_error("DDPlanarDigi: Inconsistent number of resolutions given for U and V coordinate");
   }
 
-  m_histograms[hu].reset(new Gaudi::Accumulators::Histogram<1>{this, "hu", "smearing u", {50, -5., +5.}});
-  m_histograms[hv].reset(new Gaudi::Accumulators::Histogram<1>{this, "hv", "smearing v", {50, -5., +5.}});
-  m_histograms[hT].reset(new Gaudi::Accumulators::Histogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
+  m_histograms[hu].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hu", "smearing u", {50, -5., +5.}});
+  m_histograms[hv].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hv", "smearing v", {50, -5., +5.}});
+  m_histograms[hT].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
 
-  m_histograms[diffu].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffu", "diff u", {1000, -.1, +.1}});
-  m_histograms[diffv].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffv", "diff v", {1000, -.1, +.1}});
-  m_histograms[diffT].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
+  m_histograms[diffu].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffu", "diff u", {1000, -.1, +.1}});
+  m_histograms[diffv].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffv", "diff v", {1000, -.1, +.1}});
+  m_histograms[diffT].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
 
-  m_histograms[hitE].reset(new Gaudi::Accumulators::Histogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});
+  m_histograms[hitE].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});
   m_histograms[hitsAccepted].reset(
-      new Gaudi::Accumulators::Histogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
-
-  m_sigma[hu].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hv].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hT].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[diffu].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[diffv].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[diffT].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hitE].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hitsAccepted].reset(
-      new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+      new Gaudi::Accumulators::RootHistogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
 
   m_geoSvc = serviceLocator()->service(m_geoSvcName);
 }
@@ -127,7 +117,6 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 
   for (const auto& hit : simTrackerHits) {
     ++(*m_histograms[hitE])[hit.getEDep() * (dd4hep::GeV / dd4hep::keV)];
-    *m_sigma[hitE] += hit.getEDep() * (dd4hep::GeV / dd4hep::keV);
 
     if (hit.getEDep() < m_minEnergy) {
       debug() << "Hit with insufficient energy " << hit.getEDep() * (dd4hep::GeV / dd4hep::keV) << " keV" << endmsg;
@@ -183,9 +172,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 
       double tSmear = resT > 0 ? dist(m_engine) * resT : 0;
       ++(*m_histograms[hT])[resT > 0 ? tSmear / resT : 0];
-      *m_sigma[hT] += resT > 0 ? tSmear / resT : 0;
       ++(*m_histograms[diffT])[tSmear];
-      *m_sigma[diffT] += tSmear;
 
       hitT += tSmear;
       debug() << "smeared hit at T: " << hit.getTime() << " ns to T: " << hitT
@@ -264,14 +251,10 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
         newPos    = newPosTmp;
 
         ++(*m_histograms[hu])[uSmear / resU];
-        *m_sigma[hu] += uSmear / resU;
         ++(*m_histograms[hv])[vSmear / resV];
-        *m_sigma[hv] += vSmear / resV;
 
         ++(*m_histograms[diffu])[uSmear];
-        *m_sigma[diffu] += uSmear;
         ++(*m_histograms[diffv])[vSmear];
-        *m_sigma[diffv] += vSmear;
 
         break;
       }
@@ -335,7 +318,6 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
   // Filling the fraction of accepted hits in the event
   float accFraction = nSimHits > 0 ? float(nCreatedHits) / float(nSimHits) * 100.0 : 0.0;
   ++(*m_histograms[hitsAccepted])[accFraction];
-  *m_sigma[hitsAccepted] += accFraction;
 
   debug() << "Created " << nCreatedHits << " hits, " << nDismissedHits << " hits  dismissed" << endmsg;
 
@@ -344,11 +326,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 
 StatusCode DDPlanarDigi::finalize() {
   auto file  = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
-  auto names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
-  for (int i = 0; i < hSize; ++i) {
-    info() << "Standard deviation: " << m_sigma[i].get()->standard_deviation() << endmsg;
-    info() << "Mean: " << m_sigma[i].get()->mean() << endmsg;
-  }
+  std::vector<const char*> names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
   auto it = names.begin();
   for (auto& h : m_histograms) {
     std::string name = "";

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -19,8 +19,7 @@
 #ifndef DDPLANARDIGI_H
 #define DDPLANARDIGI_H
 
-#include <Gaudi/Accumulators.h>
-#include "Gaudi/Accumulators/Histogram.h"
+#include "Gaudi/Accumulators/RootHistogram.h"
 #include "Gaudi/Property.h"
 #include "GaudiAlg/Transformer.h"
 #include "GaudiKernel/RndmGenerators.h"
@@ -67,8 +66,6 @@
  * (default value lcio::ILDDetID::VXD) <br>
  * <br>
  * 
- * @author F.Gaede CERN/DESY, S. Aplin DESY
- * @date Dec 2014
  * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigi.h
  */
 
@@ -124,10 +121,7 @@ struct DDPlanarDigi final
                                                 "Output file name for the histograms"};
 
   const dd4hep::rec::SurfaceMap*                                        surfaceMap;
-  std::array<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>, hSize> m_histograms;
-  std::array<std::unique_ptr<Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>>,
-             hSize>
-              m_sigma;
+  std::array<std::unique_ptr<Gaudi::Accumulators::RootHistogram<1>>, hSize> m_histograms;
   std::string m_collName;
 
   inline static thread_local std::mt19937 m_engine;

diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
@@ -17,7 +17,8 @@
 # limitations under the License.
 #
 from Gaudi.Configuration import INFO
-from Configurables import ApplicationMgr, k4DataSvc, PodioInput, PodioOutput
+from k4FWCore import ApplicationMgr, IOSvc
+from Configurables import EventDataSvc
 from Configurables import DDPlanarDigi
 from Configurables import GeoSvc
 from Configurables import UniqueIDGenSvc
@@ -30,30 +31,28 @@
 geoservice.OutputLevel = INFO
 geoservice.EnableGeant4Geo = False
 
-processor = DDPlanarDigi()
-processor.SubDetectorName = "Vertex"
-processor.IsStrip = False
-processor.ResolutionU = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
-processor.ResolutionV = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
-processor.SimTrackerHitCollectionName = "VertexBarrelCollection"
-processor.SimTrkHitRelCollection = "VXDTrackerHitRelations"
-processor.TrackerHitCollectionName = "VXDTrackerHits"
+digi = DDPlanarDigi()
+digi.SubDetectorName = "Vertex"
+digi.IsStrip = False
+digi.ResolutionU = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
+digi.ResolutionV = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
+digi.SimTrackerHitCollectionName = "VertexBarrelCollection"
+digi.SimTrkHitRelCollection = "VXDTrackerHitRelations"
+digi.TrackerHitCollectionName = "VXDTrackerHits"
+digi.OutputFileName = "planar_digi_histograms.root"
 
-data_svc = k4DataSvc("EventDataSvc")
-data_svc.input = "input.root"
+iosvc = IOSvc()
+iosvc.input = "input.root"
+iosvc.output = "output_digi.root"
 
-inp = PodioInput()
-inp.collections = [
-    "VertexBarrelCollection",
-    "EventHeader",
-]
+# inp.collections = [
+#     "VertexBarrelCollection",
+#     "EventHeader",
+# ]
 
-out = PodioOutput("out")
-out.filename = "planar_digi_histograms.root"
-
-ApplicationMgr(TopAlg=[inp, processor, out],
+ApplicationMgr(TopAlg=[digi],
                EvtSel="NONE",
                EvtMax=-1,
-               ExtSvc=[data_svc],
+               ExtSvc=[EventDataSvc("EventDataSvc")],
                OutputLevel=INFO,
                )