src/plugins/DDCMSDetElementCreator.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013 //
0014 // Specialized generic detector constructor
0015 //
0016 //==========================================================================
0017 #ifndef DD4HEP_DDCMSDETELEMENTCREATOR_H
0018 #define DD4HEP_DDCMSDETELEMENTCREATOR_H
0019
0020 // Framework include files
0021 #include "DD4hep/VolumeProcessor.h"
0022
0023 namespace dd4hep {
0024
0025   /// DD4hep DetElement creator for the CMS geometry.
0026   /*  Heuristically assign DetElement structures to the sensitive volume pathes.
0027    *
0028    *  \author  M.Frank
0029    *  \version 1.0
0030    *  \ingroup DD4HEP_CORE
0031    */
0032   class DDCMSDetElementCreator : public PlacedVolumeProcessor  {
0033     Detector&    description;
0034     Atom         silicon;
0035     struct Data {
0036       PlacedVolume pv {0};
0037       DetElement   element {};
0038       bool         sensitive = false;
0039       bool         has_sensitive = false;
0040       int          vol_count = 0;
0041       int          daughter_count = 0;
0042       int          sensitive_count = 0;
0043
0044       Data() = default;
0045       Data(PlacedVolume v) : pv(v) {}
0046       Data(const Data& d) = default;
0047       Data& operator=(const Data& d) = default;
0048     };
0049     struct Count {
0050       int elements = 0;
0051       int volumes = 0;
0052       int sensitives = 0;
0053       Count() = default;
0054       Count(const Count&) = default;
0055       Count& operator=(const Count&) = default;
0056     };
0057     typedef std::vector<Data> VolumeStack;
0058     typedef std::map<std::string,DetElement> Detectors;
0059     typedef std::map<DetElement,Count> Counters;
0060     typedef std::map<std::pair<DetElement,int>, std::pair<int,int> > LeafCount;
0061
0062     Counters          counters;
0063     LeafCount         leafCount;
0064     VolumeStack       stack;
0065     Detectors         subdetectors;
0066     DetElement        tracker, current_detector;
0067     SensitiveDetector current_sensitive;
0068     std::map<PlacedVolume, std::pair<int,int> > all_placements;
0069
0070     /// Add new subdetector to the detector description
0071     DetElement addSubdetector(const std::string& nam, PlacedVolume pv, bool volid);
0072     /// Create a new detector element
0073     DetElement createElement(const char* debug_tag, PlacedVolume pv, int id);
0074     /// Create the top level detectors
0075     void createTopLevelDetectors(PlacedVolume pv);
0076
0077   public:
0078     /// Initializing constructor
0079     DDCMSDetElementCreator(Detector& desc);
0080     /// Default destructor
0081     virtual ~DDCMSDetElementCreator();
0082     /// Callback to output PlacedVolume information of an single Placement
0083     virtual int operator()(PlacedVolume pv, int level);
0084     /// Callback to output PlacedVolume information of an entire Placement
0085     virtual int process(PlacedVolume pv, int level, bool recursive);
0086   };
0087 }
0088 #endif   /* DD4HEP_DDCMSDETELEMENTCREATOR_H  */
0089
0090 #include "DD4hep/detail/DetectorInterna.h"
0091 #include "DD4hep/DetFactoryHelper.h"
0092 #include "DD4hep/DetectorHelper.h"
0093 #include "DD4hep/Printout.h"
0094 #include "DDCMS/DDCMS.h"
0095
0096 #include <sstream>
0097
0098 using namespace std;
0099 using namespace dd4hep;
0100
0101 /// Initializing constructor
0102 DDCMSDetElementCreator::DDCMSDetElementCreator(Detector& desc)
0103   : description(desc)
0104 {
0105   DetectorHelper helper(description);
0106   silicon = helper.element("materials:Silicon");
0107   if ( !silicon.isValid() )   {
0108     except("DDCMSDetElementCreator",
0109            "++ Failed to extract SILICON from the element table.");
0110   }
0111   stack.reserve(32);
0112 }
0113
0114 /// Default destructor
0115 DDCMSDetElementCreator::~DDCMSDetElementCreator()   {
0116   Count total;
0117   stringstream str, id_str;
0118
0119   printout(INFO,"DDCMSDetElementCreator","+++++++++++++++ Summary of sensitve elements  ++++++++++++++++++++++++");
0120   for ( const auto& c : counters )  {
0121     printout(INFO,"DDCMSDetElementCreator","++ Summary: SD: %-24s %7d DetElements %7d sensitives out of %7d volumes",
0122              (c.first.name()+string(":")).c_str(), c.second.elements, c.second.sensitives, c.second.volumes);
0123     total.elements   += c.second.elements;
0124     total.sensitives += c.second.sensitives;
0125     total.volumes    += c.second.volumes;
0126   }
0127   printout(INFO,"DDCMSDetElementCreator",  "++ Summary:     %-24s %7d DetElements %7d sensitives out of %7d volumes",
0128            "Grand Total:",total.elements,total.sensitives,total.volumes);
0129   printout(INFO,"DDCMSDetElementCreator","+++++++++++++++ Summary of geometry depth analysis  ++++++++++++++++++");
0130   int total_cnt = 0;
0131   map<DetElement, vector<pair<int,int> > > fields;
0132   for ( const auto& l : leafCount )  {
0133     DetElement de = l.first.first;
0134     printout(INFO,"DDCMSDetElementCreator","++ Summary: SD: %-24s system:%04X Lvl:%3d Sensitives: %6d [Max: %6d].",
0135              (de.name()+string(":")).c_str(), de.id(),
0136              l.first.second, l.second.second, l.second.first);
0137     fields[de].push_back(make_pair(l.first.second,l.second.first));
0138     ++total_cnt;
0139   }
0140   printout(INFO,"DDCMSDetElementCreator","++ Summary:     %-24s  %d.","Total DetElements:",total_cnt);
0141   printout(INFO,"DDCMSDetElementCreator","+++++++++++++++ Readout structure generation  ++++++++++++++++++++++++");
0142   str << endl;
0143   for( const auto& f : fields )   {
0144     string ro_name = f.first.name() + string("Hits");
0145     int num_bits = 8;
0146     id_str.str("");
0147     id_str << "system:" << num_bits;
0148     for( const auto& q : f.second )   {
0149       int bits = 0;
0150       if      ( q.second < 1<<0  ) bits = 1;
0151       else if ( q.second < 1<<1  ) bits = 1;
0152       else if ( q.second < 1<<2  ) bits = 2;
0153       else if ( q.second < 1<<3  ) bits = 3;
0154       else if ( q.second < 1<<4  ) bits = 4;
0155       else if ( q.second < 1<<5  ) bits = 5;
0156       else if ( q.second < 1<<6  ) bits = 6;
0157       else if ( q.second < 1<<7  ) bits = 7;
0158       else if ( q.second < 1<<8  ) bits = 8;
0159       else if ( q.second < 1<<9  ) bits = 9;
0160       else if ( q.second < 1<<10 ) bits = 10;
0161       else if ( q.second < 1<<11 ) bits = 11;
0162       else if ( q.second < 1<<12 ) bits = 12;
0163       else if ( q.second < 1<<13 ) bits = 13;
0164       else if ( q.second < 1<<14 ) bits = 14;
0165       else if ( q.second < 1<<15 ) bits = 15;
0166       bits += 1;
0167       id_str << ",Lv" << q.first << ":" << bits;
0168       num_bits += bits;
0169     }
0170     string idspec = id_str.str();
0171     str << "<readout name=\"" << ro_name << "\">" << endl
0172         << "\t<id>"
0173         << idspec
0174         << "</id>  <!-- Number of bits: " << num_bits << " -->" << endl
0175         << "</readout>" << endl;
0176     /// Create ID Descriptors and readout configurations
0177     IDDescriptor dsc(ro_name,idspec);
0178     description.addIDSpecification(dsc);
0179     Readout ro(ro_name);
0180     ro.setIDDescriptor(dsc);
0181     description.addReadout(ro);
0182     SensitiveDetector sd = description.sensitiveDetector(f.first.name());
0183     sd.setHitsCollection(ro.name());
0184     sd.setReadout(ro);
0185     printout(INFO,"DDCMSDetElementCreator",
0186              "++ Setting up readout for subdetector:%-24s id:%04X",
0187              f.first.name(), f.first.id());
0188   }
0189   printout(INFO,"DDCMSDetElementCreator","+++++++++++++++ ID Descriptor generation  ++++++++++++++++++++++++++++");
0190   printout(INFO,"",str.str().c_str());
0191   char volid[32];
0192   for(auto& p : all_placements )  {
0193     PlacedVolume place = p.first;
0194     Volume vol = place.volume();
0195     ::snprintf(volid,sizeof(volid),"Lv%d", p.second.first);
0196     printout(DEBUG,"DDCMSDetElementCreator",
0197              "++ Set volid (%-24s): %-6s = %3d  -> %s  (%p)",
0198              vol.isSensitive() ? vol.sensitiveDetector().name() : "Not Sensitive",
0199              volid, p.second.second, place.name(), place.ptr());
0200     place.addPhysVolID(volid, p.second.second);
0201   }
0202   printout(ALWAYS,"DDCMSDetElementCreator",
0203            "++ Instrumented %ld subdetectors with %d DetElements %d sensitives out of %d volumes and %ld sensitive placements.",
0204            fields.size(),total.elements,total.sensitives,total.volumes,all_placements.size());
0205 }
0206
0207 /// Create a new detector element
0208 DetElement DDCMSDetElementCreator::createElement(const char* /* debug_tag */, PlacedVolume pv, int id) {
0209   string     name = cms::detElementName(pv);
0210   DetElement det(name, id);
0211   det.setPlacement(pv);
0212   /*
0213     printout(INFO,"DDCMSDetElementCreator","++ Created detector element [%s]: %s (%s)  %p",
0214     debug_tag, det.name(), name.c_str(), det.ptr());
0215   */
0216   return det;
0217 }
0218
0219 /// Create the top level detectors
0220 void DDCMSDetElementCreator::createTopLevelDetectors(PlacedVolume pv)   {
0221   auto& data = stack.back();
0222   if ( stack.size() == 2 )    {      // Main subssystem: tracker:Tracker
0223     data.element = tracker = addSubdetector(cms::detElementName(pv), pv, false);
0224     tracker->SetTitle("compound");
0225   }
0226   else if ( stack.size() == 3 )    { // Main subsystem detector: TIB, TEC, ....
0227     data.element = current_detector = addSubdetector(cms::detElementName(pv), pv, true);
0228   }
0229 }
0230
0231 /// Add new subdetector to the detector description
0232 DetElement DDCMSDetElementCreator::addSubdetector(const std::string& nam, PlacedVolume pv, bool volid)  {
0233   Detectors::iterator idet = subdetectors.find(nam);
0234   if ( idet == subdetectors.end() )   {
0235     DetElement det(nam, subdetectors.size()+1);
0236     det.setPlacement(pv);
0237     if ( volid )  {
0238       det.placement().addPhysVolID("system",det.id());
0239     }
0240     idet = subdetectors.insert(make_pair(nam,det)).first;
0241     description.add(det);
0242   }
0243   return idet->second;
0244 }
0245
0246 /// Callback to output PlacedVolume information of an single Placement
0247 int DDCMSDetElementCreator::operator()(PlacedVolume pv, int vol_level)   {
0248   double frac_si = pv.volume().material().fraction(silicon);
0249   if ( frac_si > 90e-2 )  {
0250     Data& data = stack.back();
0251     data.sensitive     = true;
0252     data.has_sensitive = true;
0253     ++data.vol_count;
0254     int   idx   = pv->GetMotherVolume()->GetIndex(pv.ptr())+1;
0255     auto& cnt   = leafCount[make_pair(current_detector,vol_level)];
0256     cnt.first   = std::max(cnt.first,idx);
0257     ++cnt.second;
0258     all_placements[pv] = make_pair(vol_level,idx);
0259     return 1;
0260   }
0261   return 0;
0262 }
0263
0264 /// Callback to output PlacedVolume information of an entire Placement
0265 int DDCMSDetElementCreator::process(PlacedVolume pv, int level, bool recursive)   {
0266   stack.push_back(Data(pv));
0267   if ( stack.size() <= 3 )   {
0268     createTopLevelDetectors(pv);
0269   }
0270   int ret = PlacedVolumeProcessor::process(pv,level,recursive);
0271
0272   /// Complete structures if the stack size is > 3!
0273   if ( stack.size() > 3 )   {
0274     // Note: short-cuts to entries in the stack MUST be local and
0275     // initialized AFTER the call to "process"! The vector may be resized!
0276     auto& data = stack.back();
0277     auto& parent = stack[stack.size()-2];
0278     auto& counts = counters[current_detector];
0279     if ( data.sensitive )   {
0280       /// If this volume is sensitve, we must attach a sensitive detector handle
0281       if ( !current_sensitive.isValid() )  {
0282         SensitiveDetector sd = description.sensitiveDetector(current_detector.name());
0283         if ( !sd.isValid() )  {
0284           sd = SensitiveDetector(current_detector.name(),"tracker");
0285           current_detector->flag |= DetElement::Object::HAVE_SENSITIVE_DETECTOR;
0286           description.add(sd);
0287         }
0288         current_sensitive = sd;
0289       }
0290       pv.volume().setSensitiveDetector(current_sensitive);
0291       ++counts.sensitives;
0292     }
0293     ++counts.volumes;
0294     bool added = false;
0295     if ( data.vol_count > 0 )   {
0296       parent.daughter_count  += data.vol_count;
0297       parent.daughter_count  += data.daughter_count;
0298       data.has_sensitive      = true;
0299     }
0300     else   {
0301       parent.daughter_count  += data.daughter_count;
0302       data.has_sensitive      = (data.daughter_count>0);
0303     }
0304
0305     if ( data.has_sensitive )  {
0306       // If we have sensitive elements at this level or below,
0307       // we must complete the DetElement hierarchy
0308       if ( !data.element.isValid() )  {
0309         data.element = createElement("Element", data.pv, current_detector.id());
0310         ++counts.elements;
0311       }
0312       if ( !parent.element.isValid() )  {
0313         parent.element = createElement("Parent ", parent.pv, current_detector.id());
0314         ++counts.elements;
0315       }
0316       printout(DEBUG,"DDCMSDetElementCreator",
0317                "++ Assign detector element: %s (%p, %ld children) to %s (%p) with %ld vols",
0318                data.element.name(), data.element.ptr(), data.element.children().size(),
0319                parent.element.name(), parent.element.ptr(),
0320                data.vol_count);
0321
0322       // Trickle up the tree only for sensitive pathes. Forget the passive rest
0323       // This should automatically omit non-sensitive pathes
0324       parent.has_sensitive = true;
0325       parent.element.add(data.element);
0326       added = true;
0327       // It is simpler to collect the volumes and later assign the volids
0328       // rather than checking if the volid already exists.
0329       int vol_level = level;
0330       int idx = data.pv->GetMotherVolume()->GetIndex(data.pv.ptr())+1;
0331       all_placements[data.pv] = make_pair(vol_level,idx); // 1...n
0332       // Update counters
0333       auto& cnt_det   = leafCount[make_pair(current_detector,vol_level)];
0334       cnt_det.first   = std::max(cnt_det.first,idx);
0335       cnt_det.second += 1;
0336     }
0337     if ( !added && data.element.isValid() )  {
0338       printout(WARNING,"MEMORY-LEAK","Level:%3d Orpahaned DetElement:%s Daugthers:%d Parent:%s",
0339                int(stack.size()), data.element.name(), data.vol_count, parent.pv.name());
0340     }
0341   }
0342   /// Now the cleanup kicks in....
0343   if ( stack.size() == 3 )  {
0344     current_sensitive = SensitiveDetector();
0345     current_detector = DetElement();
0346   }
0347   stack.pop_back();
0348   return ret;
0349 }
0350
0351 static void* create_object(Detector& description, int /* argc */, char** /* argv */)   {
0352   PlacedVolumeProcessor* proc = new DDCMSDetElementCreator(description);
0353   return (void*)proc;
0354 }
0355
0356 // first argument is the type from the xml file
0357 DECLARE_DD4HEP_CONSTRUCTOR(DDCMS_DetElementCreator,create_object)
0358