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