SimulationDrawer.cxx

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////////
/// \file    SimulationDrawer.cxx
/// \brief   Render the objects from the Simulation package
/// \author  messier@indiana.edu
/// \version $Id: SimulationDrawer.cxx,v 1.33 2012-12-05 22:14:05 bckhouse Exp $
/////////////////////////////////////////////////////////////////////////////

#include "EventDisplay/SimulationDrawer.h"

#include "TDatabasePDG.h"
#include "TLatex.h"
#include "TLine.h"
#include "TParticle.h"
#include "TPolyMarker.h"
#include "TPolyMarker3D.h"
#include "TPolyLine.h"
#include "TPolyLine3D.h"

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "EventDisplayServices/SimulationDrawingOptions.h"
#include "EventDisplay/Style.h"
#include "EventDisplayServices/SliceNavigator.h"
#include "nuevdb/EventDisplayBase/Colors.h"
#include "nuevdb/EventDisplayBase/View2D.h"
#include "nuevdb/EventDisplayBase/View3D.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/OfflineChan.h"
#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "MCCheater/BackTracker.h"

namespace {
  // Utility function to make uniform error messages.
  void writeErrMsg(const char* fcn,
                   cet::exception const& e)
  {
    mf::LogVerbatim("SimulationDrawer")<< "RawDataDrawer::" << fcn
              << " failed with message:\n"
              << e;
  }
}

namespace evd
{
  //..................................................................
  SimulationDrawer::SimulationDrawer() 
  {
  }

  //......................................................................

  SimulationDrawer::~SimulationDrawer()
  {
  }

  //......................................................................

  void SimulationDrawer::MCTruthShortText(const art::Event& evt,
                                          evdb::View2D*     view) 
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fText & 
          evd::SimulationDrawingOptions::kTEXT_SHORT) == 0) {
      return;
    }

    art::ServiceHandle<evdb::Colors> colors;

    std::vector<simb::MCTruth> mctruth;
    this->GetMCTruth(evt, mctruth);

    // Letting text get too large is a performance problem, and eventually
    // segfaults (don't know why).  Cap at a reasonable size.  There's nothing
    // special about this precise size.
    const unsigned int mctextmax = 2000;
    bool cappedtext = false;

    evd::Style style;
    std::string mctext;
    for (unsigned int i = 0; i < mctruth.size(); ++i) {
      if(i > 0) 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: "; }
      for (int j=0; j<mctruth[i].NParticles(); ++j) {

        const simb::MCParticle& p = mctruth[i].GetParticle(j);
        // Print text in color if there are other MCTruth objects drawn
        int pcolor = ( drawopt->fDraw == 0 ) ? colors->Foreground(0) : style.ColorFromPDG(p.PdgCode() );

        const unsigned int bufsize = 1024;
        char buff[bufsize];
        if (p.P()>0.05) {
          snprintf(buff, bufsize, "#color[%d]{%s #scale[0.75]{[%.1f#scale[0.5]{GeV/c}]}}",
                  pcolor,
                  Style::LatexName(p.PdgCode()), 
                  p.P());
        }
        else {
          snprintf(buff, bufsize, "#color[%d]{%s}",
                  pcolor,
                  Style::LatexName(p.PdgCode()));
        }
        if (p.StatusCode()==0) {
          if(mctext.size() + incoming.size() + 3 + strlen(buff) > mctextmax){
            cappedtext = true;
            break;
          }
          if (!firstin) incoming += " + ";
          incoming += buff;
          firstin = false;
        }
        if (p.StatusCode()==1) {
          if(mctext.size() + outgoing.size() + 3 + strlen(buff) > mctextmax){
            cappedtext = true;
            break;
          }
          if (!firstout) outgoing += " + ";
          outgoing += buff;
          firstout = false;
        }
      } // loop on j particles
      if (origin=="" && incoming=="") {
        if(mctext.size() + outgoing.size() > mctextmax){
          cappedtext = true;
          break;
        }
        mctext += outgoing;
      }
      else {
        const int intType = mctruth[i].GetNeutrino().InteractionType();
        const std::string suffix = ShortInteractionSuffix(intType);
        const std::string newtext = origin+incoming+" #rightarrow "+outgoing+suffix;
        if(mctext.size() + newtext.size() > mctextmax){
          cappedtext = true;
          break;
        }
        mctext += newtext;
      }
      if(cappedtext) break; // if capped inside particle loop
    } // loop on i mctruth objects

    if(cappedtext){
      if(mctext.empty()) mctext = "Too many particles"; // shouldn't happen
      else mctext += "... and more!";
    }

    TLatex& latex = view->AddLatex(0.01, 0.055, mctext.c_str());
    latex.SetTextColor(colors->Foreground(0));
    latex.SetTextSize(0.15);
  }

  //......................................................................

  // 'prefix' is what we need to add to our children in addition to what we draw
  void SimulationDrawer::PrintParticleAndOffspring(std::stringstream& ss,
                                                   const sim::Particle* part,
                                                   std::string prefix,
                                                   int depth) const
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    const int depthLimit = drawopt->fTextDepthLimit;
    const bool showVtxs  = drawopt->fTextIncludeVertex;
    const bool showDirs  = drawopt->fTextIncludeDirections;

    art::ServiceHandle<cheat::BackTracker> bt;
    const sim::ParticleNavigator& pnav = bt->ParticleNavigator();

    // Find out what the particle is called and print it. Fallback is PDG code
    const TDatabasePDG* db = TDatabasePDG::Instance();
    const TParticlePDG* defn = db->GetParticle(part->PdgCode());
    if(defn){
      ss << defn->GetName();
    }
    else if(part->PdgCode() > 1000000000){
      const int strange = (part->PdgCode()/10000000)%10;
      const int element = (part->PdgCode()/10000)%1000;
      const int A       = (part->PdgCode()/10)%1000;
      const int I       =  part->PdgCode()%10;
      switch(element){
      case  1: ss << "H"; break;
      case  2: ss << "He"; break;
      case  3: ss << "Li"; break;
      case  4: ss << "Be"; break;
      case  5: ss << "B"; break;
      case  6: ss << "C"; break;
      case  7: ss << "N"; break;
      case  8: ss << "O"; break;
      case  9: ss << "F"; break;
      case 10: ss << "Ne"; break;
      case 11: ss << "Na"; break;
      case 12: ss << "Mg"; break;
      case 13: ss << "Al"; break;
      case 14: ss << "Si"; break;
      case 15: ss << "P"; break;
      case 16: ss << "S"; break;
      case 17: ss << "Cl"; break;
      case 18: ss << "Ar"; break;
      case 19: ss << "K"; break;
      case 20: ss << "Ca"; break;
      case 21: ss << "Sc"; break;
      case 22: ss << "Ti"; break;
      case 23: ss << "V"; break;
      case 24: ss << "Cr"; break;
      case 25: ss << "Mn"; break;
      case 26: ss << "Fe"; break;
      case 27: ss << "Co"; break;
      case 49: ss << "In"; break;
      case 50: ss << "Sn"; break;
      case 51: ss << "Sb"; break;
      case 52: ss << "Te"; break;
      case 53: ss << "I"; break;
      case 54: ss << "Xe"; break;
      case 55: ss << "Cs"; break;
      case 56: ss << "Ba"; break;
      case 57: ss << "La"; break;
      default: ss << element << "-";
      }

      ss << A;

      // Excited nuclear state.  If it became useful, we could distinguish
      // between several different ones.
      if(I != 0) ss << "*";

      const char* Lambda  = "\u039b";
      for(int i = 0; i < strange; i++) ss << Lambda;
    }
    else{
      // Anti-nuclei will land here, among perhaps some other obscure particles
      ss << part->PdgCode();
    }

    // Print the kinetic energy, and the process that created this particle
    const double ke = part->E()-part->Mass();
    char Estr[1024];
    if(ke > 1)
      sprintf(Estr, "%.3lf GeV", ke);
    else if(ke > 1e-3)
      sprintf(Estr, "%.3lf MeV", 1e3*ke);
    else
      sprintf(Estr, "%.3lf keV", 1e6*ke);

    ss << "\t(" << Estr << ")\t[" << part->Process() << "]";
    if(showVtxs){
      ss << " at (" << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << "; " << part->T() << ") ";
    }
    if(showDirs){
      double ptot=part->P();
      if (showVtxs) ss<<",";
      if (ptot>1e-6) {
         ss << " p-hat (" << part->Px()/ptot << ", " << part->Py()/ptot << ", " << part->Pz()/ptot << ") ";
      } else {
         ss << " p=0 ";
      }
    }
    ss << "\t" << part->TrackId() << std::endl;

    // Reached depth limit, shouldn't print daughters
    if(depth >= depthLimit && depthLimit > 0) return; 
    // Unicode "light" box drawing characters
    const char* kHoriz  = "\u2500";
    const char* kVert   = "\u2502";
    const char* kCorner = "\u2514";
    const char* kTee    = "\u251c";

    // Print our daughters too, indented
    const int N = part->NumberDaughters();
    for(int n = 0; n < N; ++n){
      sim::ParticleNavigator::const_iterator childIt =  pnav.find(part->Daughter(n));
    
      // Should always be there because we created the daughter links ourselves
      // but better safe than sorry...
      if(childIt == pnav.end()) continue;

      if(n == N-1){
        ss << prefix << kCorner << kHoriz << kHoriz;
        PrintParticleAndOffspring(ss, childIt->second,
                                  prefix+"     ", depth+1);
      }
      else{
        ss << prefix << kTee << kHoriz << kHoriz;
        PrintParticleAndOffspring(ss, childIt->second,
                                  prefix+kVert+"  ", depth+1);
      }
    }

  }

  //......................................................................

  void SimulationDrawer::MCTruthLongText(const art::Event& evt,
                                         evdb::View2D* /*view*/) 
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    // Skip drawing if option is turned off
    if ((drawopt->fText & 
        evd::SimulationDrawingOptions::kTEXT_LONG)==0) {
      return;
    }

    art::ServiceHandle<cheat::BackTracker> bt;
    const sim::ParticleNavigator& pnav = bt->ParticleNavigator();

    const bool showVtxs  = drawopt->fTextIncludeVertex;
    const bool showDirs  = drawopt->fTextIncludeDirections;

    std::stringstream ss;
    ss<< "Long MCTruth text ( from sim::Particle )"<< std::endl;
    ss<<" Particle       ( KineticEnergy )    [ Process ] ";
    if (showVtxs) ss<<" at  (Vx, Vy, Vz, T)";
    if (showDirs) ss<<" with direction  (Px, Py, Pz)";
    ss<<" TrackId"<<std::endl;

    // loop over all the primary particles
    for(int n = 0; n < pnav.NumberOfPrimaries(); ++n){
      ss << "\u2500\u2500";
      PrintParticleAndOffspring(ss, pnav.Primary(n), "  ", 1);
    }
    mf::LogVerbatim("SimulationDrawer") << ss.str();
  }

  //......................................................................
  
  void SimulationDrawer::MCTruthVertices2D(const art::Event& evt,
                                           evdb::View2D*     xzview,
                                           evdb::View2D*     yzview) 
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw & 
          evd::SimulationDrawingOptions::kDRAW_VERTEX)==0) {
      return;
    }
    std::vector<double> x;
    std::vector<double> y;
    std::vector<double> z;
    this->GetVertexPoints(evt, x, y, z);
    unsigned int i;
    int    gold = kOrange-4;
    int    star = 30;
    double sz   = 2.0;
    for (i=0; i<x.size(); ++i) {
      if(xzview) xzview->AddMarker(z[i], x[i], gold, star, sz);
      if(yzview) yzview->AddMarker(z[i], y[i], gold, star, sz);
    }
 }
  
  //......................................................................

  void SimulationDrawer::MCTruthVertices3D(const art::Event& evt,
                                           evdb::View3D*     view)
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw & 
          evd::SimulationDrawingOptions::kDRAW_VERTEX)==0) {
      return;
    }
    std::vector<double> x;
    std::vector<double> y;
    std::vector<double> z;
    this->GetVertexPoints(evt, x, y, z);

    int    gold = kOrange-4;
    int    star = 30;
    double sz   = 2.0;
    TPolyMarker3D& pm = view->AddPolyMarker3D(x.size(), gold, star, sz);
    
    unsigned int i;
    for (i=0; i<x.size(); ++i) {
      pm.SetPoint(i, x[i], y[i], z[i]);
    }
  }

  //......................................................................
  
  void SimulationDrawer::MCTruthVectors2D(const art::Event& evt,
                                          evdb::View2D*     xzview,
                                          evdb::View2D*     yzview) 
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw & 
          evd::SimulationDrawingOptions::kDRAW_VECTORS)==0) {
      return;
    }

    // Unpack and draw the MC vectors
    std::vector<simb::MCTruth> mctruth;
    this->GetMCTruth(evt, mctruth);

    evd::Style style;
    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);
      
        // Skip all but incoming and out-going particles
        if (!(p.StatusCode()==0 || p.StatusCode()==1)) continue;
      
        double r  = p.P()*100.0;          // Scale length so 1 m = 1 GeV/c
        if (r<0.1) continue;              // Skip very short particles
        if (p.StatusCode()==0) r = -r;    // Flip for incoming particles
      
        double x0 = p.Vx();
        double y0 = p.Vy();
        double z0 = p.Vz();
        double x1 = x0 + r * p.Px()/p.P();
        double y1 = y0 + r * p.Py()/p.P();
        double z1 = z0 + r * p.Pz()/p.P();
        if (xzview!=0) {
          TLine& l = xzview->AddLine(z0, x0, z1, x1);
          style.FromPDG(&l,p.PdgCode());
        }
        if (yzview!=0) {
          TLine& l = yzview->AddLine(z0, y0, z1, y1);
          style.FromPDG(&l,p.PdgCode());
        }
      } // loop on j particles
    } // loop on i mctruth objects
  }

  //......................................................................

  void SimulationDrawer::MCTruthVectors3D(const art::Event& evt,
                                          evdb::View3D*     view)
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw &
          evd::SimulationDrawingOptions::kDRAW_VECTORS)==0) {
      return;
    }
    
    // Unpack and draw the MC vectors
    std::vector<simb::MCTruth> mctruth;
    this->GetMCTruth(evt, mctruth);
    
    evd::Style style;
    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);
        
        // Skip all but incoming and out-going particles
        if (!(p.StatusCode()==0 || p.StatusCode()==1)) continue;
      
        double r  = p.P()*100.0;          // Scale length so 1 m = 1 GeV/c
        if (r<0.1) continue;              // Skip very short particles
        if (p.StatusCode()==0) r = -r; // Flip for incoming particles
      
        double x0 = p.Vx();
        double y0 = p.Vy();
        double z0 = p.Vz();
        double x1 = x0 + r * p.Px()/p.P();
        double y1 = y0 + r * p.Py()/p.P();
        double z1 = z0 + r * p.Pz()/p.P();
      
        int c = style.ColorFromPDG(p.PdgCode());
        int s = style.LineStyleFromPDG(p.PdgCode());
        int w = style.LineWidthFromPDG(p.PdgCode());
        TPolyLine3D& l = view->AddPolyLine3D(2,c,w,s);
        l.SetPoint(0, x0, y0, z0);
        l.SetPoint(1, x1, y1, z1);
      } // loop on j particles
    } // loop on i mctruth objects
  }

  //......................................................................

  void SimulationDrawer::MCTruthTrajectoriesAnyD(const art::Event& evt,
                                                 evdb::View2D*     xzview,
                                                 evdb::View2D*     yzview,
                                                 evdb::View3D*     view3D)
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw &
          evd::SimulationDrawingOptions::kDRAW_TRAJECT)==0) {
      return;
    }
    
    // TODO: Seems this call is required for GetParticleList to work properly?
    std::vector<simb::MCTruth> mctruth;
    this->GetMCTruth(evt, mctruth);

    art::ServiceHandle<cheat::BackTracker> bt;
    const sim::ParticleNavigator& pnav = bt->ParticleNavigator();

    evd::Style style;
    for(sim::ParticleNavigator::const_iterator j=pnav.begin(); 
        j != pnav.end();
        ++j) {
      const sim::Particle part = *((*j).second);
      
      const int pdg = part.PdgCode();
      
      // Decay neutrinos are just messy
      if (abs(pdg) == 12 || abs(pdg) == 14 || abs(pdg) == 16) continue;
      
      const int N = part.NumberTrajectoryPoints();
      const int c = style.ColorFromPDG(pdg);
      const int w = style.LineWidthFromPDG(pdg);
      const int s = style.LineStyleFromPDG(pdg);
      
      // This drawing mode disables neutral particles
      bool draw_neutrals = (drawopt->fDraw &
                            evd::SimulationDrawingOptions::kDRAW_NEUTRALS);
      //      if (!draw_neutrals && s!=kSolid) continue;
      bool draw_gammas   = (drawopt->fDraw &
                            evd::SimulationDrawingOptions::kDRAW_GAMMAS);
      if (!draw_neutrals && s==kDotted) continue;
      if (!draw_gammas && s==kDashed) continue;
      
      TPolyLine* xz = xzview ? &xzview->AddPolyLine(N, c, w, s) : 0;
      TPolyLine* yz = yzview ? &yzview->AddPolyLine(N, c, w, s) : 0;
      TPolyLine3D* xyz = view3D ? &view3D->AddPolyLine3D(N, c, w, s) : 0;
      
      for(int n = 0; n < N; ++n){
        const TLorentzVector& vec = part.Position(n);
        if(xz) xz->SetPoint(n, vec.Z(), vec.X());
        if(yz) yz->SetPoint(n, vec.Z(), vec.Y());
        if(xyz) xyz->SetPoint(n, vec.X(), vec.Y(), vec.Z());
      }
      
      if(xz) xz->SetBit(kCannotPick);
      if(yz) yz->SetBit(kCannotPick);
      if(xyz) xyz->SetBit(kCannotPick);
    } // end for j

  }

  void SimulationDrawer::MCTruthTrajectories2D(const art::Event& evt,
                                               evdb::View2D*     xzview,
                                               evdb::View2D*     yzview)
  {
    MCTruthTrajectoriesAnyD(evt, xzview, yzview, 0);
  }

  void SimulationDrawer::MCTruthTrajectories3D(const art::Event& evt,
                                               evdb::View3D*     view)
  {
    MCTruthTrajectoriesAnyD(evt, 0, 0, view);
  }

  //......................................................................

  void SimulationDrawer::FLSHit2D(const art::Event& evt,
                                  evdb::View2D*     xzview,
                                  evdb::View2D*     yzview)
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw & 
          evd::SimulationDrawingOptions::kDRAW_HITS)==0) {
      return;
    }

    art::ServiceHandle<geo::Geometry> geo;

    // Otherwise, unpack and show the FLS hits
    for (unsigned int imod=0; imod<drawopt->fFLSHitListModules.size(); ++imod) {
      const char* which = drawopt->fFLSHitListModules[imod].c_str();
    
    std::vector<sim::FLSHitList> flshits;
    this->GetFLSHits(evt, which, flshits);

    // Indexed by color
    std::map<int, TPolyMarker*> pms_x, pms_y;
    evd::Style style;

    for (unsigned int i=0; i<flshits.size(); ++i) {
      for (unsigned int j=0; j<flshits[i].fHits.size(); ++j) {
        const sim::FLSHit& hit = flshits[i].fHits[j];
        if (hit.GetEdep()<0.001*drawopt->fFLSHitThresh) continue;

        switch (drawopt->fFLSHitStyle) {
        case evd::SimulationDrawingOptions::kFLSHIT_AS_LINE:
          for(int view = 0; view < 2; ++view){
            evdb::View2D* v2d = (view == 0) ? xzview : yzview;
            if(!v2d) continue;

            // Dummy styling, overridden in next line
            TPolyLine& l = v2d->AddPolyLine(hit.GetNPoints(), kBlack, 0, kSolid);
            style.FromPDG(&l, hit.GetPDG());
            for(int i = 0; i < hit.GetNPoints(); ++i){
              double xyzloc[3] = {hit.GetX(i), hit.GetY(i), hit.GetZ(i)};
              double xyzwor[3];
              geo->Plane(hit.GetPlaneID())->Cell(hit.GetCellID())->LocalToWorld(xyzloc, xyzwor);
              l.SetPoint(i, xyzwor[2], xyzwor[view]);
            }
          }
          break;
        default:
          double xyzloc[3] = {hit.GetXAverage(), hit.GetYAverage(), hit.GetZAverage()};
          double xyzwor[3];
          geo->Plane(hit.GetPlaneID())->Cell(hit.GetCellID())->LocalToWorld(xyzloc, xyzwor);

          const int c = style.ColorFromPDG(hit.GetPDG());
          const double sz = 0.5;
          
          if(xzview){
            TPolyMarker* pm = pms_x[c];
            if(!pm){
              pm = &xzview->AddPolyMarker(0, c, kFullDotMedium, sz);
              pms_x[c] = pm;
            }
            pm->SetNextPoint(xyzwor[2], xyzwor[0]);
          }
          if(yzview){
            TPolyMarker* pm = pms_y[c];
            if(!pm){
              pm = &yzview->AddPolyMarker(0, c, kFullDotMedium, sz);
              pms_y[c] = pm;
            }
            pm->SetNextPoint(xyzwor[2], xyzwor[1]);
          }
          
          if(fHighlite[hit.GetTrackID()]){
            if(xzview) xzview->AddMarker(xyzwor[2],xyzwor[0],c,kOpenCircle,sz);
            if(yzview) yzview->AddMarker(xyzwor[2],xyzwor[1],c,kOpenCircle,sz);
          }

        } // FLSHit style
      } // loop on j hits in list
    } // loop on i hit lists
  } // loop o imod folders
}

  //......................................................................

  void SimulationDrawer::FLSHit3D(const art::Event& evt,
                                  evdb::View3D*     view)
  {
    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;
    if ( (drawopt->fDraw &
          evd::SimulationDrawingOptions::kDRAW_HITS)==0) {
      return;
    }

    art::ServiceHandle<geo::Geometry> geo;

    evd::Style style;

    // Otherwise, unpack and show the FLS hits
    for (unsigned int imod=0; imod<drawopt->fFLSHitListModules.size(); ++imod) {
      const char* which = drawopt->fFLSHitListModules[imod].c_str();
    
      std::vector<sim::FLSHitList> flshits;
      this->GetFLSHits(evt, which, flshits);
      for (unsigned int i=0; i<flshits.size(); ++i) {
        for (unsigned int j=0; j<flshits[i].fHits.size(); ++j) {
          const sim::FLSHit& hit = flshits[i].fHits[j];

          if (hit.GetEdep()<0.001*drawopt->fFLSHitThresh) continue;

          double xyzloc[3]; // Local cell xyz coordinates
          double xyzwor[3]; // World coordinates
          xyzloc[0] = 0.5*(hit.GetEntryX() + hit.GetExitX());
          xyzloc[1] = 0.5*(hit.GetEntryY() + hit.GetExitY());
          xyzloc[2] = 0.5*(hit.GetEntryZ() + hit.GetExitZ());
        
          geo->Plane(hit.GetPlaneID())->Cell(hit.GetCellID())->LocalToWorld(xyzloc,xyzwor);
        
          int   st = 7;
          int   c  = style.ColorFromPDG(hit.GetPDG());
          float sz = 0.5;

          TPolyMarker3D& p = view->AddPolyMarker3D(1, c, st, sz);
          p.SetPoint(0, xyzwor[0], xyzwor[1], xyzwor[2]);

          if (fHighlite[hit.GetTrackID()]) {
            TPolyMarker3D& p = view->AddPolyMarker3D(1, c, 24, 2*sz);
            p.SetPoint(0, xyzwor[0], xyzwor[1], xyzwor[2]);       
          }

        } // loop on j hits in list
      } // loop on i hit lists
    } // loop o imod folders
  }


  //......................................................................

  int SimulationDrawer::GetMCTruth(const art::Event& evt,
                                   std::vector<simb::MCTruth>& mctruth)
  {
    mctruth.clear();

    art::ServiceHandle<evd::SimulationDrawingOptions> opt;

    for(const std::string& mod: opt->fMCTruthModules){
      art::ServiceHandle<evd::SliceNavigator> nav;
      std::vector<const simb::MCTruth*> truths;
      nav->GetProducts(evt, mod, "", truths);

      for(const simb::MCTruth* t: truths) mctruth.push_back(*t);
    }

    return mctruth.size();
  }

  //......................................................................

  int SimulationDrawer::GetFLSHits(const art::Event& evt,
                                   const char* which,
                                   std::vector<sim::FLSHitList>& flshits) 
  {

    flshits.clear();
    art::Handle< std::vector<sim::FLSHitList> > fcol;

    try{
      evt.getByLabel(which, fcol);
      std::vector<sim::FLSHitList> temp(fcol->size());
      for(size_t f = 0; f < fcol->size(); ++f){
        flshits.push_back(fcol->at(f));
      }
    }
    catch(cet::exception& e){
      writeErrMsg("GetFLSHits", e);
    }

    return flshits.size();
  }

  //......................................................................
  
  int SimulationDrawer::GetVertexPoints(const art::Event&    evt,
                                        std::vector<double>& x,
                                        std::vector<double>& y,
                                        std::vector<double>& z)
  {
    x.clear();
    y.clear();
    z.clear();
    //
    // The list of vertexs are attached to individual particles. To
    // get a unique list of vertexs, loop over the particles and keep
    // only the unique vertex points.
    //
    std::vector<simb::MCTruth> mctruth;
    this->GetMCTruth(evt, mctruth);
    for (unsigned int i=0; i<mctruth.size(); ++i) {
      double x0, y0, z0;
      bool   in_list;
      for (int j=0; j<mctruth[i].NParticles(); ++j) {
        const simb::MCParticle& p = mctruth[i].GetParticle(j);
        //
        // Skip all but incoming and out-going particles
        //
        if (!(p.StatusCode()==0 || p.StatusCode()==1)) continue;
        x0 = p.Vx();
        y0 = p.Vy();
        z0 = p.Vz();
        //
        // Check if this vertex is already in the list. 
        //
        in_list = false;
        for (unsigned int k=0; k<x.size(); ++k) {
          in_list = ((x0 == x[k]) && (y0 == y[k]) && (z0 == z[k]));
          if (in_list) break;
        }
        if (in_list==false) {
          x.push_back(x0);
          y.push_back(y0);
          z.push_back(z0);
        }
      }
    }
    return x.size();
  }
  
  //......................................................................

  std::string SimulationDrawer::ShortInteractionSuffix(int iType) const
  {
    switch(iType){
    case simb::kQE:
    case simb::kCCQE:
      return " (QE)";

    case simb::kRes: 
    case simb::kResCCNuProtonPiPlus:
    case simb::kResCCNuNeutronPi0:
    case simb::kResCCNuNeutronPiPlus:
    case simb::kResCCNuBarNeutronPiMinus:
    case simb::kResCCNuBarProtonPi0:
    case simb::kResCCNuBarProtonPiMinus:
      return " (RES)";

    case simb::kDIS: 
    case simb::kCCDIS:
      return " (DIS)";

    case simb::kCoh: 
    case simb::kCCCOH:
      return " (COH)";

    case simb::kNCQE:
    case simb::kNCDIS:  
    case simb::kResNCNuProtonPi0:
    case simb::kResNCNuProtonPiPlus:
    case simb::kResNCNuNeutronPi0:
    case simb::kResNCNuNeutronPiMinus:
    case simb::kResNCNuBarProtonPi0:
    case simb::kResNCNuBarProtonPiPlus:
    case simb::kResNCNuBarNeutronPi0:
    case simb::kResNCNuBarNeutronPiMinus:
      return " (NC)";

    case simb::kElectronScattering:
      return " (ES)";

    case simb::kInverseBetaDecay:
      return " (IBD)";

    default: // simb::kNuElectronElastic and simb::kInverseMuDecay
      return "";
    }
  }

  //......................................................................

  void SimulationDrawer::HiLite(int trkId, bool dohilite)
  {
    fHighlite[trkId] = dohilite;
  }

  //......................................................................

  void SimulationDrawer::GetLimits(const art::Event* evt,
                                   double& xmin, double& xmax,
                                   double& ymin, double& ymax,
                                   double& zmin, double& zmax,
                                   const std::set<geo::OfflineChan>& hmap)
  {
    assert(evt);

    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;

    art::ServiceHandle<geo::Geometry> geom;

    xmin = +1e10; ymin = +1e10; zmin = +1e10;
    xmax = -1e10; ymax = -1e10; zmax = -1e10;

    for(unsigned int imod = 0; imod < drawopt->fFLSHitListModules.size(); ++imod){
      std::vector<sim::FLSHitList> lists;
      GetFLSHits(*evt, drawopt->fFLSHitListModules[imod].c_str(), lists);

      for(unsigned int i = 0; i < lists.size(); ++i){
        const std::vector<sim::FLSHit>& hits = lists[i].fHits;

        for(unsigned int j = 0; j < hits.size(); ++j){
          const sim::FLSHit& hit = hits[j];

          if(!hmap.empty() &&
             hmap.find(geo::OfflineChan(hit.GetPlaneID(),
                                        hit.GetCellID())) == hmap.end()) continue;

          geo::View_t view;
          double xyz[3];
          geom->CellInfo(hit.GetPlaneID(), hit.GetCellID(), &view, xyz);

          if(view == geo::kX && xyz[0] < xmin) xmin = xyz[0];
          if(view == geo::kX && xyz[0] > xmax) xmax = xyz[0];
          if(view == geo::kY && xyz[1] < ymin) ymin = xyz[1];
          if(view == geo::kY && xyz[1] > ymax) ymax = xyz[1];
          if(xyz[2] < zmin) zmin = xyz[2];
          if(xyz[2] > zmax) zmax = xyz[2];
        }
      }
    } // end for imod
  }

  //......................................................................

  void SimulationDrawer::GetTimeLimits(const art::Event* evt,
                                       double& tmin, double& tmax)
  {
    tmin = FLT_MAX;
    tmax = -FLT_MAX;

    assert(evt);

    art::ServiceHandle<evd::SimulationDrawingOptions> drawopt;

    for(unsigned int imod = 0; imod < drawopt->fFLSHitListModules.size(); ++imod){
      std::vector<sim::FLSHitList> lists;
      GetFLSHits(*evt, drawopt->fFLSHitListModules[imod].c_str(), lists);

      for(unsigned int i = 0; i < lists.size(); ++i){
        const std::vector<sim::FLSHit>& hits = lists[i].fHits;

        for(unsigned int j = 0; j < hits.size(); ++j){
          const sim::FLSHit& hit = hits[j];

          if(hit.GetEntryT() < tmin) tmin = hit.GetEntryT();
          if(hit.GetExitT()  > tmax) tmax = hit.GetExitT();
        }
      }
    }
  }

} //end namespace evd
////////////////////////////////////////////////////////////////////////