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 spr::propagatedTrackDirection propagateCALO(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const CaloGeometry* geo, const MagneticField* bField, bool debug) {
00303
00304 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00305 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00306 const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00307
00308 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
00309 spr::propagatedTrackDirection trkD;
00310 trkD.ok = trk.ok;
00311 trkD.detIdECAL = DetId(0);
00312 trkD.detIdHCAL = DetId(0);
00313 trkD.detIdEHCAL= DetId(0);
00314 if (debug) std::cout << "Propagate track " << thisTrk << " charge " << trk.charge << " position " << trk.position << " p " << trk.momentum << " Flag " << trkD.ok << std::endl;
00315
00316 if (trkD.ok) {
00317 spr::propagatedTrack info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, 319.2, 129.4, 1.479, debug);
00318 GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
00319 trkD.okECAL = info.ok;
00320 trkD.pointECAL = point;
00321 trkD.directionECAL = info.direction;
00322 if (trkD.okECAL) {
00323 if (std::abs(info.point.eta())<1.479) {
00324 trkD.detIdECAL = barrelGeom->getClosestCell(point);
00325 } else {
00326 trkD.detIdECAL = endcapGeom->getClosestCell(point);
00327 }
00328 trkD.detIdEHCAL = gHB->getClosestCell(point);
00329 }
00330
00331 info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, 402.7, 180.7, 1.392, debug);
00332 point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
00333 trkD.okHCAL = info.ok;
00334 trkD.pointHCAL = point;
00335 trkD.directionHCAL = info.direction;
00336 if (trkD.okHCAL) {
00337 trkD.detIdHCAL = gHB->getClosestCell(point);
00338 }
00339 }
00340
00341 if (debug) {
00342 std::cout << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL (" << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << std::endl;
00343 if (trkD.okECAL) {
00344 std::cout << "ECAL point " << trkD.pointECAL << " direction "
00345 << trkD.directionECAL << " ";
00346 if (trkD.detIdECAL.subdetId() == EcalBarrel) {
00347 std::cout << (EBDetId)(trkD.detIdECAL);
00348 } else {
00349 std::cout << (EEDetId)(trkD.detIdECAL);
00350 }
00351 }
00352 if (trkD.okHCAL) {
00353 std::cout << " HCAL point " << trkD.pointHCAL << " direction "
00354 << trkD.directionHCAL << " " << (HcalDetId)(trkD.detIdHCAL);
00355 }
00356 if (trkD.okECAL) std::cout << " Or " << (HcalDetId)(trkD.detIdEHCAL);
00357 std::cout << std::endl;
00358 }
00359
00360 return trkD;
00361 }
00362
00363 propagatedTrack propagateTrackToECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00364 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00365 GlobalVector momentum (track->px(), track->py(), track->pz());
00366 int charge (track->charge());
00367 return spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
00368 }
00369
00370 propagatedTrack propagateTrackToECAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
00371
00372 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
00373 spr::propagatedTrack ptrk;
00374 if (trk.ok)
00375 ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 319.2, 129.4, 1.479, debug);
00376 return ptrk;
00377 }
00378
00379 std::pair<math::XYZPoint,bool> propagateECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00380 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00381 GlobalVector momentum (track->px(), track->py(), track->pz());
00382 int charge (track->charge());
00383 return spr::propagateECAL (vertex, momentum, charge, bfield, debug);
00384 }
00385
00386 std::pair<math::XYZPoint,bool> propagateECAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
00387 spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
00388 return std::pair<math::XYZPoint,bool>(track.point,track.ok);
00389 }
00390
00391 spr::propagatedTrack propagateTrackToHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00392 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00393 GlobalVector momentum (track->px(), track->py(), track->pz());
00394 int charge (track->charge());
00395 return spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
00396 }
00397
00398 spr::propagatedTrack propagateTrackToHCAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
00399 spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
00400 spr::propagatedTrack ptrk;
00401 if (trk.ok)
00402 ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 402.7, 180.7, 1.392, debug);
00403 return ptrk;
00404 }
00405
00406 std::pair<math::XYZPoint,bool> propagateHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
00407 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00408 GlobalVector momentum (track->px(), track->py(), track->pz());
00409 int charge (track->charge());
00410 return spr::propagateHCAL (vertex, momentum, charge, bfield, debug);
00411 }
00412
00413 std::pair<math::XYZPoint,bool> propagateHCAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
00414 spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
00415 return std::pair<math::XYZPoint,bool>(track.point,track.ok);
00416 }
00417
00418 std::pair<math::XYZPoint,bool> propagateTracker(const reco::Track *track, const MagneticField* bfield, bool debug) {
00419 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00420 GlobalVector momentum (track->px(), track->py(), track->pz());
00421 int charge (track->charge());
00422 spr::propagatedTrack track1 = spr::propagateCalo (vertex, momentum, charge, bfield, 290.0, 109.0, 1.705, debug);
00423 return std::pair<math::XYZPoint,bool>(track1.point,track1.ok);
00424 }
00425
00426 std::pair<math::XYZPoint,double> propagateTrackerEnd(const reco::Track *track, const MagneticField* bField, bool debug) {
00427
00428 GlobalPoint vertex (track->vx(), track->vy(), track->vz());
00429 GlobalVector momentum (track->px(), track->py(), track->pz());
00430 int charge (track->charge());
00431 float radius = track->outerPosition().Rho();
00432 float zdist = track->outerPosition().Z();
00433 if (debug) std::cout << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " Radius " << radius << " Z " << zdist << std::endl;
00434
00435 FreeTrajectoryState fts (vertex, momentum, charge, bField);
00436 Plane::PlanePointer endcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
00437 Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
00438
00439 AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
00440
00441 TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
00442 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
00443
00444 math::XYZPoint point(-999.,-999.,-999.);
00445 bool ok=false;
00446 GlobalVector direction(0,0,1);
00447 if (tsosb.isValid() && std::abs(zdist) < 110) {
00448 point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00449 direction = tsosb.globalDirection();
00450 ok = true;
00451 } else if (tsose.isValid()) {
00452 point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00453 direction = tsose.globalDirection();
00454 ok = true;
00455 }
00456
00457 double length = -1;
00458 if (ok) {
00459 math::XYZPoint vDiff(point.x()-vertex.x(), point.y()-vertex.y(), point.z()-vertex.z());
00460 double dphi = direction.phi()-momentum.phi();
00461 double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
00462 double rat = 0.5*dphi/std::sin(0.5*dphi);
00463 double dZ = vDiff.z();
00464 double dS = rdist*rat;
00465 length = std::sqrt(dS*dS+dZ*dZ);
00466 if (debug)
00467 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;
00468 }
00469
00470 return std::pair<math::XYZPoint,double>(point,length);
00471 }
00472
00473 spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex, const GlobalVector& tpMomentum, int tpCharge, const MagneticField* bField, float zdist, float radius, float corner, bool debug) {
00474
00475 spr::propagatedTrack track;
00476 if (debug) std::cout << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge " << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner << std::endl;
00477 FreeTrajectoryState fts (tpVertex, tpMomentum, tpCharge, bField);
00478
00479 Plane::PlanePointer lendcap = Plane::build(Plane::PositionType (0, 0, -zdist), Plane::RotationType());
00480 Plane::PlanePointer rendcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
00481
00482 Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
00483
00484 AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
00485
00486 TrajectoryStateOnSurface tsose;
00487 if (tpMomentum.eta() < 0) {
00488 tsose = myAP.propagate(fts, *lendcap);
00489 } else {
00490 tsose = myAP.propagate(fts, *rendcap);
00491 }
00492
00493 TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
00494
00495 track.ok=true;
00496 if (tsose.isValid() && tsosb.isValid()) {
00497 float absEta = std::abs(tsosb.globalPosition().eta());
00498 if (absEta < corner) {
00499 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00500 track.direction = tsosb.globalDirection();
00501 } else {
00502 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00503 track.direction = tsose.globalDirection();
00504 }
00505 } else if (tsose.isValid()) {
00506 track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
00507 track.direction = tsose.globalDirection();
00508 } else if (tsosb.isValid()) {
00509 track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
00510 track.direction = tsosb.globalDirection();
00511 } else {
00512 track.point.SetXYZ(-999., -999., -999.);
00513 track.direction = GlobalVector(0,0,1);
00514 track.ok = false;
00515 }
00516 if (debug) {
00517 std::cout << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << track.ok << " Point " << track.point << " Direction " << track.direction << std::endl;
00518 if (track.ok) {
00519 math::XYZPoint vDiff(track.point.x()-tpVertex.x(), track.point.y()-tpVertex.y(), track.point.z()-tpVertex.z());
00520 double dphi = track.direction.phi()-tpMomentum.phi();
00521 double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
00522 double pt = tpMomentum.perp();
00523 double rat = 0.5*dphi/std::sin(0.5*dphi);
00524 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;
00525 }
00526 }
00527 return track;
00528 }
00529
00530 spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, bool debug) {
00531
00532 spr::trackAtOrigin trk;
00533
00534 edm::SimTrackContainer::const_iterator itr = SimTk->end();
00535 for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr!= SimTk->end(); simTrkItr++) {
00536 if ( simTrkItr->trackId() == thisTrk ) {
00537 if (debug) std::cout << "matched trackId (maximum occurance) " << thisTrk << " type " << simTrkItr->type() << std::endl;
00538 itr = simTrkItr;
00539 break;
00540 }
00541 }
00542
00543 if (itr != SimTk->end()) {
00544 int vertIndex = itr->vertIndex();
00545 if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
00546 edm::SimVertexContainer::const_iterator simVtxItr= SimVtx->begin();
00547 for (int iv=0; iv<vertIndex; iv++) simVtxItr++;
00548 const math::XYZTLorentzVectorD pos = simVtxItr->position();
00549 const math::XYZTLorentzVectorD mom = itr->momentum();
00550 trk.ok = true;
00551 trk.charge = (int)(itr->charge());
00552 trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
00553 trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
00554 }
00555 }
00556 if (debug) std::cout << "Track flag " << trk.ok << " Position " << trk.position << " Momentum " << trk.momentum << std::endl;;
00557 return trk;
00558 }
00559
00560 }
1.7.3