00001 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00002 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00003 #include "DataFormats/GeometrySurface/interface/Plane.h"
00004 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
00005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00006 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00008
00009 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
00010 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
00011
00012 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
00013
00014 #include <iostream>
00015
00016 namespace spr{
00017
00018 std::vector<spr::propagatedTrackID> propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, bool debug) {
00019
00020 std::vector<spr::propagatedTrackID> vdets;
00021 spr::propagateCALO(trkCollection,geo,bField,theTrackQuality, vdets, debug);
00022 return vdets;
00023 }
00024
00025 void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, std::vector<spr::propagatedTrackID>& vdets, bool debug) {
00026
00027 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00028 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00029 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00030 reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
00031
00032 unsigned indx;
00033 reco::TrackCollection::const_iterator trkItr;
00034 for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
00035 const reco::Track* pTrack = &(*trkItr);
00036 propagatedTrackID vdet;
00037 vdet.trkItr = trkItr;
00038 vdet.ok = (pTrack->quality(trackQuality_));
00039 vdet.detIdECAL = DetId(0);
00040 vdet.detIdHCAL = DetId(0);
00041 vdet.detIdEHCAL= DetId(0);
00042 if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << vdet.ok << std::endl;
00043
00044 std::pair<math::XYZPoint,bool> info = spr::propagateECAL (pTrack, bField, debug);
00045 vdet.okECAL = info.second;
00046 if (vdet.okECAL) {
00047 const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
00048 vdet.etaECAL = point.eta();
00049 vdet.phiECAL = point.phi();
00050 if (std::abs(point.eta())<1.479) {
00051 vdet.detIdECAL = barrelGeom->getClosestCell(point);
00052 } else {
00053 vdet.detIdECAL = endcapGeom->getClosestCell(point);
00054 }
00055 vdet.detIdEHCAL = gHB->getClosestCell(point);
00056 }
00057 info = spr::propagateHCAL (pTrack, bField, debug);
00058 vdet.okHCAL = info.second;
00059 if (vdet.okHCAL) {
00060 const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
00061 vdet.etaHCAL = point.eta();
00062 vdet.phiHCAL = point.phi();
00063 vdet.detIdHCAL = gHB->getClosestCell(point);
00064 }
00065
00066 vdets.push_back(vdet);
00067 }
00068
00069 if (debug) {
00070 std::cout << "propagateCALO:: for " << vdets.size() << " tracks" << std::endl;
00071 for (unsigned int i=0; i<vdets.size(); ++i) {
00072 std::cout << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL << ") ";
00073 if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
00074 std::cout << (EBDetId)(vdets[i].detIdECAL);
00075 } else {
00076 std::cout << (EEDetId)(vdets[i].detIdECAL);
00077 }
00078 std::cout << " HCAL (" << vdets[i].okHCAL << ") " << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL) << std::endl;
00079 }
00080 }
00081 }
00082
00083 void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, std::vector<spr::propagatedTrackDirection>& trkDir, bool debug) {
00084
00085 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00086 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00087 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00088 reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
00089
00090 unsigned indx;
00091 reco::TrackCollection::const_iterator trkItr;
00092 for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
00093 const reco::Track* pTrack = &(*trkItr);
00094 propagatedTrackDirection trkD;
00095 trkD.trkItr = trkItr;
00096 trkD.ok = (pTrack->quality(trackQuality_));
00097 trkD.detIdECAL = DetId(0);
00098 trkD.detIdHCAL = DetId(0);
00099 trkD.detIdEHCAL= DetId(0);
00100 if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << trkD.ok << std::endl;
00101
00102 spr::propagatedTrack info = spr::propagateTrackToECAL (pTrack, bField, debug);
00103 GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
00104 trkD.okECAL = info.ok;
00105 trkD.pointECAL = point;
00106 trkD.directionECAL = info.direction;
00107 if (trkD.okECAL) {
00108 if (std::abs(info.point.eta())<1.479) {
00109 trkD.detIdECAL = barrelGeom->getClosestCell(point);
00110 } else {
00111 trkD.detIdECAL = endcapGeom->getClosestCell(point);
00112 }
00113 trkD.detIdEHCAL = gHB->getClosestCell(point);
00114 }
00115 info = spr::propagateTrackToHCAL (pTrack, bField, debug);
00116 point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
00117 trkD.okHCAL = info.ok;
00118 trkD.pointHCAL = point;
00119 trkD.directionHCAL = info.direction;
00120 if (trkD.okHCAL) {
00121 trkD.detIdHCAL = gHB->getClosestCell(point);
00122 }
00123 trkDir.push_back(trkD);
00124 }
00125
00126 if (debug) {
00127 std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
00128 for (unsigned int i=0; i<trkDir.size(); ++i) {
00129 std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
00130 if (trkDir[i].okECAL) {
00131 std::cout << " point " << trkDir[i].pointECAL << " direction "
00132 << trkDir[i].directionECAL << " ";
00133 if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
00134 std::cout << (EBDetId)(trkDir[i].detIdECAL);
00135 } else {
00136 std::cout << (EEDetId)(trkDir[i].detIdECAL);
00137 }
00138 }
00139 std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
00140 if (trkDir[i].okHCAL) {
00141 std::cout << " point " << trkDir[i].pointHCAL << " direction "
00142 << trkDir[i].directionHCAL << " "
00143 << (HcalDetId)(trkDir[i].detIdHCAL);
00144 }
00145 std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
00146 }
00147 }
00148 }
00149
00150 propagatedTrack propagateTrackToECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00151 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00152 GlobalVector momentum (track->px(), track->py(), track->pz());
00153 int charge (track->charge());
00154 return spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
00155 }
00156
00157 std::pair<math::XYZPoint,bool> propagateECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00158 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00159 GlobalVector momentum (track->px(), track->py(), track->pz());
00160 int charge (track->charge());
00161 return spr::propagateECAL (vertex, momentum, charge, bfield, debug);
00162 }
00163
00164 std::pair<math::XYZPoint,bool> propagateECAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
00165 spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
00166 return std::pair<math::XYZPoint,bool>(track.point,track.ok);
00167 }
00168
00169 propagatedTrack propagateTrackToHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00170 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00171 GlobalVector momentum (track->px(), track->py(), track->pz());
00172 int charge (track->charge());
00173 return spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
00174 }
00175
00176 std::pair<math::XYZPoint,bool> propagateHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00177 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00178 GlobalVector momentum (track->px(), track->py(), track->pz());
00179 int charge (track->charge());
00180 return spr::propagateHCAL (vertex, momentum, charge, bfield, debug);
00181 }
00182
00183 std::pair<math::XYZPoint,bool> propagateHCAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
00184 spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
00185 return std::pair<math::XYZPoint,bool>(track.point,track.ok);
00186 }
00187
00188 std::pair<math::XYZPoint,bool> propagateTracker(const reco::Track *track, const MagneticField* bfield, bool debug) {
00189 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00190 GlobalVector momentum (track->px(), track->py(), track->pz());
00191 int charge (track->charge());
00192 spr::propagatedTrack track1 = spr::propagateCalo (vertex, momentum, charge, bfield, 290.0, 109.0, 1.705, debug);
00193 return std::pair<math::XYZPoint,bool>(track1.point,track1.ok);
00194 }
00195
00196 std::pair<math::XYZPoint,double> propagateTrackerEnd(const reco::Track *track, const MagneticField* bField, bool debug) {
00197
00198 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00199 GlobalVector momentum (track->px(), track->py(), track->pz());
00200 int charge (track->charge());
00201 float radius = track->outerPosition().Rho();
00202 float zdist = track->outerPosition().Z();
00203 if (debug) std::cout << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " Radius " << radius << " Z " << zdist << std::endl;
00204
00205 FreeTrajectoryState fts (vertex, momentum, charge, bField);
00206 Plane::PlanePointer endcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
00207 Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
00208
00209 AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
00210
00211 TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
00212 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
00213
00214 math::XYZPoint point(-999.,-999.,-999.);
00215 bool ok=false;
00216 GlobalVector direction(0,0,1);
00217 if (tsosb.isValid() && std::abs(zdist) < 110) {
00218 point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00219 direction = tsosb.globalDirection();
00220 ok = true;
00221 } else if (tsose.isValid()) {
00222 point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00223 direction = tsose.globalDirection();
00224 ok = true;
00225 }
00226
00227 double length = -1;
00228 if (ok) {
00229 math::XYZPoint vDiff(point.x()-vertex.x(), point.y()-vertex.y(), point.z()-vertex.z());
00230 double dphi = direction.phi()-momentum.phi();
00231 double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
00232 double rat = 0.5*dphi/std::sin(0.5*dphi);
00233 double dZ = vDiff.z();
00234 double dS = rdist*rat;
00235 length = std::sqrt(dS*dS+dZ*dZ);
00236 if (debug)
00237 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;
00238 }
00239
00240 return std::pair<math::XYZPoint,double>(point,length);
00241 }
00242
00243 spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex, const GlobalVector& tpMomentum, int tpCharge, const MagneticField* bField, float zdist, float radius, float corner, bool debug) {
00244
00245 spr::propagatedTrack track;
00246 if (debug) std::cout << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge " << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner << std::endl;
00247 FreeTrajectoryState fts (tpVertex, tpMomentum, tpCharge, bField);
00248
00249 Plane::PlanePointer lendcap = Plane::build(Plane::PositionType (0, 0, -zdist), Plane::RotationType());
00250 Plane::PlanePointer rendcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
00251
00252 Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
00253
00254 AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
00255
00256 TrajectoryStateOnSurface tsose;
00257 if (tpMomentum.eta() < 0) {
00258 tsose = myAP.propagate(fts, *lendcap);
00259 } else {
00260 tsose = myAP.propagate(fts, *rendcap);
00261 }
00262
00263 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
00264
00265 track.ok=true;
00266 if (tsose.isValid() && tsosb.isValid()) {
00267 float absEta = std::abs(tsosb.globalPosition().eta());
00268 if (absEta < corner) {
00269 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00270 track.direction = tsosb.globalDirection();
00271 } else {
00272 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00273 track.direction = tsose.globalDirection();
00274 }
00275 } else if (tsose.isValid()) {
00276 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00277 track.direction = tsose.globalDirection();
00278 } else if (tsosb.isValid()) {
00279 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00280 track.direction = tsosb.globalDirection();
00281 } else {
00282 track.point.SetXYZ(-999., -999., -999.);
00283 track.direction = GlobalVector(0,0,1);
00284 track.ok = false;
00285 }
00286 if (debug) {
00287 std::cout << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << track.ok << " Point " << track.point << " Direction " << track.direction << std::endl;
00288 if (track.ok) {
00289 math::XYZPoint vDiff(track.point.x()-tpVertex.x(), track.point.y()-tpVertex.y(), track.point.z()-tpVertex.z());
00290 double dphi = track.direction.phi()-tpMomentum.phi();
00291 double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
00292 double pt = tpMomentum.perp();
00293 double rat = 0.5*dphi/std::sin(0.5*dphi);
00294 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;
00295 }
00296 }
00297 return track;
00298 }
00299
00300 }
1.7.3