TrackAssociation/src/CosmicParametersDefinerForTP.cc

0001 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0002 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "SimTracker/TrackAssociation/interface/CosmicParametersDefinerForTP.h"
0006 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0007 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0008 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0009 #include <DataFormats/GeometrySurface/interface/GloballyPositioned.h>
0010 #include <DataFormats/GeometrySurface/interface/Surface.h>
0011 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0012
0013 class TrajectoryStateClosestToBeamLineBuilder;
0014
0015 CosmicParametersDefinerForTP::CosmicParametersDefinerForTP(edm::ConsumesCollector iC)
0016     : ParametersDefinerForTP(edm::InputTag("offlineBeamSpot"), iC), geometryToken_(iC.esConsumes()) {}
0017 CosmicParametersDefinerForTP::~CosmicParametersDefinerForTP() = default;
0018
0019 TrackingParticle::Vector CosmicParametersDefinerForTP::momentum(const edm::Event &iEvent,
0020                                                                 const edm::EventSetup &iSetup,
0021                                                                 const TrackingParticleRef &tpr) const {
0022   // to add a new implementation for cosmic. For the moment, it is just as for
0023   // the base class:
0024   using namespace edm;
0025   using namespace std;
0026   using namespace reco;
0027
0028   auto const &bs = iEvent.get(bsToken_);
0029   auto const &geometry = iSetup.getData(geometryToken_);
0030   auto const &mf = iSetup.getData(mfToken_);
0031
0032   GlobalVector finalGV(0, 0, 0);
0033   GlobalPoint finalGP(0, 0, 0);
0034   double radius(9999);
0035   bool found(false);
0036   TrackingParticle::Vector momentum(0, 0, 0);
0037
0038   edm::LogVerbatim("CosmicParametersDefinerForTP") << "\t in CosmicParametersDefinerForTP::momentum";
0039   edm::LogVerbatim("CosmicParametersDefinerForTP")
0040       << "\t \t Original TP state:          pt = " << tpr->pt() << ", pz = " << tpr->pz();
0041
0042   if (simHitsTPAssoc.isValid() == 0) {
0043     LogError("TrackAssociation") << "Invalid handle!";
0044     return momentum;
0045   }
0046   std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(
0047       tpr,
0048       TrackPSimHitRef());  // SimHit is dummy: for simHitTPAssociationListGreater
0049                            // sorting only the cluster is needed
0050   auto range = std::equal_range(simHitsTPAssoc->begin(),
0051                                 simHitsTPAssoc->end(),
0052                                 clusterTPpairWithDummyTP,
0053                                 SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0054   for (auto ip = range.first; ip != range.second; ++ip) {
0055     TrackPSimHitRef it = ip->second;
0056     const GeomDet *tmpDet = geometry.idToDet(DetId(it->detUnitId()));
0057     if (!tmpDet) {
0058       edm::LogVerbatim("CosmicParametersDefinerForTP")
0059           << "***WARNING in CosmicParametersDefinerForTP::momentum: no GeomDet "
0060              "for: "
0061           << it->detUnitId() << ". Skipping it."
0062           << "\n";
0063       continue;
0064     }
0065
0066     LocalVector lv = it->momentumAtEntry();
0067     Local3DPoint lp = it->localPosition();
0068     GlobalVector gv = tmpDet->surface().toGlobal(lv);
0069     GlobalPoint gp = tmpDet->surface().toGlobal(lp);
0070
0071     // discard hits related to low energy debris from the primary particle
0072     if (it->processType() != 0)
0073       continue;
0074
0075     if (gp.perp() < radius) {
0076       found = true;
0077       radius = gp.perp();
0078       finalGV = gv;
0079       finalGP = gp;
0080     }
0081   }
0082
0083   edm::LogVerbatim("CosmicParametersDefinerForTP")
0084       //   <<"\t FINAL State at InnerMost Hit: Radius = "<< finalGP.perp() << ",
0085       //   z = "<< finalGP.z()
0086       //  <<", pt = "<< finalGV.perp() << ", pz = "<< finalGV.z();
0087       << "\t \t FINAL State at InnerMost Hit:   pt = " << finalGV.perp() << ", pz = " << finalGV.z();
0088
0089   if (found) {
0090     FreeTrajectoryState ftsAtProduction(finalGP, finalGV, TrackCharge(tpr->charge()), &mf);
0091     TSCBLBuilderNoMaterial tscblBuilder;
0092     TrajectoryStateClosestToBeamLine tsAtClosestApproach =
0093         tscblBuilder(ftsAtProduction, bs);  // as in TrackProducerAlgorithm
0094
0095     if (tsAtClosestApproach.isValid()) {
0096       GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
0097       momentum = TrackingParticle::Vector(p.x(), p.y(), p.z());
0098     } else {
0099       edm::LogVerbatim("CosmicParametersDefinerForTP") << "*** WARNING in CosmicParametersDefinerForTP::momentum: "
0100                                                           "tsAtClosestApproach is not valid."
0101                                                        << "\n";
0102     }
0103
0104     edm::LogVerbatim("CosmicParametersDefinerForTP")
0105         << "\t \t FINAL State extrap. at PCA: pt = " << sqrt(momentum.x() * momentum.x() + momentum.y() * momentum.y())
0106         << ", pz = " << momentum.z() << "\n";
0107
0108     return momentum;
0109   }
0110
0111   edm::LogVerbatim("CosmicParametersDefinerForTP")
0112       << "*** WARNING in CosmicParametersDefinerForTP::momentum: NOT found the "
0113          "innermost TP point"
0114       << "\n";
0115   edm::LogVerbatim("CosmicParametersDefinerForTP")
0116       << "*** FINAL Reference MOMENTUM TP (px,py,pz) = " << momentum.x() << momentum.y() << momentum.z() << "\n";
0117   return momentum;
0118 }
0119
0120 TrackingParticle::Point CosmicParametersDefinerForTP::vertex(const edm::Event &iEvent,
0121                                                              const edm::EventSetup &iSetup,
0122                                                              const TrackingParticleRef &tpr) const {
0123   using namespace edm;
0124   using namespace std;
0125   using namespace reco;
0126
0127   auto const &bs = iEvent.get(bsToken_);
0128   auto const &geometry = iSetup.getData(geometryToken_);
0129   auto const &mf = iSetup.getData(mfToken_);
0130
0131   GlobalVector finalGV(0, 0, 0);
0132   GlobalPoint finalGP(0, 0, 0);
0133   double radius(9999);
0134   bool found(false);
0135   TrackingParticle::Point vertex(0, 0, 0);
0136
0137   edm::LogVerbatim("CosmicParametersDefinerForTP") << "\t in CosmicParametersDefinerForTP::vertex";
0138   edm::LogVerbatim("CosmicParametersDefinerForTP")
0139       << "\t \t Original TP state:          radius = "
0140       << sqrt(tpr->vertex().x() * tpr->vertex().x() + tpr->vertex().y() * tpr->vertex().y())
0141       << ", z = " << tpr->vertex().z();
0142
0143   if (simHitsTPAssoc.isValid() == 0) {
0144     LogError("TrackAssociation") << "Invalid handle!";
0145     return vertex;
0146   }
0147   std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(
0148       tpr,
0149       TrackPSimHitRef());  // SimHit is dummy: for simHitTPAssociationListGreater
0150                            // sorting only the cluster is needed
0151   auto range = std::equal_range(simHitsTPAssoc->begin(),
0152                                 simHitsTPAssoc->end(),
0153                                 clusterTPpairWithDummyTP,
0154                                 SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0155   for (auto ip = range.first; ip != range.second; ++ip) {
0156     TrackPSimHitRef it = ip->second;
0157     const GeomDet *tmpDet = geometry.idToDet(DetId(it->detUnitId()));
0158     if (!tmpDet) {
0159       edm::LogVerbatim("CosmicParametersDefinerForTP")
0160           << "***WARNING in CosmicParametersDefinerForTP::vertex: no GeomDet "
0161              "for: "
0162           << it->detUnitId() << ". Skipping it."
0163           << "\n";
0164       continue;
0165     }
0166
0167     LocalVector lv = it->momentumAtEntry();
0168     Local3DPoint lp = it->localPosition();
0169     GlobalVector gv = tmpDet->surface().toGlobal(lv);
0170     GlobalPoint gp = tmpDet->surface().toGlobal(lp);
0171
0172     // discard hits related to low energy debris from the primary particle
0173     if (it->processType() != 0)
0174       continue;
0175
0176     if (gp.perp() < radius) {
0177       found = true;
0178       radius = gp.perp();
0179       finalGV = gv;
0180       finalGP = gp;
0181     }
0182   }
0183   edm::LogVerbatim("CosmicParametersDefinerForTP")
0184       << "\t \t FINAL State at InnerMost Hit:   radius = " << finalGP.perp() << ", z = " << finalGP.z();
0185
0186   if (found) {
0187     FreeTrajectoryState ftsAtProduction(finalGP, finalGV, TrackCharge(tpr->charge()), &mf);
0188     TSCBLBuilderNoMaterial tscblBuilder;
0189     TrajectoryStateClosestToBeamLine tsAtClosestApproach =
0190         tscblBuilder(ftsAtProduction, bs);  // as in TrackProducerAlgorithm
0191
0192     if (tsAtClosestApproach.isValid()) {
0193       GlobalPoint v = tsAtClosestApproach.trackStateAtPCA().position();
0194       vertex = TrackingParticle::Point(v.x(), v.y(), v.z());
0195     } else {
0196       // to preserve old behaviour
0197       // would be better to flag this somehow to allow ignoring in downstream
0198       vertex = TrackingParticle::Point(bs.x0(), bs.y0(), bs.z0());
0199       edm::LogVerbatim("CosmicParametersDefinerForTP") << "*** WARNING in CosmicParametersDefinerForTP::vertex: "
0200                                                           "tsAtClosestApproach is not valid."
0201                                                        << "\n";
0202     }
0203     edm::LogVerbatim("CosmicParametersDefinerForTP")
0204         << "\t \t FINAL State extrap. at PCA: radius = " << sqrt(vertex.x() * vertex.x() + vertex.y() * vertex.y())
0205         << ", z = " << vertex.z() << "\n";
0206
0207     return vertex;
0208   }
0209
0210   edm::LogVerbatim("CosmicParametersDefinerForTP")
0211       << "*** WARNING in CosmicParametersDefinerForTP::vertex: NOT found the "
0212          "innermost TP point"
0213       << "\n";
0214   edm::LogVerbatim("CosmicParametersDefinerForTP")
0215       << "*** FINAL Reference VERTEX TP   V(x,y,z) = " << vertex.x() << vertex.y() << vertex.z() << "\n";
0216
0217   return vertex;
0218 }