DDRec/src/MaterialManager.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     : F.Gaede
0011 //
0012 //==========================================================================
0013 #include "DDRec/MaterialManager.h"
0014 #include "DD4hep/Exceptions.h"
0015 #include "DD4hep/Detector.h"
0016
0017 #include "TGeoVolume.h"
0018 #include "TGeoManager.h"
0019 #include "TGeoNode.h"
0020 #include "TVirtualGeoTrack.h"
0021
0022 #define MINSTEP 1.e-5
0023
0024 namespace dd4hep {
0025   namespace rec {
0026
0027     MaterialManager::MaterialManager(Volume world) : _mV(0), _m( Material() ), _p0(),_p1(),_pos() {
0028       _tgeoMgr = world->GetGeoManager();
0029     }
0030
0031     MaterialManager::~MaterialManager(){
0032
0033     }
0034
0035     const PlacementVec& MaterialManager::placementsBetween(const Vector3D& p0, const Vector3D& p1 , double eps) {
0036       materialsBetween(p0,p1,eps);
0037       return _placeV;
0038     }
0039
0040     const MaterialManager::ScanData MaterialManager::entriesBetween(const Vector3D& p0,
0041                                                                     const Vector3D& p1,
0042                                                                     double eps)  {
0043       const auto& materials = this->materialsBetween(p0, p1, eps);
0044       return { materials, this->_placeV };
0045     }
0046
0047     const MaterialVec& MaterialManager::materialsBetween(const Vector3D& p0, const Vector3D& p1 , double eps) {
0048       if( ( p0 != _p0 ) || ( p1 != _p1 ) ) {
0049         // A backup is needed to restore the state of the navigator after the track is done
0050         // see https://github.com/AIDASoft/DD4hep/issues/1413
0051         _tgeoMgr->DoBackupState();
0052         //---------------------------------------
0053         _mV.clear() ;
0054         _placeV.clear();
0055         //
0056         // algorithm copied from TGeoGearDistanceProperties.cc (A.Munnich):
0057         //
0058
0059         double startpoint[3], endpoint[3], direction[3];
0060         double L=0;
0061         for(unsigned int i=0; i<3; i++) {
0062           startpoint[i] = p0[i];
0063           endpoint[i]   = p1[i];
0064           direction[i] = endpoint[i] - startpoint[i];
0065           L+=direction[i]*direction[i];
0066         }
0067         double totDist = sqrt( L ) ;
0068
0069         //normalize direction
0070         for(unsigned int i=0; i<3; i++)
0071           direction[i]=direction[i]/totDist;
0072
0073         _tgeoMgr->AddTrack(0, 12 ) ; // electron neutrino
0074
0075         TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
0076
0077         //check if there is a node at startpoint
0078         if(!node1)
0079           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.");
0080
0081         while ( !_tgeoMgr->IsOutside() )  {
0082
0083           // TGeoNode *node2;
0084           // TVirtualGeoTrack *track;
0085
0086           // step to (and over) the next Boundary
0087           TGeoNode * node2 = _tgeoMgr->FindNextBoundaryAndStep( 500, 1) ;
0088
0089           if( !node2 || _tgeoMgr->IsOutside() )
0090             break;
0091
0092           const double *position    =  _tgeoMgr->GetCurrentPoint();
0093           const double *previouspos =  _tgeoMgr->GetLastPoint();
0094
0095           double length = _tgeoMgr->GetStep();
0096
0097           TVirtualGeoTrack *track = _tgeoMgr->GetLastTrack();
0098
0099           //protection against infinitive loop in root which should not happen, but well it does...
0100           //work around until solution within root can be found when the step gets very small e.g. 1e-10
0101           //and the next boundary is never reached
0102
0103 #if 1   //fg: is this still needed ?
0104           if( length < MINSTEP ) {
0105
0106             _tgeoMgr->SetCurrentPoint( position[0] + MINSTEP * direction[0],
0107                                        position[1] + MINSTEP * direction[1],
0108                                        position[2] + MINSTEP * direction[2] );
0109
0110             length = _tgeoMgr->GetStep();
0111             node2  = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
0112
0113             position    = _tgeoMgr->GetCurrentPoint();
0114             previouspos = _tgeoMgr->GetLastPoint();
0115           }
0116 #endif
0117           //    printf( " --  step length :  %1.8e %1.8e   %1.8e   %1.8e   %1.8e   %1.8e   %1.8e   - %s \n" , length ,
0118           //        position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
0119
0120           Vector3D posV( position ) ;
0121
0122           double currDistance = ( posV - p0 ).r() ;
0123
0124           // //if the next boundary is further than end point
0125           //  if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
0126           //  || fabs(position[2])>fabs(endpoint[2]))
0127
0128           //if we travelled too far:
0129           if( currDistance > totDist  ) {
0130
0131             length = sqrt( pow(endpoint[0]-previouspos[0],2) +
0132                            pow(endpoint[1]-previouspos[1],2) +
0133                            pow(endpoint[2]-previouspos[2],2)   );
0134
0135             track->AddPoint( endpoint[0], endpoint[1], endpoint[2], 0. );
0136
0137
0138             if( length > eps )   {
0139               _mV.emplace_back(node1->GetMedium(), length );
0140               _placeV.emplace_back(node1,length);
0141             }
0142             break;
0143           }
0144
0145           track->AddPoint( position[0], position[1], position[2], 0.);
0146
0147           if( length > eps )   {
0148             _mV.emplace_back(node1->GetMedium(), length);
0149             _placeV.emplace_back(node1,length);
0150           }
0151           node1 = node2;
0152         }
0153
0154
0155         //fg: protect against empty list:
0156         if( _mV.empty() ){
0157           _mV.emplace_back(node1->GetMedium(), totDist);
0158           _placeV.emplace_back(node1,totDist);
0159         }
0160
0161
0162         _tgeoMgr->ClearTracks();
0163
0164         _tgeoMgr->CleanGarbage();
0165
0166         _tgeoMgr->DoRestoreState();
0167
0168         //---------------------------------------
0169
0170         _p0 = p0 ;
0171         _p1 = p1 ;
0172       }
0173
0174       return _mV ;
0175     }
0176
0177
0178     const Material& MaterialManager::materialAt(const Vector3D& pos )   {
0179       if( pos != _pos ) {
0180         TGeoNode *node = _tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
0181         if( ! node ) {
0182           std::stringstream err ;
0183           err << " MaterialManager::material: No geometry node found at location: " << pos ;
0184           throw std::runtime_error( err.str() );
0185         }
0186         _m = Material( node->GetMedium() );
0187         _pv = node;
0188         _pos = pos ;
0189       }
0190       return _m ;
0191     }
0192
0193     PlacedVolume MaterialManager::placementAt(const Vector3D& pos )   {
0194       if( pos != _pos ) {
0195         TGeoNode *node = _tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
0196         if( ! node ) {
0197           std::stringstream err ;
0198           err << " MaterialManager::material: No geometry node found at location: " << pos ;
0199           throw std::runtime_error( err.str() );
0200         }
0201         _m = Material( node->GetMedium() );
0202         _pv = node;
0203         _pos = pos;
0204       }
0205       return _pv;
0206     }
0207
0208     MaterialData MaterialManager::createAveragedMaterial( const MaterialVec& materials ) {
0209
0210       std::stringstream sstr ;
0211
0212       double sum_l = 0 ;
0213       double sum_rho_l = 0 ;
0214       double sum_rho_l_over_A = 0 ;
0215       double sum_rho_l_Z_over_A = 0 ;
0216       //double sum_rho_l_over_x = 0 ;
0217       double sum_l_over_x = 0 ;
0218       //double sum_rho_l_over_lambda = 0 ;
0219       double sum_l_over_lambda = 0 ;
0220
0221       for(unsigned i=0,n=materials.size(); i<n ; ++i){
0222
0223         Material mat = materials[i].first ;
0224         double   l   = materials[i].second ;
0225
0226         if( i != 0 ) sstr << "_" ;
0227         sstr << mat.name() << "_" << l ;
0228
0229         double rho      = mat.density() ;
0230         double  A       = mat.A() ;
0231         double  Z       = mat.Z() ;
0232         double  x       = mat.radLength() ;
0233         double  lambda  = mat.intLength() ;
0234
0235         sum_l                 +=   l ;
0236         sum_rho_l             +=   rho * l  ;
0237         sum_rho_l_over_A      +=   rho * l / A ;
0238         sum_rho_l_Z_over_A    +=   rho * l * Z / A ;
0239         sum_l_over_x          +=   l / x ;
0240         sum_l_over_lambda     +=   l / lambda ;
0241         // sum_rho_l_over_x      +=   rho * l / x ;
0242         // sum_rho_l_over_lambda +=   rho * l / lambda ;
0243       }
0244
0245       double rho      =  sum_rho_l / sum_l ;
0246
0247       double  A       =  sum_rho_l / sum_rho_l_over_A ;
0248       double  Z       =  sum_rho_l_Z_over_A / sum_rho_l_over_A ;
0249
0250       // radiation and interaction lengths already given in cm - average by length
0251
0252       // double  x       =  sum_rho_l / sum_rho_l_over_x ;
0253       double  x       =  sum_l / sum_l_over_x ;
0254
0255       //     double  lambda  =  sum_rho_l / sum_rho_l_over_lambda ;
0256       double  lambda  =  sum_l / sum_l_over_lambda ;
0257
0258
0259       return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
0260
0261     }
0262
0263   } /* namespace rec */
0264 } /* namespace dd4hep */