src/gdml/DetElementCreator.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_DETELEMENTCREATOR_H
0018 #define DD4HEP_DETELEMENTCREATOR_H
0019
0020 // Framework include files
0021 #include <DD4hep/VolumeProcessor.h>
0022 #include <DD4hep/Printout.h>
0023
0024 namespace dd4hep {
0025
0026   /// DD4hep DetElement creator for the CMS geometry.
0027   /*  Heuristically assign DetElement structures to the sensitive volume pathes.
0028    *
0029    *  \author  M.Frank
0030    *  \version 1.0
0031    *  \ingroup DD4HEP_CORE
0032    */
0033   class DetElementCreator : public PlacedVolumeProcessor  {
0034     struct Data {
0035       PlacedVolume pv {0};
0036       DetElement   element {};
0037       bool         sensitive = false;
0038       bool         has_sensitive = false;
0039       int          level = 0;
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     typedef std::map<PlacedVolume, std::pair<int,int> > AllPlacements;
0062     Detector&         description;
0063     Material          sensitive_material;
0064     Counters          counters;
0065     LeafCount         leafCount;
0066     VolumeStack       stack;
0067     Detectors         subdetectors;
0068     DetElement        current_detector;
0069     std::string       detector;
0070     std::string       sensitive_material_name;
0071     std::string       sensitive_type;
0072     std::string       detector_volume_match;
0073     std::string       detector_volume_veto;
0074     size_t            detector_volume_level = 0;
0075     int               max_volume_level = 9999;
0076     PrintLevel        printLevel = INFO;
0077     SensitiveDetector current_sensitive;
0078     AllPlacements     all_placements;
0079
0080     /// Add new subdetector to the detector description
0081     DetElement addSubdetector(const std::string& nam, PlacedVolume pv, bool volid);
0082     /// Create a new detector element
0083     DetElement createElement(const char* debug_tag, PlacedVolume pv, int id);
0084     /// Create the top level detectors
0085     void createTopLevelDetectors(PlacedVolume pv);
0086     /// Generate the name of the DetElement object from the placed volume
0087     std::string detElementName(PlacedVolume pv)  const;
0088   public:
0089     /// Initializing constructor
0090     DetElementCreator(Detector& desc,
0091                       const std::string& detector,
0092                       const std::string& sd_type,
0093                       const std::string& sd_match,
0094                       const std::string& sd_veto,
0095                       const std::string& sd_mat,
0096                       int sd_lvl,
0097                       PrintLevel p);
0098     /// Default destructor
0099     virtual ~DetElementCreator() noexcept(false);
0100     /// Callback to output PlacedVolume information of an single Placement
0101     virtual int operator()(PlacedVolume pv, int level);
0102     /// Callback to output PlacedVolume information of an entire Placement
0103     virtual int process(PlacedVolume pv, int level, bool recursive);
0104   };
0105 }
0106 #endif   /* DD4HEP_DETELEMENTCREATOR_H  */
0107
0108 // Framework include files
0109 #include <DD4hep/detail/DetectorInterna.h>
0110 #include <DD4hep/DetFactoryHelper.h>
0111 #include <DD4hep/DetectorHelper.h>
0112 #include <DD4hep/Printout.h>
0113
0114 // C/C++ include files
0115 #include <sstream>
0116
0117 using namespace std;
0118 using namespace dd4hep;
0119
0120 /// Initializing constructor
0121 DetElementCreator::DetElementCreator(Detector& desc,
0122                                      const std::string& det,
0123                                      const string& sd_match,
0124                                      const string& sd_veto,
0125                                      const string& sd_type,
0126                                      const string& sd_mat, int sd_lvl,
0127                                      PrintLevel p)
0128   : description(desc), detector(det), sensitive_material_name(sd_mat),
0129     sensitive_type(sd_type), detector_volume_match(sd_match),
0130     detector_volume_veto(sd_veto), max_volume_level(sd_lvl), printLevel(p)
0131 {
0132   DetectorHelper helper(description);
0133   sensitive_material = desc.material(sensitive_material_name);
0134   if ( !sensitive_material.isValid() )   {
0135     except("DetElementCreator",
0136            "++ Failed to extract MATERIAL from the element table.");
0137   }
0138   stack.reserve(32);
0139   detector_volume_level = 0;
0140 }
0141
0142 /// Default destructor
0143 DetElementCreator::~DetElementCreator() noexcept(false)  {
0144   Count total;
0145   stringstream str, id_str;
0146   const char* pref = detector_volume_match.c_str();
0147   printout(INFO,pref,"DetElementCreator: +++++++++++++++ Summary of sensitve elements  ++++++++++++++++++++++++");
0148   for ( const auto& c : counters )  {
0149     printout(INFO,pref,"DetElementCreator: ++ Summary: SD: %-24s %7d DetElements %7d sensitives out of %7d volumes",
0150              (c.first.name()+string(":")).c_str(), c.second.elements, c.second.sensitives, c.second.volumes);
0151     total.elements   += c.second.elements;
0152     total.sensitives += c.second.sensitives;
0153     total.volumes    += c.second.volumes;
0154   }
0155   printout(INFO,pref,"DetElementCreator: ++ Summary:     %-24s %7d DetElements %7d sensitives out of %7d volumes",
0156            "Grand Total:",total.elements,total.sensitives,total.volumes);
0157   printout(INFO,pref,"DetElementCreator: +++++++++++++++ Summary of geometry depth analysis  ++++++++++++++++++");
0158   int total_cnt = 0, total_depth = 0;
0159   map<DetElement, vector<pair<int,int> > > fields;
0160   for ( const auto& l : leafCount )  {
0161     DetElement de = l.first.first;
0162     printout(INFO,pref,"DetElementCreator: ++ Summary: SD: %-24s system:%04X Lvl:%3d Sensitives: %6d [Max: %6d].",
0163              (de.name()+string(":")).c_str(), de.id(),
0164              l.first.second, l.second.second, l.second.first);
0165     fields[de].emplace_back(l.first.second,l.second.first);
0166     total_depth += l.second.second;
0167     ++total_cnt;
0168   }
0169   if ( 0 == total_depth )  { }
0170   printout(INFO, pref, "DetElementCreator: ++ Summary:     %-24s  %d.","Total DetElements:", total_cnt);
0171   printout(INFO, pref, "DetElementCreator: +++++++++++++++ Readout structure generation  ++++++++++++++++++++++++");
0172   str << endl;
0173   for( const auto& f : fields )   {
0174     string ro_name = f.first.name() + string("Hits");
0175     int num_bits = 8;
0176     id_str.str("");
0177     id_str << "system:" << num_bits;
0178     for( const auto& q : f.second )   {
0179       int bits = 0;
0180       if      ( q.second < 1<<0  ) bits = 1;
0181       else if ( q.second < 1<<1  ) bits = 1;
0182       else if ( q.second < 1<<2  ) bits = 2;
0183       else if ( q.second < 1<<3  ) bits = 3;
0184       else if ( q.second < 1<<4  ) bits = 4;
0185       else if ( q.second < 1<<5  ) bits = 5;
0186       else if ( q.second < 1<<6  ) bits = 6;
0187       else if ( q.second < 1<<7  ) bits = 7;
0188       else if ( q.second < 1<<8  ) bits = 8;
0189       else if ( q.second < 1<<9  ) bits = 9;
0190       else if ( q.second < 1<<10 ) bits = 10;
0191       else if ( q.second < 1<<11 ) bits = 11;
0192       else if ( q.second < 1<<12 ) bits = 12;
0193       else if ( q.second < 1<<13 ) bits = 13;
0194       else if ( q.second < 1<<14 ) bits = 14;
0195       else if ( q.second < 1<<15 ) bits = 15;
0196       bits += 1;
0197       id_str << ",Lv" << q.first << ":" << bits;
0198       num_bits += bits;
0199     }
0200     string idspec = id_str.str();
0201     str << "<readout name=\"" << ro_name << "\">" << endl
0202         << "\t<id>"
0203         << idspec
0204         << "</id>  <!-- Number of bits: " << num_bits << " -->" << endl
0205         << "</readout>" << endl;
0206
0207     /// Create ID Descriptors and readout configurations
0208     try   {
0209       IDDescriptor dsc(ro_name,idspec);
0210       description.addIDSpecification(dsc);
0211       Readout ro(ro_name);
0212       ro.setIDDescriptor(dsc);
0213       description.addReadout(ro);
0214       SensitiveDetector sd = description.sensitiveDetector(f.first.name());
0215       sd.setHitsCollection(ro.name());
0216       sd.setReadout(ro);
0217       printout(INFO,pref,"DetElementCreator: ++ Setting up readout for subdetector:%-24s id:%04X",
0218                f.first.name(), f.first.id());
0219     }
0220     catch(std::exception& e)    {
0221       printout(ERROR,pref,"DetElementCreator: ++ FAILED to setup readout for subdetector:%-24s id:%04X [%s]",
0222                f.first.name(), f.first.id(), e.what());
0223     }
0224   }
0225   printout(INFO,pref,"DetElementCreator: "
0226            "+++++++++++++++ ID Descriptor generation  ++++++++++++++++++++++++++++");
0227   printout(INFO,"",str.str().c_str());
0228   char volid[32];
0229   for(auto& p : all_placements )  {
0230     try  {
0231       PlacedVolume place = p.first;
0232       Volume vol = place.volume();
0233       ::snprintf(volid,sizeof(volid),"Lv%d", p.second.first);
0234       printout(DEBUG,pref, "DetElementCreator: ++ Set volid (%-24s): %-6s = %3d  -> %s  (%p)",
0235            vol.isSensitive() ? vol.sensitiveDetector().name() : "Not Sensitive",
0236            volid, p.second.second, place.name(), place.ptr());
0237       place.addPhysVolID(volid, p.second.second);
0238     }
0239     catch(const exception& e)   {
0240       except(pref, "DetElementCreator: Exception on destruction: %s", e.what());
0241     }
0242     catch(...)   {
0243       except(pref, "DetElementCreator: UNKNOWN Exception on destruction.");
0244     }
0245   }
0246   printout(ALWAYS, pref, "DetElementCreator: ++ Instrumented %ld subdetectors with %d "
0247            "DetElements %d sensitives out of %d volumes and %ld sensitive placements.",
0248            fields.size(),total.elements,total.sensitives,total.volumes,all_placements.size());
0249 }
0250
0251 /// Generate the name of the DetElement object from the placed volume
0252 string DetElementCreator::detElementName(PlacedVolume pv)  const    {
0253   if ( pv.isValid() )  {
0254     string nam = pv.name();
0255     size_t idx = string::npos; // nam.rfind('_');
0256     string nnam = nam.substr(0, idx);
0257     return nnam;
0258   }
0259   except("DetElementCreator","++ Cannot deduce name from invalid PlacedVolume handle!");
0260   return string();
0261 }
0262
0263 /// Create a new detector element
0264 DetElement DetElementCreator::createElement(const char* /* debug_tag */, PlacedVolume pv, int id) {
0265   string     name = detElementName(pv);
0266   DetElement det(name, id);
0267   det.setPlacement(pv);
0268   /*
0269     printout(INFO,"DetElementCreator","++ Created detector element [%s]: %s (%s)  %p",
0270     debug_tag, det.name(), name.c_str(), det.ptr());
0271   */
0272   return det;
0273 }
0274
0275 /// Create the top level detectors
0276 void DetElementCreator::createTopLevelDetectors(PlacedVolume pv)   {
0277   auto& data = stack.back();
0278   data.element = current_detector = addSubdetector(detElementName(pv), pv, true);
0279 }
0280
0281 /// Add new subdetector to the detector description
0282 DetElement DetElementCreator::addSubdetector(const std::string& nam, PlacedVolume pv, bool volid)  {
0283   Detectors::iterator idet = subdetectors.find(nam);
0284   if ( idet == subdetectors.end() )   {
0285     DetElement det(nam, description.detectors().size()+1);
0286     det.setPlacement(pv);
0287     if ( volid )  {
0288       det.placement().addPhysVolID("system",det.id());
0289     }
0290     idet = subdetectors.emplace(nam,det).first;
0291     description.add(det);
0292     printout(printLevel,"DetElementCreator","++ Added sub-detector element: %s",det.path().c_str());
0293   }
0294   return idet->second;
0295 }
0296
0297 /// Callback to output PlacedVolume information of an single Placement
0298 int DetElementCreator::operator()(PlacedVolume pv, int vol_level)   {
0299   if ( detector_volume_level > 0 )   {
0300     Material mat = pv.volume().material();
0301     if ( mat == sensitive_material )  {
0302       Data& data = stack.back();
0303       data.sensitive     = true;
0304       data.has_sensitive = true;
0305       ++data.vol_count;
0306       int   idx   = pv->GetMotherVolume()->GetIndex(pv.ptr())+1;
0307       auto& cnt   = leafCount[make_pair(current_detector,vol_level)];
0308       cnt.first   = std::max(cnt.first,idx);
0309       ++cnt.second;
0310       all_placements[pv] = make_pair(vol_level,idx);
0311       return 1;
0312     }
0313   }
0314   return 0;
0315 }
0316
0317 /// Callback to output PlacedVolume information of an entire Placement
0318 int DetElementCreator::process(PlacedVolume pv, int lvl, bool recursive)   {
0319   int ret = 1;
0320   string pv_nam = pv.name();
0321   if ( detector_volume_level > 0 ||
0322        ( (!detector_volume_match.empty() &&
0323           pv_nam.find(detector_volume_match) != string::npos) &&
0324          (detector_volume_veto.empty() ||
0325           pv_nam.find(detector_volume_veto)  == string::npos)  )  )
0326   {
0327     stack.emplace_back(Data(pv));
0328     if ( 0 == detector_volume_level )   {
0329       detector_volume_level = stack.size();
0330       createTopLevelDetectors(pv);
0331     }
0332     ret = PlacedVolumeProcessor::process(pv,lvl,recursive);
0333     /// Complete structures if the stack size is > 3!
0334     if ( stack.size() > detector_volume_level )   {
0335       // Note: short-cuts to entries in the stack MUST be local and
0336       // initialized AFTER the call to "process"! The vector may be resized!
0337       auto& data   = stack.back();
0338       auto& parent = stack[stack.size()-2];
0339       auto& counts = counters[current_detector];
0340       if ( data.sensitive )   {
0341         /// If this volume is sensitve, we must attach a sensitive detector handle
0342         if ( !current_sensitive.isValid() )  {
0343           SensitiveDetector sd = description.sensitiveDetector(current_detector.name());
0344           if ( !sd.isValid() )  {
0345             sd = SensitiveDetector(current_detector.name(), sensitive_type);
0346             current_detector->flag |= DetElement::Object::HAVE_SENSITIVE_DETECTOR;
0347             description.add(sd);
0348
0349           }
0350           current_sensitive = sd;
0351         }
0352         pv.volume().setSensitiveDetector(current_sensitive);
0353         ++counts.sensitives;
0354       }
0355       ++counts.volumes;
0356       bool added = false;
0357       if ( data.vol_count > 0 )    {
0358         parent.daughter_count  += data.vol_count;
0359         parent.daughter_count  += data.daughter_count;
0360         data.has_sensitive      = true;
0361       }
0362       else   {
0363         parent.daughter_count  += data.daughter_count;
0364         data.has_sensitive      = (data.daughter_count>0);
0365       }
0366
0367       if ( data.has_sensitive )    {
0368         // If we have sensitive elements at this level or below,
0369         // we must complete the DetElement hierarchy
0370         if ( data.pv.volIDs().empty() )   {
0371           char text[32];
0372           ::snprintf(text, sizeof(text), "Lv%d", lvl);
0373           data.pv.addPhysVolID(text, data.pv->GetMotherVolume()->GetIndex(data.pv.ptr())+1);
0374         }
0375         else   {
0376           AllPlacements::const_iterator e = all_placements.find(data.pv);
0377           if ( e != all_placements.end() && (*e).second.first != lvl)  {
0378             printout(ERROR,"DetElementCreator","PLacement VOLID error: %d <> %d",lvl,(*e).second.first);
0379           }
0380         }
0381
0382         for( size_t i=1; i<stack.size(); ++i )   {
0383           auto& d = stack[i];
0384           auto& p = stack[i-1];
0385           if ( !d.element.isValid() )    {
0386             d.element = createElement("Element", d.pv, current_detector.id());
0387             (i==1 ? current_detector : p.element).add(d.element);
0388             ++counts.elements;
0389           }
0390           p.has_sensitive = true;
0391         }
0392         printout(printLevel,"DetElementCreator",
0393                  "++ Assign detector element: %s (%p, %ld children) to %s (%p) with %ld vols",
0394                  data.element.name(), data.element.ptr(), data.element.children().size(),
0395                  parent.element.name(), parent.element.ptr(), data.vol_count);
0396         added = true;
0397         // It is simpler to collect the volumes and later assign the volids
0398         // rather than checking if the volid already exists.
0399         int vol_level = lvl;
0400         int idx = data.pv->GetMotherVolume()->GetIndex(data.pv.ptr())+1;
0401         all_placements[data.pv] = make_pair(vol_level,idx); // 1...n
0402         // Update counters
0403         auto& cnt_det   = leafCount[make_pair(current_detector,vol_level)];
0404         cnt_det.first   = std::max(cnt_det.first,idx);
0405         cnt_det.second += 1;
0406         printout(printLevel,"DetElementCreator","++ [%ld] Added element: %s",
0407                  stack.size(), data.element.path().c_str());
0408       }
0409       if ( !added && data.element.isValid() )  {
0410         printout(WARNING,"MEMORY-LEAK","Level:%3d Orpahaned DetElement:%s Daugthers:%d Parent:%s",
0411                  int(stack.size()), data.element.name(), data.vol_count, parent.pv.name());
0412       }
0413     }
0414     /// Now the cleanup kicks in....
0415     if ( stack.size() == detector_volume_level )  {
0416       current_sensitive = SensitiveDetector();
0417       current_detector = DetElement();
0418       detector_volume_level = 0;
0419       ret = 0;
0420     }
0421     stack.pop_back();
0422   }
0423   else if ( lvl < max_volume_level )  {
0424     //printout(printLevel, "", "+++ Skip volume %s", pv_nam.c_str());
0425     ret = PlacedVolumeProcessor::process(pv,lvl,recursive);
0426   }
0427   return ret;
0428 }
0429
0430 static void* create_object(Detector& description, int argc, char** argv)   {
0431   PrintLevel prt  = DEBUG;
0432   size_t sd_level = 99999;
0433   string sd_mat, sd_match, sd_veto, sd_type, detector;
0434   for(int i = 0; i < argc && argv[i]; ++i)  {
0435     if ( 0 == ::strncmp("-material",argv[i],5) )
0436       sd_mat   = argv[++i];
0437     else if ( 0 == ::strncmp("-match",argv[i],5) )
0438       sd_match = argv[++i];
0439     else if ( 0 == ::strncmp("-detector",argv[i],5) )
0440       detector = argv[++i];
0441     else if ( 0 == ::strncmp("-veto",argv[i],5) )
0442       sd_veto  = argv[++i];
0443     else if ( 0 == ::strncmp("-type",argv[i],5) )
0444       sd_type  = argv[++i];
0445     else if ( 0 == ::strncmp("-level",argv[i],5) )
0446       sd_level = ::atol(argv[++i]);
0447     else if ( 0 == ::strncmp("-print",argv[i],5) )
0448       prt = decodePrintLevel(argv[++i]);
0449     else
0450       break;
0451   }
0452   if ( sd_mat.empty() || sd_match.empty() || sd_type.empty() )   {
0453     cout <<
0454       "Usage: -plugin <name> -arg [-arg]                                            \n"
0455       "     name:  factory name DD4hep_ROOTGDMLParse                           \n"
0456       "     -material <string>  Sensitive material name (identifier)           \n"
0457       "     -match    <string>  Matching string for subdetector identification \n"
0458       "\tArguments given: " << arguments(argc,argv) << endl << flush;
0459     ::exit(EINVAL);
0460   }
0461   PlacedVolumeProcessor* proc = new DetElementCreator(description, detector, sd_match, sd_veto, sd_type, sd_mat, sd_level, prt);
0462   return (void*)proc;
0463 }
0464
0465 // first argument is the type from the xml file
0466 DECLARE_DD4HEP_CONSTRUCTOR(DD4hep_DetElementCreator,create_object)
0467