diff --git a/RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h b/RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h index c2fd20c3bd957..210582a46c161 100644 --- a/RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h +++ b/RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h @@ -60,6 +60,8 @@ class TrackProducerAlgorithm { conf_.getParameter( "GeometricInnerState" ) : true); if (conf_.exists("reMatchSplitHits")) reMatchSplitHits_=conf_.getParameter("reMatchSplitHits"); + if (conf_.exists("usePropagatorForPCA")) + usePropagatorForPCA_ = conf_.getParameter("usePropagatorForPCA"); } /// Destructor @@ -136,6 +138,7 @@ class TrackProducerAlgorithm { reco::TrackBase::TrackAlgorithm algo_; bool reMatchSplitHits_; bool geometricInnerState_; + bool usePropagatorForPCA_; TrajectoryStateOnSurface getInitialState(const T * theT, TransientTrackingRecHit::RecHitContainer& hits, diff --git a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc index 65a9643a98299..9e957c27c3f39 100644 --- a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc +++ b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc @@ -20,6 +20,7 @@ #include "DataFormats/TrackerRecHit2D/interface/TrackingRecHitLessFromGlobalPosition.h" #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h" +#include "TrackingTools/PatternTools/interface/TSCBLBuilderWithPropagator.h" #include "FWCore/Utilities/interface/Exception.h" #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit2DPosConstraint.h" @@ -175,8 +176,19 @@ std::cout << algo_ << ": " << hits.size() <<'|' <measurements().size( LogDebug("TrackProducer") << "stateForProjectionToBeamLine=" << stateForProjectionToBeamLine; - TSCBLBuilderNoMaterial tscblBuilder; - TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); +// TSCBLBuilderNoMaterial tscblBuilder; +// TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + + TrajectoryStateClosestToBeamLine tscbl; + if (usePropagatorForPCA_){ + //std::cout << "PROPAGATOR FOR PCA" << std::endl; + TSCBLBuilderWithPropagator tscblBuilder(*thePropagator); + tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + } else { + TSCBLBuilderNoMaterial tscblBuilder; + tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + } + if unlikely(!tscbl.isValid()) { delete theTraj; @@ -293,8 +305,18 @@ TrackProducerAlgorithm::buildTrack (const TrajectoryFitter * the LogDebug("GsfTrackProducer") << "stateForProjectionToBeamLine=" << stateForProjectionToBeamLine; - TSCBLBuilderNoMaterial tscblBuilder; - TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); +// TSCBLBuilderNoMaterial tscblBuilder; +// TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + + TrajectoryStateClosestToBeamLine tscbl; + if (usePropagatorForPCA_){ + TSCBLBuilderWithPropagator tscblBuilder(*thePropagator); + tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + } else { + TSCBLBuilderNoMaterial tscblBuilder; + tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + } + if unlikely(tscbl.isValid()==false) { delete theTraj; diff --git a/SimG4Core/GeometryProducer/interface/GeometryProducer.h b/SimG4Core/GeometryProducer/interface/GeometryProducer.h index 128c347b4c1ab..d56434d3faf36 100644 --- a/SimG4Core/GeometryProducer/interface/GeometryProducer.h +++ b/SimG4Core/GeometryProducer/interface/GeometryProducer.h @@ -30,15 +30,22 @@ class GeometryProducer : public edm::EDProducer public: typedef std::vector > Producers; explicit GeometryProducer(edm::ParameterSet const & p); - virtual ~GeometryProducer(); - virtual void beginJob(); - virtual void endJob(); - virtual void produce(edm::Event & e, const edm::EventSetup & c); +// virtual ~GeometryProducer(); +// virtual void beginJob(); +// virtual void endJob(); +// virtual void produce(edm::Event & e, const edm::EventSetup & c); + virtual ~GeometryProducer() override; + void produce(edm::Event & e, const edm::EventSetup & c); + void beginLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&); + std::vector > producers() const { return m_producers; } std::vector& sensTkDetectors() { return m_sensTkDets; } std::vector& sensCaloDetectors() { return m_sensCaloDets; } private: + + void updateMagneticField( edm::EventSetup const& es ); + G4RunManagerKernel * m_kernel; bool m_pUseMagneticField; edm::ParameterSet m_pField; @@ -52,6 +59,7 @@ class GeometryProducer : public edm::EDProducer std::vector m_sensTkDets; std::vector m_sensCaloDets; edm::ParameterSet m_p; + bool m_firstRun; }; #endif diff --git a/SimG4Core/GeometryProducer/src/GeometryProducer.cc b/SimG4Core/GeometryProducer/src/GeometryProducer.cc index 688b6f41a05f6..c582fbdb03bb1 100644 --- a/SimG4Core/GeometryProducer/src/GeometryProducer.cc +++ b/SimG4Core/GeometryProducer/src/GeometryProducer.cc @@ -59,7 +59,8 @@ GeometryProducer::GeometryProducer(edm::ParameterSet const & p) : m_pUseMagneticField(p.getParameter("UseMagneticField")), m_pField(p.getParameter("MagneticField")), m_pUseSensitiveDetectors(p.getParameter("UseSensitiveDetectors")), - m_attach(0), m_p(p) +// m_attach(0), m_p(p) + m_attach(0), m_p(p), m_firstRun ( true ) { //Look for an outside SimActivityRegistry //this is used by the visualization code @@ -75,17 +76,48 @@ GeometryProducer::~GeometryProducer() if (m_kernel!=0) delete m_kernel; } -void GeometryProducer::beginJob(){ +//void GeometryProducer::beginJob(){ +void GeometryProducer::updateMagneticField( edm::EventSetup const& es) { + if (m_pUseMagneticField) + { + // setup the magnetic field + edm::ESHandle pMF; + es.get().get(pMF); + const GlobalPoint g(0.,0.,0.); + LogDebug("GeometryProducer") << "B-field(T) at (0,0,0)(cm): " << pMF->inTesla(g) << std::endl; + + m_fieldBuilder = std::auto_ptr(new sim::FieldBuilder(&(*pMF), m_pField)); + + G4TransportationManager * tM = G4TransportationManager::GetTransportationManager(); + // update field here ... + m_fieldBuilder->build( tM->GetFieldManager(),tM->GetPropagatorInField()); + + LogDebug("GeometryProducer") << "Magentic field updated" << std::endl; + } + } -void GeometryProducer::endJob() -{ std::cout << " GeometryProducer terminating " << std::endl; } - +//void GeometryProducer::endJob() +//{ std::cout << " GeometryProducer terminating " << std::endl; } + +void GeometryProducer::beginLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const& es) { + // mag field can change in new lumi section + updateMagneticField( es ); +} + void GeometryProducer::produce(edm::Event & e, const edm::EventSetup & es) { - m_kernel = G4RunManagerKernel::GetRunManagerKernel(); +// m_kernel = G4RunManagerKernel::GetRunManagerKernel(); + if ( !m_firstRun ) + return; + m_firstRun = false; + + LogDebug("GeometryProducer") << "Producing G4 Geom" << std::endl; + + m_kernel = G4RunManagerKernel::GetRunManagerKernel(); if (m_kernel==0) m_kernel = new G4RunManagerKernel(); std::cout << " GeometryProducer initializing " << std::endl; + LogDebug("GeometryProducer") << " GeometryProducer initializing " << std::endl; // DDDWorld: get the DDCV from the ES and use it to build the World edm::ESTransientHandle pDD; es.get().get(pDD); @@ -95,41 +127,33 @@ void GeometryProducer::produce(edm::Event & e, const edm::EventSetup & es) const DDDWorld * world = new DDDWorld(&(*pDD), map_, catalog_, false); m_registry.dddWorldSignal_(world); - if (m_pUseMagneticField) - { - // setup the magnetic field - edm::ESHandle pMF; - es.get().get(pMF); - const GlobalPoint g(0.,0.,0.); - std::cout << "B-field(T) at (0,0,0)(cm): " << pMF->inTesla(g) << std::endl; - - m_fieldBuilder = std::auto_ptr - (new sim::FieldBuilder(&(*pMF), m_pField)); - // G4TransportationManager * tM = G4TransportationManager::GetTransportationManager(); - // m_fieldBuilder->configure("MagneticFieldType",tM->GetFieldManager(),tM->GetPropagatorInField()); - } + updateMagneticField( es ); if (m_pUseSensitiveDetectors) { - std::cout << " instantiating sensitive detectors " << std::endl; - // instantiate and attach the sensitive detectors - m_trackManager = std::auto_ptr(new SimTrackManager); - if (m_attach==0) m_attach = new AttachSD; - { - std::pair< std::vector, - std::vector > - sensDets = m_attach->create(*world,(*pDD),catalog_,m_p,m_trackManager.get(),m_registry); - - m_sensTkDets.swap(sensDets.first); - m_sensCaloDets.swap(sensDets.second); - } - - std::cout << " Sensitive Detector building finished; found " << m_sensTkDets.size() - << " Tk type Producers, and " << m_sensCaloDets.size() - << " Calo type producers " << std::endl; + LogDebug("GeometryProducer") << " instantiating sensitive detectors " << std::endl; + // instantiate and attach the sensitive detectors + m_trackManager = std::auto_ptr(new SimTrackManager); + if (m_attach==0) m_attach = new AttachSD; + { + std::pair< std::vector, + std::vector > + sensDets = m_attach->create(*world,(*pDD),catalog_,m_p,m_trackManager.get(),m_registry); + + m_sensTkDets.swap(sensDets.first); + m_sensCaloDets.swap(sensDets.second); + } + + LogDebug("GeometryProducer") << " Sensitive Detector building finished; found " << m_sensTkDets.size() + << " Tk type Producers, and " << m_sensCaloDets.size() + << " Calo type producers " << std::endl; + } for(Producers::iterator itProd = m_producers.begin();itProd != m_producers.end(); - ++itProd) { (*itProd)->produce(e,es); } + ++itProd) { + (*itProd)->produce(e,es); + } + } DEFINE_FWK_MODULE(GeometryProducer); diff --git a/TrackingTools/PatternTools/BuildFile.xml b/TrackingTools/PatternTools/BuildFile.xml index 05ad6fe641029..653bd7bac64c1 100644 --- a/TrackingTools/PatternTools/BuildFile.xml +++ b/TrackingTools/PatternTools/BuildFile.xml @@ -12,6 +12,7 @@ + diff --git a/TrackingTools/PatternTools/src/TSCBLBuilderWithPropagator.cc b/TrackingTools/PatternTools/src/TSCBLBuilderWithPropagator.cc index c0884a4964956..cbe97460f3087 100644 --- a/TrackingTools/PatternTools/src/TSCBLBuilderWithPropagator.cc +++ b/TrackingTools/PatternTools/src/TSCBLBuilderWithPropagator.cc @@ -29,6 +29,9 @@ TSCBLBuilderWithPropagator::operator() TrajectoryStateOnSurface tsosfinal = tetl.extrapolate(originalFTS,bsline,*thePropagator); + if (!tsosfinal.isValid()) + return TrajectoryStateClosestToBeamLine(); + //Compute point on beamline of closest approach GlobalPoint tp = tsosfinal.globalPosition(); //position of trajectory closest approach GlobalVector hyp(tp.x() - bspos.x(),tp.y() - bspos.y(),tp.z() - bspos.z()); //difference between traj and beamline reference diff --git a/TrackingTools/PatternTools/src/TrajectoryExtrapolatorToLine.cc b/TrackingTools/PatternTools/src/TrajectoryExtrapolatorToLine.cc index 74d7986887ef5..30e79b76a2181 100644 --- a/TrackingTools/PatternTools/src/TrajectoryExtrapolatorToLine.cc +++ b/TrackingTools/PatternTools/src/TrajectoryExtrapolatorToLine.cc @@ -13,7 +13,8 @@ TrajectoryStateOnSurface TrajectoryExtrapolatorToLine::extrapolate(const FreeTra DeepCopyPointerByClone p(aPropagator.clone()); p->setPropagationDirection(anyDirection); - FreeTrajectoryState fastFts(fts.parameters()); +// FreeTrajectoryState fastFts(fts.parameters()); + FreeTrajectoryState fastFts = fts; GlobalVector T1 = fastFts.momentum().unit(); GlobalPoint T0 = fastFts.position(); double distance = 9999999.9; @@ -21,9 +22,9 @@ TrajectoryStateOnSurface TrajectoryExtrapolatorToLine::extrapolate(const FreeTra int n_iter = 0; bool refining = true; - + LogDebug("TrajectoryExtrapolatorToLine") << "START REFINING"; while (refining) { - + LogDebug("TrajectoryExtrapolatorToLine") << "Refining cycle..."; // describe orientation of target surface on basis of track parameters n_iter++; Line T(T0,T1); @@ -37,6 +38,11 @@ TrajectoryStateOnSurface TrajectoryExtrapolatorToLine::extrapolate(const FreeTra GlobalVector YY(T1.cross(XX)); Surface::RotationType rot(XX,YY); ReferenceCountingPointer surface = Plane::build(pos, rot); + LogDebug("TrajectoryExtrapolatorToLine") << "Current plane position: " << surface->toGlobal(LocalPoint(0.,0.,0.)); + LogDebug("TrajectoryExtrapolatorToLine") << "Current plane normal: " << surface->toGlobal(LocalVector(0,0,1)); + LogDebug("TrajectoryExtrapolatorToLine") << "Current momentum: " << T1; + + // extrapolate fastFts to target surface TrajectoryStateOnSurface tsos = p->propagate(fastFts, *surface);