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   printout(INFO,pref,"DetElementCreator: ++ Summary:     %-24s  %d.","Total DetElements:",total_cnt);
0170   printout(INFO,pref,"DetElementCreator: +++++++++++++++ Readout structure generation  ++++++++++++++++++++++++");
0171   str << endl;
0172   for( const auto& f : fields )   {
0173     string ro_name = f.first.name() + string("Hits");
0174     int num_bits = 8;
0175     id_str.str("");
0176     id_str << "system:" << num_bits;
0177     for( const auto& q : f.second )   {
0178       int bits = 0;
0179       if      ( q.second < 1<<0  ) bits = 1;
0180       else if ( q.second < 1<<1  ) bits = 1;
0181       else if ( q.second < 1<<2  ) bits = 2;
0182       else if ( q.second < 1<<3  ) bits = 3;
0183       else if ( q.second < 1<<4  ) bits = 4;
0184       else if ( q.second < 1<<5  ) bits = 5;
0185       else if ( q.second < 1<<6  ) bits = 6;
0186       else if ( q.second < 1<<7  ) bits = 7;
0187       else if ( q.second < 1<<8  ) bits = 8;
0188       else if ( q.second < 1<<9  ) bits = 9;
0189       else if ( q.second < 1<<10 ) bits = 10;
0190       else if ( q.second < 1<<11 ) bits = 11;
0191       else if ( q.second < 1<<12 ) bits = 12;
0192       else if ( q.second < 1<<13 ) bits = 13;
0193       else if ( q.second < 1<<14 ) bits = 14;
0194       else if ( q.second < 1<<15 ) bits = 15;
0195       bits += 1;
0196       id_str << ",Lv" << q.first << ":" << bits;
0197       num_bits += bits;
0198     }
0199     string idspec = id_str.str();
0200     str << "<readout name=\"" << ro_name << "\">" << endl
0201         << "\t<id>"
0202         << idspec
0203         << "</id>  <!-- Number of bits: " << num_bits << " -->" << endl
0204         << "</readout>" << endl;
0205
0206     /// Create ID Descriptors and readout configurations
0207     try   {
0208       IDDescriptor dsc(ro_name,idspec);
0209       description.addIDSpecification(dsc);
0210       Readout ro(ro_name);
0211       ro.setIDDescriptor(dsc);
0212       description.addReadout(ro);
0213       SensitiveDetector sd = description.sensitiveDetector(f.first.name());
0214       sd.setHitsCollection(ro.name());
0215       sd.setReadout(ro);
0216       printout(INFO,pref,"DetElementCreator: ++ Setting up readout for subdetector:%-24s id:%04X",
0217                f.first.name(), f.first.id());
0218     }
0219     catch(std::exception& e)    {
0220       printout(ERROR,pref,"DetElementCreator: ++ FAILED to setup readout for subdetector:%-24s id:%04X [%s]",
0221                f.first.name(), f.first.id(), e.what());
0222     }
0223   }
0224   printout(INFO,pref,"DetElementCreator: "
0225            "+++++++++++++++ ID Descriptor generation  ++++++++++++++++++++++++++++");
0226   printout(INFO,"",str.str().c_str());
0227   char volid[32];
0228   for(auto& p : all_placements )  {
0229     try  {
0230       PlacedVolume place = p.first;
0231       Volume vol = place.volume();
0232       ::snprintf(volid,sizeof(volid),"Lv%d", p.second.first);
0233       printout(DEBUG,pref, "DetElementCreator: ++ Set volid (%-24s): %-6s = %3d  -> %s  (%p)",
0234            vol.isSensitive() ? vol.sensitiveDetector().name() : "Not Sensitive",
0235            volid, p.second.second, place.name(), place.ptr());
0236       place.addPhysVolID(volid, p.second.second);
0237     }
0238     catch(const exception& e)   {
0239       except(pref, "DetElementCreator: Exception on destruction: %s", e.what());
0240     }
0241     catch(...)   {
0242       except(pref, "DetElementCreator: UNKNOWN Exception on destruction.");
0243     }
0244   }
0245   printout(ALWAYS, pref, "DetElementCreator: ++ Instrumented %ld subdetectors with %d "
0246            "DetElements %d sensitives out of %d volumes and %ld sensitive placements.",
0247            fields.size(),total.elements,total.sensitives,total.volumes,all_placements.size());
0248 }
0249
0250 /// Generate the name of the DetElement object from the placed volume
0251 string DetElementCreator::detElementName(PlacedVolume pv)  const    {
0252   if ( pv.isValid() )  {
0253     string nam = pv.name();
0254     size_t idx = string::npos; // nam.rfind('_');
0255     string nnam = nam.substr(0, idx);
0256     return nnam;
0257   }
0258   except("DetElementCreator","++ Cannot deduce name from invalid PlacedVolume handle!");
0259   return string();
0260 }
0261
0262 /// Create a new detector element
0263 DetElement DetElementCreator::createElement(const char* /* debug_tag */, PlacedVolume pv, int id) {
0264   string     name = detElementName(pv);
0265   DetElement det(name, id);
0266   det.setPlacement(pv);
0267   /*
0268     printout(INFO,"DetElementCreator","++ Created detector element [%s]: %s (%s)  %p",
0269     debug_tag, det.name(), name.c_str(), det.ptr());
0270   */
0271   return det;
0272 }
0273
0274 /// Create the top level detectors
0275 void DetElementCreator::createTopLevelDetectors(PlacedVolume pv)   {
0276   auto& data = stack.back();
0277   data.element = current_detector = addSubdetector(detElementName(pv), pv, true);
0278 }
0279
0280 /// Add new subdetector to the detector description
0281 DetElement DetElementCreator::addSubdetector(const std::string& nam, PlacedVolume pv, bool volid)  {
0282   Detectors::iterator idet = subdetectors.find(nam);
0283   if ( idet == subdetectors.end() )   {
0284     DetElement det(nam, description.detectors().size()+1);
0285     det.setPlacement(pv);
0286     if ( volid )  {
0287       det.placement().addPhysVolID("system",det.id());
0288     }
0289     idet = subdetectors.emplace(nam,det).first;
0290     description.add(det);
0291     printout(printLevel,"DetElementCreator","++ Added sub-detector element: %s",det.path().c_str());
0292   }
0293   return idet->second;
0294 }
0295
0296 /// Callback to output PlacedVolume information of an single Placement
0297 int DetElementCreator::operator()(PlacedVolume pv, int vol_level)   {
0298   if ( detector_volume_level > 0 )   {
0299     Material mat = pv.volume().material();
0300     if ( mat == sensitive_material )  {
0301       Data& data = stack.back();
0302       data.sensitive     = true;
0303       data.has_sensitive = true;
0304       ++data.vol_count;
0305       int   idx   = pv->GetMotherVolume()->GetIndex(pv.ptr())+1;
0306       auto& cnt   = leafCount[make_pair(current_detector,vol_level)];
0307       cnt.first   = std::max(cnt.first,idx);
0308       ++cnt.second;
0309       all_placements[pv] = make_pair(vol_level,idx);
0310       return 1;
0311     }
0312   }
0313   return 0;
0314 }
0315
0316 /// Callback to output PlacedVolume information of an entire Placement
0317 int DetElementCreator::process(PlacedVolume pv, int lvl, bool recursive)   {
0318   int ret = 1;
0319   string pv_nam = pv.name();
0320   if ( detector_volume_level > 0 ||
0321        ( (!detector_volume_match.empty() &&
0322           pv_nam.find(detector_volume_match) != string::npos) &&
0323          (detector_volume_veto.empty() ||
0324           pv_nam.find(detector_volume_veto)  == string::npos)  )  )
0325   {
0326     stack.emplace_back(Data(pv));
0327     if ( 0 == detector_volume_level )   {
0328       detector_volume_level = stack.size();
0329       createTopLevelDetectors(pv);
0330     }
0331     ret = PlacedVolumeProcessor::process(pv,lvl,recursive);
0332     /// Complete structures if the stack size is > 3!
0333     if ( stack.size() > detector_volume_level )   {
0334       // Note: short-cuts to entries in the stack MUST be local and
0335       // initialized AFTER the call to "process"! The vector may be resized!
0336       auto& data   = stack.back();
0337       auto& parent = stack[stack.size()-2];
0338       auto& counts = counters[current_detector];
0339       if ( data.sensitive )   {
0340         /// If this volume is sensitve, we must attach a sensitive detector handle
0341         if ( !current_sensitive.isValid() )  {
0342           SensitiveDetector sd = description.sensitiveDetector(current_detector.name());
0343           if ( !sd.isValid() )  {
0344             sd = SensitiveDetector(current_detector.name(), sensitive_type);
0345             current_detector->flag |= DetElement::Object::HAVE_SENSITIVE_DETECTOR;
0346             description.add(sd);
0347
0348           }
0349           current_sensitive = sd;
0350         }
0351         pv.volume().setSensitiveDetector(current_sensitive);
0352         ++counts.sensitives;
0353       }
0354       ++counts.volumes;
0355       bool added = false;
0356       if ( data.vol_count > 0 )    {
0357         parent.daughter_count  += data.vol_count;
0358         parent.daughter_count  += data.daughter_count;
0359         data.has_sensitive      = true;
0360       }
0361       else   {
0362         parent.daughter_count  += data.daughter_count;
0363         data.has_sensitive      = (data.daughter_count>0);
0364       }
0365
0366       if ( data.has_sensitive )    {
0367         // If we have sensitive elements at this level or below,
0368         // we must complete the DetElement hierarchy
0369         if ( data.pv.volIDs().empty() )   {
0370           char text[32];
0371           ::snprintf(text, sizeof(text), "Lv%d", lvl);
0372           data.pv.addPhysVolID(text, data.pv->GetMotherVolume()->GetIndex(data.pv.ptr())+1);
0373         }
0374         else   {
0375           AllPlacements::const_iterator e = all_placements.find(data.pv);
0376           if ( e != all_placements.end() && (*e).second.first != lvl)  {
0377             printout(ERROR,"DetElementCreator","PLacement VOLID error: %d <> %d",lvl,(*e).second.first);
0378           }
0379         }
0380
0381         for(size_t i=1; i<stack.size(); ++i)   {
0382           auto& d = stack[i];
0383           auto& p = stack[i-1];
0384           if ( !d.element.isValid() )    {
0385             d.element = createElement("Element", d.pv, current_detector.id());
0386             (i==1 ? current_detector : p.element).add(d.element);
0387             ++counts.elements;
0388           }
0389           p.has_sensitive = true;
0390         }
0391         printout(printLevel,"DetElementCreator",
0392                  "++ Assign detector element: %s (%p, %ld children) to %s (%p) with %ld vols",
0393                  data.element.name(), data.element.ptr(), data.element.children().size(),
0394                  parent.element.name(), parent.element.ptr(), data.vol_count);
0395         added = true;
0396         // It is simpler to collect the volumes and later assign the volids
0397         // rather than checking if the volid already exists.
0398         int vol_level = lvl;
0399         int idx = data.pv->GetMotherVolume()->GetIndex(data.pv.ptr())+1;
0400         all_placements[data.pv] = make_pair(vol_level,idx); // 1...n
0401         // Update counters
0402         auto& cnt_det   = leafCount[make_pair(current_detector,vol_level)];
0403         cnt_det.first   = std::max(cnt_det.first,idx);
0404         cnt_det.second += 1;
0405         printout(printLevel,"DetElementCreator","++ [%ld] Added element: %s",
0406                  stack.size(), data.element.path().c_str());
0407       }
0408       if ( !added && data.element.isValid() )  {
0409         printout(WARNING,"MEMORY-LEAK","Level:%3d Orpahaned DetElement:%s Daugthers:%d Parent:%s",
0410                  int(stack.size()), data.element.name(), data.vol_count, parent.pv.name());
0411       }
0412     }
0413     /// Now the cleanup kicks in....
0414     if ( stack.size() == detector_volume_level )  {
0415       current_sensitive = SensitiveDetector();
0416       current_detector = DetElement();
0417       detector_volume_level = 0;
0418       ret = 0;
0419     }
0420     stack.pop_back();
0421   }
0422   else if ( lvl < max_volume_level )  {
0423     //printout(printLevel, "", "+++ Skip volume %s", pv_nam.c_str());
0424     ret = PlacedVolumeProcessor::process(pv,lvl,recursive);
0425   }
0426   return ret;
0427 }
0428
0429 static void* create_object(Detector& description, int argc, char** argv)   {
0430   PrintLevel prt  = DEBUG;
0431   size_t sd_level = 99999;
0432   string sd_mat, sd_match, sd_veto, sd_type, detector;
0433   for(int i = 0; i < argc && argv[i]; ++i)  {
0434     if ( 0 == ::strncmp("-material",argv[i],5) )
0435       sd_mat   = argv[++i];
0436     else if ( 0 == ::strncmp("-match",argv[i],5) )
0437       sd_match = argv[++i];
0438     else if ( 0 == ::strncmp("-detector",argv[i],5) )
0439       detector = argv[++i];
0440     else if ( 0 == ::strncmp("-veto",argv[i],5) )
0441       sd_veto  = argv[++i];
0442     else if ( 0 == ::strncmp("-type",argv[i],5) )
0443       sd_type  = argv[++i];
0444     else if ( 0 == ::strncmp("-level",argv[i],5) )
0445       sd_level = ::atol(argv[++i]);
0446     else if ( 0 == ::strncmp("-print",argv[i],5) )
0447       prt = decodePrintLevel(argv[++i]);
0448     else
0449       break;
0450   }
0451   if ( sd_mat.empty() || sd_match.empty() || sd_type.empty() )   {
0452     cout <<
0453       "Usage: -plugin <name> -arg [-arg]                                            \n"
0454       "     name:  factory name DD4hep_ROOTGDMLParse                           \n"
0455       "     -material <string>  Sensitive material name (identifier)           \n"
0456       "     -match    <string>  Matching string for subdetector identification \n"
0457       "\tArguments given: " << arguments(argc,argv) << endl << flush;
0458     ::exit(EINVAL);
0459   }
0460   PlacedVolumeProcessor* proc = new DetElementCreator(description, detector, sd_match, sd_veto, sd_type, sd_mat, sd_level, prt);
0461   return (void*)proc;
0462 }
0463
0464 // first argument is the type from the xml file
0465 DECLARE_DD4HEP_CONSTRUCTOR(DD4hep_DetElementCreator,create_object)
0466