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 int 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 spr::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 int 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 spr::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 std::vector<spr::propagatedGenTrackID> propagateCALO(const HepMC::GenEvent * genEvent, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {
00151
00152 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00153 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00154 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00155
00156 std::vector<spr::propagatedGenTrackID> trkDir;
00157 unsigned int indx;
00158 HepMC::GenEvent::particle_const_iterator p;
00159 for (p=genEvent->particles_begin(),indx=0; p != genEvent->particles_end(); ++p,++indx) {
00160 spr::propagatedGenTrackID trkD;
00161 trkD.trkItr = p;
00162 trkD.detIdECAL = DetId(0);
00163 trkD.detIdHCAL = DetId(0);
00164 trkD.detIdEHCAL= DetId(0);
00165 trkD.pdgId = ((*p)->pdg_id());
00166 trkD.charge = ((pdt->particle(trkD.pdgId))->ID().threeCharge())/3;
00167 GlobalVector momentum = GlobalVector((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());
00168 if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
00169
00170
00171 if ( (*p)->status()==1 && std::abs((*p)->momentum().eta()) < etaMax ) {
00172 GlobalPoint vertex = GlobalPoint(0.1*(*p)->production_vertex()->position().x(),
00173 0.1*(*p)->production_vertex()->position().y(),
00174 0.1*(*p)->production_vertex()->position().z());
00175 trkD.ok = true;
00176 spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 319.2, 129.4, 1.479, debug);
00177 GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
00178 trkD.okECAL = info.ok;
00179 trkD.pointECAL = point;
00180 trkD.directionECAL = info.direction;
00181 if (trkD.okECAL) {
00182 if (std::abs(info.point.eta())<1.479) {
00183 trkD.detIdECAL = barrelGeom->getClosestCell(point);
00184 } else {
00185 trkD.detIdECAL = endcapGeom->getClosestCell(point);
00186 }
00187 trkD.detIdEHCAL = gHB->getClosestCell(point);
00188 }
00189
00190 info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 402.7, 180.7, 1.392, debug);
00191 point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
00192 trkD.okHCAL = info.ok;
00193 trkD.pointHCAL = point;
00194 trkD.directionHCAL = info.direction;
00195 if (trkD.okHCAL) {
00196 trkD.detIdHCAL = gHB->getClosestCell(point);
00197 }
00198 }
00199 trkDir.push_back(trkD);
00200 }
00201
00202 if (debug) {
00203 std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
00204 for (unsigned int i=0; i<trkDir.size(); ++i) {
00205 if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
00206 if (trkDir[i].okECAL) {
00207 std::cout << " point " << trkDir[i].pointECAL << " direction "
00208 << trkDir[i].directionECAL << " ";
00209 if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
00210 std::cout << (EBDetId)(trkDir[i].detIdECAL);
00211 } else {
00212 std::cout << (EEDetId)(trkDir[i].detIdECAL);
00213 }
00214 }
00215 if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
00216 if (trkDir[i].okHCAL) {
00217 std::cout << " point " << trkDir[i].pointHCAL << " direction "
00218 << trkDir[i].directionHCAL << " "
00219 << (HcalDetId)(trkDir[i].detIdHCAL);
00220 }
00221 if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
00222 }
00223 }
00224 return trkDir;
00225 }
00226
00227 std::vector<spr::propagatedGenParticleID> propagateCALO(edm::Handle<reco::GenParticleCollection>& genParticles, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {
00228
00229 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00230 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00231 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00232
00233 std::vector<spr::propagatedGenParticleID> trkDir;
00234 unsigned int indx;
00235 reco::GenParticleCollection::const_iterator p;
00236 for (p=genParticles->begin(),indx=0; p != genParticles->end(); ++p,++indx) {
00237 spr::propagatedGenParticleID trkD;
00238 trkD.trkItr = p;
00239 trkD.detIdECAL = DetId(0);
00240 trkD.detIdHCAL = DetId(0);
00241 trkD.detIdEHCAL= DetId(0);
00242 trkD.pdgId = (p->pdgId());
00243 trkD.charge = p->charge();
00244 GlobalVector momentum = GlobalVector(p->momentum().x(), p->momentum().y(), p->momentum().z());
00245 if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
00246
00247
00248 if ( p->status()==1 && std::abs(momentum.eta()) < etaMax ) {
00249 GlobalPoint vertex = GlobalPoint(p->vertex().x(), p->vertex().y(), p->vertex().z());
00250 trkD.ok = true;
00251 spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 319.2, 129.4, 1.479, debug);
00252 GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
00253 trkD.okECAL = info.ok;
00254 trkD.pointECAL = point;
00255 trkD.directionECAL = info.direction;
00256 if (trkD.okECAL) {
00257 if (std::abs(info.point.eta())<1.479) {
00258 trkD.detIdECAL = barrelGeom->getClosestCell(point);
00259 } else {
00260 trkD.detIdECAL = endcapGeom->getClosestCell(point);
00261 }
00262 trkD.detIdEHCAL = gHB->getClosestCell(point);
00263 }
00264
00265 info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 402.7, 180.7, 1.392, debug);
00266 point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
00267 trkD.okHCAL = info.ok;
00268 trkD.pointHCAL = point;
00269 trkD.directionHCAL = info.direction;
00270 if (trkD.okHCAL) {
00271 trkD.detIdHCAL = gHB->getClosestCell(point);
00272 }
00273 }
00274 trkDir.push_back(trkD);
00275 }
00276
00277 if (debug) {
00278 std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
00279 for (unsigned int i=0; i<trkDir.size(); ++i) {
00280 if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
00281 if (trkDir[i].okECAL) {
00282 std::cout << " point " << trkDir[i].pointECAL << " direction "
00283 << trkDir[i].directionECAL << " ";
00284 if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
00285 std::cout << (EBDetId)(trkDir[i].detIdECAL);
00286 } else {
00287 std::cout << (EEDetId)(trkDir[i].detIdECAL);
00288 }
00289 }
00290 if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
00291 if (trkDir[i].okHCAL) {
00292 std::cout << " point " << trkDir[i].pointHCAL << " direction "
00293 << trkDir[i].directionHCAL << " "
00294 << (HcalDetId)(trkDir[i].detIdHCAL);
00295 }
00296 if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
00297 }
00298 }
00299 return trkDir;
00300 }
00301
00302 propagatedTrack propagateTrackToECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00303 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00304 GlobalVector momentum (track->px(), track->py(), track->pz());
00305 int charge (track->charge());
00306 return spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
00307 }
00308
00309 propagatedTrack propagateTrackToECAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
00310 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
00311 spr::propagatedTrack ptrk;
00312 if (trk.ok)
00313 ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 319.2, 129.4, 1.479, debug);
00314 return ptrk;
00315 }
00316
00317 std::pair<math::XYZPoint,bool> propagateECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00318 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00319 GlobalVector momentum (track->px(), track->py(), track->pz());
00320 int charge (track->charge());
00321 return spr::propagateECAL (vertex, momentum, charge, bfield, debug);
00322 }
00323
00324 std::pair<math::XYZPoint,bool> propagateECAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
00325 spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
00326 return std::pair<math::XYZPoint,bool>(track.point,track.ok);
00327 }
00328
00329 spr::propagatedTrack propagateTrackToHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00330 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00331 GlobalVector momentum (track->px(), track->py(), track->pz());
00332 int charge (track->charge());
00333 return spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
00334 }
00335
00336 spr::propagatedTrack propagateTrackToHCAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
00337 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
00338 spr::propagatedTrack ptrk;
00339 if (trk.ok)
00340 ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 402.7, 180.7, 1.392, debug);
00341 return ptrk;
00342 }
00343
00344 std::pair<math::XYZPoint,bool> propagateHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00345 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00346 GlobalVector momentum (track->px(), track->py(), track->pz());
00347 int charge (track->charge());
00348 return spr::propagateHCAL (vertex, momentum, charge, bfield, debug);
00349 }
00350
00351 std::pair<math::XYZPoint,bool> propagateHCAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
00352 spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
00353 return std::pair<math::XYZPoint,bool>(track.point,track.ok);
00354 }
00355
00356 std::pair<math::XYZPoint,bool> propagateTracker(const reco::Track *track, const MagneticField* bfield, bool debug) {
00357 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00358 GlobalVector momentum (track->px(), track->py(), track->pz());
00359 int charge (track->charge());
00360 spr::propagatedTrack track1 = spr::propagateCalo (vertex, momentum, charge, bfield, 290.0, 109.0, 1.705, debug);
00361 return std::pair<math::XYZPoint,bool>(track1.point,track1.ok);
00362 }
00363
00364 std::pair<math::XYZPoint,double> propagateTrackerEnd(const reco::Track *track, const MagneticField* bField, bool debug) {
00365
00366 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00367 GlobalVector momentum (track->px(), track->py(), track->pz());
00368 int charge (track->charge());
00369 float radius = track->outerPosition().Rho();
00370 float zdist = track->outerPosition().Z();
00371 if (debug) std::cout << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " Radius " << radius << " Z " << zdist << std::endl;
00372
00373 FreeTrajectoryState fts (vertex, momentum, charge, bField);
00374 Plane::PlanePointer endcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
00375 Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
00376
00377 AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
00378
00379 TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
00380 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
00381
00382 math::XYZPoint point(-999.,-999.,-999.);
00383 bool ok=false;
00384 GlobalVector direction(0,0,1);
00385 if (tsosb.isValid() && std::abs(zdist) < 110) {
00386 point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00387 direction = tsosb.globalDirection();
00388 ok = true;
00389 } else if (tsose.isValid()) {
00390 point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00391 direction = tsose.globalDirection();
00392 ok = true;
00393 }
00394
00395 double length = -1;
00396 if (ok) {
00397 math::XYZPoint vDiff(point.x()-vertex.x(), point.y()-vertex.y(), point.z()-vertex.z());
00398 double dphi = direction.phi()-momentum.phi();
00399 double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
00400 double rat = 0.5*dphi/std::sin(0.5*dphi);
00401 double dZ = vDiff.z();
00402 double dS = rdist*rat;
00403 length = std::sqrt(dS*dS+dZ*dZ);
00404 if (debug)
00405 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;
00406 }
00407
00408 return std::pair<math::XYZPoint,double>(point,length);
00409 }
00410
00411 spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex, const GlobalVector& tpMomentum, int tpCharge, const MagneticField* bField, float zdist, float radius, float corner, bool debug) {
00412
00413 spr::propagatedTrack track;
00414 if (debug) std::cout << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge " << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner << std::endl;
00415 FreeTrajectoryState fts (tpVertex, tpMomentum, tpCharge, bField);
00416
00417 Plane::PlanePointer lendcap = Plane::build(Plane::PositionType (0, 0, -zdist), Plane::RotationType());
00418 Plane::PlanePointer rendcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
00419
00420 Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
00421
00422 AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
00423
00424 TrajectoryStateOnSurface tsose;
00425 if (tpMomentum.eta() < 0) {
00426 tsose = myAP.propagate(fts, *lendcap);
00427 } else {
00428 tsose = myAP.propagate(fts, *rendcap);
00429 }
00430
00431 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
00432
00433 track.ok=true;
00434 if (tsose.isValid() && tsosb.isValid()) {
00435 float absEta = std::abs(tsosb.globalPosition().eta());
00436 if (absEta < corner) {
00437 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00438 track.direction = tsosb.globalDirection();
00439 } else {
00440 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00441 track.direction = tsose.globalDirection();
00442 }
00443 } else if (tsose.isValid()) {
00444 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00445 track.direction = tsose.globalDirection();
00446 } else if (tsosb.isValid()) {
00447 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00448 track.direction = tsosb.globalDirection();
00449 } else {
00450 track.point.SetXYZ(-999., -999., -999.);
00451 track.direction = GlobalVector(0,0,1);
00452 track.ok = false;
00453 }
00454 if (debug) {
00455 std::cout << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << track.ok << " Point " << track.point << " Direction " << track.direction << std::endl;
00456 if (track.ok) {
00457 math::XYZPoint vDiff(track.point.x()-tpVertex.x(), track.point.y()-tpVertex.y(), track.point.z()-tpVertex.z());
00458 double dphi = track.direction.phi()-tpMomentum.phi();
00459 double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
00460 double pt = tpMomentum.perp();
00461 double rat = 0.5*dphi/std::sin(0.5*dphi);
00462 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;
00463 }
00464 }
00465 return track;
00466 }
00467
00468 spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, bool debug) {
00469
00470 spr::trackAtOrigin trk;
00471
00472 edm::SimTrackContainer::const_iterator itr = SimTk->end();
00473 for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr!= SimTk->end(); simTrkItr++) {
00474 if ( simTrkItr->trackId() == thisTrk ) {
00475 if (debug) std::cout << "matched trackId (maximum occurance) " << thisTrk << " type " << simTrkItr->type() << std::endl;
00476 itr = simTrkItr;
00477 break;
00478 }
00479 }
00480
00481 if (itr != SimTk->end()) {
00482 int vertIndex = itr->vertIndex();
00483 if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
00484 edm::SimVertexContainer::const_iterator simVtxItr= SimVtx->begin();
00485 for (int iv=0; iv<vertIndex; iv++) simVtxItr++;
00486 const math::XYZTLorentzVectorD pos = simVtxItr->position();
00487 const math::XYZTLorentzVectorD mom = itr->momentum();
00488 trk.ok = true;
00489 trk.charge = (int)(itr->charge());
00490 trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
00491 trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
00492 }
00493 }
00494 if (debug) std::cout << "Track flag " << trk.ok << " Position " << trk.position << " Momentum " << trk.momentum << std::endl;;
00495 return trk;
00496 }
00497
00498 }
1.7.3