/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoTracker/TkMSParametrization/src/MSLayer.cc

Go to the documentation of this file.

#include "RecoTracker/TkMSParametrization/interface/MSLayer.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
#include "MSLayersKeeper.h"

#include<iostream>


using namespace GeomDetEnumerators;
using namespace std;
template <class T> T sqr( T t) {return t*t;}

//----------------------------------------------------------------------
ostream& operator<<( ostream& s, const MSLayer & l) 
{
  s <<" face: "<<l.face() 
    <<" pos:"<<l.position()<<", "
    <<" range:"<<l.range()<<", "
    <<l.theX0Data;
  return s;
}
//----------------------------------------------------------------------
ostream& operator<<( ostream& s, const MSLayer::DataX0 & d)
{
  if (d.hasX0)  s << "x0="<<d.x0 <<" sumX0D="<<d.sumX0D;  
  else if (d.allLayers)  s << "x0 by MSLayersKeeper"; 
  else  s <<"empty DataX0"; 
  return s; 
}
//----------------------------------------------------------------------
//MP
MSLayer::MSLayer(const DetLayer* layer, DataX0 dataX0)
  : theFace(layer->location()), theSeqNum(layer->seqNum()), theX0Data(dataX0)
{
  const BarrelDetLayer* bl; const ForwardDetLayer * fl;
  theHalfThickness = layer->surface().bounds().thickness()/2;  

  switch (theFace) {
  case barrel : 
    bl = static_cast<const BarrelDetLayer* >(layer);
    thePosition = bl->specificSurface().radius();
    theRange = Range(-bl->surface().bounds().length()/2,
                     bl->surface().bounds().length()/2);
    break;
  case endcap : 
    fl = static_cast<const ForwardDetLayer* >(layer);
    thePosition = fl->position().z(); 
    theRange = Range(fl->specificSurface().innerRadius(),
                     fl->specificSurface().outerRadius());
    break;
  default:
    // should throw or simimal
    cout << " ** MSLayer ** unknown part - will not work!" <<endl; 
    break;
  } 
}
//----------------------------------------------------------------------
MSLayer::MSLayer(Location part, float position, Range range, float halfThickness,
                 DataX0 dataX0)
  : theFace(part), 
    thePosition(position), 
    theRange(range), 
    theHalfThickness(halfThickness),
    theSeqNum(-1),
    theX0Data(dataX0)
  { }


//----------------------------------------------------------------------
bool MSLayer::operator== (const MSLayer &o) const
{
  return  theFace == o.theFace && std::abs(thePosition-o.thePosition) < 1.e-3f;
}
//----------------------------------------------------------------------
bool MSLayer::operator< (const MSLayer & o) const
{

  if (theFace==barrel && o.theFace==barrel) 
    return thePosition < o.thePosition;
  else if (theFace==barrel && o.theFace==endcap)
    return thePosition < o.range().max(); 
  else if (theFace==endcap && o.theFace==endcap ) 
    return std::abs(thePosition) < std::abs(o.thePosition);
  else 
    return range().max() < o.thePosition; 
}

//----------------------------------------------------------------------
pair<PixelRecoPointRZ,bool> MSLayer::crossing( const PixelRecoLineRZ & line) const { 
  const float eps = 1.e-5;
  bool  inLayer = true;
  float value = (theFace==barrel) ? line.zAtR(thePosition) : line.rAtZ(thePosition);
  if (value > theRange.max()) { 
    value = theRange.max()-eps;
    inLayer = false;
  } else if (value < theRange.min() ) { 
    value = theRange.min()+eps;
    inLayer = false;
  }
  float z = thePosition;
  if (theFace==barrel) std::swap(z,value); // if barrel value is z
  return make_pair( PixelRecoPointRZ( value, z), inLayer);
}
pair<PixelRecoPointRZ,bool> MSLayer::crossing( const SimpleLineRZ & line) const { 
  const float eps = 1.e-5;
  bool  inLayer = true;
  float value = (theFace==barrel) ? line.zAtR(thePosition) : line.rAtZ(thePosition);
  if (value > theRange.max()) { 
    value = theRange.max()-eps;
    inLayer = false;
  } else if (value < theRange.min() ) { 
    value = theRange.min()+eps;
    inLayer = false;
  }
  float z = thePosition;
  if (theFace==barrel) std::swap(z,value); // if barrel value is z
  return make_pair( PixelRecoPointRZ( value, z), inLayer);
}

//----------------------------------------------------------------------
float MSLayer::distance2(const PixelRecoPointRZ & point) const
{
  float u = (theFace==barrel) ? point.r() : point.z();
  float v = (theFace==barrel) ? point.z() : point.r();
 
  float  du = std::abs(u-thePosition);
  if (theRange.inside(v))  return (du < theHalfThickness) ? 0.f : du*du;

  float dv = (v > theRange.max()) ?
            v-theRange.max() : theRange.min() - v;
  return sqr(du)+sqr(dv);
}


//----------------------------------------------------------------------
float MSLayer::x0(float cotTheta) const
{
  if likely(theX0Data.hasX0) {
    float OverSinTheta = std::sqrt(1.f+cotTheta*cotTheta);
    return (theFace==barrel) ? theX0Data.x0*OverSinTheta :
                               theX0Data.x0*OverSinTheta/std::abs(cotTheta);
  } else if (theX0Data.allLayers) {
    const MSLayer * dataLayer =
       theX0Data.allLayers->layers(cotTheta).findLayer(*this);
    if (dataLayer) return  dataLayer->x0(cotTheta);
  } 
  return 0.;
}

//----------------------------------------------------------------------
float MSLayer::sumX0D(float cotTheta) const
{
if likely(theX0Data.hasX0) {
    switch(theFace) {
    case barrel:  
      return theX0Data.sumX0D
        *std::sqrt( std::sqrt(   (1.f+cotTheta*cotTheta)
                                 / (1.f+theX0Data.cotTheta*theX0Data.cotTheta)
                                 )
                    );
    case endcap: 
      return (theX0Data.hasFSlope) ?  
         theX0Data.sumX0D 
             + theX0Data.slopeSumX0D * (1.f/cotTheta-1.f/theX0Data.cotTheta)
       : theX0Data.sumX0D;
    case invalidLoc: break;// make gcc happy
    }
  } else if (theX0Data.allLayers) {
    const MSLayer* dataLayer =
       theX0Data.allLayers->layers(cotTheta).findLayer(*this);
    if (dataLayer) return  dataLayer->sumX0D(cotTheta);
  }
  return 0.;

}