Skip to content

Commit

Permalink
Geant4e
Browse files Browse the repository at this point in the history
  • Loading branch information
@lviliani
lviliani committed Feb 2, 2015
1 parent 4143a49 commit 87fba30
Show file tree
Hide file tree
Showing 7 changed files with 114 additions and 47 deletions.
3 changes: 3 additions & 0 deletions RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,8 @@ class TrackProducerAlgorithm {
conf_.getParameter<bool>( "GeometricInnerState" ) : true);
if (conf_.exists("reMatchSplitHits"))
reMatchSplitHits_=conf_.getParameter<bool>("reMatchSplitHits");
if (conf_.exists("usePropagatorForPCA"))
usePropagatorForPCA_ = conf_.getParameter<bool>("usePropagatorForPCA");
}

/// Destructor
Expand Down Expand Up @@ -136,6 +138,7 @@ class TrackProducerAlgorithm {
reco::TrackBase::TrackAlgorithm algo_;
bool reMatchSplitHits_;
bool geometricInnerState_;
bool usePropagatorForPCA_;

TrajectoryStateOnSurface getInitialState(const T * theT,
TransientTrackingRecHit::RecHitContainer& hits,
Expand Down
30 changes: 26 additions & 4 deletions RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -175,8 +176,19 @@ std::cout << algo_ << ": " << hits.size() <<'|' <<theTraj->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;
Expand Down Expand Up @@ -293,8 +305,18 @@ TrackProducerAlgorithm<reco::GsfTrack>::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;
Expand Down
16 changes: 12 additions & 4 deletions SimG4Core/GeometryProducer/interface/GeometryProducer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,22 @@ class GeometryProducer : public edm::EDProducer
public:
typedef std::vector<std::shared_ptr<SimProducer> > 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<std::shared_ptr<SimProducer> > producers() const
{ return m_producers; }
std::vector<SensitiveTkDetector*>& sensTkDetectors() { return m_sensTkDets; }
std::vector<SensitiveCaloDetector*>& sensCaloDetectors() { return m_sensCaloDets; }
private:

void updateMagneticField( edm::EventSetup const& es );

G4RunManagerKernel * m_kernel;
bool m_pUseMagneticField;
edm::ParameterSet m_pField;
Expand All @@ -52,6 +59,7 @@ class GeometryProducer : public edm::EDProducer
std::vector<SensitiveTkDetector*> m_sensTkDets;
std::vector<SensitiveCaloDetector*> m_sensCaloDets;
edm::ParameterSet m_p;
bool m_firstRun;
};

#endif
96 changes: 60 additions & 36 deletions SimG4Core/GeometryProducer/src/GeometryProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,8 @@ GeometryProducer::GeometryProducer(edm::ParameterSet const & p) :
m_pUseMagneticField(p.getParameter<bool>("UseMagneticField")),
m_pField(p.getParameter<edm::ParameterSet>("MagneticField")),
m_pUseSensitiveDetectors(p.getParameter<bool>("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
Expand All @@ -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<MagneticField> pMF;
es.get<IdealMagneticFieldRecord>().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<sim::FieldBuilder>(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<DDCompactView> pDD;
es.get<IdealGeometryRecord>().get(pDD);
Expand All @@ -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<MagneticField> pMF;
es.get<IdealMagneticFieldRecord>().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<sim::FieldBuilder>
(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<SimTrackManager>(new SimTrackManager);
if (m_attach==0) m_attach = new AttachSD;
{
std::pair< std::vector<SensitiveTkDetector*>,
std::vector<SensitiveCaloDetector*> >
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<SimTrackManager>(new SimTrackManager);
if (m_attach==0) m_attach = new AttachSD;
{
std::pair< std::vector<SensitiveTkDetector*>,
std::vector<SensitiveCaloDetector*> >
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);
1 change: 1 addition & 0 deletions TrackingTools/PatternTools/BuildFile.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
<use name="TrackingTools/TrajectoryParametrization"/>
<use name="TrackingTools/TrajectoryState"/>
<use name="rootrflx"/>
<Flags CXXFLAGS="-g -D=EDM_ML_DEBUG"/>
<export>
<lib name="1"/>
</export>
3 changes: 3 additions & 0 deletions TrackingTools/PatternTools/src/TSCBLBuilderWithPropagator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 9 additions & 3 deletions TrackingTools/PatternTools/src/TrajectoryExtrapolatorToLine.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,17 +13,18 @@ TrajectoryStateOnSurface TrajectoryExtrapolatorToLine::extrapolate(const FreeTra
DeepCopyPointerByClone<Propagator> 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;
double old_distance;
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);
Expand All @@ -37,6 +38,11 @@ TrajectoryStateOnSurface TrajectoryExtrapolatorToLine::extrapolate(const FreeTra
GlobalVector YY(T1.cross(XX));
Surface::RotationType rot(XX,YY);
ReferenceCountingPointer<Plane> 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);
Expand Down

0 comments on commit 87fba30

Please sign in to comment.