ForwardDiskSectorBuilderFromDet.cc

Go to the documentation of this file.

#include "ForwardDiskSectorBuilderFromDet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
#include "DataFormats/GeometryVector/interface/VectorUtil.h"
#include "TrackingTools/DetLayers/interface/DetLayerException.h"

using namespace std;

// Warning, remember to assign this pointer to a ReferenceCountingPointer!
BoundDiskSector* 
ForwardDiskSectorBuilderFromDet::operator()( const vector<const GeomDet*>& dets) const
{
  // check that the dets are all at about the same radius and z 
  float rcheck = dets.front()->surface().position().perp();
  float zcheck = dets.front()->surface().position().z();
  for ( vector<const GeomDet*>::const_iterator i = dets.begin(); i != dets.end(); i++){
    float rdiff = (**i).surface().position().perp()-rcheck;
    if ( std::abs(rdiff) > 1.) 
      edm::LogError("TkDetLayers") << " ForwardDiskSectorBuilderFromDet: Trying to build Petal Wedge from " 
                   << "Dets at different radii !! Delta_r = " << rdiff ;
    float zdiff = zcheck - (**i).surface().position().z();
    if ( std::abs(zdiff) > 0.8) 
      edm::LogError("TkDetLayers") << " ForwardDiskSectorBuilderFromDet: Trying to build Petal Wedge from " 
                   << "Dets at different z positions !! Delta_z = " << zdiff ;
  }

  auto bo = computeBounds( dets );

  Surface::PositionType pos( bo.second.x(), bo.second.y(), bo.second.z() );
  Surface::RotationType rot = computeRotation( dets, pos);
  return new BoundDiskSector( pos, rot, bo.first);
}

pair<DiskSectorBounds*, GlobalVector>
ForwardDiskSectorBuilderFromDet::computeBounds( const vector<const GeomDet*>& dets) const
{
  // go over all corners and compute maximum deviations 
  float rmin = (**(dets.begin())).surface().position().perp();
  float rmax(rmin); 
  float zmin((**(dets.begin())).surface().position().z());
  float zmax(zmin);
  float phimin((**(dets.begin())).surface().position().phi());
  float phimax(phimin);


  for (vector<const GeomDet*>::const_iterator idet=dets.begin();
       idet != dets.end(); idet++) {
    vector<const GeomDet*> detUnits = (**idet).components();
    if( !detUnits.empty() ){
      for (vector<const GeomDet*>::const_iterator detu=detUnits.begin();
       detu!=detUnits.end(); detu++) {
    // edm::LogInfo(TkDetLayers) << " Builder: Position of detUnit :"<< (**detu).position() ;      
    vector<GlobalPoint> corners = computeTrapezoidalCorners(*detu) ;
    for (vector<GlobalPoint>::const_iterator i=corners.begin();
         i!=corners.end(); i++) {
      float r = i->perp();
      float z = i->z();
      float phi = i->phi();
      rmin = min( rmin, r);
      rmax = max( rmax, r);
      zmin = min( zmin, z);
      zmax = max( zmax, z);
      if ( Geom::phiLess( phi, phimin)) phimin = phi;
      if ( Geom::phiLess( phimax, phi)) phimax = phi;
    }
    // in addition to the corners we have to check the middle of the 
    // det +/- length/2, since the min (max) radius for typical fw
    // dets is reached there
    float rdet  = (**detu).position().perp();
    float len   = (**detu).surface().bounds().length();
    float width = (**detu).surface().bounds().width();
    
    GlobalVector xAxis = (**detu).toGlobal(LocalVector(1,0,0));
    GlobalVector yAxis = (**detu).toGlobal(LocalVector(0,1,0));
    GlobalVector perpDir = GlobalVector( (**detu).position() - GlobalPoint(0,0,(**detu).position().z()) );
    
    double xAxisCos = xAxis.unit().dot(perpDir.unit());
    double yAxisCos = yAxis.unit().dot(perpDir.unit());
    
    if( std::abs(xAxisCos) > std::abs(yAxisCos) ) {
      rmin = min( rmin, rdet-width/2.F);
      rmax = max( rmax, rdet+width/2.F);
    }else{
      rmin = min( rmin, rdet-len/2.F);
      rmax = max( rmax, rdet+len/2.F);
    }     
      }
    }else{
      vector<GlobalPoint> corners = computeTrapezoidalCorners(*idet) ;
      for (vector<GlobalPoint>::const_iterator i=corners.begin();
       i!=corners.end(); i++) {
    float r = i->perp();
    float z = i->z();
    float phi = i->phi();
    rmin = min( rmin, r);
    rmax = max( rmax, r);
    zmin = min( zmin, z);
    zmax = max( zmax, z);
    if ( Geom::phiLess( phi, phimin)) phimin = phi;
    if ( Geom::phiLess( phimax, phi)) phimax = phi;
      }
      // in addition to the corners we have to check the middle of the 
      // det +/- length/2, since the min (max) radius for typical fw
      // dets is reached there

      float rdet  = (**idet).position().perp();
      float len   = (**idet).surface().bounds().length();
      float width = (**idet).surface().bounds().width();

      GlobalVector xAxis = (**idet).toGlobal(LocalVector(1,0,0));
      GlobalVector yAxis = (**idet).toGlobal(LocalVector(0,1,0));
      GlobalVector perpDir = GlobalVector( (**idet).position() - GlobalPoint(0,0,(**idet).position().z()) );

      double xAxisCos = xAxis.unit().dot(perpDir.unit());
      double yAxisCos = yAxis.unit().dot(perpDir.unit());
      
      if( std::abs(xAxisCos) > std::abs(yAxisCos) ) {
    rmin = min( rmin, rdet-width/2.F);
    rmax = max( rmax, rdet+width/2.F);
      }else{
    rmin = min( rmin, rdet-len/2.F);
    rmax = max( rmax, rdet+len/2.F);
      }     
    }
  }
  
  
  if (!Geom::phiLess(phimin, phimax)) 
    edm::LogError("TkDetLayers") << " ForwardDiskSectorBuilderFromDet : " 
                 << "Something went wrong with Phi Sorting !" ;
  float zPos = (zmax+zmin)/2.;
  float phiWin = phimax - phimin;
  float phiPos = (phimax+phimin)/2.;
  float rmed = (rmin+rmax)/2.;
  if ( phiWin < 0. ) {
    if ( (phimin < Geom::pi() / 2.) || (phimax > -Geom::pi()/2.) ){
      edm::LogError("TkDetLayers") << " Debug: something strange going on, please check " ;
    }
    // edm::LogInfo(TkDetLayers) << " Wedge at pi: phi " << phimin << " " << phimax << " " << phiWin 
    //   << " " << 2.*Geom::pi()+phiWin << " " ;
    phiWin += 2.*Geom::pi();
    phiPos += Geom::pi(); 
  }
  GlobalVector pos( rmed*cos(phiPos), rmed*sin(phiPos), zPos);
  return make_pair(new DiskSectorBounds(rmin,rmax,zmin-zPos,zmax-zPos,phiWin), pos);
}

Surface::RotationType 
ForwardDiskSectorBuilderFromDet::computeRotation( const vector<const GeomDet*>& dets,
                          Surface::PositionType pos) const {
  
  GlobalVector yAxis = ( GlobalVector( pos.x(), pos.y(), 0.)).unit();
  
  GlobalVector zAxis( 0., 0., 1.);
  GlobalVector xAxis = yAxis.cross( zAxis);
  
  return Surface::RotationType( xAxis, yAxis);
}


vector<GlobalPoint> 
ForwardDiskSectorBuilderFromDet::computeTrapezoidalCorners( const GeomDet* det) const {


  const Plane& plane( det->specificSurface());
  
  const TrapezoidalPlaneBounds* myBounds( static_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
  
  /*
  if (myBounds == 0) {
    string errmsg="ForwardDiskSectorBuilderFromDet: problems with dynamic cast to trapezoidal bounds for DetUnits";
    throw DetLayerException(errmsg);
    edm::LogError("TkDetLayers") << errmsg ;
  }
  */

  // array<const float, 4> 
  auto const & parameters = (*myBounds).parameters();

  if ( parameters[0] == 0 ) {
    edm::LogError("TkDetLayers") << "ForwardDiskSectorBuilder: something weird going on !" ;
    edm::LogError("TkDetLayers") << " Trapezoidal parameters of GeomDet (L2/L1/T/H): " ;
    for (int i = 0; i < 4; i++ )     edm::LogError("TkDetLayers") << "  " 
                                  << 2.*parameters[i] 
                                  << "\n";  
  }


  float hbotedge = parameters[0];
  float htopedge = parameters[1];
  float hapothem = parameters[3];
  float hthick   = parameters[2];

  vector<GlobalPoint> corners;

  corners.push_back( plane.toGlobal( LocalPoint( -htopedge, hapothem, hthick)));
  corners.push_back( plane.toGlobal( LocalPoint( -htopedge, hapothem, -hthick)));
  corners.push_back( plane.toGlobal( LocalPoint(  htopedge, hapothem, hthick)));
  corners.push_back( plane.toGlobal( LocalPoint(  htopedge, hapothem, -hthick)));
  corners.push_back( plane.toGlobal( LocalPoint(  hbotedge, -hapothem, hthick)));
  corners.push_back( plane.toGlobal( LocalPoint(  hbotedge, -hapothem, -hthick)));
  corners.push_back( plane.toGlobal( LocalPoint( -hbotedge, -hapothem, hthick)));
  corners.push_back( plane.toGlobal( LocalPoint( -hbotedge, -hapothem, -hthick)));

  return corners;
}