DDRec/src/MaterialScan.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 //  Simple program to print all the materials in a detector on
0015 //  a straight line between two given points
0016 //
0017 //  Author     : M.Frank, CERN
0018 //
0019 //==========================================================================
0020 // Framework include files
0021 #include <DDRec/MaterialScan.h>
0022 #include <DD4hep/DD4hepUnits.h>
0023 #include <DD4hep/Detector.h>
0024 #include <DD4hep/Printout.h>
0025 #include <DD4hep/Shapes.h>
0026
0027 /// C/C++ include files
0028 #include <cstdio>
0029
0030 using namespace dd4hep;
0031 using namespace dd4hep::rec;
0032
0033 /// Default constructor
0034 MaterialScan::MaterialScan()
0035   : m_detector(Detector::getInstance())
0036 {
0037   m_materialMgr.reset(new MaterialManager(m_detector.world().volume()));
0038 }
0039
0040 /// Default constructor
0041 MaterialScan::MaterialScan(Detector& description)
0042   : m_detector(description)
0043 {
0044   m_materialMgr.reset(new MaterialManager(m_detector.world().volume()));
0045 }
0046
0047 /// Default destructor
0048 MaterialScan::~MaterialScan()    {
0049 }
0050
0051
0052 /// Set a specific region to limit the scan (resets the subdetector selection)
0053 void MaterialScan::setRegion(const char* region)   {
0054   Region reg;
0055   if ( region )   {
0056     reg = m_detector.region(region);
0057   }
0058   setRegion(reg);
0059 }
0060
0061 /// Set a specific region to limit the scan (resets the subdetector selection)
0062 void MaterialScan::setRegion(Region region)   {
0063   struct PvCollector  {
0064     Region rg;
0065     std::set<const TGeoNode*>& cont;
0066     PvCollector(Region r, std::set<const TGeoNode*>& c) : rg(r), cont(c) {}
0067     void collect(TGeoNode* pv)    {
0068       cont.insert(pv);
0069       for (Int_t idau = 0, ndau = pv->GetNdaughters(); idau < ndau; ++idau)
0070         collect(pv->GetDaughter(idau));
0071     }
0072     void operator()(TGeoNode* pv)    {
0073       Volume v = pv->GetVolume();
0074       Region r = v.region();
0075       if ( r.isValid() )  {
0076         collect(pv);
0077         return;
0078       }
0079       for (Int_t idau = 0, ndau = pv->GetNdaughters(); idau < ndau; ++idau)
0080         (*this)(pv->GetDaughter(idau));
0081     }
0082   };
0083   m_placements.clear();
0084   if ( region.isValid() )   {
0085     PvCollector coll(region, m_placements);
0086     coll(m_detector.world().placement().ptr());
0087     printout(ALWAYS,"MaterialScan","+++ Set new scanning region to: %s [%ld placements]",
0088              region.name(), m_placements.size());
0089   }
0090   else   {
0091     printout(ALWAYS,"MaterialScan","+++ No region restrictions set. Back to full scanning mode.");
0092   }
0093 }
0094
0095 /// Set a specific detector volume to limit the scan
0096 void MaterialScan::setDetector(const char* detector)   {
0097   DetElement det;
0098   if ( detector )   {
0099     det = m_detector.detector(detector);
0100   }
0101   setDetector(det);
0102 }
0103
0104 /// Set a specific detector volume to limit the scan
0105 void MaterialScan::setDetector(DetElement detector)   {
0106   struct PvCollector  {
0107     std::set<const TGeoNode*>& cont;
0108     PvCollector(std::set<const TGeoNode*>& c) : cont(c) {}
0109     void operator()(TGeoNode* pv)    {
0110       cont.insert(pv);
0111       for (Int_t idau = 0, ndau = pv->GetNdaughters(); idau < ndau; ++idau) {
0112         TGeoNode* daughter = pv->GetDaughter(idau);
0113         (*this)(daughter);
0114       }
0115     }
0116   };
0117   if ( detector.isValid() )   {
0118     PlacedVolume pv = detector.placement();
0119     m_placements.clear();
0120     if ( pv.isValid() )  {
0121       PvCollector coll(m_placements);
0122       coll(pv.ptr());
0123     }
0124     printout(ALWAYS,"MaterialScan","+++ Set new scanning volume to: %s [%ld placements]",
0125              detector.path().c_str(), m_placements.size());
0126   }
0127   else   {
0128     printout(ALWAYS,"MaterialScan","+++ No subdetector restrictions set. Back to full scanning mode.");
0129     m_placements.clear();
0130   }
0131 }
0132
0133 /// Set a specific volume material to limit the scan
0134 void MaterialScan::setMaterial(const char* material)   {
0135   Material mat;
0136   if ( material )   {
0137     mat = m_detector.material(material);
0138   }
0139   setMaterial(mat);
0140 }
0141
0142 /// Set a specific volume material to limit the scan
0143 void MaterialScan::setMaterial(Material material)   {
0144   struct PvCollector  {
0145     Material material;
0146     std::set<const TGeoNode*>& cont;
0147     PvCollector(Material mat, std::set<const TGeoNode*>& c) : material(mat), cont(c) {}
0148     void operator()(TGeoNode* pv)    {
0149       Volume   vol = pv->GetVolume();
0150       Material mat = vol.material();
0151       if ( mat.ptr() == material.ptr() )  {
0152         cont.insert(pv);
0153       }
0154       for (Int_t idau = 0, ndau = pv->GetNdaughters(); idau < ndau; ++idau)
0155         (*this)(pv->GetDaughter(idau));
0156     }
0157   };
0158   m_placements.clear();
0159   if ( material.isValid() )   {
0160     PvCollector coll(material, m_placements);
0161     coll(m_detector.world().placement().ptr());
0162     printout(ALWAYS,"MaterialScan","+++ Set new scanning material to: %s [%ld placements]",
0163              material.name(), m_placements.size());
0164   }
0165   else   {
0166     printout(ALWAYS,"MaterialScan","+++ No material restrictions set. Back to full scanning mode.");
0167   }
0168 }
0169
0170 /// Scan along a line and store the matrials internally
0171 const MaterialVec& MaterialScan::scan(double x0, double y0, double z0, double x1, double y1, double z1, double epsilon)  const  {
0172   Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
0173   return m_materialMgr->materialsBetween(p0, p1, epsilon);
0174 }
0175
0176 /// Scan along a line and print the materials traversed
0177 void MaterialScan::print(const Vector3D& p0, const Vector3D& p1, double epsilon)  const    {
0178   const auto& placements = m_materialMgr->placementsBetween(p0, p1, epsilon);
0179   auto& matMgr = *m_materialMgr;
0180   Vector3D end, direction;
0181   double sum_x0 = 0;
0182   double sum_lambda = 0;
0183   double path_length = 0, total_length = 0;
0184   const char* fmt1 = " | %7s %-25s %3.0f %8.3f %8.4f %11.4f  %11.4f %10.3f %8.2f %11.6f %11.6f  (%7.2f,%7.2f,%7.2f)\n";
0185   const char* fmt2 = " | %7s %-25s %3.0f %8.3f %8.4f %11.6g  %11.6g %10.3f %8.2f %11.6f %11.6f  (%7.2f,%7.2f,%7.2f)\n";
0186   const char* line = " +--------------------------------------------------------------------------------------------------------------------------------------------------\n";
0187
0188   direction = (p1-p0).unit();
0189
0190   ::printf("%s + Material scan between: x_0 = (%7.2f,%7.2f,%7.2f) [cm] and x_1 = (%7.2f,%7.2f,%7.2f) [cm] : \n%s",
0191            line,p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],line);
0192   ::printf(" |       \\   %-16s          Atomic               Radiation   Interaction               Path   Integrated  Integrated    Material\n","Material");
0193   ::printf(" | Num.   \\  %-16s   Number/Z   Mass/A  Density    Length       Length    Thickness   Length      X0        Lambda      Endpoint/Startpoint\n","Name");
0194   ::printf(" | Layer   \\ %-16s            [g/mole]  [g/cm3]     [cm]        [cm]          [cm]      [cm]     [cm]        [cm]     (     cm,     cm,     cm)\n","");
0195   ::printf("%s",line);
0196   MaterialVec materials;
0197   std::set<Solid> tessellated_solids;
0198   for( unsigned i=0, n=placements.size(); i<n; ++i){
0199     TGeoNode*  pv  = placements[i].first.ptr();
0200     double length  = placements[i].second;
0201     total_length  += length;
0202     end = p0 + total_length * direction;
0203     if ( !m_placements.empty() && m_placements.find(pv) == m_placements.end() )  {
0204 #if 0
0205       ::printf("%p %s  %s  %s\n",
0206                placements[i].first.ptr(),
0207                placements[i].first->GetName(),
0208                placements[i].first->GetVolume()->GetName(),
0209                placements[i].first->GetMedium()->GetName());
0210 #endif
0211       continue;
0212     }
0213     TGeoMaterial* curr_mat = placements[i].first->GetMedium()->GetMaterial();
0214     TGeoMaterial* next_mat = (i+1)<n ? placements[i+1].first->GetMedium()->GetMaterial() : nullptr;
0215     double nx0     = length / curr_mat->GetRadLen();
0216     double nLambda = length / curr_mat->GetIntLen();
0217
0218     sum_x0        += nx0;
0219     sum_lambda    += nLambda;
0220     path_length   += length;
0221     materials.emplace_back(placements[i].first->GetMedium(),length);
0222     const char* fmt = curr_mat->GetRadLen() >= 1e5 ? fmt2 : fmt1;
0223     std::string mname = curr_mat->GetName();
0224     if ( next_mat && (i+1)<n )  {
0225       mname += " -> ";
0226       mname += next_mat->GetName();
0227     }
0228     Volume vol(placements[i].first->GetVolume());
0229     Solid  shape(vol.solid());
0230     if ( shape->IsA() == TGeoTessellated::Class() )  {
0231       tessellated_solids.insert(shape);
0232     }
0233     if ( 0 == i )  {
0234       ::printf(fmt, "(start)" , curr_mat->GetName(), curr_mat->GetZ(), curr_mat->GetA(),
0235            curr_mat->GetDensity(), curr_mat->GetRadLen()/dd4hep::cm, curr_mat->GetIntLen()/dd4hep::cm,
0236            0e0, 0e0, 0e0, 0e0,
0237            p0[0]/dd4hep::cm, p0[1]/dd4hep::cm, p0[2]/dd4hep::cm);
0238     }
0239     // No else here!
0240     if ( (n-1) == i )  {
0241       ::printf(fmt, "(end)", curr_mat->GetName(), curr_mat->GetZ(), curr_mat->GetA(),
0242            curr_mat->GetDensity(), curr_mat->GetRadLen()/dd4hep::cm, curr_mat->GetIntLen()/dd4hep::cm,
0243            0e0, 0e0, 0e0, 0e0,
0244            p0[0]/dd4hep::cm, p0[1]/dd4hep::cm, p0[2]/dd4hep::cm);
0245     }
0246     else   {
0247       ::printf(fmt, std::to_string(i+1).c_str(), mname.c_str(), next_mat->GetZ(), next_mat->GetA(),
0248            next_mat->GetDensity(), next_mat->GetRadLen()/dd4hep::cm, next_mat->GetIntLen()/dd4hep::cm,
0249            length/dd4hep::cm, path_length/dd4hep::cm, sum_x0, sum_lambda,
0250            end[0]/dd4hep::cm, end[1]/dd4hep::cm, end[2]/dd4hep::cm);
0251     }
0252   }
0253   ::printf("%s",line);
0254   const MaterialData& avg = matMgr.createAveragedMaterial(materials);
0255   const char* fmt = avg.radiationLength() >= 1e5 ? fmt2 : fmt1;
0256   ::printf(fmt,"","Average Material",avg.Z(),avg.A(),avg.density(),
0257            avg.radiationLength()/dd4hep::cm, avg.interactionLength()/dd4hep::cm,
0258            path_length/dd4hep::cm, path_length/dd4hep::cm,
0259            path_length/avg.radiationLength(),
0260            path_length/avg.interactionLength(),
0261            end[0]/dd4hep::cm, end[1]/dd4hep::cm, end[2]/dd4hep::cm);
0262   ::printf("%s",line);
0263
0264   if ( !tessellated_solids.empty() )  {
0265     ::printf(" |  WARNING: %ld tessellated shape were encountered during the volume traversal:\n",
0266          tessellated_solids.size());
0267     for(auto shape : tessellated_solids )
0268       ::printf(" |           \t\t %s\n", shape.name());
0269     ::printf(" |  WARNING: The results of this material scan are unreliable!\n");
0270     ::printf("%s",line);
0271   }
0272 }
0273
0274 /// Scan along a line and print the materials traversed
0275 void MaterialScan::print(double x0, double y0, double z0, double x1, double y1, double z1, double epsilon)  const  {
0276   Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
0277   this->print(p0, p1, epsilon);
0278 }