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, total_depth = 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_depth += l.second.second;
0139     ++total_cnt;
0140   }
0141   printout(INFO,"DDCMSDetElementCreator","++ Summary:     %-24s  %d.","Total DetElements:",total_cnt);
0142   printout(INFO,"DDCMSDetElementCreator","+++++++++++++++ Readout structure generation  ++++++++++++++++++++++++");
0143   str << endl;
0144   for( const auto& f : fields )   {
0145     string ro_name = f.first.name() + string("Hits");
0146     int num_bits = 8;
0147     id_str.str("");
0148     id_str << "system:" << num_bits;
0149     for( const auto& q : f.second )   {
0150       int bits = 0;
0151       if      ( q.second < 1<<0  ) bits = 1;
0152       else if ( q.second < 1<<1  ) bits = 1;
0153       else if ( q.second < 1<<2  ) bits = 2;
0154       else if ( q.second < 1<<3  ) bits = 3;
0155       else if ( q.second < 1<<4  ) bits = 4;
0156       else if ( q.second < 1<<5  ) bits = 5;
0157       else if ( q.second < 1<<6  ) bits = 6;
0158       else if ( q.second < 1<<7  ) bits = 7;
0159       else if ( q.second < 1<<8  ) bits = 8;
0160       else if ( q.second < 1<<9  ) bits = 9;
0161       else if ( q.second < 1<<10 ) bits = 10;
0162       else if ( q.second < 1<<11 ) bits = 11;
0163       else if ( q.second < 1<<12 ) bits = 12;
0164       else if ( q.second < 1<<13 ) bits = 13;
0165       else if ( q.second < 1<<14 ) bits = 14;
0166       else if ( q.second < 1<<15 ) bits = 15;
0167       bits += 1;
0168       id_str << ",Lv" << q.first << ":" << bits;
0169       num_bits += bits;
0170     }
0171     string idspec = id_str.str();
0172     str << "<readout name=\"" << ro_name << "\">" << endl
0173         << "\t<id>"
0174         << idspec
0175         << "</id>  <!-- Number of bits: " << num_bits << " -->" << endl
0176         << "</readout>" << endl;
0177     /// Create ID Descriptors and readout configurations
0178     IDDescriptor dsc(ro_name,idspec);
0179     description.addIDSpecification(dsc);
0180     Readout ro(ro_name);
0181     ro.setIDDescriptor(dsc);
0182     description.addReadout(ro);
0183     SensitiveDetector sd = description.sensitiveDetector(f.first.name());
0184     sd.setHitsCollection(ro.name());
0185     sd.setReadout(ro);
0186     printout(INFO,"DDCMSDetElementCreator",
0187              "++ Setting up readout for subdetector:%-24s id:%04X",
0188              f.first.name(), f.first.id());
0189   }
0190   printout(INFO,"DDCMSDetElementCreator","+++++++++++++++ ID Descriptor generation  ++++++++++++++++++++++++++++");
0191   printout(INFO,"",str.str().c_str());
0192   char volid[32];
0193   for(auto& p : all_placements )  {
0194     PlacedVolume place = p.first;
0195     Volume vol = place.volume();
0196     ::snprintf(volid,sizeof(volid),"Lv%d", p.second.first);
0197     printout(DEBUG,"DDCMSDetElementCreator",
0198              "++ Set volid (%-24s): %-6s = %3d  -> %s  (%p)",
0199              vol.isSensitive() ? vol.sensitiveDetector().name() : "Not Sensitive",
0200              volid, p.second.second, place.name(), place.ptr());
0201     place.addPhysVolID(volid, p.second.second);
0202   }
0203   printout(ALWAYS,"DDCMSDetElementCreator",
0204            "++ Instrumented %ld subdetectors with %d DetElements %d sensitives out of %d volumes and %ld sensitive placements.",
0205            fields.size(),total.elements,total.sensitives,total.volumes,all_placements.size());
0206 }
0207
0208 /// Create a new detector element
0209 DetElement DDCMSDetElementCreator::createElement(const char* /* debug_tag */, PlacedVolume pv, int id) {
0210   string     name = cms::detElementName(pv);
0211   DetElement det(name, id);
0212   det.setPlacement(pv);
0213   /*
0214     printout(INFO,"DDCMSDetElementCreator","++ Created detector element [%s]: %s (%s)  %p",
0215     debug_tag, det.name(), name.c_str(), det.ptr());
0216   */
0217   return det;
0218 }
0219
0220 /// Create the top level detectors
0221 void DDCMSDetElementCreator::createTopLevelDetectors(PlacedVolume pv)   {
0222   auto& data = stack.back();
0223   if ( stack.size() == 2 )    {      // Main subssystem: tracker:Tracker
0224     data.element = tracker = addSubdetector(cms::detElementName(pv), pv, false);
0225     tracker->SetTitle("compound");
0226   }
0227   else if ( stack.size() == 3 )    { // Main subsystem detector: TIB, TEC, ....
0228     data.element = current_detector = addSubdetector(cms::detElementName(pv), pv, true);
0229   }
0230 }
0231
0232 /// Add new subdetector to the detector description
0233 DetElement DDCMSDetElementCreator::addSubdetector(const std::string& nam, PlacedVolume pv, bool volid)  {
0234   Detectors::iterator idet = subdetectors.find(nam);
0235   if ( idet == subdetectors.end() )   {
0236     DetElement det(nam, subdetectors.size()+1);
0237     det.setPlacement(pv);
0238     if ( volid )  {
0239       det.placement().addPhysVolID("system",det.id());
0240     }
0241     idet = subdetectors.insert(make_pair(nam,det)).first;
0242     description.add(det);
0243   }
0244   return idet->second;
0245 }
0246
0247 /// Callback to output PlacedVolume information of an single Placement
0248 int DDCMSDetElementCreator::operator()(PlacedVolume pv, int vol_level)   {
0249   double frac_si = pv.volume().material().fraction(silicon);
0250   if ( frac_si > 90e-2 )  {
0251     Data& data = stack.back();
0252     data.sensitive     = true;
0253     data.has_sensitive = true;
0254     ++data.vol_count;
0255     int   idx   = pv->GetMotherVolume()->GetIndex(pv.ptr())+1;
0256     auto& cnt   = leafCount[make_pair(current_detector,vol_level)];
0257     cnt.first   = std::max(cnt.first,idx);
0258     ++cnt.second;
0259     all_placements[pv] = make_pair(vol_level,idx);
0260     return 1;
0261   }
0262   return 0;
0263 }
0264
0265 /// Callback to output PlacedVolume information of an entire Placement
0266 int DDCMSDetElementCreator::process(PlacedVolume pv, int level, bool recursive)   {
0267   stack.push_back(Data(pv));
0268   if ( stack.size() <= 3 )   {
0269     createTopLevelDetectors(pv);
0270   }
0271   int ret = PlacedVolumeProcessor::process(pv,level,recursive);
0272
0273   /// Complete structures if the stack size is > 3!
0274   if ( stack.size() > 3 )   {
0275     // Note: short-cuts to entries in the stack MUST be local and
0276     // initialized AFTER the call to "process"! The vector may be resized!
0277     auto& data = stack.back();
0278     auto& parent = stack[stack.size()-2];
0279     auto& counts = counters[current_detector];
0280     if ( data.sensitive )   {
0281       /// If this volume is sensitve, we must attach a sensitive detector handle
0282       if ( !current_sensitive.isValid() )  {
0283         SensitiveDetector sd = description.sensitiveDetector(current_detector.name());
0284         if ( !sd.isValid() )  {
0285           sd = SensitiveDetector(current_detector.name(),"tracker");
0286           current_detector->flag |= DetElement::Object::HAVE_SENSITIVE_DETECTOR;
0287           description.add(sd);
0288         }
0289         current_sensitive = sd;
0290       }
0291       pv.volume().setSensitiveDetector(current_sensitive);
0292       ++counts.sensitives;
0293     }
0294     ++counts.volumes;
0295     bool added = false;
0296     if ( data.vol_count > 0 )   {
0297       parent.daughter_count  += data.vol_count;
0298       parent.daughter_count  += data.daughter_count;
0299       data.has_sensitive      = true;
0300     }
0301     else   {
0302       parent.daughter_count  += data.daughter_count;
0303       data.has_sensitive      = (data.daughter_count>0);
0304     }
0305
0306     if ( data.has_sensitive )  {
0307       // If we have sensitive elements at this level or below,
0308       // we must complete the DetElement hierarchy
0309       if ( !data.element.isValid() )  {
0310         data.element = createElement("Element", data.pv, current_detector.id());
0311         ++counts.elements;
0312       }
0313       if ( !parent.element.isValid() )  {
0314         parent.element = createElement("Parent ", parent.pv, current_detector.id());
0315         ++counts.elements;
0316       }
0317       printout(DEBUG,"DDCMSDetElementCreator",
0318                "++ Assign detector element: %s (%p, %ld children) to %s (%p) with %ld vols",
0319                data.element.name(), data.element.ptr(), data.element.children().size(),
0320                parent.element.name(), parent.element.ptr(),
0321                data.vol_count);
0322
0323       // Trickle up the tree only for sensitive pathes. Forget the passive rest
0324       // This should automatically omit non-sensitive pathes
0325       parent.has_sensitive = true;
0326       parent.element.add(data.element);
0327       added = true;
0328       // It is simpler to collect the volumes and later assign the volids
0329       // rather than checking if the volid already exists.
0330       int vol_level = level;
0331       int idx = data.pv->GetMotherVolume()->GetIndex(data.pv.ptr())+1;
0332       all_placements[data.pv] = make_pair(vol_level,idx); // 1...n
0333       // Update counters
0334       auto& cnt_det   = leafCount[make_pair(current_detector,vol_level)];
0335       cnt_det.first   = std::max(cnt_det.first,idx);
0336       cnt_det.second += 1;
0337     }
0338     if ( !added && data.element.isValid() )  {
0339       printout(WARNING,"MEMORY-LEAK","Level:%3d Orpahaned DetElement:%s Daugthers:%d Parent:%s",
0340                int(stack.size()), data.element.name(), data.vol_count, parent.pv.name());
0341     }
0342   }
0343   /// Now the cleanup kicks in....
0344   if ( stack.size() == 3 )  {
0345     current_sensitive = SensitiveDetector();
0346     current_detector = DetElement();
0347   }
0348   stack.pop_back();
0349   return ret;
0350 }
0351
0352 static void* create_object(Detector& description, int /* argc */, char** /* argv */)   {
0353   PlacedVolumeProcessor* proc = new DDCMSDetElementCreator(description);
0354   return (void*)proc;
0355 }
0356
0357 // first argument is the type from the xml file
0358 DECLARE_DD4HEP_CONSTRUCTOR(DDCMS_DetElementCreator,create_object)
0359