Commit 3a0eb874 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Fixed tracking units.

- The system of units is still a bit rocky.
- This might have something todo with the dd4pod output action. Need to
double check.
- Track finder works too well if it is seeded with initial track
parameters that include the exact thrown momentum.
parent f4eb7281
......@@ -58,13 +58,15 @@ namespace Jug {
//std::cout << ahit << "\n";
if(cell_hit_map.count(ahit.cellID()) == 0) {
cell_hit_map[ahit.cellID()] = rawhits->size();
eic::RawTrackerHit rawhit((long long)ahit.cellID(), ahit.truth().time*1000, (int)ahit.energyDeposit() * 10000);
eic::RawTrackerHit rawhit((long long)ahit.cellID(),
ahit.truth().time*1e6+m_gaussDist()*1e3, // ns->fs
int(ahit.energyDeposit() * 1e6));
rawhits->push_back(rawhit);
} else {
auto hit = (*rawhits)[cell_hit_map[ahit.cellID()]];
hit.time(ahit.truth().time*1e6+m_gaussDist()*1e3);
auto ch = hit.charge();
hit.charge( ch + (int)ahit.energyDeposit() * 10000);
hit.time(ahit.truth().time*1000+m_gaussDist());
hit.charge( ch + int(ahit.energyDeposit() * 1e6));
}
}
return StatusCode::SUCCESS;
......
......@@ -63,31 +63,21 @@ namespace Jug {
//debug() << trackstate.hasPredicted() << endmsg;
//debug() << trackstate.predicted() << endmsg;
auto params = trackstate.predicted() ;//<< endmsg;
debug() << 1.0/params[Acts::eBoundQOverP] << " GeV" << endmsg;
if ( std::abs(1.0 / params[Acts::eBoundQOverP]) > 10) {
double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
debug() << "track predicted p = " << p0 << " GeV" << endmsg;
if ( std::abs(p0) > 500) {
debug() << "skipping" << endmsg;
return;
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / params[Acts::eBoundQOverP], 0.105},
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / params[Acts::eBoundQOverP], 0.000511},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
(long long)13 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
debug() << p << endmsg;
(long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
//debug() << p << endmsg;
rec_parts->push_back(p);
});
//auto pos = m_geoSvc->cellIDPositionConverter()->position(ahit.cellID());
//auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(ahit.cellID());
//debug() << " dim size : " << std::size(dim) << endmsg;
//for(const auto& s : dim ) {
// debug() << "a size : " << s << endmsg;
//}
////std::array<double,3> posarr; pos.GetCoordinates(posarr);
////std::array<double,3> dimarr; dim.GetCoordinates(posarr);
////eic::TrackerHit hit;
//eic::TrackerHit hit((long long)ahit.cellID(), (long long)ahit.cellID(), (long long)ahit.time(),
// (float)ahit.charge() / 10000.0, (float)0.0, {{pos.x(), pos.y(),pos.z()}},{{dim[0],dim[1],0.0}});
//rec_hits->push_back(hit);
}
return StatusCode::SUCCESS;
}
......
......@@ -70,7 +70,7 @@ namespace Jug::Reco {
}
m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3D{0.0, 0.0, m_geoSvc->centralMagneticField()});
m_fieldctx = BFieldVariant(m_BField);
m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {100, 100}},
m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {15, 10}},
};
findTracks = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField);
......@@ -120,6 +120,7 @@ namespace Jug::Reco {
trajectories->emplace_back(std::move(trackFindingOutput.fittedStates), std::move(trackFindingOutput.trackTips),
std::move(trackFindingOutput.fittedParameters));
} else {
debug() << "Track finding failed for truth seed " << iseed << endmsg;
ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << result.error());
// Track finding failed, but still create an empty SimMultiTrajectory
trajectories->push_back(SimMultiTrajectory());
......
#include <cmath>
// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
......@@ -14,15 +15,31 @@
#include "eicd/TrackerHitCollection.h"
#include "dd4pod/Geant4ParticleCollection.h"
///// (Reconstructed) track parameters e.g. close to the vertex.
//using TrackParameters = Acts::CurvilinearTrackParameters;
///// Container of reconstructed track states for multiple tracks.
//using TrackParametersContainer = std::vector<TrackParameters>;
///// MultiTrajectory definition
//using Trajectory = Acts::MultiTrajectory<SourceLink>;
///// Container for the truth fitting/finding track(s)
//using TrajectoryContainer = std::vector<SimMultiTrajectory>;
namespace Jug::Reco {
/** Initial Track parameters from MC truth.
*
* TrackParmetersContainer
*/
class TrackParamTruthInit : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10};
Rndm::Numbers m_gaussDist;
DataHandle<dd4pod::Geant4ParticleCollection> m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader, this};
DataHandle<TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters", Gaudi::DataHandle::Writer, this};
DataHandle<dd4pod::Geant4ParticleCollection> m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader,
this};
DataHandle<TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
Gaudi::DataHandle::Writer, this};
public:
TrackParamTruthInit(const std::string& name, ISvcLocator* svcLoc)
......@@ -52,35 +69,37 @@ namespace Jug::Reco {
if(part.genStatus() != 1 ) {
continue;
}
using Acts::UnitConstants::MeV;
using Acts::UnitConstants::GeV;
using Acts::UnitConstants::mm;
double p = std::hypot( part.psx() * MeV, part.psy() * MeV, part.psz() * MeV);
// build some track cov matrix
Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 0.1 * mm*0.1 * mm;
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 0.1 * mm*0.1 * mm;
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (0.3 * GeV* 0.3 * GeV);
cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 0.98 / (p*p);
cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
init_trk_params->emplace_back(Acts::Vector4D(part.vsx() * mm, part.vsy() * mm, part.vsz() * mm, part.time() * Acts::UnitConstants::ns),
Acts::Vector3D(part.psx() * GeV, part.psy() * GeV, part.psz() * GeV),
std::sqrt(part.psx() *part.psx() + part.psy() * part.psy() + part.psz() * part.psz())*GeV,
Acts::Vector3D(part.psx() * MeV, part.psy() * MeV, part.psz() * MeV),
p+200*MeV,
((part.pdgID() > 0) ? -1 : 1),
std::make_optional(std::move(cov))
);
//part .charge()
debug() << "Invoke track finding seeded by truth particle " << part << endmsg;
debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV << " GeV" << endmsg;
//Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
//cov(Acts::eLOC_0, Acts::eLOC_0) = ahit.covMatrix(0)*ahit.covMatrix(0);
//cov(Acts::eLOC_1, Acts::eLOC_1) = ahit.covMatrix(1)*ahit.covMatrix(1);
}
return StatusCode::SUCCESS;
}
};
};
DECLARE_COMPONENT(TrackParamTruthInit)
} // namespace Jug::reco
......
......@@ -21,6 +21,8 @@
#include "eicd/TrackerHitCollection.h"
#include "JugReco/SourceLinks.h"
#include "DD4hep/DD4hepUnits.h"
namespace Jug {
namespace Reco {
......@@ -68,18 +70,22 @@ namespace Jug {
// Create output collections
auto rec_hits = m_outputHitCollection.createAndPut();
for(const auto& ahit : *rawhits) {
debug() << "cell ID : " << ahit.cellID() << endmsg;
//debug() << "cell ID : " << ahit.cellID() << endmsg;
auto pos = m_geoSvc->cellIDPositionConverter()->position(ahit.cellID());
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(ahit.cellID());
debug() << " dim size : " << std::size(dim) << endmsg;
for(const auto& s : dim ) {
debug() << "a size : " << s << endmsg;
}
//debug() << " dim size : " << std::size(dim) << endmsg;
//for(const auto& s : dim ) {
// debug() << "a size : " << s << endmsg;
//}
//std::array<double,3> posarr; pos.GetCoordinates(posarr);
//std::array<double,3> dimarr; dim.GetCoordinates(posarr);
//eic::TrackerHit hit;
eic::TrackerHit hit((long long)ahit.cellID(), (long long)ahit.time(),
(float)ahit.charge() / 10000.0, (float)0.0, {{pos.x(), pos.y(),pos.z()}},{{dim[0],dim[1],0.0}});
eic::TrackerHit hit((long long)ahit.cellID(),
(long long)ahit.time()/1000, // ps
(float)ahit.charge() / 1000.0, // MeV
(float)0.0,
{{pos.x()/dd4hep::mm, pos.y()/dd4hep::mm,pos.z()/dd4hep::mm}},
{{dim[0]/dd4hep::mm,dim[1]/dd4hep::mm,0.0}});
rec_hits->push_back(hit);
}
return StatusCode::SUCCESS;
......
......@@ -24,6 +24,9 @@
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DD4hep/Volumes.h"
#include "DD4hep/DD4hepUnits.h"
#include "Acts/Utilities/Units.hpp"
#include "Acts/Utilities/Definitions.hpp"
......@@ -38,8 +41,6 @@ namespace Jug::Reco {
class TrackerSourceLinker : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10};
Rndm::Numbers m_gaussDist;
DataHandle<eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<SourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
/// Pointer to the geometry service
......@@ -65,11 +66,6 @@ namespace Jug::Reco {
<< endmsg;
return StatusCode::FAILURE;
}
IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, m_timeResolution.value()));
if (!sc.isSuccess()) {
return StatusCode::FAILURE;
}
debug() << "visiting all the surfaces " << endmsg;
m_geoSvc->trackingGeometry()->visitSurfaces([this](const Acts::Surface* surface) {
// for now we just require a valid surface
......@@ -95,34 +91,40 @@ namespace Jug::Reco {
// setup local covariance
// TODO add support for per volume/layer/module settings
debug() << (*hits).size() << " hits " << endmsg;
for(const auto& ahit : *hits) {
Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = ahit.covMatrix(0)*Acts::UnitConstants::mm;//*ahit.covMatrix(0);
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = ahit.covMatrix(1)*Acts::UnitConstants::mm;//*ahit.covMatrix(1);
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = ahit.covMatrix(0)*Acts::UnitConstants::mm*ahit.covMatrix(0)*Acts::UnitConstants::mm;
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = ahit.covMatrix(1)*Acts::UnitConstants::mm*ahit.covMatrix(1)*Acts::UnitConstants::mm;
debug() << "cell ID : " << ahit.cellID() << endmsg;
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
debug() << " vol_id : " << vol_id << endmsg;
debug() << " hit : " << ahit << endmsg;
debug() << " hit : \n" << ahit << endmsg;
//debug() << "cell ID : " << ahit.cellID() << endmsg;
//debug() << " position : (" << ahit.position(0) << ", " << ahit.position(1) << ", "<< ahit.position(2) << ") " << endmsg;
debug() << " vol_id : " << vol_id << endmsg;
debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg;
const auto is = m_surfaces.find(vol_id);
if (is == m_surfaces.end()) {
debug() << " vol_id (" << vol_id << ") not found in m_surfaces." <<endmsg;
continue;
}
const Acts::Surface* surface = is->second;
debug() << " surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
// transform global position into local coordinates
Acts::Vector2D pos(0, 0);
// geometry context contains nothing here
surface->globalToLocal(Acts::GeometryContext(), {ahit.position(0), ahit.position(1), ahit.position(2)},
{0, 0, 0});//, pos);
pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.position(0), ahit.position(1), ahit.position(2)}, {0, 0, 0}).value();//, pos);
//// smear truth to create local measurement
Acts::BoundVector loc = Acts::BoundVector::Zero();
//loc[Acts::eLOC_0] = pos[0] + m_cfg.sigmaLoc0 * stdNormal(rng);
//loc[Acts::eLOC_1] = pos[1] + m_cfg.sigmaLoc1 * stdNormal(rng);
loc[Acts::eBoundLoc0] = pos[0] ;//+ m_cfg.sigmaLoc0 * stdNormal(rng);
loc[Acts::eBoundLoc0] = pos[1] ;//+ m_cfg.sigmaLoc1 * stdNormal(rng);
debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
// create source link at the end of the container
//auto it = source_links->emplace_hint(source_links->end(), *surface, hit, 2, loc, cov);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment