CaloSegment.cc

Go to the documentation of this file.

//FAMOS headers
#include "FastSimulation/CaloGeometryTools/interface/CaloSegment.h"
#include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
#include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer1Properties.h"
#include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer2Properties.h"
#include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
#include "FastSimulation/CalorimeterProperties/interface/ECALProperties.h"

#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

CaloSegment::CaloSegment(const CaloPoint& in,
             const CaloPoint& out,
             double si,
             double siX0, 
             double siL0,
             Material mat,
             const CaloGeometryHelper* myCalorimeter):
  entrance_(in),
  exit_(out),
  sentrance_(si),
  sX0entrance_(siX0),
  sL0entrance_(siL0),
  material_(mat)

{

  sexit_= sentrance_+std::sqrt((exit_-entrance_).mag2());
  // Change this. CaloProperties from FamosShower should be used instead 
  double radLenIncm=999999;
  double intLenIncm=999999;
  detector_=in.whichDetector();
  if(detector_!=out.whichDetector()&&mat!=CRACK&&mat!=GAP&&mat!=ECALHCALGAP) 
    {
      std::cout << " Problem in the segments " << detector_ << " " << out.whichDetector() <<std::endl;
    }
  switch (mat)
    {
    case PbWO4:
      {

    int det = 0;
    if (in.whichSubDetector()==EcalBarrel) det = 1;
    if (in.whichSubDetector()==EcalEndcap) det = 2;

    radLenIncm =
      myCalorimeter->ecalProperties(det)->radLenIncm();
    intLenIncm =
      myCalorimeter->ecalProperties(det)->interactionLength();
      }
      break;
    case CRACK:
      {
    radLenIncm=8.9;//cracks : Al
    intLenIncm=35.4;
      }
      break;
    case PS:
      {
    radLenIncm = 
      myCalorimeter->layer1Properties(1)->radLenIncm();
    intLenIncm = 
      myCalorimeter->layer1Properties(1)->interactionLength();
      }
      break;
    case HCAL:
      {
    radLenIncm = 
      myCalorimeter->hcalProperties(1)->radLenIncm();
    intLenIncm = 
      myCalorimeter->hcalProperties(1)->interactionLength();
      }
      break;
    case ECALHCALGAP:
      {
    // From Olga's & Patrick's talk PRS/JetMET 21 Sept 2004 
    radLenIncm = 22.3;
    intLenIncm = 140;
      }
      break;
    case PSEEGAP:
      {
    // according to Sunanda 0.19 X0 (0.08X0 of polyethylene), support (0.06X0 of aluminium)  + other stuff
    // in the geometry 12 cm between layer and entrance of EE. Polyethylene is rather 48 and Al 8.9 (PDG)
    // for the inLen, just rescale according to PDG (85cm)
    radLenIncm = myCalorimeter->layer2Properties(1)->pseeRadLenIncm();
    intLenIncm = myCalorimeter->layer2Properties(1)->pseeIntLenIncm();
      }
      break;
    default:
      radLenIncm=999999;
    }
  sX0exit_ = sX0entrance_+(sexit_-sentrance_)/radLenIncm;
  sL0exit_ = sL0entrance_+(sexit_-sentrance_)/intLenIncm;
  if(mat==GAP) 
    {
      sX0exit_=sX0entrance_;
      sL0exit_=sL0entrance_;
    }
  length_ = sexit_-sentrance_;
  X0length_ = sX0exit_-sX0entrance_;
  L0length_ = sL0exit_-sL0entrance_;
}

CaloSegment::XYZPoint 
CaloSegment::positionAtDepthincm(double depth) const
{
  if (depth<sentrance_||depth>sexit_) return XYZPoint();
  return XYZPoint(entrance_+((depth-sentrance_)/(sexit_-sentrance_)*(exit_-entrance_)));
}

CaloSegment::XYZPoint 
CaloSegment::positionAtDepthinX0(double depth) const
{
  if (depth<sX0entrance_||depth>sX0exit_) return XYZPoint();
  return XYZPoint(entrance_+((depth-sX0entrance_)/(sX0exit_-sX0entrance_)*(exit_-entrance_)));
}

CaloSegment::XYZPoint 
CaloSegment::positionAtDepthinL0(double depth) const
{
  if (depth<sL0entrance_||depth>sL0exit_) return XYZPoint();
  return XYZPoint(entrance_+((depth-sL0entrance_)/(sL0exit_-sL0entrance_)*(exit_-entrance_)));
}

double 
CaloSegment::x0FromCm(double cm) const{
  return sX0entrance_+cm/length_*X0length_;
}

std::ostream & operator<<(std::ostream& ost ,const CaloSegment& seg)
{
  ost << " DetId " ;
  if(!seg.entrance().getDetId().null()) 
    ost << seg.entrance().getDetId()() ;
  else 
    {
      ost << seg.entrance().whichDetector() ;
      //  ost<< " Entrance side " << seg.entrance().getSide()
      ost << " Point " << (math::XYZVector)seg.entrance() << std::endl;
    }
  ost  << "DetId " ;
  if(!seg.exit().getDetId().null()) 
    ost << seg.exit().getDetId()() ;
  else
    ost << seg.exit().whichDetector() ;

  //  ost << " Exit side " << seg.exit().getSide() 
  ost << " Point " << (math::XYZVector)seg.exit() << " " 
      << seg.length() << " cm " 
      << seg.X0length() << " X0 " 
      <<  seg.L0length() << " Lambda0 " ;
  switch (seg.material())
    {
    case CaloSegment::PbWO4:
      ost << "PbWO4 " ;
      break;
    case CaloSegment::CRACK:
      ost << "CRACK ";
      break;
    case CaloSegment::PS:
      ost << "PS ";
      break;
    case CaloSegment::HCAL:
      ost << "HCAL ";
      break;
    case CaloSegment::ECALHCALGAP:
      ost << "ECAL-HCAL GAP ";
      break;
    case CaloSegment::PSEEGAP:
      ost  << "PS-ECAL GAP";
      break;
    default:
      ost << "GAP " ;
    }
  return ost;
}