DDCMSDetElementCreator.cc

Go to the documentation of this file.

#include "DD4hep/VolumeProcessor.h"
#include "DD4hep/detail/DetectorInterna.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/DetectorHelper.h"
#include "DD4hep/Printout.h"
#include "DetectorDescription/DDCMS/interface/DDAlgoArguments.h"

#include <sstream>

using namespace std;
using namespace cms;
using namespace dd4hep;

namespace cms {
  
  // Heuristically assign DetElement structures
  // to the sensitive volume pathes
  //
  class DDCMSDetElementCreator : public dd4hep::PlacedVolumeProcessor
  {
  public:
    
    DDCMSDetElementCreator( dd4hep::Detector& );
    ~DDCMSDetElementCreator() override;
    
    int operator()( dd4hep::PlacedVolume volume, int level ) override;
    int process( dd4hep::PlacedVolume volume, int level, bool recursive ) override;

  private:

    dd4hep::DetElement addSubdetector( const std::string& nam, dd4hep::PlacedVolume volume, bool valid );
    dd4hep::DetElement createElement( const char* debugTag, dd4hep::PlacedVolume volume, int id );
    void createTopLevelDetectors( dd4hep::PlacedVolume volume );    

    struct Data
    {
      Data() = default;
      Data( dd4hep::PlacedVolume v ) : volume( v ) {}
      Data( const Data& d ) = default;
      Data& operator=( const Data& d ) = default;

      dd4hep::PlacedVolume volume {nullptr};
      dd4hep::DetElement   element {};
      bool         sensitive = false;
      bool         hasSensitive = false;
      int          volumeCount = 0;
      int          daughterCount = 0;
      int          sensitiveCount = 0;      
    };
    
    struct Count
    {
      Count() = default;
      Count( const Count& ) = default;
      Count& operator=( const Count& ) = default;

      int elements = 0;
      int volumes = 0;
      int sensitives = 0;
    };

    using Detectors = std::map<std::string, dd4hep::DetElement>;
    using Counters = std::map<dd4hep::DetElement, Count>;
    using LeafCount = std::map<std::pair<dd4hep::DetElement, int>, std::pair<int, int> >;
    using VolumeStack = std::vector<Data>;
    
    std::map<dd4hep::PlacedVolume, std::pair<int,int> > m_allPlacements;
        
    Counters          m_counters;
    LeafCount         m_leafCount;
    VolumeStack       m_stack;
    Detectors         m_subdetectors;
    dd4hep::DetElement m_tracker, m_currentDetector;
    dd4hep::SensitiveDetector m_currentSensitive;
    dd4hep::Detector&    m_description;
    dd4hep::Atom         m_silicon;
  };
 
  std::string detElementName( dd4hep::PlacedVolume volume );
}

std::string
cms::detElementName( dd4hep::PlacedVolume volume ) {
  if( volume.isValid()) {
    std::string name = volume.name();
    std::string nnam = name.substr( name.find('_') + 1 );
    return nnam;
  }
  except("MyDDCMS","++ Cannot deduce name from invalid PlacedVolume handle!");
  return std::string();
}

DDCMSDetElementCreator::DDCMSDetElementCreator( dd4hep::Detector& desc )
  : m_description( desc )
{
  dd4hep::DetectorHelper helper( m_description );
  m_silicon = helper.element( "SI" );
  if( !m_silicon.isValid()) {
    except( "DDCMSDetElementCreator",
        "++ Failed to extract SILICON from the element table." );
  }
  m_stack.reserve(32);
}

DDCMSDetElementCreator::~DDCMSDetElementCreator()
{
  Count total;
  stringstream str, id_str;

  printout( INFO, "DDCMSDetElementCreator", "+++++++++++++++ Summary of sensitve elements  ++++++++++++++++++++++++" );
  for( const auto& c : m_counters ) {
    printout( INFO, "DDCMSDetElementCreator", "++ Summary: SD: %-24s %7d DetElements %7d sensitives out of %7d volumes",
          ( c.first.name() + string(":")).c_str(), c.second.elements, c.second.sensitives, c.second.volumes );
    total.elements   += c.second.elements;
    total.sensitives += c.second.sensitives;
    total.volumes    += c.second.volumes;
  }
  printout( INFO, "DDCMSDetElementCreator", "++ Summary:     %-24s %7d DetElements %7d sensitives out of %7d volumes",
        "Grand Total:", total.elements, total.sensitives, total.volumes );
  printout( INFO, "DDCMSDetElementCreator", "+++++++++++++++ Summary of geometry depth analysis  ++++++++++++++++++" );
  int totalCount = 0, totalDepth = 0;
  map<dd4hep::DetElement, vector<pair<int,int> > > fields;
  for( const auto& l : m_leafCount ) {
    dd4hep::DetElement de = l.first.first;
    printout( INFO,"DDCMSDetElementCreator", "++ Summary: SD: %-24s system:%04X Lvl:%3d Sensitives: %6d [Max: %6d].",
          ( de.name() + string(":")).c_str(), de.id(),
          l.first.second, l.second.second, l.second.first );
    fields[de].push_back( make_pair( l.first.second, l.second.first ));
    totalDepth += l.second.second;
    ++totalCount;
  }
  printout( INFO, "DDCMSDetElementCreator", "++ Summary:     %-24s  %d.","Total DetElements:", totalCount );
  printout( INFO, "DDCMSDetElementCreator", "+++++++++++++++ Readout structure generation  ++++++++++++++++++++++++" );
  str << endl;
  for( const auto& f : fields ) {
    string roName = f.first.name() + string("Hits");
    int num_bits = 8;
    id_str.str("");
    id_str << "system:" << num_bits;
    for( const auto& q : f.second )   {
      int bits = 0;
      if      ( q.second < 1<<0  ) bits = 1;
      else if ( q.second < 1<<1  ) bits = 1;
      else if ( q.second < 1<<2  ) bits = 2;
      else if ( q.second < 1<<3  ) bits = 3;
      else if ( q.second < 1<<4  ) bits = 4;
      else if ( q.second < 1<<5  ) bits = 5;
      else if ( q.second < 1<<6  ) bits = 6;
      else if ( q.second < 1<<7  ) bits = 7;
      else if ( q.second < 1<<8  ) bits = 8;
      else if ( q.second < 1<<9  ) bits = 9;
      else if ( q.second < 1<<10 ) bits = 10;
      else if ( q.second < 1<<11 ) bits = 11;
      else if ( q.second < 1<<12 ) bits = 12;
      else if ( q.second < 1<<13 ) bits = 13;
      else if ( q.second < 1<<14 ) bits = 14;
      else if ( q.second < 1<<15 ) bits = 15;
      bits += 1;
      id_str << ",Lv" << q.first << ":" << bits;
      num_bits += bits;
    }
    string idspec = id_str.str();
    str << "<readout name=\"" << roName << "\">" << endl
        << "\t<id>"
        << idspec
        << "</id>  <!-- Number of bits: " << num_bits << " -->" << endl
        << "</readout>" << endl;

    IDDescriptor dsc( roName, idspec );
    m_description.addIDSpecification( dsc );
    Readout ro( roName );
    ro.setIDDescriptor( dsc );
    m_description.addReadout( ro );
    dd4hep::SensitiveDetector sd = m_description.sensitiveDetector( f.first.name());
    sd.setHitsCollection( ro.name());
    sd.setReadout( ro );
    printout( INFO, "DDCMSDetElementCreator",
          "++ Setting up readout for subdetector:%-24s id:%04X",
          f.first.name(), f.first.id());
  }
  printout( INFO, "DDCMSDetElementCreator","+++++++++++++++ ID Descriptor generation  ++++++++++++++++++++++++++++");
  printout( INFO, "", str.str().c_str());
  char volId[32];
  for(auto& p : m_allPlacements ) {
    dd4hep::PlacedVolume place = p.first;
    dd4hep::Volume volume = place.volume();
    ::snprintf( volId, sizeof(volId), "Lv%d", p.second.first );
    printout( DEBUG,"DDCMSDetElementCreator",
          "++ Set volid (%-24s): %-6s = %3d  -> %s  (%p)",
          volume.isSensitive() ? volume.sensitiveDetector().name() : "Not Sensitive",
          volId, p.second.second, place.name(), place.ptr());
    place.addPhysVolID( volId, p.second.second);
  }
  printout( ALWAYS, "DDCMSDetElementCreator",
        "++ Instrumented %ld subdetectors with %d DetElements %d sensitives out of %d volumes and %ld sensitive placements.",
        fields.size(), total.elements, total.sensitives, total.volumes, m_allPlacements.size());
}

dd4hep::DetElement
DDCMSDetElementCreator::createElement( const char*, PlacedVolume volume, int id )
{
  string name = detElementName( volume );
  dd4hep::DetElement det( name, id );
  det.setPlacement( volume );
  return det;
}
void
DDCMSDetElementCreator::createTopLevelDetectors( PlacedVolume volume )
{
  auto& data = m_stack.back();
  if( m_stack.size() == 2 ) {      // Main subssystem: tracker:Tracker
    data.element = m_tracker = addSubdetector( cms::detElementName(volume), volume, false );
    m_tracker->SetTitle( "compound" );
  }
  else if( m_stack.size() == 3 ) { // Main subsystem detector: TIB, TEC, ....
    data.element = m_currentDetector = addSubdetector( cms::detElementName(volume), volume, true );
  }
}

dd4hep::DetElement
DDCMSDetElementCreator::addSubdetector( const std::string& nam, dd4hep::PlacedVolume volume, bool valid )
{
  auto idet = m_subdetectors.find(nam);
  if( idet == m_subdetectors.end()) {
    dd4hep::DetElement det( nam, m_subdetectors.size() + 1 );
    det.setPlacement( volume );
    if( valid ) {
      det.placement().addPhysVolID( "system", det.id());
    }
    idet = m_subdetectors.insert( make_pair( nam, det)).first;
    m_description.add( det );
  }
  return idet->second;
}

int
DDCMSDetElementCreator::operator()( dd4hep::PlacedVolume volume, int volumeLevel )
{
  double fracSi = volume.volume().material().fraction( m_silicon );
  if( fracSi > 90e-2 ) {
    Data& data = m_stack.back();
    data.sensitive     = true;
    data.hasSensitive = true;
    ++data.volumeCount;
    int   idx   = volume->GetMotherVolume()->GetIndex( volume.ptr())+1;
    auto& cnt   = m_leafCount[make_pair( m_currentDetector, volumeLevel )];
    cnt.first   = std::max( cnt.first, idx );
    ++cnt.second;
    m_allPlacements[volume] = make_pair( volumeLevel, idx );
    return 1;
  }
  return 0;
}

    int
DDCMSDetElementCreator::process( dd4hep::PlacedVolume volume, int level,
                 bool recursive )
    {
  m_stack.push_back(Data(volume));
  if ( m_stack.size() <= 3 )   {
    createTopLevelDetectors(volume);
  }
  int ret = dd4hep::PlacedVolumeProcessor::process( volume, level, recursive );

  if( m_stack.size() > 3 ) {
    // Note: short-cuts to entries in the m_stack MUST be local and
    // initialized AFTER the call to "process"! The vector may be resized!
    auto& data = m_stack.back();
    auto& parent = m_stack[ m_stack.size() - 2 ];
    auto& counts = m_counters[ m_currentDetector ];
    if( data.sensitive ) {
      if ( !m_currentSensitive.isValid()) {
        dd4hep::SensitiveDetector sd = m_description.sensitiveDetector( m_currentDetector.name());
        if( !sd.isValid()) {
          sd = dd4hep::SensitiveDetector( m_currentDetector.name(), "tracker" );
          m_currentDetector->flag |= DetElement::Object::HAVE_SENSITIVE_DETECTOR;
          m_description.add( sd );
        }
        m_currentSensitive = sd;
      }
      volume.volume().setSensitiveDetector( m_currentSensitive );
      ++counts.sensitives;
    }
    ++counts.volumes;
    bool added = false;
    if( data.volumeCount > 0 ) {
      parent.daughterCount  += data.volumeCount;
      parent.daughterCount  += data.daughterCount;
      data.hasSensitive      = true;
    }
    else {
      parent.daughterCount  += data.daughterCount;
      data.hasSensitive      = ( data.daughterCount > 0 );
    }

    if( data.hasSensitive ) {
      // If we have sensitive elements at this level or below,
      // we must complete the DetElement hierarchy
      if( !data.element.isValid()) {
        data.element = createElement( "Element", data.volume, m_currentDetector.id());
        ++counts.elements;
      }
      if( !parent.element.isValid()) {
        parent.element = createElement( "Parent ", parent.volume, m_currentDetector.id());
        ++counts.elements;
      }
      printout( DEBUG, "DDCMSDetElementCreator",
        "++ Assign detector element: %s (%p, %ld children) to %s (%p) with %ld vols",
        data.element.name(), data.element.ptr(), data.element.children().size(),
        parent.element.name(), parent.element.ptr(),
        data.volumeCount );

      // Trickle up the tree only for sensitive pathes. Forget the passive rest
      // This should automatically omit non-sensitive pathes
      parent.hasSensitive = true;
      parent.element.add( data.element );
      added = true;
      // It is simpler to collect the volumes and later assign the volids
      // rather than checking if the volid already exists.
      int volumeLevel = level;
      int idx = data.volume->GetMotherVolume()->GetIndex( data.volume.ptr())+1;
      m_allPlacements[data.volume] = make_pair( volumeLevel, idx ); // 1...n
      // Update counters
      auto& cnt_det = m_leafCount[make_pair( m_currentDetector, volumeLevel )];
      cnt_det.first = std::max( cnt_det.first, idx );
      cnt_det.second += 1;
    }
    if ( !added && data.element.isValid() )  {
      printout(WARNING,"MEMORY-LEAK","Level:%3d Orpahaned DetElement:%s Daugthers:%d Parent:%s",
               int(m_stack.size()), data.element.name(), data.volumeCount, parent.volume.name());
    }
  }
  if( m_stack.size() == 3 ) {
    m_currentSensitive = SensitiveDetector();
    m_currentDetector = DetElement();
    ret = 0;
  }
  m_stack.pop_back();
  return ret;
}

static void*
createObject( dd4hep::Detector& description, int /* argc */, char** /* argv */) {
  dd4hep::PlacedVolumeProcessor* proc = new DDCMSDetElementCreator( description );
  return (void*)proc;
}

// first argument is the type from the xml file
DECLARE_DD4HEP_CONSTRUCTOR( DDCMS_DetElementCreator, createObject);