25 #include "cetlib_except/exception.h" 28 #include "TDatabasePDG.h" 29 #include "TPolyLine3D.h" 30 #include "TPolyMarker3D.h" 68 if (!mcParticleHandle.
isValid())
return;
73 int neutrinoColor(38);
80 <<
"Starting loop over " << mcParticleHandle->size() <<
" McParticles, voxel list size is " 81 << voxels.
size() << std::endl;
93 std::map<int, const simb::MCParticle*> trackToMcParticleMap;
96 double minPartEnergy(0.01);
98 for (
size_t p = 0; p < mcParticleHandle->size(); ++p) {
101 trackToMcParticleMap[mcParticle->
TrackId()] = mcParticle.
get();
108 int pdgCode(mcParticle->
PdgCode());
110 TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
111 double partCharge = partPDG ? partPDG->Charge() : 0.;
112 double partEnergy = mcParticle->
E();
116 if (!mcTraj.
empty() && partEnergy > minPartEnergy && mcParticle->
TrackId() < 100000000) {
117 double g4Ticks(clockData.TPCG4Time2Tick(mcParticle->
T()) -
trigger_offset(clockData));
119 double xPosMinTick = 0.;
120 double xPosMaxTick = std::numeric_limits<double>::max();
123 int numTrajPoints = mcTraj.
size();
125 std::unique_ptr<double[]> hitPositions(
new double[3 * numTrajPoints]);
128 for (
int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
129 double xPos = mcTraj.
X(hitIdx);
130 double yPos = mcTraj.
Y(hitIdx);
131 double zPos = mcTraj.
Z(hitIdx);
142 xPosMinTick = detProp.ConvertTicksToX(0, planeID);
143 xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
144 xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
146 if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
156 if (xPos > xPosMinTick && xPos < xPosMaxTick) {
157 hitPositions[3 * hitCount] = xPos;
158 hitPositions[3 * hitCount + 1] = yPos;
159 hitPositions[3 * hitCount + 2] = zPos;
167 if (partCharge == 0.) {
168 pl.SetLineColor(neutralColor);
172 pl.SetPolyLine(hitCount, hitPositions.get(),
"");
178 std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
181 for (vxitr = voxels.
begin(); vxitr != voxels.
end(); vxitr++) {
186 int trackId = vxd.
TrackID(partIdx);
191 partToPosMap[mcPart].push_back(std::vector<double>(3));
193 partToPosMap[mcPart].back()[0] = vxd.
VoxelID().
X();
194 partToPosMap[mcPart].back()[1] = vxd.
VoxelID().
Y();
195 partToPosMap[mcPart].back()[2] = vxd.
VoxelID().
Z();
202 std::map<const simb::MCParticle*, std::vector<std::vector<double>>>
::iterator partToPosMapItr;
204 for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
210 if (!mcPart || partToPosMapItr->second.empty())
continue;
212 double g4Ticks(clockData.TPCG4Time2Tick(mcPart->
T()) -
trigger_offset(clockData));
214 double xPosMinTick = 0.;
215 double xPosMaxTick = std::numeric_limits<double>::max();
218 int markerIdx(kFullDotSmall);
222 colorIdx = grayedColor;
227 std::unique_ptr<double[]> hitPositions(
new double[3 * partToPosMapItr->second.size()]);
231 for (
size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
232 const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
235 geo::Point_t hitLocation(posVec[0], posVec[1], posVec[2]);
241 xPosMinTick = detProp.ConvertTicksToX(0, planeID);
242 xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
243 xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
245 if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
251 double xCoord = posVec[0] + xOffset;
255 if (xCoord > xPosMinTick && xCoord < xPosMaxTick) {
256 hitPositions[3 * hitCount] = xCoord;
257 hitPositions[3 * hitCount + 1] = posVec[1];
258 hitPositions[3 * hitCount + 2] = posVec[2];
263 TPolyMarker3D& pm = view->
AddPolyMarker3D(1, colorIdx, markerIdx, markerSize);
264 pm.SetPolyMarker(hitCount, hitPositions.get(), markerIdx);
268 std::vector<const simb::MCTruth*> mctruth;
272 for (
unsigned int idx = 0; idx < mctruth.size(); idx++) {
274 for (
int particleIdx = 0; particleIdx < mctruth[idx]->NParticles(); particleIdx++) {
279 mf::LogDebug(
"SimulationDrawer") << mcPart << std::endl;
282 TVector3 particlePosition(mcPart.
Vx(), mcPart.
Vy(), mcPart.
Vz());
285 TVector3 oppPartDir(-mcPart.
Px(), -mcPart.
Py(), -mcPart.
Pz());
287 if (oppPartDir.Mag2() > 0.) oppPartDir.SetMag(1.);
289 double arcLenToDraw = -particlePosition.Z() / oppPartDir.CosTheta();
292 if (arcLenToDraw > 0.) {
294 TPolyLine3D& pl(view->
AddPolyLine3D(2, neutrinoColor, 1, 2));
296 pl.SetPoint(0, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
298 particlePosition += std::min(arcLenToDraw + 10., 1000.) * oppPartDir;
300 pl.SetPoint(1, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
314 std::vector<const simb::MCTruth*>& mcvec)
const 320 std::vector<const simb::MCTruth*> temp;
322 std::vector<art::Handle<std::vector<simb::MCTruth>>> mctcol;
327 mctcol = evt.
getMany<std::vector<simb::MCTruth>>();
329 for (
size_t mctc = 0; mctc < mctcol.size(); ++mctc) {
332 for (
size_t i = 0; i < mclistHandle->size(); ++i) {
333 temp.push_back(&(mclistHandle->at(i)));
340 <<
" failed with message:\n" double E(const int i=0) const
double Z(const size_type i) const
double X(const size_type i) const
double Py(const int i=0) const
Container of LAr voxel information.
const simb::MCTrajectory & Trajectory() const
The data type to uniquely identify a Plane.
double Px(const int i=0) const
static sim::LArVoxelList GetLArVoxelList(const art::Event &evt, std::string moduleLabel)
art::InputTag fG4ModuleLabel
module label producing sim::SimChannel objects
list_type::const_iterator const_iterator
bool isValid() const noexcept
TPolyLine3D & AddPolyLine3D(int n, int c, int w, int s)
double Y(const size_type i) const
Encapsulates the information we want store for a voxel.
static int ColorFromPDG(int pdgcode)
double T(const int i=0) const
double fMinEnergyDeposition
bool fShowMCTruthTrajectories
size_type NumberParticles() const
The data type to uniquely identify a TPC.
double Vx(const int i=0) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
mapped_type Energy() const
art::InputTag fSimChannelLabel
SimChannels may be independent of MC stuff.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double Pz(const int i=0) const
double Vz(const int i=0) const
int trigger_offset(DetectorClocksData const &data)
sim::LArVoxelID VoxelID() const
A collection of 3D drawable objects.
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
double Vy(const int i=0) const
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
TPolyMarker3D & AddPolyMarker3D(int n, int c, int st, double sz)
bool fShowMCTruthFullSize
const key_type & TrackID(const size_type) const
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const