Skip to content

Commit

Permalink
Merge pull request #38906 from mmusich/FixGeomProd
Browse files Browse the repository at this point in the history
[12.4.X] fix GeometryProducer and G4e examples
  • Loading branch information
cmsbuild committed Jul 31, 2022
2 parents 9ac63d8 + cb639a8 commit e313621
Show file tree
Hide file tree
Showing 10 changed files with 351 additions and 111 deletions.
2 changes: 2 additions & 0 deletions SimG4Core/GeometryProducer/interface/GeometryProducer.h
Expand Up @@ -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;
Expand Down
19 changes: 13 additions & 6 deletions SimG4Core/GeometryProducer/src/GeometryProducer.cc
Expand Up @@ -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>();
Expand Down Expand Up @@ -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);
Expand All @@ -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";

Expand Down Expand Up @@ -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);
33 changes: 10 additions & 23 deletions TrackPropagation/Geant4e/python/geantRefit_cff.py
Expand Up @@ -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 *
Expand Down
4 changes: 3 additions & 1 deletion TrackPropagation/Geant4e/test/BuildFile.xml
Expand Up @@ -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"/>
12 changes: 8 additions & 4 deletions TrackPropagation/Geant4e/test/Geant4ePropagatorAnalyzer.cc
Expand Up @@ -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");

Expand Down Expand Up @@ -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();

Expand Down
27 changes: 15 additions & 12 deletions 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)))
Expand All @@ -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'
),
)

Expand All @@ -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 )

Expand Down
137 changes: 137 additions & 0 deletions 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);

0 comments on commit e313621

Please sign in to comment.