CaloPropagateTrack.cc

Go to the documentation of this file.

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "DataFormats/GeometrySurface/interface/Cylinder.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"

#include "Calibration/IsolatedParticles/interface/CaloConstants.h"
#include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"

#include <iostream>

//#define EDM_ML_DEBUG

namespace spr{

  std::vector<spr::propagatedTrackID> propagateCosmicCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, const std::string & theTrackQuality, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB =        geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
    reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
    std::vector<spr::propagatedTrackID> vdets;

    unsigned int indx;
    reco::TrackCollection::const_iterator trkItr;
    for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
      const reco::Track* pTrack = &(*trkItr);
      spr::propagatedTrackID vdet;
      vdet.trkItr = trkItr;
      vdet.ok     = (trackQuality_ != reco::TrackBase::undefQuality) ?
    (pTrack->quality(trackQuality_)) : true;
      vdet.detIdECAL = DetId(0);
      vdet.detIdHCAL = DetId(0);
      vdet.detIdEHCAL= DetId(0);
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << vdet.ok << std::endl;
#endif
      GlobalPoint  vertex;
      GlobalVector momentum;
      int charge (pTrack->charge());
      if (((pTrack->innerPosition()).Perp2()) < 
      ((pTrack->outerPosition()).Perp2())) {
    vertex = GlobalPoint(((pTrack->innerPosition()).X()),
                 ((pTrack->innerPosition()).Y()),
                 ((pTrack->innerPosition()).Z()));
    momentum = GlobalVector(((pTrack->innerMomentum()).X()),
                ((pTrack->innerMomentum()).Y()),
                ((pTrack->innerMomentum()).Z()));
      } else {
    vertex = GlobalPoint(((pTrack->outerPosition()).X()),
                 ((pTrack->outerPosition()).Y()),
                 ((pTrack->outerPosition()).Z()));
    momentum = GlobalVector(((pTrack->outerMomentum()).X()),
                ((pTrack->outerMomentum()).Y()),
                ((pTrack->outerMomentum()).Z()));
      }
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Track charge " << charge << " p " << momentum << " position " << vertex << std::endl;
#endif
      std::pair<math::XYZPoint,bool> info = 
    spr::propagateECAL (vertex, momentum, charge, bField, debug);
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate to ECAL " << info.second << " at (" << info.first.x() << ", "<< info.first.y() << ", " << info.first.z() << ")\n";
#endif

      vdet.okECAL = info.second;
      if (vdet.okECAL) {
    const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
    vdet.etaECAL = point.eta();
    vdet.phiECAL = point.phi();
    if (std::abs(point.eta())<spr::etaBEEcal) {
      vdet.detIdECAL = barrelGeom->getClosestCell(point);
    } else {
      if (endcapGeom) 
        vdet.detIdECAL = endcapGeom->getClosestCell(point);
      else
        vdet.okECAL    = false;
    }
    vdet.detIdEHCAL = gHB->getClosestCell(point);
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "Point at ECAL (" << vdet.etaECAL << ", " << vdet.phiECAL << " ";
      if (std::abs(point.eta())<spr::etaBEEcal)
        std::cout << EBDetId(vdet.detIdECAL);
      else
        std::cout << EEDetId(vdet.detIdECAL);
      std::cout << " " << HcalDetId(vdet.detIdEHCAL) << std::endl;
    }
#endif
      }
      info = spr::propagateHCAL (vertex, momentum, charge, bField, debug);
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate to HCAL " << info.second << " at (" << info.first.x() << ", "<< info.first.y() << ", " << info.first.z() << ")\n";
#endif
      vdet.okHCAL = info.second;
      if (vdet.okHCAL) {
    const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
    vdet.etaHCAL = point.eta();
    vdet.phiHCAL = point.phi();
    vdet.detIdHCAL = gHB->getClosestCell(point);
      }
#ifdef EDM_ML_DEBUG
      if (debug) {
    std::cout << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL (" 
          << vdet.okECAL << ") ";
    if (vdet.detIdECAL.subdetId() == EcalBarrel) 
      std::cout << (EBDetId)(vdet.detIdECAL);
    else 
      std::cout << (EEDetId)(vdet.detIdECAL); 
    std::cout << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or " << (HcalDetId)(vdet.detIdEHCAL) << std::endl;
      }
#endif
      vdets.push_back(vdet);
    }
    
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for " << vdets.size() << " tracks" << std::endl;
      for (unsigned int i=0; i<vdets.size(); ++i) {
    std::cout << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL << ") ";
    if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
      std::cout << (EBDetId)(vdets[i].detIdECAL);
    } else {
      std::cout << (EEDetId)(vdets[i].detIdECAL); 
    }
    std::cout << " HCAL (" << vdets[i].okHCAL << ") " << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL) << std::endl;
      }
    }
#endif
    return vdets;
  }
 
  std::vector<spr::propagatedTrackID> propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, const std::string & theTrackQuality, bool debug) {

    std::vector<spr::propagatedTrackID> vdets;
    spr::propagateCALO(trkCollection,geo,bField,theTrackQuality, vdets, debug);
    return vdets;
  }

  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, const std::string & theTrackQuality, std::vector<spr::propagatedTrackID>& vdets, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB        = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
    reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);

    unsigned int indx;
    reco::TrackCollection::const_iterator trkItr;
    for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
      const reco::Track* pTrack = &(*trkItr);
      spr::propagatedTrackID vdet;
      vdet.trkItr = trkItr;
      vdet.ok     = (trackQuality_ != reco::TrackBase::undefQuality) ?
    (pTrack->quality(trackQuality_)) : true;
      vdet.detIdECAL = DetId(0);
      vdet.detIdHCAL = DetId(0);
      vdet.detIdEHCAL= DetId(0);
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << vdet.ok << std::endl;
#endif
      std::pair<math::XYZPoint,bool> info = spr::propagateECAL (pTrack, bField, debug);
      vdet.okECAL = info.second;
      if (vdet.okECAL) {
    const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
    vdet.etaECAL = point.eta();
    vdet.phiECAL = point.phi();
    if (std::abs(point.eta())<spr::etaBEEcal) {
      vdet.detIdECAL = barrelGeom->getClosestCell(point);
    } else {
      if (endcapGeom) 
        vdet.detIdECAL = endcapGeom->getClosestCell(point);
      else 
        vdet.okECAL    = false;
    }
    vdet.detIdEHCAL = gHB->getClosestCell(point);
      }
      info = spr::propagateHCAL (pTrack, bField, debug);
      vdet.okHCAL = info.second;
      if (vdet.okHCAL) {
    const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
    vdet.etaHCAL = point.eta();
    vdet.phiHCAL = point.phi();
    vdet.detIdHCAL = gHB->getClosestCell(point);
      }
#ifdef EDM_ML_DEBUG
      if (debug) {
    std::cout << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL (" 
          << vdet.okECAL << ") ";
    if (vdet.detIdECAL.subdetId() == EcalBarrel) 
      std::cout << (EBDetId)(vdet.detIdECAL);
    else 
      std::cout << (EEDetId)(vdet.detIdECAL); 
    std::cout << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or " << (HcalDetId)(vdet.detIdEHCAL) << std::endl;
      }
#endif
      vdets.push_back(vdet);
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for " << vdets.size() << " tracks" << std::endl;
      for (unsigned int i=0; i<vdets.size(); ++i) {
    std::cout << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL << ") ";
    if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
      std::cout << (EBDetId)(vdets[i].detIdECAL);
    } else {
      std::cout << (EEDetId)(vdets[i].detIdECAL); 
    }
    std::cout << " HCAL (" << vdets[i].okHCAL << ") " << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL) << std::endl;
      }
    }
#endif
  }

  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, const std::string & theTrackQuality, std::vector<spr::propagatedTrackDirection>& trkDir, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB =        geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
    reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);

    unsigned int indx;
    reco::TrackCollection::const_iterator trkItr;
    for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
      const reco::Track* pTrack = &(*trkItr);
      spr::propagatedTrackDirection trkD;
      trkD.trkItr = trkItr;
      trkD.ok     = (trackQuality_ != reco::TrackBase::undefQuality) ?
    (pTrack->quality(trackQuality_)) : true;
      trkD.detIdECAL = DetId(0);
      trkD.detIdHCAL = DetId(0);
      trkD.detIdEHCAL= DetId(0);
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << trkD.ok << std::endl;
#endif
      spr::propagatedTrack info = spr::propagateTrackToECAL (pTrack, bField, debug);
      GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
      trkD.okECAL        = info.ok;
      trkD.pointECAL     = point;
      trkD.directionECAL = info.direction;
      if (trkD.okECAL) {
    if (std::abs(info.point.eta())<spr::etaBEEcal) {
      trkD.detIdECAL = barrelGeom->getClosestCell(point);
    } else {
      if (endcapGeom) 
        trkD.detIdECAL = endcapGeom->getClosestCell(point);
      else 
        trkD.okECAL    = false;
    }
    trkD.detIdEHCAL = gHB->getClosestCell(point);
      }
      info = spr::propagateTrackToHCAL (pTrack, bField, debug);
      point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
      trkD.okHCAL        = info.ok;
      trkD.pointHCAL     = point;
      trkD.directionHCAL = info.direction;
      if (trkD.okHCAL) {
    trkD.detIdHCAL = gHB->getClosestCell(point);
      }
      trkDir.push_back(trkD);
    }    
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
      for (unsigned int i=0; i<trkDir.size(); ++i) {
    std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
    if (trkDir[i].okECAL) {
      std::cout << " point " << trkDir[i].pointECAL << " direction "
            << trkDir[i].directionECAL << " "; 
      if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
        std::cout << (EBDetId)(trkDir[i].detIdECAL);
      } else {
        std::cout << (EEDetId)(trkDir[i].detIdECAL); 
      }
    }
    std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
    if (trkDir[i].okHCAL) {
      std::cout << " point " << trkDir[i].pointHCAL << " direction "
            << trkDir[i].directionHCAL << " " 
            << (HcalDetId)(trkDir[i].detIdHCAL); 
    }
    std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
      }
    }
#endif
  }

  spr::propagatedTrackID propagateCALO(const reco::Track* pTrack, const CaloGeometry* geo, const MagneticField* bField, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB =        geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);

    spr::propagatedTrackID vdet;
    vdet.ok        = true;
    vdet.detIdECAL = DetId(0);
    vdet.detIdHCAL = DetId(0);
    vdet.detIdEHCAL= DetId(0);
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "Propagate track:  p " << pTrack->p() << " eta " << pTrack->eta() << " phi " << pTrack->phi() << " Flag " << vdet.ok << std::endl;
#endif
    std::pair<math::XYZPoint,bool> info = spr::propagateECAL (pTrack, bField, debug);
    vdet.okECAL = info.second;
    if (vdet.okECAL) {
      const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
      vdet.etaECAL = point.eta();
      vdet.phiECAL = point.phi();
      if (std::abs(point.eta())<spr::etaBEEcal) {
    vdet.detIdECAL = barrelGeom->getClosestCell(point);
      } else {
    if (endcapGeom) 
      vdet.detIdECAL = endcapGeom->getClosestCell(point);
    else
      vdet.okECAL    = false;
      }
      vdet.detIdEHCAL = gHB->getClosestCell(point);
    }
    info = spr::propagateHCAL (pTrack, bField, debug);
    vdet.okHCAL = info.second;
    if (vdet.okHCAL) {
      const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
      vdet.etaHCAL = point.eta();
      vdet.phiHCAL = point.phi();
      vdet.detIdHCAL = gHB->getClosestCell(point);
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for 1 track" << std::endl;
      std::cout << "Track [0] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") ";
      if (vdet.detIdECAL.subdetId() == EcalBarrel) {
    std::cout << (EBDetId)(vdet.detIdECAL);
      } else {
    std::cout << (EEDetId)(vdet.detIdECAL); 
      }
      std::cout << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or " << (HcalDetId)(vdet.detIdEHCAL) << std::endl;
    }
#endif
    return vdet;
  }

  std::vector<spr::propagatedGenTrackID> propagateCALO(const HepMC::GenEvent * genEvent, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB =        geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);

    std::vector<spr::propagatedGenTrackID> trkDir;
    unsigned int indx;
    HepMC::GenEvent::particle_const_iterator p;
    for (p=genEvent->particles_begin(),indx=0;   p != genEvent->particles_end(); ++p,++indx) {
      spr::propagatedGenTrackID trkD;
      trkD.trkItr    = p;
      trkD.detIdECAL = DetId(0);
      trkD.detIdHCAL = DetId(0);
      trkD.detIdEHCAL= DetId(0);
      trkD.pdgId  = ((*p)->pdg_id());
      trkD.charge = ((pdt->particle(trkD.pdgId))->ID().threeCharge())/3;
      const GlobalVector momentum = GlobalVector((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
#endif      
      // consider stable particles
      if ( (*p)->status()==1 && std::abs((*p)->momentum().eta()) < etaMax ) { 
    const GlobalPoint vertex = GlobalPoint(0.1*(*p)->production_vertex()->position().x(), 
                           0.1*(*p)->production_vertex()->position().y(), 
                           0.1*(*p)->production_vertex()->position().z());
    trkD.ok = true;
    spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
    GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
    trkD.okECAL        = info.ok;
    trkD.pointECAL     = point;
    trkD.directionECAL = info.direction;
    if (trkD.okECAL) {
      if (std::abs(info.point.eta())<spr::etaBEEcal) {
        trkD.detIdECAL = barrelGeom->getClosestCell(point);
      } else {
        if (endcapGeom) 
          trkD.detIdECAL = endcapGeom->getClosestCell(point);
        else 
          trkD.okECAL    = false;
      }
      trkD.detIdEHCAL = gHB->getClosestCell(point);
    }

    info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
    point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
    trkD.okHCAL        = info.ok;
    trkD.pointHCAL     = point;
    trkD.directionHCAL = info.direction;
    if (trkD.okHCAL) {
      trkD.detIdHCAL = gHB->getClosestCell(point);
    }
      }
      trkDir.push_back(trkD);
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
      for (unsigned int i=0; i<trkDir.size(); ++i) {
    if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
    if (trkDir[i].okECAL) {
      std::cout << " point " << trkDir[i].pointECAL << " direction "
            << trkDir[i].directionECAL << " "; 
      if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
        std::cout << (EBDetId)(trkDir[i].detIdECAL);
      } else {
        std::cout << (EEDetId)(trkDir[i].detIdECAL); 
      }
    }
    if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
    if (trkDir[i].okHCAL) {
      std::cout << " point " << trkDir[i].pointHCAL << " direction "
            << trkDir[i].directionHCAL << " " 
            << (HcalDetId)(trkDir[i].detIdHCAL); 
    }
    if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
      }
    }
#endif
    return trkDir;
  }

  std::vector<spr::propagatedGenParticleID> propagateCALO(edm::Handle<reco::GenParticleCollection>& genParticles, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB =        geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);

    std::vector<spr::propagatedGenParticleID> trkDir;
    unsigned int indx;
    reco::GenParticleCollection::const_iterator p;
    for (p=genParticles->begin(),indx=0;   p != genParticles->end(); ++p,++indx) {
      spr::propagatedGenParticleID trkD;
      trkD.trkItr    = p;
      trkD.detIdECAL = DetId(0);
      trkD.detIdHCAL = DetId(0);
      trkD.detIdEHCAL= DetId(0);
      trkD.pdgId     = (p->pdgId());
      trkD.charge    = p->charge();
      const GlobalVector momentum = GlobalVector(p->momentum().x(), p->momentum().y(), p->momentum().z());
#ifdef EDM_ML_DEBUG
      if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
#endif      
      // consider stable particles
      if ( p->status()==1 && std::abs(momentum.eta()) < etaMax ) { 
    const GlobalPoint vertex = GlobalPoint(p->vertex().x(), p->vertex().y(), p->vertex().z());
    trkD.ok = true;
    spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
    GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
    trkD.okECAL        = info.ok;
    trkD.pointECAL     = point;
    trkD.directionECAL = info.direction;
    if (trkD.okECAL) {
      if (std::abs(info.point.eta())<spr::etaBEEcal) {
        trkD.detIdECAL = barrelGeom->getClosestCell(point);
      } else {
        if (endcapGeom) 
          trkD.detIdECAL = endcapGeom->getClosestCell(point);
        else 
          trkD.okECAL    = false;
      }
      trkD.detIdEHCAL = gHB->getClosestCell(point);
    }

    info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
    point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
    trkD.okHCAL        = info.ok;
    trkD.pointHCAL     = point;
    trkD.directionHCAL = info.direction;
    if (trkD.okHCAL) {
      trkD.detIdHCAL = gHB->getClosestCell(point);
    }
      }
      trkDir.push_back(trkD);
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
      for (unsigned int i=0; i<trkDir.size(); ++i) {
    if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
    if (trkDir[i].okECAL) {
      std::cout << " point " << trkDir[i].pointECAL << " direction "
            << trkDir[i].directionECAL << " "; 
      if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
        std::cout << (EBDetId)(trkDir[i].detIdECAL);
      } else {
        std::cout << (EEDetId)(trkDir[i].detIdECAL); 
      }
    }
    if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
    if (trkDir[i].okHCAL) {
      std::cout << " point " << trkDir[i].pointHCAL << " direction "
            << trkDir[i].directionHCAL << " " 
            << (HcalDetId)(trkDir[i].detIdHCAL); 
    }
    if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
      }
    }
#endif
    return trkDir;
  }

  spr::propagatedTrackDirection propagateCALO(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const CaloGeometry* geo, const MagneticField* bField, bool debug) {

    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    const CaloSubdetectorGeometry* gHB =        geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);

    spr::trackAtOrigin   trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
    spr::propagatedTrackDirection trkD;
    trkD.ok     = trk.ok;
    trkD.detIdECAL = DetId(0);
    trkD.detIdHCAL = DetId(0);
    trkD.detIdEHCAL= DetId(0);
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "Propagate track " << thisTrk << " charge " << trk.charge << " position " << trk.position << " p " << trk.momentum << " Flag " << trkD.ok << std::endl;
#endif
    if (trkD.ok) {
      spr::propagatedTrack info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
      GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
      trkD.okECAL        = info.ok;
      trkD.pointECAL     = point;
      trkD.directionECAL = info.direction;
      if (trkD.okECAL) {
    if (std::abs(info.point.eta())<spr::etaBEEcal) {
      trkD.detIdECAL = barrelGeom->getClosestCell(point);
    } else {
      if (endcapGeom) 
        trkD.detIdECAL = endcapGeom->getClosestCell(point);
      else 
        trkD.okECAL    = false;
    }
    trkD.detIdEHCAL = gHB->getClosestCell(point);
      }

      info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
      point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
      trkD.okHCAL        = info.ok;
      trkD.pointHCAL     = point;
      trkD.directionHCAL = info.direction;
      if (trkD.okHCAL) {
    trkD.detIdHCAL = gHB->getClosestCell(point);
      }
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL (" << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << std::endl;
      if (trkD.okECAL) {
    std::cout << "ECAL point " << trkD.pointECAL << " direction "
          << trkD.directionECAL << " "; 
    if (trkD.detIdECAL.subdetId() == EcalBarrel) {
      std::cout << (EBDetId)(trkD.detIdECAL);
    } else {
      std::cout << (EEDetId)(trkD.detIdECAL); 
    }
      }
      if (trkD.okHCAL) {
    std::cout << " HCAL point " << trkD.pointHCAL << " direction "
          << trkD.directionHCAL << " " << (HcalDetId)(trkD.detIdHCAL); 
      }
      if (trkD.okECAL) std::cout << " Or " << (HcalDetId)(trkD.detIdEHCAL);
      std::cout << std::endl;
    }
#endif
    return trkD;
  }

  spr::propagatedTrackDirection propagateHCALBack(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const CaloGeometry* geo, const MagneticField* bField, bool debug) {

    const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
    spr::trackAtOrigin   trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
    spr::propagatedTrackDirection trkD;
    trkD.ok     = trk.ok;
    trkD.detIdECAL = DetId(0);
    trkD.detIdHCAL = DetId(0);
    trkD.detIdEHCAL= DetId(0);
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "Propagate track " << thisTrk << " charge " << trk.charge << " position " << trk.position << " p " << trk.momentum << " Flag " << trkD.ok << std::endl;
#endif
    if (trkD.ok) {
      spr::propagatedTrack info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, spr::zBackHE, spr::rBackHB, spr::etaBEHcal, debug);
      const GlobalPoint point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
      trkD.okHCAL        = info.ok;
      trkD.pointHCAL     = point;
      trkD.directionHCAL = info.direction;
      if (trkD.okHCAL) {
    trkD.detIdHCAL = gHB->getClosestCell(point);
      }
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL (" << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << std::endl;
            if (trkD.okHCAL) {
    std::cout << " HCAL point " << trkD.pointHCAL << " direction "
          << trkD.directionHCAL << " " << (HcalDetId)(trkD.detIdHCAL); 
      }
    }
#endif
    return trkD;
  }


  std::pair<bool,HcalDetId> propagateHCALBack(const reco::Track* track, const CaloGeometry* geo, const MagneticField* bField, bool debug) {
    const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    spr::propagatedTrack info = spr::propagateCalo(vertex, momentum, charge, bField, spr::zBackHE, spr::rBackHB, spr::etaBEHcal, debug);
    if (info.ok) {
      const GlobalPoint point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
      return std::pair<bool,HcalDetId>(true,HcalDetId(gHB->getClosestCell(point)));
    } else {
      return std::pair<bool,HcalDetId>(false,HcalDetId());
    }
  }

  propagatedTrack propagateTrackToECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    return spr::propagateCalo (vertex, momentum, charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
  }

  propagatedTrack propagateTrackToECAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {

    spr::trackAtOrigin   trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
    spr::propagatedTrack ptrk;
    if (trk.ok) 
      ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
    return ptrk;
  }

  std::pair<math::XYZPoint,bool> propagateECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {    
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    return spr::propagateECAL (vertex, momentum, charge, bfield, debug);
  }

  std::pair<DetId,bool> propagateIdECAL(const HcalDetId& id, const CaloGeometry* geo, const MagneticField* bField, bool debug) {

    const HcalGeometry* gHB = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel));
    const GlobalPoint vertex(0,0,0);
    const GlobalPoint hit(gHB->getPosition(id));
    const GlobalVector momentum = GlobalVector(hit.x(),hit.y(),hit.z());
    std::pair<math::XYZPoint,bool> info = propagateECAL(vertex,momentum,0,bField,debug);
    DetId eId(0);
    if (info.second) {
      const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
      if (std::abs(point.eta())<spr::etaBEEcal) {
    const CaloSubdetectorGeometry *barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    eId = barrelGeom->getClosestCell(point);
      } else {
    const CaloSubdetectorGeometry *endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
    if (endcapGeom) 
       eId = endcapGeom->getClosestCell(point);
      else 
        info.second = false;
      }
    }
    return std::pair<DetId,bool>(eId,info.second);
  }

  std::pair<math::XYZPoint,bool> propagateECAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
    spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
    return std::pair<math::XYZPoint,bool>(track.point,track.ok);
  }

  spr::propagatedTrack propagateTrackToHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    return spr::propagateCalo (vertex, momentum, charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
  }

  spr::propagatedTrack propagateTrackToHCAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
    spr::trackAtOrigin   trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
    spr::propagatedTrack ptrk;
    if (trk.ok) 
      ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
    return ptrk;
  }

  std::pair<math::XYZPoint,bool> propagateHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    return spr::propagateHCAL (vertex, momentum, charge, bfield, debug);
  }

  std::pair<math::XYZPoint,bool> propagateHCAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
    spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
    return std::pair<math::XYZPoint,bool>(track.point,track.ok);
  }

  std::pair<math::XYZPoint,bool> propagateTracker(const reco::Track *track, const MagneticField* bfield, bool debug) {
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    spr::propagatedTrack track1 = spr::propagateCalo (vertex, momentum, charge, bfield, spr::zBackTE, spr::rBackTB, spr::etaBETrak, debug);
    return std::pair<math::XYZPoint,bool>(track1.point,track1.ok);
  }

  std::pair<math::XYZPoint,double> propagateTrackerEnd(const reco::Track *track,
                               const MagneticField* bField, bool
#ifdef EDM_ML_DEBUG
                               debug
#endif
                               ) {

    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
    float radius = track->outerPosition().Rho();
    float zdist  = track->outerPosition().Z();
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " Radius " << radius << " Z " << zdist << std::endl;
#endif
    FreeTrajectoryState fts (vertex, momentum, charge, bField);
    Plane::PlanePointer endcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
    Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);

    AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);

    TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
    TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);

    math::XYZPoint point(-999.,-999.,-999.);
    bool ok=false;
    GlobalVector direction(0,0,1);
    if (tsosb.isValid() && std::abs(zdist) < spr::zFrontTE) {
      point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
      direction = tsosb.globalDirection();
      ok = true;
    } else if (tsose.isValid()) {
      point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
      direction = tsose.globalDirection();
      ok = true;
    }

    double length = -1;
    if (ok) {
      math::XYZPoint vDiff(point.x()-vertex.x(), point.y()-vertex.y(), point.z()-vertex.z());
      double dphi  = direction.phi()-momentum.phi();
      double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
      double rat   = 0.5*dphi/std::sin(0.5*dphi);
      double dZ    = vDiff.z();
      double dS    = rdist*rat; //dZ*momentum.z()/momentum.perp();
      length       = std::sqrt(dS*dS+dZ*dZ);
#ifdef EDM_ML_DEBUG
      if (debug) 
    std::cout << "propagateTracker:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << ok << " Point " << point << " RDist " << rdist << " dS " << dS << " dS/pt " << rdist*rat/momentum.perp() << " zdist " << dZ << " dz/pz " << dZ/momentum.z() << " Length " << length << std::endl;
#endif
    }

    return std::pair<math::XYZPoint,double>(point,length);
  }

  spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex, 
                     const GlobalVector& tpMomentum, 
                     int tpCharge, const MagneticField* bField,
                     float zdist, float radius, float corner, bool
#ifdef EDM_ML_DEBUG
                     debug
#endif
                     ) {
    
    spr::propagatedTrack track;
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge " << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner << std::endl;
#endif
    FreeTrajectoryState fts (tpVertex, tpMomentum, tpCharge, bField);
    
    Plane::PlanePointer lendcap = Plane::build(Plane::PositionType (0, 0, -zdist), Plane::RotationType());
    Plane::PlanePointer rendcap = Plane::build(Plane::PositionType (0, 0,  zdist), Plane::RotationType());
    
    Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
  
    AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);

    TrajectoryStateOnSurface tsose;
    if (tpMomentum.eta() < 0) {
      tsose = myAP.propagate(fts, *lendcap);
    } else {
      tsose = myAP.propagate(fts, *rendcap);
    }

    TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);

    track.ok=true;
    if (tsose.isValid() && tsosb.isValid()) {
      float absEta = std::abs(tsosb.globalPosition().eta());
      if (absEta < corner) {
    track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
    track.direction = tsosb.globalDirection();
      } else {
    track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
    track.direction = tsose.globalDirection();
      }
    } else if (tsose.isValid()) {
      track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
      track.direction = tsose.globalDirection();
    } else if (tsosb.isValid()) {
      track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
      track.direction = tsosb.globalDirection();
    } else {
      track.point.SetXYZ(-999., -999., -999.);
      track.direction = GlobalVector(0,0,1);
      track.ok = false;
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << track.ok << " Point " << track.point << " Direction " << track.direction << std::endl;
      if (track.ok) {
    math::XYZPoint vDiff(track.point.x()-tpVertex.x(), track.point.y()-tpVertex.y(), track.point.z()-tpVertex.z());
    double dphi = track.direction.phi()-tpMomentum.phi();
    double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
    double pt    = tpMomentum.perp();
    double rat   = 0.5*dphi/std::sin(0.5*dphi);
    std::cout << "RDist " << rdist << " pt " << pt << " r/pt " << rdist*rat/pt << " zdist " << vDiff.z() << " pz " << tpMomentum.z() << " z/pz " << vDiff.z()/tpMomentum.z() << std::endl;
      }
    }
#endif
    return track;
  }

  spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, 
                      edm::Handle<edm::SimTrackContainer>& SimTk, 
                      edm::Handle<edm::SimVertexContainer>& SimVtx, bool 
#ifdef EDM_ML_DEBUG
                      debug
#endif
                      ) {

    spr::trackAtOrigin trk;

    edm::SimTrackContainer::const_iterator itr = SimTk->end();
    for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr!= SimTk->end(); simTrkItr++) {
      if ( simTrkItr->trackId() == thisTrk ) {
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "matched trackId (maximum occurance) " << thisTrk << " type " << simTrkItr->type() << std::endl;
#endif
    itr = simTrkItr;
    break;
      }
    }

    if (itr != SimTk->end()) {
      int vertIndex = itr->vertIndex();
      if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
    edm::SimVertexContainer::const_iterator simVtxItr= SimVtx->begin();
    for (int iv=0; iv<vertIndex; iv++) simVtxItr++;
    const math::XYZTLorentzVectorD pos = simVtxItr->position();
    const math::XYZTLorentzVectorD mom = itr->momentum();
    trk.ok = true;
    trk.charge   = (int)(itr->charge());
    trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
    trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
      }
    }
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "Track flag " << trk.ok << " Position " << trk.position << " Momentum " << trk.momentum << std::endl;
#endif
    return trk;
  }

  bool propagateHCAL(const reco::Track *track, const CaloGeometry* geo, const MagneticField* bField, bool typeRZ, const std::pair<double,double> rz, bool debug) {
    const GlobalPoint  vertex (track->vx(), track->vy(), track->vz());
    const GlobalVector momentum (track->px(), track->py(), track->pz());
    int charge (track->charge());
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "Propagate track with charge " << charge << " position " << vertex << " p " << momentum << std::endl;
#endif
    std::pair<HcalDetId,HcalDetId> ids = propagateHCAL(geo, bField, vertex, momentum, charge, typeRZ, rz, debug);
    bool ok = ((ids.first != HcalDetId()) && 
           (ids.first.ieta() == ids.second.ieta()) && 
           (ids.first.iphi() == ids.second.iphi()));
    return ok;
  }

  bool propagateHCAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const CaloGeometry* geo, const MagneticField* bField, bool typeRZ, const std::pair<double,double> rz, bool debug) {

    spr::trackAtOrigin   trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "Propagate track " << thisTrk << " charge " << trk.charge << " position " << trk.position << " p " << trk.momentum << std::endl;
#endif
    std::pair<HcalDetId,HcalDetId> ids = propagateHCAL(geo, bField, trk.position, trk.momentum, trk.charge, typeRZ, rz, debug);
    bool ok = ((ids.first != HcalDetId()) && 
           (ids.first.ieta() == ids.second.ieta()) && 
           (ids.first.iphi() == ids.second.iphi()));
    return ok;
  }
  
  std::pair<HcalDetId,HcalDetId> propagateHCAL(const CaloGeometry* geo,
                           const MagneticField* bField,
                           const GlobalPoint& vertex, 
                           const GlobalVector& momentum, 
                           int charge, bool typeRZ,
                           const std::pair<double,double> rz, bool
#ifdef EDM_ML_DEBUG
                           debug
#endif
                           ) {
    
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "propagateCalo:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " R/Z " << rz.first << " : " << rz.second << " Type " << typeRZ << std::endl;
#endif
    const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
    FreeTrajectoryState fts (vertex, momentum, charge, bField);
    AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
    
    HcalDetId id1, id2;
    for (int k=0; k<2; ++k) {
      TrajectoryStateOnSurface tsos;
      double rzv = (k == 0) ? rz.first : rz.second;
      if (typeRZ) {
    Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), rzv);
    tsos = myAP.propagate(fts, *barrel);
      } else {
    Plane::PlanePointer endcap = Plane::build(Plane::PositionType (0, 0, rzv), Plane::RotationType());
    tsos = myAP.propagate(fts, *endcap);
      }

      if (tsos.isValid()) {
    GlobalPoint point = tsos.globalPosition();
    if (k == 0) id1 = gHB->getClosestCell(point);
    else        id2 = gHB->getClosestCell(point);
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "Iteration " << k << " Point " << point << " ID ";
      if (k == 0) std::cout << id1;
      else        std::cout << id2;
      std::cout << std::endl;
    }
#endif
      } 
    }
#ifdef EDM_ML_DEBUG
    if (debug) std::cout << "propagateCalo:: Front " << id1 << " Back " << id2 << std::endl;
#endif
    return std::pair<HcalDetId,HcalDetId>(id1,id2);
  }
}