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