IsolatedParticles/src/CaloSimInfo.cc

0001 #include "DataFormats/DetId/interface/DetId.h"
0002 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0003 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0004 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
0005 #include "Calibration/IsolatedParticles/interface/CaloSimInfo.h"
0006 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0007
0008 #include <CLHEP/Units/PhysicalConstants.h>
0009 #include <CLHEP/Units/SystemOfUnits.h>
0010
0011 #include <sstream>
0012
0013 namespace spr {
0014
0015   double timeOfFlight(DetId id, const CaloGeometry* geo, bool debug) {
0016     GlobalPoint point = (id.det() == DetId::Hcal)
0017                             ? (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(id)))->getPosition(id)
0018                             : geo->getPosition(id);
0019     double R = point.mag();
0020     double tmp = R / CLHEP::c_light / CLHEP::ns;
0021     if (debug) {
0022       DetId::Detector det = id.det();
0023       int subdet = id.subdetId();
0024       double eta = point.eta();
0025       double theta = 2.0 * atan(exp(-std::abs(eta)));
0026       double dist = 0;
0027       if (det == DetId::Ecal) {
0028         if (subdet == static_cast<int>(EcalBarrel)) {
0029           const double rEB = spr::rFrontEB * CLHEP::cm;
0030           dist = rEB / sin(theta);
0031         } else if (subdet == static_cast<int>(EcalEndcap)) {
0032           const double zEE = spr::zFrontEE * CLHEP::cm;
0033           dist = zEE / cos(theta);
0034         } else {
0035           const double zES = spr::zFrontES * CLHEP::cm;
0036           dist = zES / cos(theta);
0037         }
0038       } else if (det == DetId::Hcal) {
0039         if (subdet == static_cast<int>(HcalBarrel)) {
0040           const double rHB = spr::rFrontHB * CLHEP::cm;
0041           dist = rHB / sin(theta);
0042         } else if (subdet == static_cast<int>(HcalEndcap)) {
0043           const double zHE = spr::zFrontHE * CLHEP::cm;
0044           dist = zHE / cos(theta);
0045         } else if (subdet == static_cast<int>(HcalOuter)) {
0046           const double rHO = spr::rFrontHO * CLHEP::cm;
0047           dist = rHO / sin(theta);
0048         } else {
0049           const double zHF = spr::zFrontHF * CLHEP::cm;
0050           dist = zHF / cos(theta);
0051         }
0052       }
0053       double tmp1 = dist / CLHEP::c_light / CLHEP::ns;
0054
0055       edm::LogVerbatim("IsoTrack") << "Detector " << det << "/" << subdet << " Eta/Theta " << eta << "/"
0056                                    << theta / CLHEP::deg << " Dist " << dist / CLHEP::cm << " R " << R << " TOF " << tmp
0057                                    << ":" << tmp1;
0058     }
0059     return tmp;
0060   }
0061
0062 }  // namespace spr