[12.4.X] fix GeometryProducer and G4e examples

SimG4Core/GeometryProducer/interface/GeometryProducer.h

@@ -59,6 +59,8 @@ class GeometryProducer : public edm::one::EDProducer<edm::one::SharedResources,
 private:
   void updateMagneticField(edm::EventSetup const &es);
 
+  void makeGeom(const edm::EventSetup &c);
+
   G4RunManagerKernel *m_kernel;
   edm::ParameterSet m_pField;
   SimActivityRegistry m_registry;

SimG4Core/GeometryProducer/src/GeometryProducer.cc

@@ -61,15 +61,18 @@ GeometryProducer::GeometryProducer(edm::ParameterSet const &p)
       m_firstRun(true),
       m_pUseMagneticField(p.getParameter<bool>("UseMagneticField")),
       m_pUseSensitiveDetectors(p.getParameter<bool>("UseSensitiveDetectors")),
-      m_pGeoFromDD4hep(false) {
+      m_pGeoFromDD4hep(p.getParameter<bool>("GeoFromDD4hep")) {
   // Look for an outside SimActivityRegistry
   // this is used by the visualization code
+
   edm::Service<SimActivityRegistry> otherRegistry;
   if (otherRegistry)
     m_registry.connect(*otherRegistry);
   createWatchers(m_p, m_registry, m_watchers, m_producers);
 
-  m_sdMakers = sim::sensitiveDetectorMakers(m_p, consumesCollector(), std::vector<std::string>());
+  if (m_pUseSensitiveDetectors)
+    m_sdMakers = sim::sensitiveDetectorMakers(m_p, consumesCollector(), std::vector<std::string>());
+
   tokMF_ = esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();
   if (m_pGeoFromDD4hep) {
     tokDD4hep_ = esConsumes<cms::DDCompactView, IdealGeometryRecord, edm::Transition::BeginRun>();
@@ -103,6 +106,7 @@ void GeometryProducer::beginLuminosityBlock(edm::LuminosityBlock &, edm::EventSe
 }
 
 void GeometryProducer::beginRun(const edm::Run &run, const edm::EventSetup &es) {
+  makeGeom(es);
   updateMagneticField(es);
   for (auto &maker : m_sdMakers) {
     maker.second->beginRun(es);
@@ -115,6 +119,13 @@ void GeometryProducer::produce(edm::Event &e, const edm::EventSetup &es) {
   if (!m_firstRun)
     return;
   m_firstRun = false;
+  for (Producers::iterator itProd = m_producers.begin(); itProd != m_producers.end(); ++itProd) {
+    (*itProd)->produce(e, es);
+  }
+}
+void GeometryProducer::makeGeom(const edm::EventSetup &es) {
+  if (!m_firstRun)
+    return;
 
   edm::LogVerbatim("GeometryProducer") << "Producing G4 Geom";
 
@@ -156,10 +167,6 @@ void GeometryProducer::produce(edm::Event &e, const edm::EventSetup &es) {
     edm::LogInfo("GeometryProducer") << " Sensitive Detector building finished; found " << m_sensTkDets.size()
                                      << " Tk type Producers, and " << m_sensCaloDets.size() << " Calo type producers ";
   }
-
-  for (Producers::iterator itProd = m_producers.begin(); itProd != m_producers.end(); ++itProd) {
-    (*itProd)->produce(e, es);
-  }
 }
 
 DEFINE_FWK_MODULE(GeometryProducer);

TrackPropagation/Geant4e/python/geantRefit_cff.py

@@ -6,33 +6,20 @@
 
 from TrackingTools.TrackRefitter.TracksToTrajectories_cff import *
 
+from SimG4Core.Application.g4SimHits_cfi import g4SimHits as _g4SimHits
+
+
 ## Set up geometry
 geopro = cms.EDProducer("GeometryProducer",
+     GeoFromDD4hep = cms.bool(False),
      UseMagneticField = cms.bool(True),
      UseSensitiveDetectors = cms.bool(False),
-     MagneticField = cms.PSet(
-         UseLocalMagFieldManager = cms.bool(False),
-         Verbosity = cms.untracked.bool(False),
-         ConfGlobalMFM = cms.PSet(
-             Volume = cms.string('OCMS'),
-             OCMS = cms.PSet(
-                 Stepper = cms.string('G4ClassicalRK4'),
-                 Type = cms.string('CMSIMField'),
-                 StepperParam = cms.PSet(
-                     MaximumEpsilonStep = cms.untracked.double(0.01), ## in mm
-                     DeltaOneStep = cms.double(1e-4),## in mm
-                     MaximumLoopCounts = cms.untracked.double(1000.0),
-                     DeltaChord = cms.double(0.001), ## in mm
-                     MinStep = cms.double(0.1), ## in mm
-                     DeltaIntersectionAndOneStep = cms.untracked.double(-1.0),
-                     DeltaIntersection = cms.double(1e-6), ## in mm
-                     MinimumEpsilonStep = cms.untracked.double(1e-05) ## in mm
-                 )
-             )
-         ),
-         delta = cms.double(1.0)
-     )
-   )
+     MagneticField =  _g4SimHits.MagneticField.clone()
+)
+
+from Configuration.ProcessModifiers.dd4hep_cff import dd4hep
+dd4hep.toModify(geopro, GeoFromDD4hep = True )
+
 
 # load this to do a track refit
 from RecoTracker.TrackProducer.TrackRefitters_cff import *

TrackPropagation/Geant4e/test/BuildFile.xml

@@ -11,6 +11,8 @@
 <use name="Geometry/RPCGeometry"/>
 <use name="Geometry/CSCGeometry"/>
 <use name="DataFormats/MuonDetId"/>
-<library file="Geant4ePropagatorAnalyzer.cc" name="Geant4ePropagatorAnalyzer">
+<library file="*.cc" name="Geant4ePropagatorTests">
   <flags EDM_PLUGIN="1"/>
 </library>
+
+<test name="testG4Refitter" command="cmsRun ${LOCALTOP}/src/TrackPropagation/Geant4e/test/Geant4e_example_cfg.py"/>

TrackPropagation/Geant4e/test/Geant4ePropagatorAnalyzer.cc

@@ -147,6 +147,8 @@ Geant4ePropagatorAnalyzer::Geant4ePropagatorAnalyzer(const edm::ParameterSet &iC
       thePropagator(nullptr),
       G4VtxSrc_(iConfig.getParameter<edm::InputTag>("G4VtxSrc")),
       G4TrkSrc_(iConfig.getParameter<edm::InputTag>("G4TrkSrc")) {
+  auto token = consumes(G4TrkSrc_);
+
   // debug_ = iConfig.getParameter<bool>("debug");
   fStudyStation = iConfig.getParameter<int>("StudyStation");
 
@@ -334,15 +336,17 @@ void Geant4ePropagatorAnalyzer::analyze(const edm::Event &iEvent, const edm::Eve
   iEvent.getByLabel<SimTrackContainer>(G4TrkSrc_, simTracks);
   if (!simTracks.isValid()) {
     LogWarning("Geant4e") << "No tracks found" << std::endl;
-    return;
+    //    return;
   }
   LogDebug("Geant4e") << "G4e -- Got simTracks of size " << simTracks->size();
+  std::cout << "Geant4e "
+            << "G4e -- Got simTracks of size " << simTracks->size() << std::endl;
 
   Handle<SimVertexContainer> simVertices;
-  iEvent.getByLabel<SimVertexContainer>(G4VtxSrc_, simVertices);
+  iEvent.getByLabel<SimVertexContainer>("g4SimHits", simVertices);
   if (!simVertices.isValid()) {
-    LogWarning("Geant4e") << "No tracks found" << std::endl;
-    return;
+    LogWarning("Geant4e") << "No verticess found" << std::endl;
+    //    return;
   }
   LogDebug("Geant4e") << "Got simVertices of size " << simVertices->size();
 

TrackPropagation/Geant4e/test/Geant4e_example_cfg.py

@@ -1,19 +1,24 @@
-import FWCore.ParameterSet.Config 
+import FWCore.ParameterSet.Config as cms
 
-process = cms.Process("G4eRefit")
+from Configuration.Eras.Era_Run3_cff import Run3
 
+process = cms.Process("G4eRefit",Run3)
 
+# import of standard configurations
 process.load('Configuration.StandardSequences.Services_cff')
-process.load('Configuration.StandardSequences.GeometryDB_cff')
-process.load("Configuration.EventContent.EventContent_cff")
-process.load("Configuration.StandardSequences.Reconstruction_cff")
-process.load('Configuration.StandardSequences.MagneticField_38T_cff')
-process.load('Configuration.StandardSequences.EndOfProcess_cff')
-process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.StandardSequences.GeometrySimDB_cff')
+process.load('Configuration.StandardSequences.GeometryRecoDB_cff')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('Configuration.StandardSequences.Reconstruction_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
 
 
 from Configuration.AlCa.GlobalTag import GlobalTag
-process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc', '')
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase1_2022_realistic', '')
 
 process.load("FWCore.MessageService.MessageLogger_cfi")
 process.MessageLogger.default = cms.untracked.PSet(ERROR = cms.untracked.PSet(limit = cms.untracked.int32(5)))
@@ -23,8 +28,7 @@
 
 process.source = cms.Source("PoolSource",
     fileNames = cms.untracked.vstring(
-       '/store/relval/CMSSW_7_4_0_pre6/RelValSingleMuPt10_UP15/GEN-SIM-RECO/MCRUN2_74_V1-v1/00000/4C6BF9BF-95A9-E411-ADA9-0025905A60FE.root',
-       '/store/relval/CMSSW_7_4_0_pre6/RelValSingleMuPt10_UP15/GEN-SIM-RECO/MCRUN2_74_V1-v1/00000/DE6B164D-9FA9-E411-9BAF-0025905B85D0.root' 
+        '/store/relval/CMSSW_12_4_0/RelValSingleMuPt10/GEN-SIM-RECO/124X_mcRun3_2022_realistic_v5-v1/2580000/7a892809-eb85-4c21-8fc2-3723b4f85c85.root'
     ),
 )
 
@@ -42,7 +46,6 @@
   fileName = cms.untracked.string( 'test.root' )
 )
 
-process.g4RefitPath = cms.Path( process.MeasurementTrackerEvent * process.geant4eTrackRefit )
 process.e = cms.EndPath(process.out)
 process.schedule = cms.Schedule( process.g4RefitPath, process.e )
 

TrackPropagation/Geant4e/test/SimpleGeant4ePropagatorTest.cc

@@ -0,0 +1,137 @@
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"  //For define_fwk_module
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+//- Timing
+//#include "Utilities/Timing/interface/TimingReport.h"
+
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+
+//- Magnetic field
+#include "MagneticField/Engine/interface/MagneticField.h"
+#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+#include "SimG4Core/MagneticField/interface/Field.h"
+#include "SimG4Core/MagneticField/interface/FieldBuilder.h"
+
+//- Propagator
+#include "TrackPropagation/Geant4e/interface/ConvertFromToCLHEP.h"
+#include "TrackPropagation/Geant4e/interface/Geant4ePropagator.h"
+#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
+#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
+
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "DataFormats/GeometryVector/interface/GlobalVector.h"
+#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
+#include "MagneticField/VolumeGeometry/interface/MagVolumeOutsideValidity.h"
+#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
+
+//- Geant4
+#include "G4TransportationManager.hh"
+
+class SimpleGeant4ePropagatorTest final : public edm::one::EDAnalyzer<> {
+public:
+  explicit SimpleGeant4ePropagatorTest(const edm::ParameterSet &);
+  ~SimpleGeant4ePropagatorTest() override {}
+
+  void analyze(const edm::Event &, const edm::EventSetup &) override;
+
+protected:
+  // tokens
+  const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken;
+
+  Propagator *thePropagator;
+};
+
+SimpleGeant4ePropagatorTest::SimpleGeant4ePropagatorTest(const edm::ParameterSet &iConfig)
+    : magFieldToken(esConsumes()), thePropagator(nullptr) {}
+
+namespace {
+  Surface::RotationType rotation(const GlobalVector &zDir) {
+    GlobalVector zAxis = zDir.unit();
+    GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
+    GlobalVector xAxis = yAxis.cross(zAxis);
+    return Surface::RotationType(xAxis, yAxis, zAxis);
+  }
+}  // namespace
+
+void SimpleGeant4ePropagatorTest::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
+  using namespace edm;
+
+  std::cout << "Starting G4e test..." << std::endl;
+
+  ///////////////////////////////////////
+  // Construct Magnetic Field
+  const ESHandle<MagneticField> bField = iSetup.getHandle(magFieldToken);
+  if (bField.isValid())
+    std::cout << "G4e -- Magnetic field is valid. Value in (0,0,0): " << bField->inTesla(GlobalPoint(0, 0, 0)).mag()
+              << " Tesla " << std::endl;
+  else
+    LogError("Geant4e") << "G4e -- NO valid Magnetic field" << std::endl;
+
+  // Initialise the propagator
+  if (!thePropagator)
+    thePropagator = new Geant4ePropagator(bField.product());
+
+  if (thePropagator)
+    std::cout << "Propagator built!" << std::endl;
+  else
+    LogError("Geant4e") << "Could not build propagator!" << std::endl;
+
+  GlobalVector p3T(10., 10., 2.);
+  std::cout << "*** Phi (rad): " << p3T.phi() << " - Phi(deg)" << p3T.phi().degrees();
+  std::cout << "Track P.: " << p3T << "\nTrack P.: PT=" << p3T.perp() << "\tEta=" << p3T.eta()
+            << "\tPhi=" << p3T.phi().degrees() << "--> Rad: Phi=" << p3T.phi() << std::endl;
+
+  GlobalPoint r3T(0., 0., 0.);
+  std::cout << "Init point: " << r3T << "\nInit point Ro=" << r3T.perp() << "\tEta=" << r3T.eta()
+            << "\tPhi=" << r3T.phi().degrees() << std::endl;
+
+  //- Charge
+  int charge = 1;
+  std::cout << "Track charge = " << charge << std::endl;
+
+  //- Initial covariance matrix is unity 10-6
+  ROOT::Math::SMatrixIdentity id;
+  AlgebraicSymMatrix55 C(id);
+  C *= 0.01;
+  CurvilinearTrajectoryError covT(C);
+
+  PlaneBuilder pb;
+  Surface::RotationType rot = rotation(p3T);
+  // Define end planes
+  for (float d = 50.; d < 700.; d += 50.) {
+    float propDistance = d;  // 100 cm
+    std::cout << "G4e -- Extrapolatation distance " << d << " cm" << std::endl;
+    GlobalPoint targetPos = r3T + propDistance * p3T.unit();
+    auto endPlane = pb.plane(targetPos, rot);
+
+    //- Build FreeTrajectoryState
+    GlobalTrajectoryParameters trackPars(r3T, p3T, charge, &*bField);
+    FreeTrajectoryState ftsTrack(trackPars, covT);
+
+    // Propagate: Need to explicetely
+    TrajectoryStateOnSurface tSOSDest = thePropagator->propagate(ftsTrack, *endPlane);
+    if (!tSOSDest.isValid()) {
+      std::cout << "TSOS not valid? Propagation failed...." << std::endl;
+      continue;
+    }
+
+    auto posExtrap = tSOSDest.freeState()->position();
+    auto momExtrap = tSOSDest.freeState()->momentum();
+    std::cout << "G4e -- Extrapolated position:" << posExtrap << " cm\n"
+              << "G4e --       (Rho, eta, phi): (" << posExtrap.perp() << " cm, " << posExtrap.eta() << ", "
+              << posExtrap.phi() << ')' << std::endl;
+    std::cout << "G4e -- Extrapolated momentum:" << momExtrap << " GeV\n"
+              << "G4e --       (pt, eta, phi): (" << momExtrap.perp() << " cm, " << momExtrap.eta() << ", "
+              << momExtrap.phi() << ')' << std::endl;
+  }
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(SimpleGeant4ePropagatorTest);