Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-09 22:24:50

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       // consider stable particles
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       // consider stable particles
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;  //dZ*momentum.z()/momentum.perp();
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 }  // namespace spr