In This Package:

DrawHistoryAlg.cc

Go to the documentation of this file.

#include "DrawHistoryAlg.h"
#include "Event/SimHeader.h"
#include "Event/SimParticleHistory.h"
#include "Event/SimTrack.h"
#include "Event/SimVertex.h"
#include "Event/SimHit.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/SimPmtHit.h"
#include "Event/SimRpcHit.h"

#include "DetDesc/IDetectorElement.h"
#include "DetDesc/GeometryInfoPlus.h"
#include "DetDesc/IPVolume.h"
#include "GaudiKernel/DeclareFactoryEntries.h"
#include "GaudiKernel/ParticleProperty.h"
#include "GaudiKernel/IParticlePropertySvc.h"

#include "CLHEP/Units/SystemOfUnits.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"

#include <fstream>
#include <sstream>
#include <string>

DECLARE_ALGORITHM_FACTORY(DrawHistoryAlg);

using namespace DayaBay;

DrawHistoryAlg::DrawHistoryAlg(const std::string& name, ISvcLocator* pSvcLocator)
    : GaudiAlgorithm(name,pSvcLocator)
{
    declareProperty("Location",m_location= DayaBay::SimHeaderLocation::Default,
                     "Location in the TES to get Event data");
    declareProperty("track_filename",m_track_filename = "tracks_%d.dot",
                    "Filename to write track structure. Use '\%d' to indicate event number.");
    declareProperty("trackandvertex_filename",m_trackandvertex_filename = "tracks_and_vertices_%d.dot",
                    "Filename to write track structure. Use '\%d' to indicate event number.");
    declareProperty("do_hits",m_do_hits = true,
                    "1 to make hits in track plot, 0 to ignore them.");
    declareProperty("Geometry",m_topGeoPath = "/dd/Structure/DayaBay",
                    "Top level geometry structure element.");
}

DrawHistoryAlg::~DrawHistoryAlg()
{
}

StatusCode DrawHistoryAlg::initialize()
{
 // service retrieval
  m_ppSvc =0;
  StatusCode sc = service("ParticlePropertySvc",m_ppSvc,true);
  if(sc.isFailure()){
    err() << "Cant' get a ParticlePropertySvc pointer." << endreq;
    return StatusCode::FAILURE;
  }  

  return StatusCode::SUCCESS;
}

StatusCode DrawHistoryAlg::finalize()
{
  return StatusCode::SUCCESS;
}

std::string DrawHistoryAlg::particleName(int pdg)
{
  ParticleProperty* p = m_ppSvc->findByStdHepID(pdg);
  if(p) return p->particle();
  std::ostringstream oss;
  oss << "[" << pdg << "]";
  return oss.str();
}

bool DrawHistoryAlg::volumeInfo(SimVertex& v, std::string& name, std::string& material)
{
    static IDetectorElement* topDet = getDet<IDetectorElement>(m_topGeoPath);
    static IDetectorElement* lastDet = 0;

    // Perfer lastDet to topDet
    IDetectorElement* start = lastDet ? lastDet : topDet;
    std::string child_de_path = start->geometry()->belongsToPath(v.position(),-1);
    // If failed with lastDet, try again with topDet
    if ("" == child_de_path && lastDet) {
        child_de_path = topDet->geometry()->belongsToPath(v.position(),-1);
    }
    IDetectorElement* child = getDet<IDetectorElement>(child_de_path);


    // If still failed then game over
    if (!child) return false;

    // Optimizing cache, next vertex likely in this same DE
    lastDet = child;

    name = child->name();
    material = child->geometry()->lvolume()->materialName();
    return true;
}

StatusCode DrawHistoryAlg::execute()
{
  using std::endl;
  
  DayaBay::SimHeader* header = 0;
  if (exist<DayaBay::SimHeader>(evtSvc(),m_location)) {
     header = get<DayaBay::SimHeader>(m_location);
  }
  if(!header) {
    warning() << "No SimHeader found." <<endreq;
    return StatusCode::SUCCESS;
  }

  const SimParticleHistory* h = header->particleHistory();
  if(!h) {
    warning() << "No SimParticleHistory found!" << endreq;
    return StatusCode::SUCCESS;
  }

  
  // Dot file of vertex and tracks.

  // sort filename
  std::string filename = m_trackandvertex_filename;
  {
    std::string::size_type dpos = filename.find("%d");
    if(dpos!=std::string::npos) {
      std::ostringstream oss; 
      oss << header->execNumber();
      filename.replace(dpos,2,oss.str());
    }
  }
  
  std::ofstream dot(filename.c_str());
  dot << "digraph G {" << endl;
  const std::list<SimTrack*>& tracks = h->tracks();
  std::list<SimTrack*>::const_iterator tit;
  for(tit = tracks.begin();tit!= tracks.end(); ++tit) {
    const SimTrack* t= *tit;
    const std::vector<SimVertex*>& vertices = t->vertices();
    dot << "// SimTrack id " << t->trackId() << endl;
    // Name the first vertex.
    dot << "v" << (long)(vertices[0]) 
               << " [label=\""
               << "SimTrack Start " << t->trackId()
               <<"\"]" << endl;

    // point back to parent vertex.
    // SimVertex* parent = t->parentVertex().vertex();
    // if(parent)  {
    //   dot << "v" << (int)(parent) <<" -> v" << (int)(vertices[0]);
    //   if(t->parentVertex().isDirect()) dot << " [ style = bold] ";
    //   else                             dot << " [ style = dotted ]";
    //   dot << ";" << endl;
    // }
    
    // Now the track chain.
    dot << "subgraph cluster" << t->trackId() <<" {" << endl;
    dot << "  label = \"SimTrack " << t->trackId()
        << "\\n" << particleName(t->particle())
        << "\\nparent=" << t->parentParticle()
        << "\\n" << "KE=" << t->vertices()[0]->kineticEnergy()/CLHEP::MeV << " MeV";
      for(SimTrack::descendants_map::const_iterator udit = t->unrecordedDescendants().begin();
          udit != t->unrecordedDescendants().end();
          ++udit) {
          dot << "\\n" << udit->second << " skipped of type " << particleName(udit->first);
      }
      dot << "\";" << endl;
    
    assert(vertices.size()>0);
         
    // The chain of vertices along the track.
    for(unsigned int iv = 1; iv<vertices.size(); iv++) {
      dot << "  v" << long(vertices[iv-1]) << "  -> v" << long(vertices[iv]);
      // Compute distance, dE,dx.
      double dx = (vertices[iv]->position()-vertices[iv-1]->position()).R();
      double dE = (vertices[iv]->totalEnergy()-vertices[iv-1]->totalEnergy());
      dot << "[label=\""
          << "dx=" << dx/CLHEP::cm << " cm"
          << "\\ndE=" << dE/CLHEP::MeV << " MeV"
          << "\"]";
      dot << ";" << endl;
    }

    dot << "}" << endl;
    // Do secondaries.
    for(unsigned int iv = 0; iv<vertices.size(); iv++) {
      const SimVertex* v = vertices[iv];
      const std::vector<SimTrackReference>& secs = v->secondaries();
      for(unsigned int isec = 0; isec < secs.size(); ++ isec){
        const SimTrack*  trk = secs[isec].track();
        if(trk) {
          dot << "v" << (long)v << " -> v" << (long)(trk->vertices()[0]);
          if(secs[isec].isIndirect())
            dot << "[ style = dotted ]";
          dot << ";" << endl;
        }
      }
    }    
  }
  
  // Just declare each vertex, see if there are wierd ones around.
  const std::list<SimVertex*>& verts = h->vertices();
  dot << endl << "// Dump of all " << verts.size() << " vertices" << endl;
  std::list<SimVertex*>::const_iterator vit;
  for(vit = verts.begin(); vit!= verts.end(); ++vit) {
    SimVertex* v = *vit;
    dot << "v" << (long)(v);  
    dot << " [label=\""; 
    dot << v->kineticEnergy()/CLHEP::MeV << " MeV";

    std::string name, material;
    if (volumeInfo(*v,name,material)) {
        dot << "\\n" << name;
        dot << "\\n" << material;
    }
    dot << "\"]" << endl;
    dot << ";" << endl;
  }
  dot << "}" << endl;
  dot.close();
  
  // Dot file of tracks, hits, and primaries.
  
  // sort filename
  filename = m_track_filename;
  {
    std::string::size_type dpos = filename.find("%d");
    if(dpos!=std::string::npos) {
      std::ostringstream oss;
      oss << header->execNumber();
      filename.replace(dpos,2,oss.str());
    }
  }
  
  dot.open(filename.c_str());
  dot << "digraph G {" << endl;

  // Create hooks from primary particles to first tracks.
  std::set<const HepMC::GenEvent*> primaryEvents;
  const std::list<const SimTrack*>& primaries = h->primaryTracks();
  std::list<const SimTrack*>::const_iterator it;
  for(it = primaries.begin();it!= primaries.end(); ++it) {
    // Get the primary.
    const HepMC::GenParticle* primaryParticle =  (*it)->primaryParticle();
    const HepMC::GenEvent* gev = primaryParticle->parent_event();
    primaryEvents.insert(gev);
    dot << "  primaryEvent" << (long)(gev) << ":genpar" << (long)(primaryParticle);
    dot << " -> track" << (long)(*it) << endl; 
  }
  
  // create record-like nodes for the primary events.
  std::set<const HepMC::GenEvent*>::iterator gevit;
  for(gevit = primaryEvents.begin() ; gevit != primaryEvents.end(); gevit++){
    const HepMC::GenEvent* gev = *gevit;
    dot << "  primaryEvent" << (long)(gev) << " [shape=record,rank=source,label=\"{";
    dot << "event " << gev->event_number() <<"\\nprocess_id=" << gev->signal_process_id() << "| {";
    
    HepMC::GenEvent::particle_const_iterator parit = gev->particles_begin();
    while(parit != gev->particles_end()) {
      const HepMC::GenParticle* gpar = *parit;
      HepMC::FourVector momentum(gpar->momentum());
      double pmag = HepMC::ThreeVector(momentum).mag();
      double ke = momentum.e() - sqrt(momentum.e()*momentum.e() - pmag*pmag);
      dot << "<" << "genpar" << (long)(gpar) << "> ";
      dot << particleName(gpar->pdg_id())
          << "\\n KE=" << ke/CLHEP::MeV << " MeV";
      parit++;
      if(parit != gev->particles_end()) dot << " | ";
    }
    
    dot << "} }\"];" << endl;    
  }
   

  // Now the actual tracks.
  dot << "subgraph cluster_of_the_tracks {" << endl;
  for(tit = tracks.begin();tit!= tracks.end(); ++tit) {
    const SimTrack* t= *tit;
    dot << "track" << (long)(t)
        << "[label=\""
        << "SimTrack " << t->trackId()
        << "\\n" << particleName(t->particle())
        //<< "\\nparent=" << t->parentParticle()
        << "\\n" << t->vertices()[0]->process().type() << " " <<  t->vertices()[0]->process().name()
        << "\\n" << "KE=" << t->vertices()[0]->kineticEnergy()/CLHEP::MeV << " MeV"
        << "\\n" << "with " << t->vertices().size() << " vertices";
      for(SimTrack::descendants_map::const_iterator udit = t->unrecordedDescendants().begin();
          udit != t->unrecordedDescendants().end();
          ++udit) {
          dot << "\\n" << udit->second << " skipped of type " << particleName(udit->first);
       }        
    dot << "\"];" << endl;

    if(t->ancestorTrack().track()) {
      dot << "track" << (long)(t->ancestorTrack().track()) << " -> track" <<  (long)(t);
      if(! t->ancestorTrack().isDirect())
        dot << " [ style = dotted label=\"" << t->ancestorTrack().indirection() << "\"]";
      dot << ";" << endl;
    }
  }
  dot << "}" << endl; // end subgraph
  
  
  // Now the hits that resulted from the tracks.
  if(m_do_hits) {
    SimHitHeader::hc_map::const_iterator hcit = header->hits()->hitCollection().begin();
    for( ; hcit != header->hits()->hitCollection().end(); hcit ++){
      const SimHitCollection* hc = hcit->second;
      for(unsigned int ihit = 0; ihit<hc->collection().size(); ihit++) {
        const SimHit* simhit = hc->collection()[ihit];
        // the simhit description.
      
        dot << "simhit" << (long)(simhit) << "[shape=triangle,rank=sink,label=\"";
        if(dynamic_cast<const SimPmtHit*>(simhit))       dot << "Pmt";
        else if(dynamic_cast<const SimRpcHit*>(simhit))  dot << "Rpc";
        else                                             dot << "Hit";      
        dot << "\\nID=" 
            << ((simhit->sensDetId()&(0xFF000000))>>24) << "."
            << ((simhit->sensDetId()&(0xFF0000))>>16) <<  "."
            << ((simhit->sensDetId()&(0xFF00))>>8) << "."
            << ((simhit->sensDetId()&(0xFF)))
            <<"\"];" << endl;
      
        // connect track to hit.
        dot << " track" << (long)(simhit->ancestor().track()) << " -> simhit" << (long)(simhit);
        if(! simhit->ancestor().isDirect())
           dot << " [ headport = n, style = dotted, label=\"" << simhit->ancestor().indirection() << "\"]";
        dot << ";" << endl;
      }
    }
  }

  dot << "}" << endl;
  dot.close();
  
  return StatusCode::SUCCESS;
}

---