00001
00002
00003
00004
00005
00006
00007 #include <functional>
00008 #include <algorithm>
00009
00011
00012 #include "GaudiKernel/AlgFactory.h"
00013 #include "GaudiKernel/IRndmGenSvc.h"
00014 #include "GaudiKernel/RndmGenerators.h"
00015 #include "GaudiKernel/SystemOfUnits.h"
00016 #include "GaudiKernel/Vector3DTypes.h"
00017
00018
00019
00020 #include "DetDesc/ILVolume.h"
00021 #include "DetDesc/IPVolume.h"
00022 #include "DetDesc/ISolid.h"
00023 #include "DetDesc/Material.h"
00024 #include "DetDesc/SolidBox.h"
00025 #include "DetDesc/VolumeIntersectionIntervals.h"
00026
00027
00028
00029 #include "VolumeCheckAlg.h"
00030
00031
00032
00033 #include "boost/progress.hpp"
00034
00035
00036
00042
00043
00044
00048
00049 DECLARE_ALGORITHM_FACTORY(VolumeCheckAlg)
00050
00051
00052
00053
00058
00059 VolumeCheckAlg::VolumeCheckAlg
00060 ( const std::string& name ,
00061 ISvcLocator* svcloc )
00062 : GaudiHistoAlg ( name , svcloc )
00063 , m_volumeName ( "Undefined Volume" )
00064 , m_volume ( 0 )
00065 , m_minx ( -10 * Gaudi::Units::m )
00066 , m_maxx ( 10 * Gaudi::Units::m )
00067 , m_miny ( -10 * Gaudi::Units::m )
00068 , m_maxy ( 10 * Gaudi::Units::m )
00069 , m_minz ( -10 * Gaudi::Units::m )
00070 , m_maxz ( 10 * Gaudi::Units::m )
00071 , m_shotsSphere ( 50000 )
00072 , m_shotsXY ( 10000 )
00073 , m_shotsYZ ( 10000 )
00074 , m_shotsZX ( 10000 )
00075 , m_vrtx ( 3 , 0 )
00076 , m_vertex ( )
00077 {
00078 declareProperty ( "Volume" , m_volumeName ) ;
00079 declareProperty ( "MinX" , m_minx ) ;
00080 declareProperty ( "MinY" , m_miny ) ;
00081 declareProperty ( "MinZ" , m_minz ) ;
00082 declareProperty ( "MaxX" , m_maxx ) ;
00083 declareProperty ( "MaxY" , m_maxy ) ;
00084 declareProperty ( "MaxZ" , m_maxz ) ;
00085 declareProperty ( "Shots3D" , m_shotsSphere ) ;
00086 declareProperty ( "ShotsXY" , m_shotsXY ) ;
00087 declareProperty ( "ShotsYZ" , m_shotsYZ ) ;
00088 declareProperty ( "ShotsZX" , m_shotsZX ) ;
00089 declareProperty ( "Null" , m_vrtx ) ;
00090 };
00091
00092
00093
00095
00096 VolumeCheckAlg::~VolumeCheckAlg() {};
00097
00098
00099
00104
00105 StatusCode VolumeCheckAlg::initialize()
00106 {
00107
00108 StatusCode sc = GaudiHistoAlg::initialize() ;
00109 if ( sc.isFailure() ) { return sc ; }
00110
00111 Assert ( 0 != randSvc () , "randSvc() points to NULL!" ) ;
00112 Assert ( 0 != detSvc () , "detSvc() points to NULL!" ) ;
00113
00114 m_volume = getDet<ILVolume>( m_volumeName ) ;
00115
00116 const std::string stars(80,'*');
00117
00118 info() << stars << endreq ;
00119 info() << m_volume << endreq ;
00120 info() << stars << endreq ;
00121 if ( !m_volume->pvolumes().empty() )
00122 {
00123 info() << " Has " << m_volume->pvolumes().size()
00124 << " daughter volumes " << endreq ;
00125 for ( ILVolume::PVolumes::const_iterator pv = m_volume->pvBegin() ;
00126 m_volume->pvEnd() != pv ; ++pv )
00127 {
00128 info() << " **** Daughter Physical Volume #"
00129 << pv - m_volume->pvBegin() << endreq ;
00130 info() << *pv << endreq ;
00131 }
00132 }
00133
00134 if ( !m_volume->isAssembly() )
00135 {
00136 const ISolid* solid = m_volume->solid() ;
00137 info() << " **** Solid " << endreq ;
00138 info() << solid << endreq ;
00139 const Material* material = m_volume->material() ;
00140 info() << " **** Material " << endreq ;
00141 info() << material << endreq ;
00142 }
00143
00144
00145 if ( m_vrtx.size() <= 3 )
00146 { while ( 3 != m_vrtx.size() ) { m_vrtx.push_back( 0.0 ); } }
00147 else
00148 {
00149 warning() << " Ignore extra fields in 'ShootingPoint' "<< endreq ;
00150 }
00151 m_vertex.SetXYZ( m_vrtx[0], m_vrtx[1], m_vrtx[2] ) ;
00152
00153 if ( !m_volume->isAssembly() && 0 != m_volume->solid() )
00154 {
00155 const ISolid* top = m_volume->solid()->coverTop();
00156 if( 0 == top )
00157 { return Error( "CoverTop* points to NULL!" ) ; }
00158 const SolidBox* box = dynamic_cast<const SolidBox*> ( top );
00159 if( 0 == box )
00160 { return Error("SolidBox* points to NULL!"); }
00161
00162 m_minx = -1 * box->xHalfLength() * 1.05 ;
00163 m_maxx = box->xHalfLength() * 1.05 ;
00164 m_miny = -1 * box->yHalfLength() * 1.05 ;
00165 m_maxy = box->yHalfLength() * 1.05 ;
00166 m_minz = -1 * box->zHalfLength() * 1.05 ;
00167 m_maxz = box->zHalfLength() * 1.05 ;
00168 }
00169
00170 return StatusCode::SUCCESS;
00171 };
00172
00173
00174
00175
00180
00181 StatusCode VolumeCheckAlg::execute()
00182 {
00183
00184 {
00185 info() << " Start 3D shots #" << m_shotsSphere << endreq ;
00186
00187 Rndm::Numbers flat( randSvc() , Rndm::Flat( -1.0 , 1.0 ) );
00188
00189 boost::progress_display progress ( m_shotsSphere ) ;
00190 for ( int shot = 0 ; shot < m_shotsSphere ; ++shot )
00191 {
00192 Gaudi::XYZVector vect;
00193
00194 double s = 0 ;
00195 do
00196 {
00197 vect.SetX( flat.shoot() ) ;
00198 vect.SetY( flat.shoot() ) ;
00199 s = vect.mag2() ;
00200 }
00201 while ( s > 1.0 || s == 0.0 );
00202
00203 vect.SetZ( -1 + 2 * s );
00204 const double a = 2. * sqrt( 1 - s ) ;
00205
00206 vect.SetX( vect.x() * a ) ;
00207 vect.SetY( vect.y() * a ) ;
00208
00209 ILVolume::Intersections intersections;
00210 m_volume->intersectLine ( m_vertex , vect , intersections , 0 );
00211
00212 const double radLength =
00213 std::accumulate
00214 ( intersections.begin() ,
00215 intersections.end () ,
00216 0.0 ,
00217 VolumeIntersectionIntervals::AccumulateIntersections() );
00218
00219 plot2D ( vect.phi () / Gaudi::Units::degree ,
00220 vect.theta () / Gaudi::Units::degree ,
00221 1 , " 3D-Material Budget for " + m_volumeName ,
00222 -180.0 , 180.0 ,
00223 0.0 , 180.0 ,
00224 50 , 50 ,
00225 radLength ) ;
00226
00227 plot2D ( vect.phi () / Gaudi::Units::degree ,
00228 vect.theta () / Gaudi::Units::degree ,
00229 2 , " 3D-Material Budget (N) for " + m_volumeName ,
00230 -180.0 , 180.0 ,
00231 0.0 , 180.0 ,
00232 50 , 50 ) ;
00233
00234
00235 ++progress ;
00236 }
00237 }
00238
00239
00240 {
00241 info() << " Start XY shots #" << m_shotsXY << endreq ;
00242
00243 Rndm::Numbers flat( randSvc() , Rndm::Flat( 0.0 , 1.0 ) );
00244
00245 boost::progress_display progress ( m_shotsXY ) ;
00246 for ( int shot = 0 ; shot < m_shotsXY ; ++shot )
00247 {
00248 const double x = m_minx + flat.shoot() * ( m_maxx - m_minx ) ;
00249 const double y = m_miny + flat.shoot() * ( m_maxy - m_miny ) ;
00250 const double z = m_minz + flat.shoot() * ( m_maxz - m_minz ) ;
00251 Gaudi::XYZPoint point( x , y , z ) ;
00252 Gaudi::XYZVector vect ( 0 , 0 , 1 ) ;
00253
00254 ILVolume::Intersections intersections;
00255 m_volume->intersectLine ( point , vect , intersections , 0 );
00256
00257 const double radLength =
00258 std::accumulate
00259 ( intersections.begin() ,
00260 intersections.end () ,
00261 0.0 ,
00262 VolumeIntersectionIntervals::AccumulateIntersections() );
00263
00264 plot2D ( point.x() , point.y() ,
00265 3 , " XY-Material Budget for " + m_volumeName ,
00266 m_minx , m_maxx ,
00267 m_miny , m_maxy ,
00268 50 , 50 ,
00269 radLength ) ;
00270
00271 plot2D ( point.x() , point.y() ,
00272 4 , " XY-Material Budget (N) for " + m_volumeName ,
00273 m_minx , m_maxx ,
00274 m_miny , m_maxy ,
00275 50 , 50 ) ;
00276
00277
00278 ++progress ;
00279 }
00280 }
00281
00282
00283 {
00284 info() << " Start YZ shots #" << m_shotsYZ << endreq ;
00285
00286 Rndm::Numbers flat( randSvc() , Rndm::Flat( 0.0 , 1.0 ) );
00287
00288 boost::progress_display progress ( m_shotsYZ ) ;
00289 for ( int shot = 0 ; shot < m_shotsYZ ; ++shot )
00290 {
00291 const double x = m_minx + flat.shoot() * ( m_maxx - m_minx ) ;
00292 const double y = m_miny + flat.shoot() * ( m_maxy - m_miny ) ;
00293 const double z = m_minz + flat.shoot() * ( m_maxz - m_minz ) ;
00294 Gaudi::XYZPoint point( x , y , z ) ;
00295 Gaudi::XYZVector vect ( 1 , 0 , 0 ) ;
00296
00297 ILVolume::Intersections intersections;
00298 m_volume->intersectLine ( point , vect , intersections , 0 );
00299
00300 const double radLength =
00301 std::accumulate
00302 ( intersections.begin() ,
00303 intersections.end () ,
00304 0.0 ,
00305 VolumeIntersectionIntervals::AccumulateIntersections() );
00306
00307 plot2D ( point.y() , point.z() ,
00308 5 , " YZ-Material Budget for " + m_volumeName ,
00309 m_miny , m_maxy ,
00310 m_minz , m_maxz ,
00311 50 , 50 ,
00312 radLength ) ;
00313
00314 plot2D ( point.y() , point.z() ,
00315 6 , " YZ-Material Budget (N) for " + m_volumeName ,
00316 m_miny , m_maxy ,
00317 m_minz , m_maxz ,
00318 50 , 50 ) ;
00319
00320
00321
00322 ++progress ;
00323
00324 }
00325 }
00326
00327
00328 {
00329 info() << " Start ZX shots #" << m_shotsZX << endreq ;
00330
00331 Rndm::Numbers flat( randSvc() , Rndm::Flat( 0.0 , 1.0 ) );
00332
00333 boost::progress_display progress ( m_shotsZX ) ;
00334 for ( int shot = 0 ; shot < m_shotsZX ; ++shot )
00335 {
00336
00337 const double x = m_minx + flat.shoot() * ( m_maxx - m_minx ) ;
00338 const double y = m_miny + flat.shoot() * ( m_maxy - m_miny ) ;
00339 const double z = m_minz + flat.shoot() * ( m_maxz - m_minz ) ;
00340 Gaudi::XYZPoint point( x , y , z ) ;
00341 Gaudi::XYZVector vect ( 0 , 1 , 0 ) ;
00342
00343 ILVolume::Intersections intersections;
00344 m_volume->intersectLine ( point , vect , intersections , 0 );
00345
00346 const double radLength =
00347 std::accumulate
00348 ( intersections.begin() ,
00349 intersections.end () ,
00350 0.0 ,
00351 VolumeIntersectionIntervals::AccumulateIntersections() );
00352
00353 plot2D ( point.y() , point.z() ,
00354 7 , " ZX-Material Budget for " + m_volumeName ,
00355 m_minz , m_maxz ,
00356 m_minx , m_maxx ,
00357 50 , 50 ,
00358 radLength ) ;
00359
00360 plot2D ( point.y() , point.z() ,
00361 8 , " ZX-Material Budget (N) for " + m_volumeName ,
00362 m_minz , m_maxz ,
00363 m_minx , m_maxx ,
00364 50 , 50 ) ;
00365
00366
00367 ++progress ;
00368
00369 }
00370 }
00371
00372 return StatusCode::SUCCESS;
00373 };
00374
00375
00376
00377
00378
00379
00380
