CaloSimInfo.cc

Go to the documentation of this file.

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "Calibration/IsolatedParticles/interface/CaloConstants.h"
#include "Calibration/IsolatedParticles/interface/CaloSimInfo.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"

#include <CLHEP/Units/PhysicalConstants.h>
#include <CLHEP/Units/SystemOfUnits.h>

#include <sstream>

namespace spr {

  double timeOfFlight(DetId id, const CaloGeometry* geo, bool debug) {
    GlobalPoint point = (id.det() == DetId::Hcal)
                            ? (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(id)))->getPosition(id)
                            : geo->getPosition(id);
    double R = point.mag();
    double tmp = R / CLHEP::c_light / CLHEP::ns;
    if (debug) {
      DetId::Detector det = id.det();
      int subdet = id.subdetId();
      double eta = point.eta();
      double theta = 2.0 * atan(exp(-std::abs(eta)));
      double dist = 0;
      if (det == DetId::Ecal) {
        if (subdet == static_cast<int>(EcalBarrel)) {
          const double rEB = spr::rFrontEB * CLHEP::cm;
          dist = rEB / sin(theta);
        } else if (subdet == static_cast<int>(EcalEndcap)) {
          const double zEE = spr::zFrontEE * CLHEP::cm;
          dist = zEE / cos(theta);
        } else {
          const double zES = spr::zFrontES * CLHEP::cm;
          dist = zES / cos(theta);
        }
      } else if (det == DetId::Hcal) {
        if (subdet == static_cast<int>(HcalBarrel)) {
          const double rHB = spr::rFrontHB * CLHEP::cm;
          dist = rHB / sin(theta);
        } else if (subdet == static_cast<int>(HcalEndcap)) {
          const double zHE = spr::zFrontHE * CLHEP::cm;
          dist = zHE / cos(theta);
        } else if (subdet == static_cast<int>(HcalOuter)) {
          const double rHO = spr::rFrontHO * CLHEP::cm;
          dist = rHO / sin(theta);
        } else {
          const double zHF = spr::zFrontHF * CLHEP::cm;
          dist = zHF / cos(theta);
        }
      }
      double tmp1 = dist / CLHEP::c_light / CLHEP::ns;

      edm::LogVerbatim("IsoTrack") << "Detector " << det << "/" << subdet << " Eta/Theta " << eta << "/"
                                   << theta / CLHEP::deg << " Dist " << dist / CLHEP::cm << " R " << R << " TOF " << tmp
                                   << ":" << tmp1;
    }
    return tmp;
  }

}  // namespace spr