#include <SimulationDrawer.h>

Public Member Functions
	SimulationDrawer ()

	~SimulationDrawer ()

void	MCTruthShortText (const art::Event &evt, evdb::View2D *view)

void	MCTruthLongText (const art::Event &evt, evdb::View2D *view)

void	MCTruth3D (const art::Event &evt, evdb::View3D *view)

void	HiLite (int trkId, bool hlt=true)

Public Attributes
double	minx

double	maxx

double	miny

double	maxy

double	minz

double	maxz

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

int	GetParticle (const art::Event &evt, std::vector< const simb::MCParticle * > &plist)

Private Attributes
std::map< int, bool >	fHighlite

Detailed Description

Definition at line 25 of file SimulationDrawer.h.

Constructor & Destructor Documentation

gar::evd::SimulationDrawer::SimulationDrawer ( )

Definition at line 51 of file SimulationDrawer.cxx.

     {
       art::ServiceHandle<geo::GeometryGAr> geom;
       minx =  1e9;
       maxx = -1e9;
       miny =  1e9;
       maxy = -1e9;
       minz =  1e9;
       maxz = -1e9;
       double world[3] = {geom->TPCXCent(),geom->TPCYCent(),geom->TPCZCent()};
 
       if (minx>world[0]-geom->TPCRadius())
         minx = world[0]-geom->TPCRadius();
       if (maxx<world[0]+geom->TPCRadius())
         maxx = world[0]+geom->TPCRadius();
       if (miny>world[1]-geom->TPCRadius())
         miny = world[1]-geom->TPCRadius();
       if (maxy<world[1]+geom->TPCRadius())
         maxy = world[1]+geom->TPCRadius();
       if (minz>world[2]-geom->TPCLength()/2.)
         minz = world[2]-geom->TPCLength()/2.;
       if (maxz<world[2]+geom->TPCLength()/2.)
         maxz = world[2]+geom->TPCLength()/2.;
     }

gar::evd::SimulationDrawer::~SimulationDrawer ( )

Definition at line 78 of file SimulationDrawer.cxx.

79 {

80 }

Member Function Documentation

int gar::evd::SimulationDrawer::GetMCTruth	(	const art::Event &	evt,
		std::vector< const simb::MCTruth * > &	mctruth
	)

private

Definition at line 439 of file SimulationDrawer.cxx.

     {
       mcvec.clear();
 
       if( evt.isRealData() ) return 0;
 
       std::vector<const simb::MCTruth*> temp;
 
       // use get by Type because there should only be one collection of these in the event
       try{
         auto 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){
         writeErrMsg("GetMCTruth", e);
       }
 
       return mcvec.size();
     }

int gar::evd::SimulationDrawer::GetParticle	(	const art::Event &	evt,
		std::vector< const simb::MCParticle * > &	plist
	)

private

Definition at line 410 of file SimulationDrawer.cxx.

     {
       plist.clear();
 
       if( evt.isRealData() ) return 0;
 
       art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
 
       std::vector<const simb::MCParticle*> temp;
 
       art::View<simb::MCParticle> plcol;
       // use get by Type because there should only be one collection of these in the event
       try{
         evt.getView(drawopt->fG4ModuleLabel, plcol);
         for(unsigned int i = 0; i < plcol.vals().size(); ++i){
           temp.push_back(plcol.vals().at(i));
         }
         temp.swap(plist);
       }
       catch(cet::exception& e){
         writeErrMsg("GetRawDigits", e);
       }
 
       return plist.size();
 
     }

void gar::evd::SimulationDrawer::HiLite	(	int	trkId,
		bool	hlt = `true`
	)

Definition at line 468 of file SimulationDrawer.cxx.

     {
       fHighlite[trkId] = dohilite;
     }

void gar::evd::SimulationDrawer::MCTruth3D	(	const art::Event &	evt,
		evdb::View3D *	view
	)

Definition at line 177 of file SimulationDrawer.cxx.

     {
       if( evt.isRealData() ) return;
 
       art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
 
       art::ServiceHandle<geo::GeometryGAr> geo;
       double xcent = geo->TPCXCent();
       double ycent = geo->TPCYCent();
       double zcent = geo->TPCZCent();
       xcent = 0;
       ycent = 0;
       zcent = 0;
 
       // get the particles from the Geant4 step
       std::vector<const simb::MCParticle*> plist;
       this->GetParticle(evt, plist);
 
       // Useful variables
       double xMinimum(-1.*(maxx-minx));
       double xMaximum( 2.*(maxx-minx));
 
       // Define a couple of colors for neutrals and if we gray it out...
       int neutralColor(12);
       int grayedColor(15);
       int neutrinoColor(38);
 
       // 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);
 
       std::cout << "MC particle listing (limit 20)" << std::endl;
       std::cout << "    ID          PDG       Px       Py       Pz    E (GeV)       origin (cm) " << std::endl;
 
       for(size_t p = 0; p < plist.size(); ++p){
 
         trackToMcParticleMap[plist[p]->TrackId()] = plist[p];
 
         const simb::MCParticle*   mcPart = plist[p];
         const simb::MCTrajectory& mcTraj = mcPart->Trajectory();
 
         int           pdgCode(mcPart->PdgCode());
         int           colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
         TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
         double        partCharge = partPDG ? partPDG->Charge() : 0.;
         double        partEnergy = mcPart->E();
 
         if (p<20)
           {
             printf("%6d %12d %8.3f %8.3f %8.3f %8.3f ",(int) p,pdgCode,mcPart->Px(),mcPart->Py(),mcPart->Pz(),partEnergy);
             if (!mcTraj.empty())
               {
                 double xPos = mcTraj.X(0) + xcent;
                 double yPos = mcTraj.Y(0) + ycent;
                 double zPos = mcTraj.Z(0) + zcent;
                 printf("%8.3f %8.3f %8.3f\n",xPos, yPos, zPos);
               }
             else
               {
                 printf("\n");
               }
           }
 
         // If the option is turned off, stop here after printing it out
         if (!drawopt->fShowMCTruthTrajectories) continue;
 
         // cut out all neutrals if asked to do so
         if ( !drawopt->fShowNeutrals && partCharge > -0.1 && partCharge < 0.1) continue;
 
         // energy cut on photons
         if (pdgCode == 22 && partEnergy < drawopt->fPhotonEnergyCut) continue;
 
         // energy cut on neutrons
         if (pdgCode == 2112 && partEnergy < drawopt->fNeutronEnergyCut) continue;
 
         // energy cut on all other neutral particles, not photons or neutrons
         if ( partEnergy < drawopt->fOtherNeutralEnergyCut  && partCharge > -0.1 && partCharge < 0.1
              && pdgCode != 22 && pdgCode != 2112) continue;
 
         if (!drawopt->fShowMCTruthColors) colorIdx = grayedColor;
 
         if (pdgCode >= 1000000000) continue;
         if (pdgCode == 111) continue; // don't draw pizeros
 
         //if ( partCharge > -0.1 && partCharge < 0.1 )
         //{
         //  std::cout << "Drawing neutral particle: " << pdgCode << " " << partEnergy << std::endl;
         //}
 
 
         if (!mcTraj.empty() && partEnergy > minPartEnergy && mcPart->TrackId() < 100000000){
           // The following is meant to get the correct offset for drawing the particle trajectory
           // In particular, the cosmic rays will not be correctly placed without this
 
 
           // 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) + xcent;
             double yPos = mcTraj.Y(hitIdx) + ycent;
             double zPos = mcTraj.Z(hitIdx) + zcent;
 
             // If the original simulated hit did not occur in the TPC volume then don't draw it
             if (xPos < minx || xPos > maxx || yPos < miny || yPos > maxy|| zPos < minz || zPos > maxz) continue;
 
             // Check fiducial limits
             // nb. coordinates are so Y is up, but ROOT thinks Z is up, so report (z,x,y) to the polyline coords
             if (xPos > xMinimum && xPos < xMaximum){
               hitPositions[3*hitCount    ] = zPos;
               hitPositions[3*hitCount + 1] = xPos;
               hitPositions[3*hitCount + 2] = yPos;
               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;
 
       // 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;
 
           // The following is meant to get the correct offset for drawing the particle trajectory
           // In particular, the cosmic rays will not be correctly placed without this
 
           int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
           int markerIdx(kFullDotSmall);
           int markerSize(2);
 
           if (!drawopt->fShowMCTruthFullSize){
             colorIdx   = grayedColor;
             markerIdx  = kDot;
             markerSize = 1;
           }
 
           TPolyMarker3D& pm = view->AddPolyMarker3D(partToPosMapItr->second.size(), colorIdx, markerIdx, markerSize);
 
           // 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];
 
             double xCoord = posVec[0];
 
             // nb. coordinates are so Y is up, but ROOT thinks Z is up, so report (z,x,y)
             if (xCoord > xMinimum && xCoord < xMaximum)
               pm.SetPoint(posIdx, posVec[2], xCoord, posVec[1]);
           }
         }
 
       // 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();
             double arcLenToDraw = 50.0;  // always draw 50 cm of neutrinos
 
             // 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));
 
               // nb. coordinates are so Y is up, but ROOT thinks Z is up, so report (z,x,y)
 
               pl.SetPoint(0,particlePosition.Z()+zcent,particlePosition.X()+xcent,particlePosition.Y()+ycent);
 
               particlePosition += std::min(arcLenToDraw + 10.,1000.) * oppPartDir;
 
               pl.SetPoint(1,particlePosition.Z()+zcent,particlePosition.X()+xcent,particlePosition.Y()+ycent);
             }
           }
           // 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;
     }

void gar::evd::SimulationDrawer::MCTruthLongText	(	const art::Event &	evt,
		evdb::View2D *	view
	)

Definition at line 144 of file SimulationDrawer.cxx.

     {
       if( evt.isRealData() ) return;
 
       art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
       // Skip drawing if option is turned off
       if (!drawopt->fShowMCTruthText) return;
 
       std::vector<const simb::MCTruth*> mctruth;
       this->GetMCTruth(evt, mctruth);
       std::cout<<"\nMCTruth Ptcl trackID            PDG      P      T   Moth  Process\n";
       for (unsigned int i=0; i<mctruth.size(); ++i) {
         for (int j=0; j<mctruth[i]->NParticles(); ++j) {
           const simb::MCParticle& p = mctruth[i]->GetParticle(j);
           if(p.StatusCode() == 0 || p.StatusCode() == 1) {
             int KE = 1000 * (p.E() - p.Mass());
             std::cout<<std::right<<std::setw(7)<<i<<std::setw(5)<<j
                      <<std::setw(8)<<p.TrackId()
                      <<" "<<std::setw(14)<<Style::LatexName(p.PdgCode())
                      <<std::setw(7)<<int(1000 * p.P())
                      <<std::setw(7)<<KE<<std::setw(7)<<p.Mother()
                      <<" "<<p.Process()
                      <<"\n";
           }
         } // loop on j particles in list
       }
       std::cout<<"Note: Momentum, P, and kinetic energy, T, in MeV/c\n";
     } // MCTruthLongText

void gar::evd::SimulationDrawer::MCTruthShortText	(	const art::Event &	evt,
		evdb::View2D *	view
	)

Definition at line 83 of file SimulationDrawer.cxx.

     {
 
       if( evt.isRealData() ) return;
 
       art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
       // Skip drawing if option is turned off
       if (!drawopt->fShowMCTruthText) return;
 
       std::vector<const simb::MCTruth*> mctruth;
       this->GetMCTruth(evt, mctruth);
 
       for (unsigned int i = 0; i<mctruth.size(); ++i) {
         std::string mctext;
         bool firstin  = true;
         bool firstout = true;
         std::string origin;
         std::string incoming;
         std::string outgoing;
         // Label cosmic rays -- others are pretty obvious
         if (mctruth[i]->Origin()==simb::kCosmicRay)  origin = "c-ray: ";
         int jmax = TMath::Min(20,mctruth[i]->NParticles());
         for (int j=0; j<jmax; ++j) {
           const simb::MCParticle& p = mctruth[i]->GetParticle(j);
           char buff[1024];
           if (p.P()>0.05) {
             sprintf(buff,"#color[%d]{%s #scale[0.75]{[%.1f GeV/c]}}",
                     Style::ColorFromPDG(p.PdgCode()),
                     Style::LatexName(p.PdgCode()),
                     p.P());
           }
           else {
             sprintf(buff,"#color[%d]{%s}",
                     Style::ColorFromPDG(p.PdgCode()),
                     Style::LatexName(p.PdgCode()));
           }
           if (p.StatusCode()==0) {
             if (firstin==false) incoming += " + ";
             incoming += buff;
             firstin = false;
           }
           if (p.StatusCode()==1) {
             if (firstout==false) outgoing += " + ";
             outgoing += buff;
             firstout = false;
           }
         } // loop on j particles
         if (origin=="" && incoming=="") {
           mctext = outgoing;
         }
         else {
           mctext = origin+incoming+" #rightarrow "+outgoing;
         }
         TLatex& latex = view->AddLatex(0.03, 0.2, mctext.c_str());
         latex.SetTextSize(0.6);
 
       } // loop on i mctruth objects
     }

Member Data Documentation

std::map<int,bool> gar::evd::SimulationDrawer::fHighlite

private

Definition at line 56 of file SimulationDrawer.h.

double gar::evd::SimulationDrawer::maxx

Definition at line 42 of file SimulationDrawer.h.

double gar::evd::SimulationDrawer::maxy

Definition at line 44 of file SimulationDrawer.h.

double gar::evd::SimulationDrawer::maxz

Definition at line 46 of file SimulationDrawer.h.

double gar::evd::SimulationDrawer::minx

Definition at line 41 of file SimulationDrawer.h.

double gar::evd::SimulationDrawer::miny

Definition at line 43 of file SimulationDrawer.h.

double gar::evd::SimulationDrawer::minz

Definition at line 45 of file SimulationDrawer.h.

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

garsoft/EventDisplay/EVD/SimulationDrawer.h
garsoft/EventDisplay/EVD/SimulationDrawer.cxx

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation