diff --git a/DDDigi/plugins/DigiTestAction.cpp b/DDDigi/plugins/DigiTestAction.cpp
index a9d404b95680c9b23159a8257f9e6344abcc36ec..6ba3c8e02a26ba6152ed5e07415ebd8a77d039b7 100644
--- a/DDDigi/plugins/DigiTestAction.cpp
+++ b/DDDigi/plugins/DigiTestAction.cpp
@@ -81,6 +81,11 @@ namespace dd4hep {
// C/C++ include files
+#ifdef __APPLE__
+static void noop(int) {}
+#define usleep(x) noop(x)
+#endif
+
using namespace std;
using namespace dd4hep::digi;
diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp
index 9de626bec8eaa5443758bb4a6984a69253b5b095..fd6a629c708b08346c33f5abbbf7ba78ae80263c 100644
--- a/DDG4/src/Geant4Converter.cpp
+++ b/DDG4/src/Geant4Converter.cpp
@@ -116,7 +116,7 @@ using namespace std;
#include "G4VPhysicalVolume.hh"
#include "G4ReflectionFactory.hh"
-static constexpr double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
+static const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
void Geant4AssemblyVolume::imprint(Geant4GeometryInfo& info,
const TGeoNode* parent,
diff --git a/DDRec/include/DDRec/MaterialScan.h b/DDRec/include/DDRec/MaterialScan.h
new file mode 100644
index 0000000000000000000000000000000000000000..03286c577870c409d2ec6934158f2633ec755f04
--- /dev/null
+++ b/DDRec/include/DDRec/MaterialScan.h
@@ -0,0 +1,72 @@
+//==========================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author : M.Frank
+//
+//==========================================================================
+#ifndef DD4HEP_DDREC_MATERIALSCAN_H
+#define DD4HEP_DDREC_MATERIALSCAN_H
+
+// Framework include files
+#include "DDRec/MaterialManager.h"
+
+#include <memory>
+
+/// Namespace for the AIDA detector description toolkit
+namespace dd4hep {
+
+ /// Forward declarations
+ class Detector;
+
+ /// Namespace for the reconstruction part of the AIDA detector description toolkit
+ namespace rec {
+
+ /// Class to perform material scans along a straight line
+ /**
+ *
+ * \author M.Frank
+ * \version 1.0
+ * \ingroup DD4HEP_REC
+ */
+ class MaterialScan {
+ private:
+
+ /// Reference to detector setup
+ Detector& m_detector;
+ /// Material manager
+ std::unique_ptr<MaterialManager> m_materialMgr; //!
+
+ /// Default constructor
+ MaterialScan();
+
+ public:
+
+ /// Standard constructor for the master instance
+ MaterialScan(Detector& description);
+
+ /// Default destructor
+ virtual ~MaterialScan();
+
+ /// Set a specific detector volume to limit the scan
+ void setDetector(DetElement detector);
+
+ /// Scan along a line and store the matrials internally
+ const MaterialVec& scan(double x0, double y0, double z0, double x1, double y1, double z1, double epsilon=1e-4) const;
+
+ /// Scan along a line and print the materials traversed
+ void print(const Vector3D& start, const Vector3D& end, double epsilon=1e-4) const;
+
+ /// Scan along a line and print the materials traversed
+ void print(double x0, double y0, double z0,
+ double x1, double y1, double z1,
+ double epsilon=1e-4) const;
+ };
+ } // End namespace rec
+} // End namespace dd4hep
+#endif // DD4HEP_DDREC_MATERIALSCAN_H
diff --git a/DDRec/src/MaterialManager.cpp b/DDRec/src/MaterialManager.cpp
index 7331880f293784250e0d0b19156da5564f5829cc..d786bd824cdeabf5d72346bbaf24caade1c25663 100644
--- a/DDRec/src/MaterialManager.cpp
+++ b/DDRec/src/MaterialManager.cpp
@@ -27,121 +27,121 @@ namespace dd4hep {
if( ( p0 != _p0 ) || ( p1 != _p1 ) ) {
- //---------------------------------------
- _mV.clear() ;
+ //---------------------------------------
+ _mV.clear() ;
- //
- // algorithm copied from TGeoGearDistanceProperties.cc (A.Munnich):
- //
+ //
+ // algorithm copied from TGeoGearDistanceProperties.cc (A.Munnich):
+ //
- double startpoint[3], endpoint[3], direction[3];
- double L=0;
- for(unsigned int i=0; i<3; i++) {
- startpoint[i] = p0[i];
- endpoint[i] = p1[i];
- direction[i] = endpoint[i] - startpoint[i];
- L+=direction[i]*direction[i];
- }
- double totDist = sqrt( L ) ;
+ double startpoint[3], endpoint[3], direction[3];
+ double L=0;
+ for(unsigned int i=0; i<3; i++) {
+ startpoint[i] = p0[i];
+ endpoint[i] = p1[i];
+ direction[i] = endpoint[i] - startpoint[i];
+ L+=direction[i]*direction[i];
+ }
+ double totDist = sqrt( L ) ;
- //normalize direction
- for(unsigned int i=0; i<3; i++)
- direction[i]=direction[i]/totDist;
+ //normalize direction
+ for(unsigned int i=0; i<3; i++)
+ direction[i]=direction[i]/totDist;
- _tgeoMgr->AddTrack(0, 12 ) ; // electron neutrino
+ _tgeoMgr->AddTrack(0, 12 ) ; // electron neutrino
- TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
+ TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
- //check if there is a node at startpoint
- if(!node1)
- throw std::runtime_error("No geometry node found at given location. Either there is no node placed here or position is outside of top volume.");
+ //check if there is a node at startpoint
+ if(!node1)
+ throw std::runtime_error("No geometry node found at given location. Either there is no node placed here or position is outside of top volume.");
- while ( !_tgeoMgr->IsOutside() ) {
+ while ( !_tgeoMgr->IsOutside() ) {
- // TGeoNode *node2;
- // TVirtualGeoTrack *track;
+ // TGeoNode *node2;
+ // TVirtualGeoTrack *track;
- // step to (and over) the next Boundary
- TGeoNode * node2 = _tgeoMgr->FindNextBoundaryAndStep( 500, 1) ;
+ // step to (and over) the next Boundary
+ TGeoNode * node2 = _tgeoMgr->FindNextBoundaryAndStep( 500, 1) ;
- if( !node2 || _tgeoMgr->IsOutside() )
- break;
+ if( !node2 || _tgeoMgr->IsOutside() )
+ break;
- const double *position = _tgeoMgr->GetCurrentPoint();
- const double *previouspos = _tgeoMgr->GetLastPoint();
+ const double *position = _tgeoMgr->GetCurrentPoint();
+ const double *previouspos = _tgeoMgr->GetLastPoint();
- double length = _tgeoMgr->GetStep();
+ double length = _tgeoMgr->GetStep();
- TVirtualGeoTrack *track = _tgeoMgr->GetLastTrack();
+ TVirtualGeoTrack *track = _tgeoMgr->GetLastTrack();
- //protection against infinitive loop in root which should not happen, but well it does...
- //work around until solution within root can be found when the step gets very small e.g. 1e-10
- //and the next boundary is never reached
+ //protection against infinitive loop in root which should not happen, but well it does...
+ //work around until solution within root can be found when the step gets very small e.g. 1e-10
+ //and the next boundary is never reached
#if 1 //fg: is this still needed ?
- if( length < MINSTEP ) {
+ if( length < MINSTEP ) {
- _tgeoMgr->SetCurrentPoint( position[0] + MINSTEP * direction[0],
- position[1] + MINSTEP * direction[1],
- position[2] + MINSTEP * direction[2] );
+ _tgeoMgr->SetCurrentPoint( position[0] + MINSTEP * direction[0],
+ position[1] + MINSTEP * direction[1],
+ position[2] + MINSTEP * direction[2] );
- length = _tgeoMgr->GetStep();
- node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
+ length = _tgeoMgr->GetStep();
+ node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
- position = _tgeoMgr->GetCurrentPoint();
- previouspos = _tgeoMgr->GetLastPoint();
- }
+ position = _tgeoMgr->GetCurrentPoint();
+ previouspos = _tgeoMgr->GetLastPoint();
+ }
#endif
- // printf( " -- step length : %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e - %s \n" , length ,
- // position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
+ // printf( " -- step length : %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e - %s \n" , length ,
+ // position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
- Vector3D posV( position ) ;
+ Vector3D posV( position ) ;
- double currDistance = ( posV - p0 ).r() ;
+ double currDistance = ( posV - p0 ).r() ;
- // //if the next boundary is further than end point
- // if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
- // || fabs(position[2])>fabs(endpoint[2]))
+ // //if the next boundary is further than end point
+ // if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
+ // || fabs(position[2])>fabs(endpoint[2]))
- //if we travelled too far:
- if( currDistance > totDist ) {
+ //if we travelled too far:
+ if( currDistance > totDist ) {
- length = sqrt( pow(endpoint[0]-previouspos[0],2) +
- pow(endpoint[1]-previouspos[1],2) +
- pow(endpoint[2]-previouspos[2],2) );
+ length = sqrt( pow(endpoint[0]-previouspos[0],2) +
+ pow(endpoint[1]-previouspos[1],2) +
+ pow(endpoint[2]-previouspos[2],2) );
- track->AddPoint( endpoint[0], endpoint[1], endpoint[2], 0. );
+ track->AddPoint( endpoint[0], endpoint[1], endpoint[2], 0. );
- if( length > epsilon )
- _mV.push_back( std::make_pair( Material( node1->GetMedium() ) , length ) ) ;
+ if( length > epsilon )
+ _mV.push_back( std::make_pair( Material( node1->GetMedium() ) , length ) ) ;
- break;
- }
+ break;
+ }
- track->AddPoint( position[0], position[1], position[2], 0.);
+ track->AddPoint( position[0], position[1], position[2], 0.);
- if( length > epsilon )
- _mV.push_back( std::make_pair( Material( node1->GetMedium() ), length ) ) ;
+ if( length > epsilon )
+ _mV.push_back( std::make_pair( Material( node1->GetMedium() ), length ) ) ;
- node1 = node2 ;
- }
+ node1 = node2 ;
+ }
- //fg: protect against empty list:
- if( _mV.empty() ){
- _mV.push_back( std::make_pair( Material( node1->GetMedium() ), totDist ) ) ;
- }
+ //fg: protect against empty list:
+ if( _mV.empty() ){
+ _mV.push_back( std::make_pair( Material( node1->GetMedium() ), totDist ) ) ;
+ }
- _tgeoMgr->ClearTracks();
+ _tgeoMgr->ClearTracks();
- _tgeoMgr->CleanGarbage();
+ _tgeoMgr->CleanGarbage();
- //---------------------------------------
+ //---------------------------------------
- _p0 = p0 ;
- _p1 = p1 ;
+ _p0 = p0 ;
+ _p1 = p1 ;
}
return _mV ; ;
@@ -152,19 +152,19 @@ namespace dd4hep {
if( pos != _pos ) {
- TGeoNode *node=_tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
+ TGeoNode *node=_tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
- if( ! node ) {
- std::stringstream err ;
- err << " MaterialManager::material: No geometry node found at location: " << pos ;
- throw std::runtime_error( err.str() );
- }
+ if( ! node ) {
+ std::stringstream err ;
+ err << " MaterialManager::material: No geometry node found at location: " << pos ;
+ throw std::runtime_error( err.str() );
+ }
- // std::cout << " @@@ MaterialManager::material @ " << pos << " found volume : " << node->GetName() << std::endl ;
+ // std::cout << " @@@ MaterialManager::material @ " << pos << " found volume : " << node->GetName() << std::endl ;
- _m = Material( node->GetMedium() ) ;
+ _m = Material( node->GetMedium() ) ;
- _pos = pos ;
+ _pos = pos ;
}
return _m ; ;
@@ -183,45 +183,45 @@ namespace dd4hep {
//double sum_rho_l_over_lambda = 0 ;
double sum_l_over_lambda = 0 ;
- for(unsigned i=0,n=materials.size(); i<n ; ++i){
+ for(unsigned i=0,n=materials.size(); i<n ; ++i){
- Material mat = materials[i].first ;
- double l = materials[i].second ;
+ Material mat = materials[i].first ;
+ double l = materials[i].second ;
- if( i != 0 ) sstr << "_" ;
- sstr << mat.name() << "_" << l ;
+ if( i != 0 ) sstr << "_" ;
+ sstr << mat.name() << "_" << l ;
- double rho = mat.density() ;
- double A = mat.A() ;
- double Z = mat.Z() ;
- double x = mat.radLength() ;
- double lambda = mat.intLength() ;
+ double rho = mat.density() ;
+ double A = mat.A() ;
+ double Z = mat.Z() ;
+ double x = mat.radLength() ;
+ double lambda = mat.intLength() ;
- sum_l += l ;
- sum_rho_l += rho * l ;
- sum_rho_l_over_A += rho * l / A ;
- sum_rho_l_Z_over_A += rho * l * Z / A ;
- sum_l_over_x += l / x ;
- sum_l_over_lambda += l / lambda ;
- // sum_rho_l_over_x += rho * l / x ;
- // sum_rho_l_over_lambda += rho * l / lambda ;
- }
-
- double rho = sum_rho_l / sum_l ;
-
- double A = sum_rho_l / sum_rho_l_over_A ;
- double Z = sum_rho_l_Z_over_A / sum_rho_l_over_A ;
-
- // radiation and interaction lengths already given in cm - average by length
+ sum_l += l ;
+ sum_rho_l += rho * l ;
+ sum_rho_l_over_A += rho * l / A ;
+ sum_rho_l_Z_over_A += rho * l * Z / A ;
+ sum_l_over_x += l / x ;
+ sum_l_over_lambda += l / lambda ;
+ // sum_rho_l_over_x += rho * l / x ;
+ // sum_rho_l_over_lambda += rho * l / lambda ;
+ }
+
+ double rho = sum_rho_l / sum_l ;
+
+ double A = sum_rho_l / sum_rho_l_over_A ;
+ double Z = sum_rho_l_Z_over_A / sum_rho_l_over_A ;
+
+ // radiation and interaction lengths already given in cm - average by length
- // double x = sum_rho_l / sum_rho_l_over_x ;
- double x = sum_l / sum_l_over_x ;
+ // double x = sum_rho_l / sum_rho_l_over_x ;
+ double x = sum_l / sum_l_over_x ;
- // double lambda = sum_rho_l / sum_rho_l_over_lambda ;
- double lambda = sum_l / sum_l_over_lambda ;
+ // double lambda = sum_rho_l / sum_rho_l_over_lambda ;
+ double lambda = sum_l / sum_l_over_lambda ;
- return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
+ return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
}
diff --git a/DDRec/src/MaterialScan.cpp b/DDRec/src/MaterialScan.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2453fff0e1d1f750264bd4638d6636259de75f6d
--- /dev/null
+++ b/DDRec/src/MaterialScan.cpp
@@ -0,0 +1,116 @@
+//==========================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author : M.Frank
+//
+//==========================================================================
+//
+// Simple program to print all the materials in a detector on
+// a straight line between two given points
+//
+// Author : M.Frank, CERN
+//
+//==========================================================================
+// Framework include files
+#include "DDRec/MaterialScan.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/Printout.h"
+
+#include <cstdio>
+
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+/// Default constructor
+MaterialScan::MaterialScan()
+ : m_detector(Detector::getInstance())
+{
+ m_materialMgr.reset(new MaterialManager(m_detector.world().volume()));
+}
+
+/// Default constructor
+MaterialScan::MaterialScan(Detector& description)
+ : m_detector(description)
+{
+ m_materialMgr.reset(new MaterialManager(m_detector.world().volume()));
+}
+
+/// Default destructor
+MaterialScan::~MaterialScan() {
+}
+
+
+/// Set a specific detector volume to limit the scan
+void MaterialScan::setDetector(DetElement detector) {
+ if ( detector.isValid() ) {
+ printout(ALWAYS,"MaterialScan","+++ Set new scanning volume to: %s",detector.path().c_str());
+ m_materialMgr.reset(new MaterialManager(detector.volume()));
+ }
+}
+
+/// Scan along a line and store the matrials internally
+const MaterialVec& MaterialScan::scan(double x0, double y0, double z0, double x1, double y1, double z1, double epsilon) const {
+ Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
+ return m_materialMgr->materialsBetween(p0, p1, epsilon);
+}
+
+/// Scan along a line and print the materials traversed
+void MaterialScan::print(const Vector3D& p0, const Vector3D& p1, double epsilon) const {
+ const auto& materials = m_materialMgr->materialsBetween(p0, p1, epsilon);
+ auto& matMgr = *m_materialMgr;
+ Vector3D end, direction;
+ double sum_x0 = 0;
+ double sum_lambda = 0;
+ double path_length = 0;
+ const char* fmt1 = " | %5d %-20s %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";
+ const char* fmt2 = " | %5d %-20s %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";
+ const char* line = " +--------------------------------------------------------------------------------------------------------------------------------------------------\n";
+
+ direction = (p1-p0).unit();
+
+ ::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",
+ line,p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],line);
+ ::printf(" | \\ %-11s Atomic Radiation Interaction Path Integrated Integrated Material\n","Material");
+ ::printf(" | Num. \\ %-11s Number/Z Mass/A Density Length Length Thickness Length X0 Lambda Endpoint \n","Name");
+ ::printf(" | Layer \\ %-11s [g/mole] [g/cm3] [cm] [cm] [cm] [cm] [cm] [cm] ( cm, cm, cm)\n","");
+ ::printf("%s",line);
+
+ for( unsigned i=0,n=materials.size();i<n;++i){
+ TGeoMaterial* mat = materials[i].first->GetMaterial();
+ double length = materials[i].second;
+ double nx0 = length / mat->GetRadLen();
+ sum_x0 += nx0;
+ double nLambda = length / mat->GetIntLen();
+ sum_lambda += nLambda;
+ path_length += length;
+ end = p0 + path_length * direction;
+ const char* fmt = mat->GetRadLen() >= 1e5 ? fmt2 : fmt1;
+ ::printf(fmt, i+1, mat->GetName(), mat->GetZ(), mat->GetA(),
+ mat->GetDensity(), mat->GetRadLen(), mat->GetIntLen(),
+ length, path_length, sum_x0, sum_lambda, end[0], end[1], end[2]);
+ //mat->Print();
+ }
+ printf("%s",line);
+ const MaterialData& avg = matMgr.createAveragedMaterial(materials);
+ const char* fmt = avg.radiationLength() >= 1e5 ? fmt2 : fmt1;
+ end = p0 + path_length * direction;
+ ::printf(fmt,0,"Average Material",avg.Z(),avg.A(),avg.density(),
+ avg.radiationLength(), avg.interactionLength(),
+ path_length, path_length,
+ path_length/avg.radiationLength(),
+ path_length/avg.interactionLength(),
+ end[0], end[1], end[2]);
+ printf("%s",line);
+}
+
+/// Scan along a line and print the materials traversed
+void MaterialScan::print(double x0, double y0, double z0, double x1, double y1, double z1, double epsilon) const {
+ Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
+ print(p0, p1, epsilon);
+}
diff --git a/DDRec/src/RecDictionary.h b/DDRec/src/RecDictionary.h
index d6866de486d3e533e9f9063827894bae153d73e9..81c4f8aae2c17eb72e0ab06d6a0b5f50ef01d2a6 100644
--- a/DDRec/src/RecDictionary.h
+++ b/DDRec/src/RecDictionary.h
@@ -19,6 +19,7 @@
#include "DDRec/DetectorData.h"
#include "DDRec/DetectorSurfaces.h"
#include "DDRec/MaterialManager.h"
+#include "DDRec/MaterialScan.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"
@@ -84,6 +85,7 @@ using namespace dd4hep::rec;
// DDRec/Material.h
#pragma link C++ class MaterialData+;
#pragma link C++ class MaterialManager+;
+#pragma link C++ class MaterialScan+;
#pragma link C++ class VolSurfaceBase+;
#pragma link C++ class VolSurface+;
#pragma link C++ class VolSurfaceList+;
diff --git a/UtilityApps/src/materialScan.cpp b/UtilityApps/src/materialScan.cpp
index b823fe3ecf18935791aeee667ab60ae90680c671..0d4bcc859f105e3981cbabc27903da55e8903178 100644
--- a/UtilityApps/src/materialScan.cpp
+++ b/UtilityApps/src/materialScan.cpp
@@ -18,11 +18,12 @@
//==========================================================================
#include "TError.h"
+#include "TInterpreter.h"
// Framework include files
#include "DD4hep/Detector.h"
#include "DD4hep/Printout.h"
-#include "DDRec/MaterialManager.h"
+#include "DDRec/MaterialScan.h"
#include "main.h"
using namespace dd4hep;
@@ -35,71 +36,38 @@ int main_wrapper(int argc, char** argv) {
if ( level > kInfo || abort ) ::printf("%s: %s\n", location, msg);
}
static void usage() {
- std::cout << " usage: materialScan compact.xml x0 y0 z0 x1 y1 z1 " << std::endl
- << " -> prints the materials on a straight line between the two given points ( unit is cm) "
+ std::cout << " usage: materialScan compact.xml x0 y0 z0 x1 y1 z1 [-interactive]" << std::endl
+ << " -> prints the materials on a straight line between the two given points ( unit is cm) " << std::endl
+ << " -interactive Load geometry once, then allow for shots from the ROOT prompt"
<< std::endl;
exit(1);
}
} _handler;
- if( argc != 8 ) Handler::usage();
-
- Vector3D p0, p1, end, direction;
+ if( argc < 8 ) Handler::usage();
+ bool interactive = argc == 9 && ::strncmp(argv[8],"-interactive",5) == 0;
+ double x0, y0, z0, x1, y1, z1;
std::string inFile = argv[1];
std::stringstream sstr;
sstr << argv[2] << " " << argv[3] << " " << argv[4] << " " << argv[5] << " " << argv[6] << " " << argv[7] << " " << "NONE";
- sstr >> p0.array()[0] >> p0.array()[1] >> p0.array()[2];
- sstr >> p1.array()[0] >> p1.array()[1] >> p1.array()[2];
+ sstr >> x0 >> y0 >> z0 >> x1 >> y1 >> z1;
if ( !sstr.good() ) Handler::usage();
setPrintLevel(WARNING);
Detector& description = Detector::getInstance();
description.fromXML(inFile);
- direction = (p1-p0).unit();
-
- MaterialManager matMgr( description.world().volume() ) ;
- //MaterialManager matMgr( description.world() ) ;
-
- const MaterialVec& materials = matMgr.materialsBetween(p0, p1);
- double sum_x0 = 0;
- double sum_lambda = 0;
- double path_length = 0;
- const char* fmt1 = " | %5d %-20s %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";
- const char* fmt2 = " | %5d %-20s %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";
- const char* line = " +--------------------------------------------------------------------------------------------------------------------------------------------------\n";
-
- ::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",
- line,p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],line);
- ::printf(" | \\ %-11s Atomic Radiation Interaction Path Integrated Integrated Material\n","Material");
- ::printf(" | Num. \\ %-11s Number/Z Mass/A Density Length Length Thickness Length X0 Lambda Endpoint \n","Name");
- ::printf(" | Layer \\ %-11s [g/mole] [g/cm3] [cm] [cm] [cm] [cm] [cm] [cm] ( cm, cm, cm)\n","");
- ::printf("%s",line);
-
- for( unsigned i=0,n=materials.size();i<n;++i){
- TGeoMaterial* mat = materials[i].first->GetMaterial();
- double length = materials[i].second;
- double nx0 = length / mat->GetRadLen();
- sum_x0 += nx0;
- double nLambda = length / mat->GetIntLen();
- sum_lambda += nLambda;
- path_length += length;
- end = p0 + path_length * direction;
- const char* fmt = mat->GetRadLen() >= 1e5 ? fmt2 : fmt1;
- ::printf(fmt, i+1, mat->GetName(), mat->GetZ(), mat->GetA(),
- mat->GetDensity(), mat->GetRadLen(), mat->GetIntLen(),
- length, path_length, sum_x0, sum_lambda, end[0], end[1], end[2]);
- //mat->Print();
+ MaterialScan scan(description);
+ scan.print(x0, y0, z0, x1, y1, z1);
+ if ( interactive ) {
+ char cmd[256];
+ description.apply("DD4hep_InteractiveUI",0,0);
+ ::snprintf(cmd,sizeof(cmd),
+ "dd4hep::rec::MaterialScan& gMaterialScan = "
+ "*(dd4hep::rec::MaterialScan*)%p",(void*)&scan);
+ gInterpreter->ProcessLine(cmd);
+ printout(ALWAYS,"materialScan","Use the ROOT interpreter variable gMaterialScan to interact with the material scanner");
+ gInterpreter->ProcessLine(".class dd4hep::rec::MaterialScan");
+ description.apply("DD4hep_Rint",0,0);
}
- printf("%s",line);
- const MaterialData& avg = matMgr.createAveragedMaterial(materials);
- const char* fmt = avg.radiationLength() >= 1e5 ? fmt2 : fmt1;
- end = p0 + path_length * direction;
- ::printf(fmt,0,"Average Material",avg.Z(),avg.A(),avg.density(),
- avg.radiationLength(), avg.interactionLength(),
- path_length, path_length,
- path_length/avg.radiationLength(),
- path_length/avg.interactionLength(),
- end[0], end[1], end[2]);
- printf("%s",line);
return 0;
}