BladeShapeBuilderFromDet.cc

Go to the documentation of this file.

#include "BladeShapeBuilderFromDet.h"

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

#include "TrackingTools/DetLayers/interface/PhiLess.h"
#include "DataFormats/GeometrySurface/interface/BoundingBox.h"

#include <iomanip>

using namespace std;


namespace {
  
  pair<DiskSectorBounds *, GlobalVector>
  computeBounds( const vector<const GeomDet*>& dets,
         const Plane& plane)
  {
    Surface::PositionType tmpPos = dets.front()->surface().position();
    
    
    float rmin(plane.toLocal(tmpPos).perp());
    float rmax(plane.toLocal(tmpPos).perp());
    float zmin(plane.toLocal(tmpPos).z());
    float zmax(plane.toLocal(tmpPos).z());
    float phimin(plane.toLocal(tmpPos).phi());
    float phimax(plane.toLocal(tmpPos).phi());
    
    for(vector<const GeomDet*>::const_iterator it=dets.begin(); it!=dets.end(); it++){
      vector<GlobalPoint> corners = BoundingBox().corners( (*it)->specificSurface() );
      
      for(vector<GlobalPoint>::const_iterator i=corners.begin();
      i!=corners.end(); i++) {
    
    float r   = plane.toLocal(*i).perp();
    float z   = plane.toLocal(*i).z();
    float phi = plane.toLocal(*i).phi();
    rmin = min( rmin, r);
    rmax = max( rmax, r);
    zmin = min( zmin, z);
    zmax = max( zmax, z);
    if ( PhiLess()( phi, phimin)) phimin = phi;
    if ( 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 = (*it)->position().perp();
      float height  = (*it)->surface().bounds().width();
      rmin = min( rmin, rdet-height/2.F);
      rmax = max( rmax, rdet+height/2.F);  
      
      
    }
    
    if (!PhiLess()(phimin, phimax)) 
      edm::LogError("TkDetLayers") << " BladeShapeBuilderFromDet : " 
                   << "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") << " 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(); 
    }
    
    LocalVector localPos( rmed*cos(phiPos), rmed*sin(phiPos), zPos);
    
    LogDebug("TkDetLayers") << "localPos in computeBounds: " << localPos << "\n"
                << "rmin:   " << rmin << "\n"
                << "rmax:   " << rmax << "\n"
                << "zmin:   " << zmin << "\n"
                << "zmax:   " << zmax << "\n"
                << "phiWin: " << phiWin ;
    
    return make_pair(new DiskSectorBounds(rmin,rmax,zmin,zmax,phiWin),
             plane.toGlobal(localPos) );
    
  }
  
  
  Surface::RotationType 
  computeRotation( const vector<const GeomDet*>& dets,
           const Surface::PositionType& meanPos)
  {
    const Plane& plane = dets.front()->surface();
    
    GlobalVector xAxis;
    GlobalVector yAxis;
    GlobalVector zAxis;
    
    GlobalVector planeXAxis    = plane.toGlobal( LocalVector( 1, 0, 0));
    GlobalPoint  planePosition = plane.position();
    
    if(planePosition.x()*planeXAxis.x()+planePosition.y()*planeXAxis.y() > 0.){
      yAxis = planeXAxis;
    }else{
      yAxis = -planeXAxis;
    }
    
    GlobalVector planeZAxis = plane.toGlobal( LocalVector( 0, 0, 1));
    if(planeZAxis.z()*planePosition.z() > 0.){
      zAxis = planeZAxis;
    }else{
      zAxis = -planeZAxis;
    }
    
    xAxis = yAxis.cross( zAxis);
    
    return Surface::RotationType( xAxis, yAxis);
  }
  
}

BoundDiskSector* 
BladeShapeBuilderFromDet::build( const vector<const GeomDet*>& dets)
{
  // find mean position
  typedef Surface::PositionType::BasicVectorType Vector;
  Vector posSum(0,0,0);
  for (vector<const GeomDet*>::const_iterator i=dets.begin(); i!=dets.end(); i++) {
    posSum += (**i).surface().position().basicVector();
  }
  Surface::PositionType meanPos( 0.,0.,posSum.z()/float(dets.size()) );

  // temporary plane - for the computation of bounds
  Surface::RotationType rotation = computeRotation( dets, meanPos);
  Plane tmpPlane( meanPos, rotation);


  auto bo = 
    computeBounds( dets,tmpPlane );
  GlobalPoint pos = meanPos+bo.second;
  //edm::LogInfo(TkDetLayers) << "global pos in operator: " << pos ;
  return new BoundDiskSector( pos, rotation, bo.first);
}