File indexing completed on 2024-04-06 11:59:19
0001 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0002 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0003 #include "DataFormats/GeometrySurface/interface/Plane.h"
0004 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0006 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0011
0012 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
0013 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0014
0015 #include <iostream>
0016 #include <sstream>
0017
0018 namespace spr {
0019
0020 std::vector<spr::propagatedTrackID> propagateCosmicCALO(edm::Handle<reco::TrackCollection>& trkCollection,
0021 const CaloGeometry* geo,
0022 const MagneticField* bField,
0023 const std::string& theTrackQuality,
0024 bool debug) {
0025 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0026 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0027 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0028 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0029 std::vector<spr::propagatedTrackID> vdets;
0030
0031 unsigned int indx;
0032 reco::TrackCollection::const_iterator trkItr;
0033 for (trkItr = trkCollection->begin(), indx = 0; trkItr != trkCollection->end(); ++trkItr, indx++) {
0034 const reco::Track* pTrack = &(*trkItr);
0035 spr::propagatedTrackID vdet;
0036 vdet.trkItr = trkItr;
0037 vdet.ok = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack->quality(trackQuality_)) : true;
0038 vdet.detIdECAL = DetId(0);
0039 vdet.detIdHCAL = DetId(0);
0040 vdet.detIdEHCAL = DetId(0);
0041 if (debug)
0042 edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta()
0043 << " phi " << trkItr->phi() << " Flag " << vdet.ok;
0044 GlobalPoint vertex;
0045 GlobalVector momentum;
0046 int charge(pTrack->charge());
0047 if (((pTrack->innerPosition()).Perp2()) < ((pTrack->outerPosition()).Perp2())) {
0048 vertex = GlobalPoint(
0049 ((pTrack->innerPosition()).X()), ((pTrack->innerPosition()).Y()), ((pTrack->innerPosition()).Z()));
0050 momentum = GlobalVector(
0051 ((pTrack->innerMomentum()).X()), ((pTrack->innerMomentum()).Y()), ((pTrack->innerMomentum()).Z()));
0052 } else {
0053 vertex = GlobalPoint(
0054 ((pTrack->outerPosition()).X()), ((pTrack->outerPosition()).Y()), ((pTrack->outerPosition()).Z()));
0055 momentum = GlobalVector(
0056 ((pTrack->outerMomentum()).X()), ((pTrack->outerMomentum()).Y()), ((pTrack->outerMomentum()).Z()));
0057 }
0058 if (debug)
0059 edm::LogVerbatim("IsoTrack") << "Track charge " << charge << " p " << momentum << " position " << vertex;
0060 std::pair<math::XYZPoint, bool> info = spr::propagateECAL(vertex, momentum, charge, bField, debug);
0061 if (debug)
0062 edm::LogVerbatim("IsoTrack") << "Propagate to ECAL " << info.second << " at (" << info.first.x() << ", "
0063 << info.first.y() << ", " << info.first.z() << ")";
0064
0065 vdet.okECAL = info.second;
0066 if (vdet.okECAL) {
0067 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0068 vdet.etaECAL = point.eta();
0069 vdet.phiECAL = point.phi();
0070 if (std::abs(point.eta()) < spr::etaBEEcal) {
0071 vdet.detIdECAL = barrelGeom->getClosestCell(point);
0072 } else {
0073 if (endcapGeom)
0074 vdet.detIdECAL = endcapGeom->getClosestCell(point);
0075 else
0076 vdet.okECAL = false;
0077 }
0078 vdet.detIdEHCAL = gHB->getClosestCell(point);
0079 if (debug) {
0080 std::ostringstream st1;
0081 if (std::abs(point.eta()) < spr::etaBEEcal)
0082 st1 << EBDetId(vdet.detIdECAL);
0083 else
0084 st1 << EEDetId(vdet.detIdECAL);
0085 edm::LogVerbatim("IsoTrack") << "Point at ECAL (" << vdet.etaECAL << ", " << vdet.phiECAL << " " << st1.str()
0086 << " " << HcalDetId(vdet.detIdEHCAL);
0087 }
0088 }
0089 info = spr::propagateHCAL(vertex, momentum, charge, bField, debug);
0090 if (debug)
0091 edm::LogVerbatim("IsoTrack") << "Propagate to HCAL " << info.second << " at (" << info.first.x() << ", "
0092 << info.first.y() << ", " << info.first.z() << ")";
0093 vdet.okHCAL = info.second;
0094 if (vdet.okHCAL) {
0095 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0096 vdet.etaHCAL = point.eta();
0097 vdet.phiHCAL = point.phi();
0098 vdet.detIdHCAL = gHB->getClosestCell(point);
0099 }
0100 if (debug) {
0101 std::ostringstream st1;
0102 if (vdet.detIdECAL.subdetId() == EcalBarrel)
0103 st1 << (EBDetId)(vdet.detIdECAL);
0104 else
0105 st1 << (EEDetId)(vdet.detIdECAL);
0106 edm::LogVerbatim("IsoTrack") << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") "
0107 << st1.str() << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL)
0108 << " Or " << (HcalDetId)(vdet.detIdEHCAL);
0109 }
0110 vdets.push_back(vdet);
0111 }
0112
0113 if (debug) {
0114 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << vdets.size() << " tracks";
0115 for (unsigned int i = 0; i < vdets.size(); ++i) {
0116 std::ostringstream st1;
0117 if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
0118 st1 << (EBDetId)(vdets[i].detIdECAL);
0119 } else {
0120 st1 << (EEDetId)(vdets[i].detIdECAL);
0121 }
0122 edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL
0123 << ") " << st1.str() << " HCAL (" << vdets[i].okHCAL << ") "
0124 << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL);
0125 }
0126 }
0127 return vdets;
0128 }
0129
0130 std::vector<spr::propagatedTrackID> propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection,
0131 const CaloGeometry* geo,
0132 const MagneticField* bField,
0133 const std::string& theTrackQuality,
0134 bool debug) {
0135 std::vector<spr::propagatedTrackID> vdets;
0136 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, vdets, debug);
0137 return vdets;
0138 }
0139
0140 void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection,
0141 const CaloGeometry* geo,
0142 const MagneticField* bField,
0143 const std::string& theTrackQuality,
0144 std::vector<spr::propagatedTrackID>& vdets,
0145 bool debug) {
0146 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0147 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0148 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0149 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0150
0151 unsigned int indx;
0152 reco::TrackCollection::const_iterator trkItr;
0153 for (trkItr = trkCollection->begin(), indx = 0; trkItr != trkCollection->end(); ++trkItr, indx++) {
0154 const reco::Track* pTrack = &(*trkItr);
0155 spr::propagatedTrackID vdet;
0156 vdet.trkItr = trkItr;
0157 vdet.ok = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack->quality(trackQuality_)) : true;
0158 vdet.detIdECAL = DetId(0);
0159 vdet.detIdHCAL = DetId(0);
0160 vdet.detIdEHCAL = DetId(0);
0161 if (debug)
0162 edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta()
0163 << " phi " << trkItr->phi() << " Flag " << vdet.ok;
0164 std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack, bField, debug);
0165 vdet.okECAL = info.second;
0166 if (vdet.okECAL) {
0167 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0168 vdet.etaECAL = point.eta();
0169 vdet.phiECAL = point.phi();
0170 if (std::abs(point.eta()) < spr::etaBEEcal) {
0171 vdet.detIdECAL = barrelGeom->getClosestCell(point);
0172 } else {
0173 if (endcapGeom)
0174 vdet.detIdECAL = endcapGeom->getClosestCell(point);
0175 else
0176 vdet.okECAL = false;
0177 }
0178 vdet.detIdEHCAL = gHB->getClosestCell(point);
0179 }
0180 info = spr::propagateHCAL(pTrack, bField, debug);
0181 vdet.okHCAL = info.second;
0182 if (vdet.okHCAL) {
0183 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0184 vdet.etaHCAL = point.eta();
0185 vdet.phiHCAL = point.phi();
0186 vdet.detIdHCAL = gHB->getClosestCell(point);
0187 }
0188 if (debug) {
0189 std::ostringstream st1;
0190 if (vdet.detIdECAL.subdetId() == EcalBarrel)
0191 st1 << (EBDetId)(vdet.detIdECAL);
0192 else
0193 st1 << (EEDetId)(vdet.detIdECAL);
0194 edm::LogVerbatim("IsoTrack") << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") "
0195 << st1.str() << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL)
0196 << " Or " << (HcalDetId)(vdet.detIdEHCAL);
0197 }
0198 vdets.push_back(vdet);
0199 }
0200 if (debug) {
0201 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << vdets.size() << " tracks";
0202 for (unsigned int i = 0; i < vdets.size(); ++i) {
0203 std::ostringstream st1;
0204 if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
0205 st1 << (EBDetId)(vdets[i].detIdECAL);
0206 } else {
0207 st1 << (EEDetId)(vdets[i].detIdECAL);
0208 }
0209 edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL
0210 << ") " << st1.str() << " HCAL (" << vdets[i].okHCAL << ") "
0211 << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL);
0212 }
0213 }
0214 }
0215
0216 void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection,
0217 const CaloGeometry* geo,
0218 const MagneticField* bField,
0219 const std::string& theTrackQuality,
0220 std::vector<spr::propagatedTrackDirection>& trkDir,
0221 bool debug) {
0222 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0223 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0224 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0225 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0226
0227 unsigned int indx;
0228 reco::TrackCollection::const_iterator trkItr;
0229 for (trkItr = trkCollection->begin(), indx = 0; trkItr != trkCollection->end(); ++trkItr, indx++) {
0230 const reco::Track* pTrack = &(*trkItr);
0231 spr::propagatedTrackDirection trkD;
0232 trkD.trkItr = trkItr;
0233 trkD.ok = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack->quality(trackQuality_)) : true;
0234 trkD.detIdECAL = DetId(0);
0235 trkD.detIdHCAL = DetId(0);
0236 trkD.detIdEHCAL = DetId(0);
0237 if (debug)
0238 edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta()
0239 << " phi " << trkItr->phi() << " Flag " << trkD.ok;
0240 spr::propagatedTrack info = spr::propagateTrackToECAL(pTrack, bField, debug);
0241 GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
0242 trkD.okECAL = info.ok;
0243 trkD.pointECAL = point;
0244 trkD.directionECAL = info.direction;
0245 if (trkD.okECAL) {
0246 if (std::abs(info.point.eta()) < spr::etaBEEcal) {
0247 trkD.detIdECAL = barrelGeom->getClosestCell(point);
0248 } else {
0249 if (endcapGeom)
0250 trkD.detIdECAL = endcapGeom->getClosestCell(point);
0251 else
0252 trkD.okECAL = false;
0253 }
0254 trkD.detIdEHCAL = gHB->getClosestCell(point);
0255 }
0256 info = spr::propagateTrackToHCAL(pTrack, bField, debug);
0257 point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
0258 trkD.okHCAL = info.ok;
0259 trkD.pointHCAL = point;
0260 trkD.directionHCAL = info.direction;
0261 if (trkD.okHCAL) {
0262 trkD.detIdHCAL = gHB->getClosestCell(point);
0263 }
0264 trkDir.push_back(trkD);
0265 }
0266 if (debug) {
0267 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << trkDir.size() << " tracks";
0268 for (unsigned int i = 0; i < trkDir.size(); ++i) {
0269 std::ostringstream st1, st2;
0270 if (trkDir[i].okECAL) {
0271 st1 << " point " << trkDir[i].pointECAL << " direction " << trkDir[i].directionECAL << " ";
0272 if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
0273 st1 << (EBDetId)(trkDir[i].detIdECAL);
0274 } else {
0275 st1 << (EEDetId)(trkDir[i].detIdECAL);
0276 }
0277 }
0278 if (trkDir[i].okHCAL) {
0279 st2 << " point " << trkDir[i].pointHCAL << " direction " << trkDir[i].directionHCAL << " "
0280 << (HcalDetId)(trkDir[i].detIdHCAL);
0281 }
0282 edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL
0283 << ")" << st1.str() << " HCAL (" << trkDir[i].okHCAL << ")" << st2.str() << " Or "
0284 << (HcalDetId)(trkDir[i].detIdEHCAL);
0285 }
0286 }
0287 }
0288
0289 spr::propagatedTrackID propagateCALO(const reco::Track* pTrack,
0290 const CaloGeometry* geo,
0291 const MagneticField* bField,
0292 bool debug) {
0293 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0294 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0295 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0296
0297 spr::propagatedTrackID vdet;
0298 vdet.ok = true;
0299 vdet.detIdECAL = DetId(0);
0300 vdet.detIdHCAL = DetId(0);
0301 vdet.detIdEHCAL = DetId(0);
0302 if (debug)
0303 edm::LogVerbatim("IsoTrack") << "Propagate track: p " << pTrack->p() << " eta " << pTrack->eta() << " phi "
0304 << pTrack->phi() << " Flag " << vdet.ok;
0305 std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack, bField, debug);
0306 vdet.okECAL = info.second;
0307 if (vdet.okECAL) {
0308 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0309 vdet.etaECAL = point.eta();
0310 vdet.phiECAL = point.phi();
0311 if (std::abs(point.eta()) < spr::etaBEEcal) {
0312 vdet.detIdECAL = barrelGeom->getClosestCell(point);
0313 } else {
0314 if (endcapGeom)
0315 vdet.detIdECAL = endcapGeom->getClosestCell(point);
0316 else
0317 vdet.okECAL = false;
0318 }
0319 vdet.detIdEHCAL = gHB->getClosestCell(point);
0320 }
0321 info = spr::propagateHCAL(pTrack, bField, debug);
0322 vdet.okHCAL = info.second;
0323 if (vdet.okHCAL) {
0324 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0325 vdet.etaHCAL = point.eta();
0326 vdet.phiHCAL = point.phi();
0327 vdet.detIdHCAL = gHB->getClosestCell(point);
0328 }
0329 if (debug) {
0330 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for 1 track";
0331 std::ostringstream st1;
0332 if (vdet.detIdECAL.subdetId() == EcalBarrel) {
0333 st1 << (EBDetId)(vdet.detIdECAL);
0334 } else {
0335 st1 << (EEDetId)(vdet.detIdECAL);
0336 }
0337 edm::LogVerbatim("IsoTrack") << "Track [0] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") " << st1.str()
0338 << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or "
0339 << (HcalDetId)(vdet.detIdEHCAL);
0340 }
0341 return vdet;
0342 }
0343
0344 std::vector<spr::propagatedGenTrackID> propagateCALO(const HepMC::GenEvent* genEvent,
0345 const ParticleDataTable* pdt,
0346 const CaloGeometry* geo,
0347 const MagneticField* bField,
0348 double etaMax,
0349 bool debug) {
0350 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0351 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0352 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0353
0354 std::vector<spr::propagatedGenTrackID> trkDir;
0355 unsigned int indx;
0356 HepMC::GenEvent::particle_const_iterator p;
0357 for (p = genEvent->particles_begin(), indx = 0; p != genEvent->particles_end(); ++p, ++indx) {
0358 spr::propagatedGenTrackID trkD;
0359 trkD.trkItr = p;
0360 trkD.detIdECAL = DetId(0);
0361 trkD.detIdHCAL = DetId(0);
0362 trkD.detIdEHCAL = DetId(0);
0363 trkD.pdgId = ((*p)->pdg_id());
0364 trkD.charge = ((pdt->particle(trkD.pdgId))->ID().threeCharge()) / 3;
0365 const GlobalVector momentum = GlobalVector((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());
0366 if (debug)
0367 edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge
0368 << " p " << momentum;
0369
0370 if ((*p)->status() == 1 && std::abs((*p)->momentum().eta()) < etaMax) {
0371 const GlobalPoint vertex = GlobalPoint(0.1 * (*p)->production_vertex()->position().x(),
0372 0.1 * (*p)->production_vertex()->position().y(),
0373 0.1 * (*p)->production_vertex()->position().z());
0374 trkD.ok = true;
0375 spr::propagatedTrack info = spr::propagateCalo(
0376 vertex, momentum, trkD.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
0377 GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
0378 trkD.okECAL = info.ok;
0379 trkD.pointECAL = point;
0380 trkD.directionECAL = info.direction;
0381 if (trkD.okECAL) {
0382 if (std::abs(info.point.eta()) < spr::etaBEEcal) {
0383 trkD.detIdECAL = barrelGeom->getClosestCell(point);
0384 } else {
0385 if (endcapGeom)
0386 trkD.detIdECAL = endcapGeom->getClosestCell(point);
0387 else
0388 trkD.okECAL = false;
0389 }
0390 trkD.detIdEHCAL = gHB->getClosestCell(point);
0391 }
0392
0393 info = spr::propagateCalo(
0394 vertex, momentum, trkD.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
0395 point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
0396 trkD.okHCAL = info.ok;
0397 trkD.pointHCAL = point;
0398 trkD.directionHCAL = info.direction;
0399 if (trkD.okHCAL) {
0400 trkD.detIdHCAL = gHB->getClosestCell(point);
0401 }
0402 }
0403 trkDir.push_back(trkD);
0404 }
0405 if (debug) {
0406 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << trkDir.size() << " tracks";
0407 for (unsigned int i = 0; i < trkDir.size(); ++i) {
0408 std::ostringstream st1, st2;
0409 if (trkDir[i].okECAL) {
0410 st1 << " point " << trkDir[i].pointECAL << " direction " << trkDir[i].directionECAL << " ";
0411 if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
0412 st1 << (EBDetId)(trkDir[i].detIdECAL);
0413 } else {
0414 st1 << (EEDetId)(trkDir[i].detIdECAL);
0415 }
0416 }
0417 st2 << " HCAL (" << trkDir[i].okHCAL << ")";
0418 if (trkDir[i].okHCAL) {
0419 st2 << " point " << trkDir[i].pointHCAL << " direction " << trkDir[i].directionHCAL << " "
0420 << (HcalDetId)(trkDir[i].detIdHCAL);
0421 }
0422 edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL
0423 << ")" << st1.str() << st2.str() << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL);
0424 }
0425 }
0426 return trkDir;
0427 }
0428
0429 std::vector<spr::propagatedGenParticleID> propagateCALO(edm::Handle<reco::GenParticleCollection>& genParticles,
0430 const ParticleDataTable* pdt,
0431 const CaloGeometry* geo,
0432 const MagneticField* bField,
0433 double etaMax,
0434 bool debug) {
0435 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0436 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0437 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0438
0439 std::vector<spr::propagatedGenParticleID> trkDir;
0440 unsigned int indx;
0441 reco::GenParticleCollection::const_iterator p;
0442 for (p = genParticles->begin(), indx = 0; p != genParticles->end(); ++p, ++indx) {
0443 spr::propagatedGenParticleID trkD;
0444 trkD.trkItr = p;
0445 trkD.detIdECAL = DetId(0);
0446 trkD.detIdHCAL = DetId(0);
0447 trkD.detIdEHCAL = DetId(0);
0448 trkD.pdgId = (p->pdgId());
0449 trkD.charge = p->charge();
0450 const GlobalVector momentum = GlobalVector(p->momentum().x(), p->momentum().y(), p->momentum().z());
0451 if (debug)
0452 edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge
0453 << " p " << momentum;
0454
0455 if (p->status() == 1 && std::abs(momentum.eta()) < etaMax) {
0456 const GlobalPoint vertex = GlobalPoint(p->vertex().x(), p->vertex().y(), p->vertex().z());
0457 trkD.ok = true;
0458 spr::propagatedTrack info = spr::propagateCalo(
0459 vertex, momentum, trkD.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
0460 GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
0461 trkD.okECAL = info.ok;
0462 trkD.pointECAL = point;
0463 trkD.directionECAL = info.direction;
0464 if (trkD.okECAL) {
0465 if (std::abs(info.point.eta()) < spr::etaBEEcal) {
0466 trkD.detIdECAL = barrelGeom->getClosestCell(point);
0467 } else {
0468 if (endcapGeom)
0469 trkD.detIdECAL = endcapGeom->getClosestCell(point);
0470 else
0471 trkD.okECAL = false;
0472 }
0473 trkD.detIdEHCAL = gHB->getClosestCell(point);
0474 }
0475
0476 info = spr::propagateCalo(
0477 vertex, momentum, trkD.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
0478 point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
0479 trkD.okHCAL = info.ok;
0480 trkD.pointHCAL = point;
0481 trkD.directionHCAL = info.direction;
0482 if (trkD.okHCAL) {
0483 trkD.detIdHCAL = gHB->getClosestCell(point);
0484 }
0485 }
0486 trkDir.push_back(trkD);
0487 }
0488 if (debug) {
0489 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
0490 for (unsigned int i = 0; i < trkDir.size(); ++i) {
0491 std::ostringstream st1, st2;
0492 if (trkDir[i].okECAL) {
0493 st1 << " point " << trkDir[i].pointECAL << " direction " << trkDir[i].directionECAL << " ";
0494 if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
0495 st1 << (EBDetId)(trkDir[i].detIdECAL);
0496 } else {
0497 st1 << (EEDetId)(trkDir[i].detIdECAL);
0498 }
0499 }
0500 st2 << " HCAL (" << trkDir[i].okHCAL << ")";
0501 if (trkDir[i].okHCAL) {
0502 st2 << " point " << trkDir[i].pointHCAL << " direction " << trkDir[i].directionHCAL << " "
0503 << (HcalDetId)(trkDir[i].detIdHCAL);
0504 }
0505 edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL
0506 << ")" << st1.str() << st2.str() << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL);
0507 }
0508 }
0509 return trkDir;
0510 }
0511
0512 spr::propagatedTrackDirection propagateCALO(unsigned int thisTrk,
0513 edm::Handle<edm::SimTrackContainer>& SimTk,
0514 edm::Handle<edm::SimVertexContainer>& SimVtx,
0515 const CaloGeometry* geo,
0516 const MagneticField* bField,
0517 bool debug) {
0518 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0519 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0520 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0521
0522 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
0523 spr::propagatedTrackDirection trkD;
0524 trkD.ok = trk.ok;
0525 trkD.detIdECAL = DetId(0);
0526 trkD.detIdHCAL = DetId(0);
0527 trkD.detIdEHCAL = DetId(0);
0528 if (debug)
0529 edm::LogVerbatim("IsoTrack") << "Propagate track " << thisTrk << " charge " << trk.charge << " position "
0530 << trk.position << " p " << trk.momentum << " Flag " << trkD.ok;
0531 if (trkD.ok) {
0532 spr::propagatedTrack info = spr::propagateCalo(
0533 trk.position, trk.momentum, trk.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
0534 GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
0535 trkD.okECAL = info.ok;
0536 trkD.pointECAL = point;
0537 trkD.directionECAL = info.direction;
0538 if (trkD.okECAL) {
0539 if (std::abs(info.point.eta()) < spr::etaBEEcal) {
0540 trkD.detIdECAL = barrelGeom->getClosestCell(point);
0541 } else {
0542 if (endcapGeom)
0543 trkD.detIdECAL = endcapGeom->getClosestCell(point);
0544 else
0545 trkD.okECAL = false;
0546 }
0547 trkD.detIdEHCAL = gHB->getClosestCell(point);
0548 }
0549
0550 info = spr::propagateCalo(
0551 trk.position, trk.momentum, trk.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
0552 point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
0553 trkD.okHCAL = info.ok;
0554 trkD.pointHCAL = point;
0555 trkD.directionHCAL = info.direction;
0556 if (trkD.okHCAL) {
0557 trkD.detIdHCAL = gHB->getClosestCell(point);
0558 }
0559 }
0560 if (debug) {
0561 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL ("
0562 << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")";
0563 std::ostringstream st1;
0564 if (trkD.okECAL) {
0565 st1 << "ECAL point " << trkD.pointECAL << " direction " << trkD.directionECAL << " ";
0566 if (trkD.detIdECAL.subdetId() == EcalBarrel) {
0567 st1 << (EBDetId)(trkD.detIdECAL);
0568 } else {
0569 st1 << (EEDetId)(trkD.detIdECAL);
0570 }
0571 }
0572 if (trkD.okHCAL) {
0573 st1 << " HCAL point " << trkD.pointHCAL << " direction " << trkD.directionHCAL << " "
0574 << (HcalDetId)(trkD.detIdHCAL);
0575 }
0576 if (trkD.okECAL)
0577 st1 << " Or " << (HcalDetId)(trkD.detIdEHCAL);
0578 edm::LogVerbatim("IsoTrack") << st1.str();
0579 }
0580 return trkD;
0581 }
0582
0583 spr::propagatedTrackDirection propagateHCALBack(unsigned int thisTrk,
0584 edm::Handle<edm::SimTrackContainer>& SimTk,
0585 edm::Handle<edm::SimVertexContainer>& SimVtx,
0586 const CaloGeometry* geo,
0587 const MagneticField* bField,
0588 bool debug) {
0589 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0590 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
0591 spr::propagatedTrackDirection trkD;
0592 trkD.ok = trk.ok;
0593 trkD.detIdECAL = DetId(0);
0594 trkD.detIdHCAL = DetId(0);
0595 trkD.detIdEHCAL = DetId(0);
0596 if (debug)
0597 edm::LogVerbatim("IsoTrack") << "Propagate track " << thisTrk << " charge " << trk.charge << " position "
0598 << trk.position << " p " << trk.momentum << " Flag " << trkD.ok;
0599 if (trkD.ok) {
0600 spr::propagatedTrack info = spr::propagateCalo(
0601 trk.position, trk.momentum, trk.charge, bField, spr::zBackHE, spr::rBackHB, spr::etaBEHcal, debug);
0602 const GlobalPoint point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
0603 trkD.okHCAL = info.ok;
0604 trkD.pointHCAL = point;
0605 trkD.directionHCAL = info.direction;
0606 if (trkD.okHCAL) {
0607 trkD.detIdHCAL = gHB->getClosestCell(point);
0608 }
0609 }
0610 if (debug) {
0611 std::ostringstream st1;
0612 if (trkD.okHCAL) {
0613 st1 << " HCAL point " << trkD.pointHCAL << " direction " << trkD.directionHCAL << " "
0614 << (HcalDetId)(trkD.detIdHCAL);
0615 }
0616 edm::LogVerbatim("IsoTrack") << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL ("
0617 << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << st1.str();
0618 }
0619 return trkD;
0620 }
0621
0622 std::pair<bool, HcalDetId> propagateHCALBack(const reco::Track* track,
0623 const CaloGeometry* geo,
0624 const MagneticField* bField,
0625 bool debug) {
0626 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0627 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0628 const GlobalVector momentum(track->px(), track->py(), track->pz());
0629 int charge(track->charge());
0630 spr::propagatedTrack info =
0631 spr::propagateCalo(vertex, momentum, charge, bField, spr::zBackHE, spr::rBackHB, spr::etaBEHcal, debug);
0632 if (info.ok) {
0633 const GlobalPoint point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
0634 return std::pair<bool, HcalDetId>(true, HcalDetId(gHB->getClosestCell(point)));
0635 } else {
0636 return std::pair<bool, HcalDetId>(false, HcalDetId());
0637 }
0638 }
0639
0640 propagatedTrack propagateTrackToECAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
0641 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0642 const GlobalVector momentum(track->px(), track->py(), track->pz());
0643 int charge(track->charge());
0644 return spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
0645 }
0646
0647 propagatedTrack propagateTrackToECAL(unsigned int thisTrk,
0648 edm::Handle<edm::SimTrackContainer>& SimTk,
0649 edm::Handle<edm::SimVertexContainer>& SimVtx,
0650 const MagneticField* bfield,
0651 bool debug) {
0652 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
0653 spr::propagatedTrack ptrk;
0654 if (trk.ok)
0655 ptrk = spr::propagateCalo(
0656 trk.position, trk.momentum, trk.charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
0657 return ptrk;
0658 }
0659
0660 std::pair<math::XYZPoint, bool> propagateECAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
0661 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0662 const GlobalVector momentum(track->px(), track->py(), track->pz());
0663 int charge(track->charge());
0664 return spr::propagateECAL(vertex, momentum, charge, bfield, debug);
0665 }
0666
0667 std::pair<DetId, bool> propagateIdECAL(const HcalDetId& id,
0668 const CaloGeometry* geo,
0669 const MagneticField* bField,
0670 bool debug) {
0671 const HcalGeometry* gHB = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0672 const GlobalPoint vertex(0, 0, 0);
0673 const GlobalPoint hit(gHB->getPosition(id));
0674 const GlobalVector momentum = GlobalVector(hit.x(), hit.y(), hit.z());
0675 std::pair<math::XYZPoint, bool> info = propagateECAL(vertex, momentum, 0, bField, debug);
0676 DetId eId(0);
0677 if (info.second) {
0678 const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
0679 if (std::abs(point.eta()) < spr::etaBEEcal) {
0680 const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0681 eId = barrelGeom->getClosestCell(point);
0682 } else {
0683 const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0684 if (endcapGeom)
0685 eId = endcapGeom->getClosestCell(point);
0686 else
0687 info.second = false;
0688 }
0689 }
0690 return std::pair<DetId, bool>(eId, info.second);
0691 }
0692
0693 std::pair<math::XYZPoint, bool> propagateECAL(
0694 const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
0695 spr::propagatedTrack track =
0696 spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
0697 return std::pair<math::XYZPoint, bool>(track.point, track.ok);
0698 }
0699
0700 spr::propagatedTrack propagateTrackToHCAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
0701 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0702 const GlobalVector momentum(track->px(), track->py(), track->pz());
0703 int charge(track->charge());
0704 return spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
0705 }
0706
0707 spr::propagatedTrack propagateTrackToHCAL(unsigned int thisTrk,
0708 edm::Handle<edm::SimTrackContainer>& SimTk,
0709 edm::Handle<edm::SimVertexContainer>& SimVtx,
0710 const MagneticField* bfield,
0711 bool debug) {
0712 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
0713 spr::propagatedTrack ptrk;
0714 if (trk.ok)
0715 ptrk = spr::propagateCalo(
0716 trk.position, trk.momentum, trk.charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
0717 return ptrk;
0718 }
0719
0720 std::pair<math::XYZPoint, bool> propagateHCAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
0721 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0722 const GlobalVector momentum(track->px(), track->py(), track->pz());
0723 int charge(track->charge());
0724 return spr::propagateHCAL(vertex, momentum, charge, bfield, debug);
0725 }
0726
0727 std::pair<math::XYZPoint, bool> propagateHCAL(
0728 const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
0729 spr::propagatedTrack track =
0730 spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
0731 return std::pair<math::XYZPoint, bool>(track.point, track.ok);
0732 }
0733
0734 std::pair<math::XYZPoint, bool> propagateTracker(const reco::Track* track, const MagneticField* bfield, bool debug) {
0735 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0736 const GlobalVector momentum(track->px(), track->py(), track->pz());
0737 int charge(track->charge());
0738 spr::propagatedTrack track1 =
0739 spr::propagateCalo(vertex, momentum, charge, bfield, spr::zBackTE, spr::rBackTB, spr::etaBETrak, debug);
0740 return std::pair<math::XYZPoint, bool>(track1.point, track1.ok);
0741 }
0742
0743 std::pair<math::XYZPoint, double> propagateTrackerEnd(const reco::Track* track,
0744 const MagneticField* bField,
0745 bool debug) {
0746 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0747 const GlobalVector momentum(track->px(), track->py(), track->pz());
0748 int charge(track->charge());
0749 float radius = track->outerPosition().Rho();
0750 float zdist = track->outerPosition().Z();
0751 if (debug)
0752 edm::LogVerbatim("IsoTrack") << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum
0753 << " Charge " << charge << " Radius " << radius << " Z " << zdist;
0754 FreeTrajectoryState fts(vertex, momentum, charge, bField);
0755 Plane::PlanePointer endcap = Plane::build(Plane::PositionType(0, 0, zdist), Plane::RotationType());
0756 Cylinder::CylinderPointer barrel =
0757 Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), radius);
0758
0759 AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
0760
0761 TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
0762 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
0763
0764 math::XYZPoint point(-999., -999., -999.);
0765 bool ok = false;
0766 GlobalVector direction(0, 0, 1);
0767 if (tsosb.isValid() && std::abs(zdist) < spr::zFrontTE) {
0768 point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
0769 direction = tsosb.globalDirection();
0770 ok = true;
0771 } else if (tsose.isValid()) {
0772 point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
0773 direction = tsose.globalDirection();
0774 ok = true;
0775 }
0776
0777 double length = -1;
0778 if (ok) {
0779 math::XYZPoint vDiff(point.x() - vertex.x(), point.y() - vertex.y(), point.z() - vertex.z());
0780 double dphi = direction.phi() - momentum.phi();
0781 double rdist = std::sqrt(vDiff.x() * vDiff.x() + vDiff.y() * vDiff.y());
0782 double rat = 0.5 * dphi / std::sin(0.5 * dphi);
0783 double dZ = vDiff.z();
0784 double dS = rdist * rat;
0785 length = std::sqrt(dS * dS + dZ * dZ);
0786 if (debug)
0787 edm::LogVerbatim("IsoTrack") << "propagateTracker:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid()
0788 << " OverAll " << ok << " Point " << point << " RDist " << rdist << " dS " << dS
0789 << " dS/pt " << rdist * rat / momentum.perp() << " zdist " << dZ << " dz/pz "
0790 << dZ / momentum.z() << " Length " << length;
0791 }
0792
0793 return std::pair<math::XYZPoint, double>(point, length);
0794 }
0795
0796 spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex,
0797 const GlobalVector& tpMomentum,
0798 int tpCharge,
0799 const MagneticField* bField,
0800 float zdist,
0801 float radius,
0802 float corner,
0803 bool debug) {
0804 spr::propagatedTrack track;
0805 if (debug)
0806 edm::LogVerbatim("IsoTrack") << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge "
0807 << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner;
0808 FreeTrajectoryState fts(tpVertex, tpMomentum, tpCharge, bField);
0809
0810 Plane::PlanePointer lendcap = Plane::build(Plane::PositionType(0, 0, -zdist), Plane::RotationType());
0811 Plane::PlanePointer rendcap = Plane::build(Plane::PositionType(0, 0, zdist), Plane::RotationType());
0812
0813 Cylinder::CylinderPointer barrel =
0814 Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), radius);
0815
0816 AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
0817
0818 TrajectoryStateOnSurface tsose;
0819 if (tpMomentum.eta() < 0) {
0820 tsose = myAP.propagate(fts, *lendcap);
0821 } else {
0822 tsose = myAP.propagate(fts, *rendcap);
0823 }
0824
0825 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
0826
0827 track.ok = true;
0828 if (tsose.isValid() && tsosb.isValid()) {
0829 float absEta = std::abs(tsosb.globalPosition().eta());
0830 if (absEta < corner) {
0831 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
0832 track.direction = tsosb.globalDirection();
0833 } else {
0834 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
0835 track.direction = tsose.globalDirection();
0836 }
0837 } else if (tsose.isValid()) {
0838 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
0839 track.direction = tsose.globalDirection();
0840 } else if (tsosb.isValid()) {
0841 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
0842 track.direction = tsosb.globalDirection();
0843 } else {
0844 track.point.SetXYZ(-999., -999., -999.);
0845 track.direction = GlobalVector(0, 0, 1);
0846 track.ok = false;
0847 }
0848 if (debug) {
0849 edm::LogVerbatim("IsoTrack") << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid()
0850 << " OverAll " << track.ok << " Point " << track.point << " Direction "
0851 << track.direction;
0852 if (track.ok) {
0853 math::XYZPoint vDiff(
0854 track.point.x() - tpVertex.x(), track.point.y() - tpVertex.y(), track.point.z() - tpVertex.z());
0855 double dphi = track.direction.phi() - tpMomentum.phi();
0856 double rdist = std::sqrt(vDiff.x() * vDiff.x() + vDiff.y() * vDiff.y());
0857 double pt = tpMomentum.perp();
0858 double rat = 0.5 * dphi / std::sin(0.5 * dphi);
0859 edm::LogVerbatim("IsoTrack") << "RDist " << rdist << " pt " << pt << " r/pt " << rdist * rat / pt << " zdist "
0860 << vDiff.z() << " pz " << tpMomentum.z() << " z/pz " << vDiff.z() / tpMomentum.z()
0861 << std::endl;
0862 }
0863 }
0864 return track;
0865 }
0866
0867 spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk,
0868 edm::Handle<edm::SimTrackContainer>& SimTk,
0869 edm::Handle<edm::SimVertexContainer>& SimVtx,
0870 bool debug) {
0871 spr::trackAtOrigin trk;
0872
0873 edm::SimTrackContainer::const_iterator itr = SimTk->end();
0874 for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0875 if (simTrkItr->trackId() == thisTrk) {
0876 if (debug)
0877 edm::LogVerbatim("IsoTrack") << "matched trackId (maximum occurance) " << thisTrk << " type "
0878 << simTrkItr->type();
0879 itr = simTrkItr;
0880 break;
0881 }
0882 }
0883
0884 if (itr != SimTk->end()) {
0885 int vertIndex = itr->vertIndex();
0886 if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
0887 edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
0888 for (int iv = 0; iv < vertIndex; iv++)
0889 simVtxItr++;
0890 const math::XYZTLorentzVectorD pos = simVtxItr->position();
0891 const math::XYZTLorentzVectorD mom = itr->momentum();
0892 trk.ok = true;
0893 trk.charge = (int)(itr->charge());
0894 trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
0895 trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
0896 }
0897 }
0898 if (debug)
0899 edm::LogVerbatim("IsoTrack") << "Track flag " << trk.ok << " Position " << trk.position << " Momentum "
0900 << trk.momentum;
0901 return trk;
0902 }
0903
0904 bool propagateHCAL(const reco::Track* track,
0905 const CaloGeometry* geo,
0906 const MagneticField* bField,
0907 bool typeRZ,
0908 const std::pair<double, double> rz,
0909 bool debug) {
0910 const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
0911 const GlobalVector momentum(track->px(), track->py(), track->pz());
0912 int charge(track->charge());
0913 if (debug)
0914 edm::LogVerbatim("IsoTrack") << "Propagate track with charge " << charge << " position " << vertex << " p "
0915 << momentum;
0916 std::pair<HcalDetId, HcalDetId> ids = propagateHCAL(geo, bField, vertex, momentum, charge, typeRZ, rz, debug);
0917 bool ok = ((ids.first != HcalDetId()) && (ids.first.ieta() == ids.second.ieta()) &&
0918 (ids.first.iphi() == ids.second.iphi()));
0919 return ok;
0920 }
0921
0922 bool propagateHCAL(unsigned int thisTrk,
0923 edm::Handle<edm::SimTrackContainer>& SimTk,
0924 edm::Handle<edm::SimVertexContainer>& SimVtx,
0925 const CaloGeometry* geo,
0926 const MagneticField* bField,
0927 bool typeRZ,
0928 const std::pair<double, double> rz,
0929 bool debug) {
0930 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
0931 if (debug)
0932 edm::LogVerbatim("IsoTrack") << "Propagate track " << thisTrk << " charge " << trk.charge << " position "
0933 << trk.position << " p " << trk.momentum;
0934 std::pair<HcalDetId, HcalDetId> ids =
0935 propagateHCAL(geo, bField, trk.position, trk.momentum, trk.charge, typeRZ, rz, debug);
0936 bool ok = ((ids.first != HcalDetId()) && (ids.first.ieta() == ids.second.ieta()) &&
0937 (ids.first.iphi() == ids.second.iphi()));
0938 return ok;
0939 }
0940
0941 std::pair<HcalDetId, HcalDetId> propagateHCAL(const CaloGeometry* geo,
0942 const MagneticField* bField,
0943 const GlobalPoint& vertex,
0944 const GlobalVector& momentum,
0945 int charge,
0946 bool typeRZ,
0947 const std::pair<double, double> rz,
0948 bool debug) {
0949 if (debug)
0950 edm::LogVerbatim("IsoTrack") << "propagateCalo:: Vertex " << vertex << " Momentum " << momentum << " Charge "
0951 << charge << " R/Z " << rz.first << " : " << rz.second << " Type " << typeRZ;
0952 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0953 FreeTrajectoryState fts(vertex, momentum, charge, bField);
0954 AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
0955
0956 HcalDetId id1, id2;
0957 for (int k = 0; k < 2; ++k) {
0958 TrajectoryStateOnSurface tsos;
0959 double rzv = (k == 0) ? rz.first : rz.second;
0960 if (typeRZ) {
0961 Cylinder::CylinderPointer barrel =
0962 Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), rzv);
0963 tsos = myAP.propagate(fts, *barrel);
0964 } else {
0965 Plane::PlanePointer endcap = Plane::build(Plane::PositionType(0, 0, rzv), Plane::RotationType());
0966 tsos = myAP.propagate(fts, *endcap);
0967 }
0968
0969 if (tsos.isValid()) {
0970 GlobalPoint point = tsos.globalPosition();
0971 if (k == 0)
0972 id1 = gHB->getClosestCell(point);
0973 else
0974 id2 = gHB->getClosestCell(point);
0975 if (debug) {
0976 std::ostringstream st1;
0977 if (k == 0)
0978 st1 << id1;
0979 else
0980 st1 << id2;
0981 edm::LogVerbatim("IsoTrack") << "Iteration " << k << " Point " << point << " ID " << st1.str();
0982 }
0983 }
0984 }
0985 if (debug)
0986 edm::LogVerbatim("IsoTrack") << "propagateCalo:: Front " << id1 << " Back " << id2;
0987 return std::pair<HcalDetId, HcalDetId>(id1, id2);
0988 }
0989 }
