Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:05:36

0001 #include <cstdlib>
0002 #include <chrono>
0003 #include <fstream>
0004 #include <map>
0005 #include <memory>
0006 #include <string>
0007 #include <utility>
0008 #include <vector>
0009 
0010 #include "DetectorDescription/Core/interface/DDRotationMatrix.h"
0011 #include "DetectorDescription/Core/interface/DDTranslation.h"
0012 #include "DetectorDescription/Core/interface/DDCompactView.h"
0013 #include "DetectorDescription/Core/interface/DDEnums.h"
0014 #include "DetectorDescription/Core/interface/DDExpandedNode.h"
0015 #include "DetectorDescription/Core/interface/DDExpandedView.h"
0016 #include "DetectorDescription/Core/interface/DDFilter.h"
0017 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0018 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0019 #include "DetectorDescription/Core/interface/DDMaterial.h"
0020 #include "DetectorDescription/Core/interface/DDName.h"
0021 #include "DetectorDescription/Core/interface/DDPartSelection.h"
0022 #include "DetectorDescription/Core/interface/DDPosData.h"
0023 #include "DetectorDescription/Core/interface/DDScope.h"
0024 #include "DetectorDescription/Core/interface/DDSolid.h"
0025 #include "DetectorDescription/Core/interface/DDTransform.h"
0026 #include "DetectorDescription/Core/interface/DDValue.h"
0027 #include "DetectorDescription/Core/interface/DDsvalues.h"
0028 #include "DataFormats/Math/interface/GeantUnits.h"
0029 #include "DataFormats/Math/interface/Graph.h"
0030 #include "DetectorDescription/Core/interface/ClhepEvaluator.h"
0031 #include "FWCore/Utilities/interface/Exception.h"
0032 
0033 using namespace geant_units::operators;
0034 
0035 namespace {
0036   class GroupFilter : public DDFilter {
0037   public:
0038     GroupFilter(std::vector<DDSpecificsFilter*>& filters) : filters_(filters) {}
0039 
0040     bool accept(const DDExpandedView& cv) const final {
0041       bool returnValue = true;
0042       for (const auto& f : filters_) {
0043         returnValue = returnValue and f->accept(cv);
0044         if (not returnValue) {
0045           break;
0046         }
0047       }
0048       return returnValue;
0049     };
0050 
0051   private:
0052     std::vector<DDSpecificsFilter*> filters_;
0053   };
0054 }  // namespace
0055 
0056 DDTranslation calc(const DDGeoHistory& aHist) {
0057   const DDGeoHistory& h = aHist;
0058   unsigned int sz = h.size();
0059   std::vector<DDRotationMatrix> vr;
0060   std::vector<DDTranslation> vt;
0061   DDRotationMatrix r;
0062   vr.emplace_back(r);
0063 
0064   if (h.size() > 1) {
0065     vt.emplace_back(h[1].posdata()->translation());
0066     unsigned int i = 1;
0067     for (; i <= sz - 2; ++i) {
0068       vr.emplace_back(vr.back() * h[i].posdata()->ddrot().rotation());
0069       vt.emplace_back(h[i + 1].posdata()->translation());
0070     }
0071   }
0072 
0073   DDTranslation t;
0074   for (unsigned int i = 0; i < vt.size(); ++i) {
0075     t += vr[i] * vt[i];
0076   }
0077   return t;
0078 }
0079 
0080 void debugHistory(const DDGeoHistory& h) {
0081   static constexpr char const c = 'a';
0082   DDGeoHistory::const_iterator it = h.begin();
0083   std::string fname("hdebug_");
0084   std::ofstream file((fname + c).c_str());
0085   std::vector<DDRotationMatrix> rmv;
0086   std::vector<DDTranslation> tv;
0087   for (; it != h.end(); ++it) {
0088   }
0089 }
0090 
0091 void goPersistent(const DDCompactView& cv, const std::string& file) {
0092   std::ofstream f(file.c_str());
0093   const auto& g = cv.graph();
0094   unsigned int node = 0;
0095   auto it = g.begin();
0096   for (; it != g.end(); ++it) {
0097     auto eit = it->begin();
0098     for (; eit != it->end(); ++eit) {
0099       unsigned int eindex = eit->first;
0100       int copyno = g.edgeData(eit->second)->copyno();
0101       double x, y, z;
0102 
0103       x = g.edgeData(eit->second)->trans().x();
0104       y = g.edgeData(eit->second)->trans().y();
0105       z = g.edgeData(eit->second)->trans().z();
0106       f << node << " " << eindex << " " << copyno << " " << x << " " << y << " " << z << " "
0107         << g.edgeData(eit->second)->ddrot().ddname().ns() << " " << g.edgeData(eit->second)->ddrot().ddname().name()
0108         << std::endl;
0109     }
0110     ++node;
0111   }
0112   f.close();
0113 }
0114 
0115 bool NEXT(DDFilteredView& fv, int& count) {
0116   bool result = false;
0117 
0118   if (fv.firstChild()) {
0119     result = true;
0120   } else if (fv.nextSibling()) {
0121     result = true;
0122   } else {
0123     while (fv.parent()) {
0124       if (fv.nextSibling()) {
0125         result = true;
0126         break;
0127       }
0128     }
0129   }
0130 
0131   if (result)
0132     ++count;
0133   return result;
0134 }
0135 
0136 void dumpHistory(const DDGeoHistory& h, bool short_dump = false) {
0137   DDGeoHistory::size_type i = 0;
0138   for (; i < h.size(); ++i) {
0139     std::cout << h[i].logicalPart().name() << "[" << h[i].copyno() << "]-";
0140     if (!short_dump) {
0141       DDAxisAngle ra(h[i].absRotation());
0142       std::cout << h[i].absTranslation() << ra.Axis() << convertRadToDeg(ra.Angle());
0143     }
0144   }
0145 }
0146 
0147 // function object to compare to ExpandedNodes
0148 struct NodeComp {
0149   bool operator()(const DDExpandedNode& n1, const DDExpandedNode& n2) {
0150     const DDTranslation& t1 = n1.absTranslation();
0151     const DDTranslation& t2 = n2.absTranslation();
0152 
0153     bool result = false;
0154 
0155     // 'alphabetical ordering' according to absolute position
0156 
0157     if (t1.z() < t2.z()) {
0158       result = true;
0159     } else if ((t1.z() == t2.z()) && (t1.y() < t2.y())) {
0160       result = true;
0161     } else if ((t1.z() == t2.z()) && (t1.y() == t2.y()) && (t1.x() < t2.x())) {
0162       result = true;
0163     }
0164 
0165     return result;
0166   }
0167 };
0168 
0169 #include <algorithm>  // sort
0170 
0171 using DDNodes = std::vector<DDExpandedNode>;
0172 
0173 void dump_nodes(DDNodes& nodes, int max = 100) {
0174   DDNodes::iterator it = nodes.begin();
0175   DDNodes::iterator ed = nodes.end();
0176 
0177   sort(it, ed, NodeComp());
0178   int nodeCount = 1;
0179   for (; it != nodes.end(); ++it) {
0180     std::cout << nodeCount << " " << it->logicalPart() << " trans=" << it->absTranslation() << " [mm]" << std::endl;
0181     std::cout << "  specifics: " << it->logicalPart().specifics().size() << std::endl;
0182     if (nodeCount == max)
0183       break;
0184     ++nodeCount;
0185   }
0186   if (nodeCount == max)
0187     std::cout << " ... truncated ... " << std::endl;
0188 }
0189 
0190 void tutorial() {
0191   std::cout << std::endl
0192             << std::endl
0193             << "   >>>      Entering DDD user-code in tutorial.cc       <<< " << std::endl
0194             << "            -------------------------------------" << std::endl
0195             << std::endl;
0196 
0197   // Initialize a CompactView.
0198   // (During XML parsing internal DDD objects are created in memory.
0199   //  These objects are then put into a handy structure called DDCompactView
0200   //  to allow hierarchical access to them)
0201   // The default constructor creates a CompactView using the content of
0202   // DDRoot::root() singleton being the root of the hierarchy.
0203   // That's why we call DDRootDef::instance().set(..) before!
0204   // The CompactView is an acyclic directed multigraph with nodes of type
0205   // DDLogicalPart and edges of type DDPosData* .
0206   // ( refere to header documentation of these interface-classes located
0207   //   in DDD/DDCore/interface/..)
0208   // DDLogicalPart provides access to material and solid definitons
0209   // DDPosData provides access to relative spacial positionings of the
0210   // current node relative to its parent node.
0211 
0212   // Parser handles it correctly now!
0213   // DDRootDef::instance().set(DDName("CMS","cms.xml"));
0214   DDCompactView cpv;
0215 
0216   // As we want to demonstrate some DDD features in the ECal, we set the
0217   // default namespace to 'ecal-endcap.xml' (which is where the endcap
0218   // geometry is defined).
0219   // (Setting a default namespace can help to shorten instantiations of
0220   // references to objects defined in the file the namespace belongs to)
0221   //typedef DDCurrentNamespace CNS;
0222   //CNS::ns()="ecal-endcap.xml";
0223 
0224   // Navigating the CompactView (Usecase for building up G4-geometries):
0225   // -------------------------------------------------------------------
0226   // In the Protoype SW the CompactView can be navigated only be taking
0227   // its internal data structure (Graph<DDLogcialPart,DDPosData*> available
0228   // as type 'graph_type').
0229   // The corresponding interfaces for navigation of graphs are described in
0230   // DDD/DDCore/interface/graph.h
0231 
0232   // First one takes the graph representation of CompactView
0233   const auto& cpvGraph = cpv.graph();
0234 
0235   // Using the Graph.h interface, some basic information is available:
0236   std::cout << "CompactView, basic information: " << std::endl
0237             << "   LogicalParts = " << cpvGraph.size() << std::endl
0238             << "   PosParts( converted to DDPosData *) = [not yet available]"  //<< cpvGraph.edgeCount()
0239             << std::endl
0240             << std::endl;
0241 
0242   // Now we navigate a bit:
0243   // The root of the graph can be obtained from DDCompactView:
0244   std::cout << "Root volume of the Detector Description is named [" << cpv.root() << "]" << std::endl << std::endl;
0245 
0246   // The same, but creating a reference to it:
0247   const DDLogicalPart& root = cpv.root();
0248   const DDLogicalPart& world = root;  //(DDName("CMS","cms"));
0249   std::cout << "The world volume is described by following solid:" << std::endl;
0250   std::cout << world.solid() << std::endl << std::endl;
0251   const DDMaterial& worldMaterial = root.material();
0252   std::cout << "The world volume is filled with following material:" << std::endl;
0253   std::cout << worldMaterial << std::endl << std::endl;
0254 
0255   // Of course, if the names are well known (from the XML) one can always
0256   // obtain direct access:
0257   // (note the implicit usage of the default namespace! Whenever reference
0258   //  variable are created without using DDName(..) but a simple std::string
0259   //  the default namespace DDCurrentNamespace::ns_ is automatically
0260   //  taken into account)
0261   DDLogicalPart endcapXtal("ECalEndcapCrystal");
0262   DDLogicalPart endcapBasket("ECalEndcapE7");
0263 
0264   // Let's switch from the CompactView to the ExpandedView
0265   // interfaced through DDExpandedView
0266   // The expanded view is a tree of DDExpandedNode s. The usual tree
0267   // navigation is supported. Please refere to the documentation
0268   // in DDD/DDCore/interface/DDExpandedView
0269   DDExpandedView exv(cpv);
0270 
0271   // Let's count all endcap-xtals in the detector
0272   // (the input geometry from ecal-endcap.xml is the XMLized version
0273   //  of OSCARs 1.3.0 endcap. When counting these xtals we get about
0274   //  1000 more than described in the references availabe from the
0275   //  Encap-Web ...)
0276   int sensVolCount(0);
0277   int overall(0);
0278   std::cout << "Start sensitive volumes counting ... " << std::endl;
0279   clock_t StartT, EndT;
0280   StartT = clock();
0281   while (exv.next()) {  // loop over the !whole! ExpandedView ...
0282     if (exv.logicalPart().category() == DDEnums::sensitive) {
0283       ++sensVolCount;
0284     }
0285     ++overall;
0286   }
0287   EndT = clock();
0288   std::cout << "Time: " << ((double)(EndT - StartT)) / double(CLOCKS_PER_SEC) << " sec" << std::endl;
0289   std::cout << "There were " << sensVolCount << " sensitive volumes counted" << std::endl
0290             << "out of " << overall << " expanded nodes! " << std::endl
0291             << std::endl;
0292 
0293   // Test the SCOPE
0294   DDCompactView ccv;
0295   DDExpandedView ev(ccv);
0296   ev.firstChild();
0297   ev.firstChild();
0298   ev.firstChild();
0299   ev.nextSibling();
0300   ev.firstChild();
0301   ev.firstChild();
0302   std::cout << "now: " << ev.logicalPart() << std::endl;
0303   DDGeoHistory h = ev.geoHistory();
0304   std::cout << "now-hist: " << h << std::endl;
0305 
0306   DDExpandedView sev(ccv);
0307   std::cout << "scope-set=" << sev.setScope(h, 1) << std::endl;
0308   std::cout << "scope-root=" << sev.logicalPart() << std::endl;
0309   std::cout << "scope-scope=" << sev.scope() << std::endl;
0310   int si = 0;
0311   while (sev.next()) {
0312     std::cout << "scope-next=" << sev.geoHistory() << std::endl;
0313     ++si;
0314   }
0315   std::cout << "counted " << si << " nodes in the scope " << std::endl << std::endl;
0316 
0317   // test the filter
0318 
0319   std::map<std::string, DDCompOp> cop;
0320   cop["=="] = DDCompOp::equals;
0321   cop["!="] = DDCompOp::not_equals;
0322   bool moreFilters = true;
0323   bool moreQueries = true;
0324   bool moreFilterCriteria = true;
0325   std::string flog, ls, p, cs, v, q;
0326   while (moreQueries) {
0327     std::vector<DDSpecificsFilter*> vecF;
0328     while (moreFilters) {
0329       DDSpecificsFilter* f = new DDSpecificsFilter();
0330       std::string flog;
0331       std::string asString;
0332       std::cout << "filter LogOp = ";
0333       std::cin >> flog;
0334       if (flog == "end")
0335         break;
0336       vecF.emplace_back(f);
0337       while (moreFilterCriteria) {
0338         std::cout << " logic   = ";
0339         std::cin >> ls;
0340         if (ls == "end")
0341           break;
0342         std::cout << " par-name= ";
0343         std::cin >> p;
0344         std::cout << " comp    = ";
0345         std::cin >> cs;
0346         std::cout << " par-val = ";
0347         std::cin >> v;
0348 
0349         double dv = 0.;
0350         try {
0351           dv = ClhepEvaluator().eval("", v);
0352         } catch (const cms::Exception& e) {
0353           dv = 0;
0354         }
0355         DDValue ddval(p, v, dv);
0356         vecF.back()->setCriteria(ddval, cop[cs]);
0357 
0358       }  //<- moreFilterCriteria
0359     }    //<- morFilters
0360 
0361     std::string ans;
0362     ans = "";
0363     DDCompactView aaaaa;
0364 
0365     std::cout << "go persistent (filename,n)?";
0366     std::cin >> ans;
0367     if (ans != "n") {
0368       goPersistent(aaaaa, ans);
0369     }
0370 
0371     std::cout << "default num-scheme? ";
0372     std::cin >> ans;
0373     if (ans == "y") {
0374       std::cout << "creating the default numbering scheme ..." << std::endl;
0375 
0376       DDExpandedView eeeee(aaaaa);
0377 
0378       std::cout << "do sibling-stack navigation?";
0379       std::cin >> ans;
0380       if (ans == "y") {
0381         bool go = true;
0382         while (go) {
0383           DDCompactView c;
0384           DDExpandedView e(c);
0385           DDExpandedView::nav_type n;
0386           std::cout << "  size ( 0 to stop )=";
0387           int s;
0388           std::cin >> s;
0389           go = (bool)s;
0390           int i = 0;
0391           for (; i < s; ++i) {
0392             int k;
0393             std::cin >> k;
0394             n.emplace_back(k);
0395           }
0396           std::cout << "input=" << n << std::endl;
0397           if (e.goTo(n)) {
0398             std::cout << "node=" << e.geoHistory() << std::endl;
0399           } else {
0400             std::cout << "no match!" << std::endl;
0401           }
0402         }
0403       }
0404 
0405       std::cout << "node calculation based on id?" << std::endl;
0406       std::cin >> ans;
0407       if (ans == "y") {
0408         bool go = true;
0409         while (go) {
0410           DDCompactView c;
0411           DDExpandedView e(c);
0412           DDExpandedView::nav_type n;
0413           std::cout << "  id ( 0 to stop)=";
0414           int s;
0415           std::cin >> s;
0416           go = (bool)s;
0417         }
0418       }
0419     }
0420     std::cout << "iterate the FilteredView (y/n)";
0421     std::cin >> ans;
0422     DDCompactView compactview;
0423 
0424     if (ans == "y") {
0425       GroupFilter gf(vecF);
0426       DDFilteredView fv(compactview, gf);
0427 
0428       int count = 0;
0429       std::cout << "The filtered-view contained " << count << " nodes." << std::endl;
0430       fv.reset();
0431       std::cout << "Now entering interactive navigation: f = (f)irstChild," << std::endl
0432                 << "                                     n = (n)extSibling," << std::endl
0433                 << "                                     p = (p)arent," << std::endl
0434                 << "                                     s = (s)tatus," << std::endl
0435                 << "                                     h = (h)istory debugging," << std::endl
0436                 << "                                     e = (e)nd" << std::endl;
0437       std::string nav = "";
0438       DDCompactView wcpv;
0439       while (nav != "e") {
0440         std::cout << "  nav = ";
0441         std::cin >> nav;
0442         char c = nav[0];
0443         typedef std::vector<std::pair<const DDPartSelection*, const DDsvalues_type*> > spectype;
0444         spectype v = fv.logicalPart().attachedSpecifics();
0445         std::vector<const DDsvalues_type*> vlp = fv.specifics();
0446         std::vector<const DDsvalues_type*> only = fv.logicalPart().specifics();
0447         DDsvalues_type merged = fv.mergedSpecifics();
0448         DDLogicalPart curlp = fv.logicalPart();
0449         bool result = false;
0450         switch (c) {
0451           case 'f':
0452             result = fv.firstChild();
0453             break;
0454           case 'n':
0455             result = fv.nextSibling();
0456             break;
0457           case 'p':
0458             result = fv.parent();
0459             break;
0460           case 'h':
0461             debugHistory(fv.geoHistory());
0462             break;
0463           case 's':
0464             fv.print();
0465             std::cout << std::endl << "specifics sets = " << v.size() << ":" << std::endl;
0466             for (const auto& o : v) {
0467               std::cout << *(o.first) << " = " << *(o.second) << std::endl;  // << std::endl;
0468             }
0469             std::cout << std::endl;
0470             std::cout << "merged-specifics:" << std::endl;
0471             std::cout << merged << std::endl;
0472 
0473             std::cout << "specifics only at logicalPart:" << std::endl;
0474             for (const auto& o : only) {
0475               std::cout << *o << std::endl;
0476             }
0477             std::cout << std::endl;
0478             std::cout << "translation: " << fv.translation() << std::endl;
0479             std::cout << "       calc: " << calc(fv.geoHistory()) << std::endl;
0480             {
0481               DDAxisAngle aa(fv.rotation());
0482               std::cout << "rotation: axis=" << aa.Axis() << " angle=" << convertRadToDeg(aa.Angle()) << std::endl
0483                         << std::endl;
0484             }
0485             std::cout << "sibling-stack=" << fv.navPos() << std::endl << std::endl;
0486             std::cout << "material=" << fv.logicalPart().material().ddname() << std::endl;
0487             std::cout << "solid=" << fv.logicalPart().solid().ddname()
0488                       << " volume[m3]=" << convertMm3ToM3(fv.logicalPart().solid().volume()) << std::endl;
0489             break;
0490           case 'e':
0491             break;
0492           default:
0493             std::cout << "  > not understood, try again < " << std::endl;
0494         }
0495         std::cout << "  node = " << fv.geoHistory().back() << " isNew=" << result << std::endl;
0496       }
0497     }
0498 
0499     std::cout << " again (y/n)?";
0500     std::cin >> ans;
0501     if (ans != "y")
0502       moreQueries = false;
0503 
0504     int fv_count = 0;
0505 
0506     clock_t Start, End;
0507     Start = clock();
0508     //while (NEXT(fv,fv_count)) ;
0509     int cc = 0;
0510     auto g = math::GraphWalker<DDLogicalPart, DDPosData*>(DDCompactView().graph(), DDCompactView().root());
0511     while (g.next())
0512       ++cc;
0513     End = clock();
0514     std::cout << "Time : " << ((double)(End - Start)) / double(CLOCKS_PER_SEC) << " sec" << std::endl;
0515     //std::cout << fv.history().back() << std::endl;
0516     std::cout << "Nodes: " << cc << std::endl;
0517     std::cout << "Using navigation the filtered-view has " << fv_count << " nodes." << std::endl;
0518   }
0519 
0520   /*
0521     std::cout << "Deleting transient store!" << std::endl;
0522     DDLogicalPart::clear();
0523     DDRotation::clear();
0524     DDAlgo::clear();
0525     DDMaterial::clear();
0526     DDSolid::clear();
0527     DDSpecifics::clear();
0528     std::cout << "Transient store deleted!" << std::endl;
0529   */
0530   return;
0531   // Or we look for any specific data attached to any nodes
0532   // while iterating the expanded-view:
0533   DDCompactView cp;
0534   DDExpandedView ex(cp);
0535   bool goon(true);  // just for truncation of lengthy output
0536   int trunc = 0;    // just for truncation of lengthy output
0537   std::cout << "retrieving all nodes with specific data attached ..." << std::endl;
0538   while (ex.next() && goon) {
0539     // ask each expanded-not for its specifics
0540     // std::vector<..>.size() will be 0 if there are no specifics
0541     std::vector<const DDsvalues_type*> spec = ex.specifics();
0542     if (!spec.empty()) {
0543       std::cout << spec.size() << " different specific-data sets found for " << std::endl;
0544       dumpHistory(ex.geoHistory(), true);
0545       std::cout << std::endl;
0546       if (trunc > 3) {
0547         std::cout << " ... truncated .. " << std::endl;
0548         goon = false;
0549       }
0550       ++trunc;
0551     }
0552   }
0553 
0554   // At last we create a graphical represention of the compact view by
0555   // iterating over the PosPart-Registry
0556 
0557 #ifdef DD_OLD_SCHEME
0558 
0559   std::ofstream file("graph.dot");
0560   std::cout << "Writing a dot-file for the graph-structure of the compact-view .." << std::endl
0561             << "File: graph.dot" << std::endl;
0562   file << "digraph G {" << std::endl;
0563   DDPosPartReg::instance_t::const_iterator PPIT = DDPosPartReg::instance().begin();
0564   for (; PPIT != DDPosPartReg::instance().end(); ++PPIT) {
0565     file << PPIT->second.myself().name() << " -> " << PPIT->second.mother().name() << " ; " << std::endl;
0566   }
0567   file << "}" << std::endl;
0568   file.close();
0569 
0570 #endif  // DD_OLD_SCHEME
0571 
0572   std::cout << std::endl
0573             << std::endl
0574             << "   >>>       End of DDD user-code in tutorial.cc       <<< " << std::endl
0575             << "            -------------------------------------" << std::endl
0576             << "                         ============            " << std::endl;
0577 
0578   std::cout << std::endl;
0579   int n_c = 500;
0580 
0581   exit(1);
0582 
0583   std::cout << "DUMPING the first " << n_c << " node-histories of the expanded-view" << std::endl << std::endl;
0584 
0585   DDCompactView scratch;
0586   DDExpandedView exx(scratch);
0587 
0588   while (exx.next() && n_c--) {
0589     dumpHistory(exx.geoHistory());
0590   }
0591 }