Inheritance diagram for evdb_tool::DrawLArVoxel3D:

Public Member Functions
	DrawLArVoxel3D (const fhicl::ParameterSet &)

void	Draw (const art::Event &, evdb::View3D *) const override

Private Member Functions
int	GetMCTruth (const art::Event &, std::vector< const simb::MCTruth * > &) const

Detailed Description

Definition at line 34 of file DrawLArVoxel3D_tool.cc.

Constructor & Destructor Documentation

evdb_tool::DrawLArVoxel3D::DrawLArVoxel3D ( const fhicl::ParameterSet & pset )

explicit

Definition at line 46 of file DrawLArVoxel3D_tool.cc.

46 {}

Member Function Documentation

void evdb_tool::DrawLArVoxel3D::Draw	(	const art::Event &	evt,
		evdb::View3D *	view
	)		const

overridevirtual

Implements evdb_tool::ISim3DDrawer.

Definition at line 48 of file DrawLArVoxel3D_tool.cc.

   {
     art::ServiceHandle<evd::SimulationDrawingOptions const> drawOpt;
 
     // If the option is turned off, there's nothing to do
     if (!drawOpt->fShowSimChannelInfo) return;
 
     // If the option is turned off, there's nothing to do
     if (!drawOpt->fShowMCTruthTrajectories) return;
 
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
     auto const detProp =
       art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
     art::ServiceHandle<geo::Geometry const> geom;
 
     // Recover a handle to the collection of MCParticles
     art::Handle<std::vector<simb::MCParticle>> mcParticleHandle;
 
     evt.getByLabel(drawOpt->fG4ModuleLabel, mcParticleHandle);
 
     if (!mcParticleHandle.isValid()) return;
 
     // Define a couple of colors for neutrals and if we gray it out...
     int neutralColor(12);
     int grayedColor(15);
     int neutrinoColor(38);
 
     // Use the LArVoxelList to get the true energy deposition locations as opposed to using MCTrajectories
     const sim::LArVoxelList voxels =
       sim::SimListUtils::GetLArVoxelList(evt, drawOpt->fSimChannelLabel.encode());
 
     mf::LogDebug("SimulationDrawer")
       << "Starting loop over " << mcParticleHandle->size() << " McParticles, voxel list size is "
       << voxels.size() << std::endl;
 
     // Using the voxel information can be slow (see previous implementation of this code).
     // In order to speed things up we have modified the strategy:
     // 1) Make one pass through the list of voxels
     // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
     //    which is done by keeping a map between the MCParticle and a vector of positions
     // 3) Then loop through the map to draw the particle trajectories.
     // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
     // more loop to make a map of track id's and MCParticles.
 
     // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
     std::map<int, const simb::MCParticle*> trackToMcParticleMap;
 
     // Should we display the trajectories too?
     double minPartEnergy(0.01);
 
     for (size_t p = 0; p < mcParticleHandle->size(); ++p) {
       art::Ptr<simb::MCParticle> mcParticle(mcParticleHandle, p);
 
       trackToMcParticleMap[mcParticle->TrackId()] = mcParticle.get();
 
       // Quick loop through to draw trajectories...
       if (drawOpt->fShowMCTruthTrajectories) {
         // Is there an associated McTrajectory?
         const simb::MCTrajectory& mcTraj = mcParticle->Trajectory();
 
         int pdgCode(mcParticle->PdgCode());
         int colorIdx(evd::Style::ColorFromPDG(mcParticle->PdgCode()));
         TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
         double partCharge = partPDG ? partPDG->Charge() : 0.;
         double partEnergy = mcParticle->E();
 
         if (!drawOpt->fShowMCTruthColors) colorIdx = grayedColor;
 
         if (!mcTraj.empty() && partEnergy > minPartEnergy && mcParticle->TrackId() < 100000000) {
           double g4Ticks(clockData.TPCG4Time2Tick(mcParticle->T()) - trigger_offset(clockData));
           double xOffset(0.);
           double xPosMinTick = 0.;
           double xPosMaxTick = std::numeric_limits<double>::max();
 
           // collect the points from this particle
           int numTrajPoints = mcTraj.size();
 
           std::unique_ptr<double[]> hitPositions(new double[3 * numTrajPoints]);
           int hitCount(0);
 
           for (int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
             double xPos = mcTraj.X(hitIdx);
             double yPos = mcTraj.Y(hitIdx);
             double zPos = mcTraj.Z(hitIdx);
 
             // If we have cosmic rays then we need to get the offset which allows translating from
             // when they were generated vs when they were tracked.
             // Note that this also explicitly checks that they are in a TPC volume
             geo::Point_t hitLocation(xPos, yPos, zPos);
 
             try {
               geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
               geo::PlaneID planeID(tpcID, 0);
 
               xPosMinTick = detProp.ConvertTicksToX(0, planeID);
               xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
               xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
 
               if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
             }
             catch (...) {
               continue;
             }
 
             // Now move the hit position to correspond to the timing
             xPos += xOffset;
 
             // Check fiducial limits
             if (xPos > xPosMinTick && xPos < xPosMaxTick) {
               hitPositions[3 * hitCount] = xPos;
               hitPositions[3 * hitCount + 1] = yPos;
               hitPositions[3 * hitCount + 2] = zPos;
               hitCount++;
             }
           }
 
           TPolyLine3D& pl(view->AddPolyLine3D(1, colorIdx, 1, 1));
 
           // Draw neutrals as a gray dotted line to help fade into background a bit...
           if (partCharge == 0.) {
             pl.SetLineColor(neutralColor);
             pl.SetLineStyle(3);
             pl.SetLineWidth(1);
           }
           pl.SetPolyLine(hitCount, hitPositions.get(), "");
         }
       }
     }
 
     // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
     std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
 
     sim::LArVoxelList::const_iterator vxitr;
     for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++) {
       const sim::LArVoxelData& vxd = (*vxitr).second;
 
       for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++) {
         if (vxd.Energy(partIdx) > drawOpt->fMinEnergyDeposition) {
           int trackId = vxd.TrackID(partIdx);
 
           // It can be in some instances that mcPart here could be zero.
           const simb::MCParticle* mcPart = trackToMcParticleMap[trackId];
 
           partToPosMap[mcPart].push_back(std::vector<double>(3));
 
           partToPosMap[mcPart].back()[0] = vxd.VoxelID().X();
           partToPosMap[mcPart].back()[1] = vxd.VoxelID().Y();
           partToPosMap[mcPart].back()[2] = vxd.VoxelID().Z();
         }
       } // end if this track id is in the current voxel
     }   // end loop over voxels
 
     // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
     // draw the trajectories
     std::map<const simb::MCParticle*, std::vector<std::vector<double>>>::iterator partToPosMapItr;
 
     for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
          partToPosMapItr++) {
       // Recover the McParticle, we'll need to access several data members so may as well dereference it
       const simb::MCParticle* mcPart = partToPosMapItr->first;
 
       // Apparently, it can happen that we get a null pointer here or maybe no points to plot
       if (!mcPart || partToPosMapItr->second.empty()) continue;
 
       double g4Ticks(clockData.TPCG4Time2Tick(mcPart->T()) - trigger_offset(clockData));
       double xOffset = 0.;
       double xPosMinTick = 0.;
       double xPosMaxTick = std::numeric_limits<double>::max();
 
       int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
       int markerIdx(kFullDotSmall);
       int markerSize(2);
 
       if (!drawOpt->fShowMCTruthFullSize) {
         colorIdx = grayedColor;
         markerIdx = kDot;
         markerSize = 1;
       }
 
       std::unique_ptr<double[]> hitPositions(new double[3 * partToPosMapItr->second.size()]);
       int hitCount(0);
 
       // Now loop over points and add to trajectory
       for (size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
         const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
 
         // Check xOffset state and set if necessary
         geo::Point_t hitLocation(posVec[0], posVec[1], posVec[2]);
 
         try {
           geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
           geo::PlaneID planeID(tpcID, 0);
 
           xPosMinTick = detProp.ConvertTicksToX(0, planeID);
           xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
           xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
 
           if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
         }
         catch (...) {
           continue;
         }
 
         double xCoord = posVec[0] + xOffset;
 
         // If a voxel records an energy deposit then must have been in the TPC
         // But because things get shifted still need to cut off if outside drift
         if (xCoord > xPosMinTick && xCoord < xPosMaxTick) {
           hitPositions[3 * hitCount] = xCoord;
           hitPositions[3 * hitCount + 1] = posVec[1];
           hitPositions[3 * hitCount + 2] = posVec[2];
           hitCount++;
         }
       }
 
       TPolyMarker3D& pm = view->AddPolyMarker3D(1, colorIdx, markerIdx, markerSize);
       pm.SetPolyMarker(hitCount, hitPositions.get(), markerIdx);
     }
 
     // Finally, let's see if we can draw the incoming particle from the MCTruth information
     std::vector<const simb::MCTruth*> mctruth;
     this->GetMCTruth(evt, mctruth);
 
     // Loop through the MCTruth vector
     for (unsigned int idx = 0; idx < mctruth.size(); idx++) {
       // Go through each MCTruth object in the list
       for (int particleIdx = 0; particleIdx < mctruth[idx]->NParticles(); particleIdx++) {
         const simb::MCParticle& mcPart = mctruth[idx]->GetParticle(particleIdx);
 
         // A negative mother id indicates the "primary" particle
         if (mcPart.Mother() == -1 && mcPart.StatusCode() == 0) {
           mf::LogDebug("SimulationDrawer") << mcPart << std::endl;
 
           // Get position vector
           TVector3 particlePosition(mcPart.Vx(), mcPart.Vy(), mcPart.Vz());
 
           // Get direction vector (in opposite direction)
           TVector3 oppPartDir(-mcPart.Px(), -mcPart.Py(), -mcPart.Pz());
 
           if (oppPartDir.Mag2() > 0.) oppPartDir.SetMag(1.);
 
           double arcLenToDraw = -particlePosition.Z() / oppPartDir.CosTheta();
 
           // No point in drawing if arc length is zero (e.g. Ar nucleus)
           if (arcLenToDraw > 0.) {
             // Draw the line, use an off color to be unique
             TPolyLine3D& pl(view->AddPolyLine3D(2, neutrinoColor, 1, 2));
 
             pl.SetPoint(0, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
 
             particlePosition += std::min(arcLenToDraw + 10., 1000.) * oppPartDir;
 
             pl.SetPoint(1, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
           }
         }
         // The particles we want to draw will be early in the list so break out if we didn't find them
         else
           break;
       } // loop on particles in list
     }
 
     return;
   }

int evdb_tool::DrawLArVoxel3D::GetMCTruth	(	const art::Event &	evt,
		std::vector< const simb::MCTruth * > &	mcvec
	)		const

private

Definition at line 313 of file DrawLArVoxel3D_tool.cc.

References DEFINE_ART_CLASS_TOOL, e, art::ProductRetriever::getMany(), and art::Event::isRealData().

Referenced by Draw().

   {
     mcvec.clear();
 
     if (evt.isRealData()) return 0;
 
     std::vector<const simb::MCTruth*> temp;
 
     std::vector<art::Handle<std::vector<simb::MCTruth>>> mctcol;
 
     // use get by Type because there should only be one collection of these in the event
     try {
       //evt.getManyByType(mctcol);
       mctcol = evt.getMany<std::vector<simb::MCTruth>>();
 
       for (size_t mctc = 0; mctc < mctcol.size(); ++mctc) {
         art::Handle<std::vector<simb::MCTruth>> mclistHandle = mctcol[mctc];
 
         for (size_t i = 0; i < mclistHandle->size(); ++i) {
           temp.push_back(&(mclistHandle->at(i)));
         }
       }
       temp.swap(mcvec);
     }
     catch (cet::exception& e) {
       mf::LogWarning("DrawLArVoxel3D") << "GetMCTruth:"
                                        << " failed with message:\n"
                                        << e;
     }
 
     return mcvec.size();
   }

The documentation for this class was generated from the following file:

lareventdisplay/v09_11_05/source/lareventdisplay/EventDisplay/SimDrawers/DrawLArVoxel3D_tool.cc

Public Member Functions

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation