packages/micromegas/CylinderGeomMicromegas.cc

0001 /*!
0002  * \file CylinderGeomMicromegas.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "CylinderGeomMicromegas.h"
0007
0008 #include <g4main/PHG4Hit.h>
0009 #include <trackbase/ActsGeometry.h>
0010
0011 #include <Acts/Surfaces/RectangleBounds.hpp>
0012 #include <Acts/Surfaces/SurfaceBounds.hpp>
0013 #include <TVector3.h>
0014
0015 #include <cassert>
0016
0017 namespace
0018 {
0019   // convenient square function
0020   template<class T>
0021     inline constexpr T square( const T& x ) { return x*x; }
0022
0023   //! radius
0024   template<class T>
0025     inline constexpr T get_r( T x, T y ) { return std::sqrt( square(x) + square(y) ); }
0026
0027   // bind angle to [-M_PI,+M_PI[. This is useful to avoid edge effects when making the difference between two angles
0028   template<class T>
0029     inline T bind_angle( const T& angle )
0030   {
0031     if( angle >= M_PI ) { return angle - 2*M_PI;
0032     } else if( angle < -M_PI ) { return angle + 2*M_PI;
0033     } else { return angle;
0034 }
0035   }
0036
0037 }
0038
0039 //________________________________________________________________________________
0040 TVector3 CylinderGeomMicromegas::get_local_from_world_coords( uint tileid, ActsGeometry* geometry, const TVector3& world_coordinates ) const
0041 {
0042   assert( tileid < m_tiles.size() );
0043   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0044   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0045
0046   const Acts::Vector3 global(
0047     world_coordinates.x()*Acts::UnitConstants::cm,
0048     world_coordinates.y()*Acts::UnitConstants::cm,
0049     world_coordinates.z()*Acts::UnitConstants::cm );
0050
0051   // convert to local
0052   /* this is equivalent to calling surface->globalToLocal but without the "on surface" check, and while returning a full Acts::Vector3 */
0053   const auto local = surface->transform(geometry->geometry().getGeoContext()).inverse()*global;
0054   return TVector3(
0055     local.x()/Acts::UnitConstants::cm,
0056     local.y()/Acts::UnitConstants::cm,
0057     local.z()/Acts::UnitConstants::cm );
0058 }
0059
0060 //________________________________________________________________________________
0061 TVector3 CylinderGeomMicromegas::get_local_from_world_vect( uint tileid, ActsGeometry* geometry, const TVector3& world_vect ) const
0062 {
0063   assert( tileid < m_tiles.size() );
0064   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0065   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0066
0067   const Acts::Vector3 global(
0068     world_vect.x()*Acts::UnitConstants::cm,
0069     world_vect.y()*Acts::UnitConstants::cm,
0070     world_vect.z()*Acts::UnitConstants::cm );
0071
0072   const Acts::Vector3 local = surface->referenceFrame(geometry->geometry().getGeoContext(), Acts::Vector3(), Acts::Vector3()).inverse()*global;
0073   return TVector3(
0074     local.x()/Acts::UnitConstants::cm,
0075     local.y()/Acts::UnitConstants::cm,
0076     local.z()/Acts::UnitConstants::cm );
0077 }
0078
0079 //________________________________________________________________________________
0080 TVector3 CylinderGeomMicromegas::get_world_from_local_coords( uint tileid, ActsGeometry* geometry, const TVector2& local_coordinates ) const
0081 {
0082   assert( tileid < m_tiles.size() );
0083   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0084   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0085
0086   const Acts::Vector2 local(
0087     local_coordinates.X()*Acts::UnitConstants::cm,
0088     local_coordinates.Y()*Acts::UnitConstants::cm );
0089
0090   // convert to global
0091   const auto global =  surface->localToGlobal(geometry->geometry().getGeoContext(), local, Acts::Vector3());
0092   return TVector3(
0093     global.x()/Acts::UnitConstants::cm,
0094     global.y()/Acts::UnitConstants::cm,
0095     global.z()/Acts::UnitConstants::cm );
0096 }
0097
0098 //________________________________________________________________________________
0099 TVector3 CylinderGeomMicromegas::get_world_from_local_coords( uint tileid, ActsGeometry* geometry, const TVector3& local_coordinates ) const
0100 {
0101   assert( tileid < m_tiles.size() );
0102   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0103   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0104
0105   const Acts::Vector3 local(
0106     local_coordinates.x()*Acts::UnitConstants::cm,
0107     local_coordinates.y()*Acts::UnitConstants::cm,
0108     local_coordinates.z()*Acts::UnitConstants::cm );
0109
0110   // convert to global
0111   /* this is equivalent to calling surface->localToGlobal but without assuming that the local point is on surface */
0112   const auto global =  surface->transform(geometry->geometry().getGeoContext())*local;
0113   return TVector3(
0114     global.x()/Acts::UnitConstants::cm,
0115     global.y()/Acts::UnitConstants::cm,
0116     global.z()/Acts::UnitConstants::cm );
0117 }
0118
0119 //________________________________________________________________________________
0120 TVector3 CylinderGeomMicromegas::get_world_from_local_vect( uint tileid, ActsGeometry* geometry, const TVector3& local_vect ) const
0121 {
0122   assert( tileid < m_tiles.size() );
0123   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0124   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0125
0126   Acts::Vector3 local(
0127     local_vect.x()*Acts::UnitConstants::cm,
0128     local_vect.y()*Acts::UnitConstants::cm,
0129     local_vect.z()*Acts::UnitConstants::cm );
0130
0131   const Acts::Vector3 global = surface->referenceFrame(geometry->geometry().getGeoContext(), Acts::Vector3(), Acts::Vector3())*local;
0132   return TVector3(
0133     global.x()/Acts::UnitConstants::cm,
0134     global.y()/Acts::UnitConstants::cm,
0135     global.z()/Acts::UnitConstants::cm );
0136 }
0137
0138 //________________________________________________________________________________
0139 int CylinderGeomMicromegas::find_tile_cylindrical( const TVector3& world_coordinates ) const
0140 {
0141   // convert to polar coordinates
0142   const auto phi = std::atan2( world_coordinates.y(), world_coordinates.x() );
0143   const auto z = world_coordinates.z();
0144
0145   for( size_t itile = 0; itile < m_tiles.size(); ++itile )
0146   {
0147     const auto& tile = m_tiles.at(itile);
0148
0149     if( std::abs( z - tile.m_centerZ ) > tile.m_sizeZ/2 ) { continue;
0150 }
0151     if( std::abs( bind_angle( phi - tile.m_centerPhi ) ) > tile.m_sizePhi/2 ) { continue;
0152 }
0153
0154     return itile;
0155   }
0156
0157   return -1;
0158 }
0159
0160 //________________________________________________________________________________
0161 int CylinderGeomMicromegas::find_strip_from_world_coords( uint tileid, ActsGeometry* geometry, const TVector3& world_coordinates ) const
0162 {
0163   // convert to local coordinates
0164   const auto local_coordinates = get_local_from_world_coords( tileid, geometry, world_coordinates );
0165
0166   // check against thickness
0167   if( std::abs(local_coordinates.z() ) >= m_thickness/2 )
0168   {
0169     std::cout << " CylinderGeomMicromegas::find_strip_from_world_coords - invalid world coordinates" << std::endl;
0170     return -1;
0171   } else {
0172     return find_strip_from_local_coords( tileid, geometry, { local_coordinates.x(), local_coordinates.y() } );
0173   }
0174 }
0175
0176 //________________________________________________________________________________
0177 int CylinderGeomMicromegas::find_strip_from_local_coords( uint tileid, ActsGeometry* geometry, const TVector2& local_coordinates ) const
0178 {
0179   assert( tileid < m_tiles.size() );
0180   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0181   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0182
0183   // get boundaries and corresponding dimension
0184   assert( surface->bounds().type() == Acts::SurfaceBounds::BoundsType::eRectangle );
0185   auto rectangle_bounds = static_cast<const Acts::RectangleBounds*>( &surface->bounds() );
0186   const auto half_length_x = rectangle_bounds->halfLengthX()/Acts::UnitConstants::cm;
0187   const auto half_length_y = rectangle_bounds->halfLengthY()/Acts::UnitConstants::cm;
0188
0189   // check azimuth
0190   if( std::abs( local_coordinates.X() ) >  half_length_x ) { return -1;
0191 }
0192
0193   // check z extend
0194   if( std::abs( local_coordinates.Y() ) > half_length_y ) { return -1;
0195 }
0196
0197   // calculate strip index, depending on segmentation
0198   switch( m_segmentation_type )
0199   {
0200     case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0201     return (int) std::floor((local_coordinates.X() + half_length_x)/m_pitch);
0202
0203     case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0204     return (int) std::floor((local_coordinates.Y() + half_length_y)/m_pitch);
0205   }
0206
0207   // unreachable
0208   return -1;
0209 }
0210
0211 //________________________________________________________________________________
0212 double CylinderGeomMicromegas::get_strip_length( uint tileid, ActsGeometry* geometry ) const
0213 {
0214   assert( tileid < m_tiles.size() );
0215   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0216   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0217
0218   // get boundaries and return corresponding dimension
0219   assert( surface->bounds().type() == Acts::SurfaceBounds::BoundsType::eRectangle );
0220   auto rectangle_bounds = static_cast<const Acts::RectangleBounds*>( &surface->bounds() );
0221   switch( m_segmentation_type )
0222   {
0223     case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0224     return 2.*rectangle_bounds->halfLengthY()/Acts::UnitConstants::cm;
0225
0226     case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0227     return 2.*rectangle_bounds->halfLengthX()/Acts::UnitConstants::cm;
0228   }
0229
0230   // unreachable
0231   return 0;
0232 }
0233
0234 //________________________________________________________________________________
0235 uint CylinderGeomMicromegas::get_strip_count( uint tileid, ActsGeometry* geometry ) const
0236 {
0237   assert( tileid < m_tiles.size() );
0238   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0239   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0240
0241   // get boundaries and corresponding dimension
0242   assert( surface->bounds().type() == Acts::SurfaceBounds::BoundsType::eRectangle );
0243   auto rectangle_bounds = static_cast<const Acts::RectangleBounds*>( &surface->bounds() );
0244   switch( m_segmentation_type )
0245   {
0246     case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0247     return std::floor( (2.*rectangle_bounds->halfLengthX()/Acts::UnitConstants::cm)/m_pitch );
0248
0249     case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0250     return std::floor( (2.*rectangle_bounds->halfLengthY()/Acts::UnitConstants::cm)/m_pitch );
0251   }
0252
0253   // unreachable
0254   return 0;
0255 }
0256
0257 //________________________________________________________________________________
0258 TVector2 CylinderGeomMicromegas::get_local_coordinates( uint tileid, ActsGeometry* geometry, uint stripnum ) const
0259 {
0260   assert( tileid < m_tiles.size() );
0261   const auto hitsetkey = MicromegasDefs::genHitSetKey(get_layer(), get_segmentation_type(), tileid);
0262   const auto surface = geometry->maps().getMMSurface(hitsetkey);
0263
0264   // get boundaries and return corresponding dimension
0265   assert( surface->bounds().type() == Acts::SurfaceBounds::BoundsType::eRectangle );
0266   auto rectangle_bounds = static_cast<const Acts::RectangleBounds*>( &surface->bounds() );
0267
0268   switch( m_segmentation_type )
0269   {
0270     case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0271     return TVector2( (0.5+stripnum)*m_pitch-rectangle_bounds->halfLengthX()/Acts::UnitConstants::cm, 0 );
0272
0273     case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0274     return TVector2( 0, (0.5+stripnum)*m_pitch-rectangle_bounds->halfLengthY()/Acts::UnitConstants::cm );
0275   }
0276
0277   // unreachable
0278   return TVector2();
0279 }
0280
0281 //________________________________________________________________________________
0282 TVector3 CylinderGeomMicromegas::get_world_coordinates( uint tileid, ActsGeometry* geometry, uint stripnum ) const
0283 { return get_world_from_local_coords( tileid, geometry, get_local_coordinates( tileid, geometry, stripnum ) ); }
0284
0285 //________________________________________________________________________________
0286 CylinderGeomMicromegas::range_t CylinderGeomMicromegas::get_phi_range(uint tileid, ActsGeometry* geometry ) const
0287 {
0288   switch( m_segmentation_type )
0289   {
0290     // for phi views we use the phi of the first and last strip, at middle of strip
0291     case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0292     {
0293       const auto strip_count = get_strip_count(tileid, geometry);
0294       const auto first_strip_center = get_world_coordinates(tileid, geometry, 0);
0295       const auto last_strip_center = get_world_coordinates(tileid, geometry, strip_count-1);
0296       return {
0297         std::atan2(first_strip_center.y(), first_strip_center.x()),
0298         std::atan2(last_strip_center.y(), last_strip_center.x())
0299       };
0300     }
0301
0302     // for z views we use the begin and end phi of the center strip
0303     case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0304     {
0305       const auto strip_count = get_strip_count(tileid, geometry);
0306       const auto strip_length = get_strip_length(tileid, geometry);
0307       const auto mid_strip_center_local = get_local_coordinates( tileid, geometry, strip_count/2 );
0308       const auto mid_strip_begin = get_world_from_local_coords(tileid, geometry, mid_strip_center_local + TVector2(-strip_length/2,0));
0309       const auto mid_strip_end = get_world_from_local_coords(tileid, geometry, mid_strip_center_local + TVector2(strip_length/2,0));
0310       return {
0311         std::atan2(mid_strip_begin.y(), mid_strip_begin.x()),
0312         std::atan2(mid_strip_end.y(), mid_strip_end.x())
0313       };
0314     }
0315
0316   }
0317
0318   // unreachable
0319   return {};
0320 }
0321
0322 //________________________________________________________________________________
0323 CylinderGeomMicromegas::range_t CylinderGeomMicromegas::get_theta_range(uint tileid, ActsGeometry* geometry ) const
0324 {
0325   switch( m_segmentation_type )
0326   {
0327     // for phi views we use the begin and end theta of the center strip
0328     case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0329     {
0330       const auto strip_count = get_strip_count(tileid, geometry);
0331       const auto strip_length = get_strip_length(tileid, geometry);
0332       const auto mid_strip_center_local = get_local_coordinates(tileid, geometry, strip_count/2);
0333       const auto mid_strip_begin = get_world_from_local_coords(tileid, geometry, mid_strip_center_local + TVector2(0,-strip_length/2));
0334       const auto mid_strip_end = get_world_from_local_coords( tileid, geometry, mid_strip_center_local + TVector2(0,strip_length/2));
0335       return {
0336         std::atan2(mid_strip_begin.z(), get_r(mid_strip_begin.x(), mid_strip_begin.y())),
0337         std::atan2(mid_strip_end.z(), get_r(mid_strip_end.x(), mid_strip_end.y()))
0338       };
0339     }
0340
0341     // for z views we use the theta of the first and last strip, at middle of strip
0342     case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0343     {
0344       const auto strip_count = get_strip_count(tileid, geometry);
0345       // need to add special care for tile 0, whose half detector is disconnected
0346       const auto first_strip_center = tileid == 0 ?
0347         get_world_coordinates(tileid, geometry, 128):
0348         get_world_coordinates(tileid, geometry, 0);
0349       const auto last_strip_center = get_world_coordinates(tileid, geometry, strip_count-1);
0350       return {
0351         std::atan2(first_strip_center.z(), get_r(first_strip_center.x(), first_strip_center.y())),
0352         std::atan2(last_strip_center.z(), get_r(last_strip_center.x(), last_strip_center.y()))
0353       };
0354     }
0355   }
0356
0357   // unreachable
0358   return {};
0359
0360 }
0361
0362 //________________________________________________________________________________
0363 void CylinderGeomMicromegas::identify( std::ostream& out ) const
0364 {
0365   out << "CylinderGeomMicromegas" << std::endl;
0366   out << "layer: " << m_layer << std::endl;
0367   out << "segmentation_type: " << (m_segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_PHI ? "SEGMENTATION_PHI":"SEGMENTATION_Z") << std::endl;
0368   out << "drift_direction: " << (m_drift_direction == MicromegasDefs::DriftDirection::INWARD ? "INWARD":"OUTWARD") << std::endl;
0369   out << "radius: " << m_radius << "cm" << std::endl;
0370   out << "thickness: " << m_thickness << "cm" << std::endl;
0371   out << "zmin: " << m_zmin << "cm" << std::endl;
0372   out << "zmax: " << m_zmax << "cm" << std::endl;
0373   out << "pitch: " << m_pitch << "cm" << std::endl;
0374   out << std::endl;
0375 }
0376
0377 //________________________________________________________________________________
0378 bool  CylinderGeomMicromegas::check_radius( const TVector3& world_coordinates ) const
0379 {
0380   const auto radius = get_r( world_coordinates.x(), world_coordinates.y() );
0381   return std::abs( radius - m_radius ) <= m_thickness/2;
0382 }