ConversionSeedGenerators/plugins/PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.cc

0001 #include "PhotonConversionTrajectorySeedProducerFromSingleLegAlgo.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003 #include "FWCore/Utilities/interface/Likely.h"
0004
0005 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "MagneticField/Engine/interface/MagneticField.h"
0008
0009 //#define debugTSPFSLA
0010
0011 namespace {
0012   inline double sqr(double a) { return a * a; }
0013 }  // namespace
0014
0015 PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::PhotonConversionTrajectorySeedProducerFromSingleLegAlgo(
0016     const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0017     : theHitsGenerator(new CombinedHitPairGeneratorForPhotonConversion(
0018           conf.getParameter<edm::ParameterSet>("OrderedHitsFactoryPSet"), iC)),
0019       theSeedCreator(new SeedForPhotonConversion1Leg(conf.getParameter<edm::ParameterSet>("SeedCreatorPSet"), iC)),
0020       theRegionProducer(
0021           new GlobalTrackingRegionProducerFromBeamSpot(conf.getParameter<edm::ParameterSet>("RegionFactoryPSet"), iC)),
0022       theClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet"), iC),
0023       theSilentOnClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet")
0024                                   .getUntrackedParameter<bool>("silentClusterCheck", false)),
0025       _vtxMinDoF(conf.getParameter<double>("vtxMinDoF")),
0026       _maxDZSigmas(conf.getParameter<double>("maxDZSigmas")),
0027       _maxNumSelVtx(conf.getParameter<uint32_t>("maxNumSelVtx")),
0028       _applyTkVtxConstraint(conf.getParameter<bool>("applyTkVtxConstraint")),
0029       _countSeedTracks(0),
0030       _primaryVtxInputTag(conf.getParameter<edm::InputTag>("primaryVerticesTag")),
0031       _beamSpotInputTag(conf.getParameter<edm::InputTag>("beamSpotInputTag")) {
0032   token_vertex = iC.consumes<reco::VertexCollection>(_primaryVtxInputTag);
0033   token_bs = iC.consumes<reco::BeamSpot>(_beamSpotInputTag);
0034   token_refitter = iC.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("TrackRefitter"));
0035   token_magField = iC.esConsumes();
0036 }
0037
0038 PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::~PhotonConversionTrajectorySeedProducerFromSingleLegAlgo() {}
0039
0040 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::find(const edm::Event& event,
0041                                                                    const edm::EventSetup& setup,
0042                                                                    TrajectorySeedCollection& output) {
0043   myEsetup = &setup;
0044   myEvent = &event;
0045
0046   assert(seedCollection == nullptr);
0047
0048   seedCollection = &output;
0049
0050   size_t clustsOrZero = theClusterCheck.tooManyClusters(event);
0051   if (clustsOrZero) {
0052     if (!theSilentOnClusterCheck)
0053       edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
0054     seedCollection = nullptr;
0055     return;
0056   }
0057
0058   const auto& magField = setup.getData(token_magField);
0059   if (UNLIKELY(magField.inTesla(GlobalPoint(0., 0., 0.)).z() < 0.01)) {
0060     seedCollection = nullptr;
0061     return;
0062   }
0063
0064   _IdealHelixParameters.setMagnField(&magField);
0065
0066   event.getByToken(token_vertex, vertexHandle);
0067   if (!vertexHandle.isValid() || vertexHandle->empty()) {
0068     edm::LogError("PhotonConversionFinderFromTracks")
0069         << "Error! Can't get the product primary Vertex Collection " << _primaryVtxInputTag << "\n";
0070     seedCollection = nullptr;
0071     return;
0072   }
0073
0074   event.getByToken(token_bs, recoBeamSpotHandle);
0075
0076   regions = theRegionProducer->regions(event, setup);
0077
0078   // find seeds
0079   loopOnTracks();
0080
0081 #ifdef debugTSPFSLA
0082   std::stringstream ss;
0083   ss.str("");
0084   ss << "\n++++++++++++++++++\n";
0085   ss << "seed collection size " << seedCollection->size();
0086   for (auto const& tjS : *seedCollection) {
0087     po.print(ss, tjS);
0088   }
0089   edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
0090   //-------------------------------------------------
0091 #endif
0092
0093   // clear memory
0094   theHitsGenerator->clearCache();
0095
0096   seedCollection = nullptr;
0097 }
0098
0099 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnTracks() {
0100   //--- Get Tracks
0101   myEvent->getByToken(token_refitter, trackCollectionH);
0102
0103   if (trackCollectionH.isValid() == 0) {
0104     edm::LogError("MissingInput")
0105         << " could not find track collecion in PhotonConversionTrajectorySeedProducerFromSingleLegAlgo";
0106     return;
0107   }
0108
0109   size_t idx = 0;
0110   _countSeedTracks = 0;
0111
0112 #ifdef debugTSPFSLA
0113   size_t sel = 0;
0114   ss.str("");
0115 #endif
0116
0117   for (reco::TrackCollection::const_iterator tr = trackCollectionH->begin(); tr != trackCollectionH->end();
0118        tr++, idx++) {
0119     if (rejectTrack(*tr))
0120       continue;
0121     std::vector<reco::Vertex> selectedPriVtxCompatibleWithTrack;
0122     if (!_applyTkVtxConstraint) {
0123       selectedPriVtxCompatibleWithTrack.push_back(*(vertexHandle->begin()));  //Same approach as before
0124     } else {
0125       if (!selectPriVtxCompatibleWithTrack(*tr, selectedPriVtxCompatibleWithTrack))
0126         continue;
0127     }
0128 #ifdef debugTSPFSLA
0129     sel++;
0130 #endif
0131     loopOnPriVtx(*tr, selectedPriVtxCompatibleWithTrack);
0132   }
0133 #ifdef debugTSPFSLA
0134   edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
0135   edm::LogInfo("debugTrajSeedFromSingleLeg") << "Inspected " << sel << " tracks over " << idx
0136                                              << " tracks. \t # tracks providing at least one seed " << _countSeedTracks;
0137 #endif
0138 }
0139
0140 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::selectPriVtxCompatibleWithTrack(
0141     const reco::Track& tk, std::vector<reco::Vertex>& selectedPriVtxCompatibleWithTrack) {
0142   std::vector<std::pair<double, short> > idx;
0143   short count = -1;
0144
0145   double cosPhi = tk.px() / tk.pt();
0146   double sinPhi = tk.py() / tk.pt();
0147   double sphi2 = tk.covariance(2, 2);
0148   double stheta2 = tk.covariance(1, 1);
0149
0150   for (const reco::Vertex& vtx : *vertexHandle) {
0151     count++;
0152     if (vtx.ndof() <= _vtxMinDoF)
0153       continue;
0154
0155     double _dz = tk.dz(vtx.position());
0156     double _dzError = tk.dzError();
0157
0158     double cotTheta = tk.pz() / tk.pt();
0159     double dx = vtx.position().x();
0160     double dy = vtx.position().y();
0161     double sx2 = vtx.covariance(0, 0);
0162     double sy2 = vtx.covariance(1, 1);
0163
0164     double sxy2 = sqr(cosPhi * cotTheta) * sx2 + sqr(sinPhi * cotTheta) * sy2 +
0165                   sqr(cotTheta * (-dx * sinPhi + dy * cosPhi)) * sphi2 +
0166                   sqr((1 + cotTheta * cotTheta) * (dx * cosPhi + dy * sinPhi)) * stheta2;
0167
0168     _dzError = sqrt(
0169         _dzError * _dzError + vtx.covariance(2, 2) +
0170         sxy2);  //there is a missing component, related to the element (vtx.x*px+vtx.y*py)/pt * pz/pt. since the tk ref point is at the point of closest approach, this scalar product should be almost zero.
0171
0172 #ifdef debugTSPFSLA
0173     ss << " primary vtx " << vtx.position() << " \tk vz " << tk.vz() << " vx " << tk.vx() << " vy " << tk.vy()
0174        << " pz/pt " << tk.pz() / tk.pt() << " \t dz " << _dz << " \t " << _dzError << " sxy2 " << sxy2
0175        << " \t dz/dzErr " << _dz / _dzError << std::endl;
0176 #endif
0177
0178     if (fabs(_dz) / _dzError > _maxDZSigmas)
0179       continue;
0180
0181     idx.push_back(std::pair<double, short>(fabs(_dz), count));
0182   }
0183   if (idx.empty()) {
0184 #ifdef debugTSPFSLA
0185     ss << "no vertex selected " << std::endl;
0186 #endif
0187     return false;
0188   }
0189
0190   std::stable_sort(
0191       idx.begin(), idx.end(), [](std::pair<double, short> a, std::pair<double, short> b) { return a.first < b.first; });
0192
0193   for (size_t i = 0; i < _maxNumSelVtx && i < idx.size(); ++i) {
0194     selectedPriVtxCompatibleWithTrack.push_back((*vertexHandle)[idx[i].second]);
0195 #ifdef debugTSPFSLA
0196     ss << "selected vtx dz " << idx[0].first << "  position" << idx[0].second << std::endl;
0197 #endif
0198   }
0199
0200   return true;
0201 }
0202
0203 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnPriVtx(
0204     const reco::Track& tk, const std::vector<reco::Vertex>& selectedPriVtxCompatibleWithTrack) {
0205   bool foundAtLeastASeedCand = false;
0206   for (auto const& vtx : selectedPriVtxCompatibleWithTrack) {
0207     math::XYZPoint primaryVertexPoint = math::XYZPoint(vtx.position());
0208
0209     for (IR ir = regions.begin(), irEnd = regions.end(); ir < irEnd; ++ir) {
0210       const TrackingRegion& region = **ir;
0211
0212 #ifdef debugTSPFSLA
0213       ss << "[PrintRegion] " << region.print() << std::endl;
0214 #endif
0215
0216       //This if is just for the _countSeedTracks. otherwise
0217       //inspectTrack(&tk,region, primaryVertexPoint);
0218       //would be enough
0219
0220       if (inspectTrack(&tk, region, primaryVertexPoint) and !foundAtLeastASeedCand) {
0221         foundAtLeastASeedCand = true;
0222         _countSeedTracks++;
0223       }
0224     }
0225   }
0226 }
0227
0228 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::rejectTrack(const reco::Track& track) {
0229   math::XYZVector beamSpot;
0230   if (recoBeamSpotHandle.isValid()) {
0231     beamSpot = math::XYZVector(recoBeamSpotHandle->position());
0232
0233     _IdealHelixParameters.setData(&track, beamSpot);
0234     if (_IdealHelixParameters.GetTangentPoint().r() == 0) {
0235       //this case means a null results on the _IdealHelixParameters side
0236       return true;
0237     }
0238
0239     float rMin = 2.;  //cm
0240     if (_IdealHelixParameters.GetTangentPoint().rho() < rMin) {
0241       //this case means a track that has the tangent point nearby the primary vertex
0242       // if the track is primary, this number tends to be the primary vertex itself
0243       //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
0244       //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
0245       return true;
0246     }
0247   }
0248
0249   //-------------------------------------------------------
0250   /*
0251   float maxPt2=64.; //Cut on pt^2 Indeed doesn't do almost nothing
0252   if(track.momentum().Perp2() > maxPt2)
0253     return true;
0254   */
0255   //-------------------------------------------------------
0256   //Cut in the barrel eta region FIXME: to be extended to endcaps
0257   /*
0258   float maxEta=1.3;
0259   if(fabs(track.eta()) > maxEta)
0260     return true;
0261   */
0262   //-------------------------------------------------------
0263   //Reject tracks that have a first valid hit in the pixel barrel/endcap layer/disk 1
0264   //assume that the hits are aligned along momentum
0265   /*
0266   const reco::HitPattern& p=track.hitPattern();
0267   for (int i=0; i<p.numberOfHits(); i++) {
0268     uint32_t hit = p.getHitPattern(i);
0269     // if the hit is valid and in pixel barrel, print out the layer
0270     if (! p.validHitFilter(hit) ) continue;
0271     if( (p.pixelBarrelHitFilter(hit) || p.pixelEndcapHitFilter(hit))
0272     &&
0273     p.getLayer(hit) == 1
0274     )
0275       return true;
0276     else
0277       break; //because the first valid hit is in a different layer
0278   }
0279   */
0280   //-------------------------------------------------------
0281
0282   return false;
0283 }
0284
0285 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::inspectTrack(const reco::Track* track,
0286                                                                            const TrackingRegion& region,
0287                                                                            math::XYZPoint& primaryVertexPoint) {
0288   _IdealHelixParameters.setData(track, primaryVertexPoint);
0289
0290   if (edm::isNotFinite(_IdealHelixParameters.GetTangentPoint().r()) ||
0291       (_IdealHelixParameters.GetTangentPoint().r() == 0)) {
0292     //this case means a null results on the _IdealHelixParameters side
0293     return false;
0294   }
0295
0296   float rMin = 3.;  //cm
0297   if (_IdealHelixParameters.GetTangentPoint().rho() < rMin) {
0298     //this case means a track that has the tangent point nearby the primary vertex
0299     // if the track is primary, this number tends to be the primary vertex itself
0300     //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
0301     //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
0302     return false;
0303   }
0304
0305   float ptmin = 0.5;
0306   float originRBound = 3;
0307   float originZBound = 3.;
0308
0309   GlobalPoint originPos;
0310   originPos = GlobalPoint(_IdealHelixParameters.GetTangentPoint().x(),
0311                           _IdealHelixParameters.GetTangentPoint().y(),
0312                           _IdealHelixParameters.GetTangentPoint().z());
0313   float cotTheta;
0314   if (std::abs(_IdealHelixParameters.GetMomentumAtTangentPoint().rho()) > 1.e-4f) {
0315     cotTheta =
0316         _IdealHelixParameters.GetMomentumAtTangentPoint().z() / _IdealHelixParameters.GetMomentumAtTangentPoint().rho();
0317   } else {
0318     if (_IdealHelixParameters.GetMomentumAtTangentPoint().z() > 0)
0319       cotTheta = 99999.f;
0320     else
0321       cotTheta = -99999.f;
0322   }
0323   GlobalVector originBounds(originRBound, originRBound, originZBound);
0324
0325   GlobalPoint pvtxPoint(primaryVertexPoint.x(), primaryVertexPoint.y(), primaryVertexPoint.z());
0326
0327   ConversionRegion convRegion(originPos, pvtxPoint, cotTheta, track->thetaError(), -1 * track->charge());
0328
0329 #ifdef debugTSPFSLA
0330   ss << "\nConversion Point " << originPos << " " << originPos.perp() << "\n";
0331 #endif
0332
0333   const OrderedSeedingHits& hitss = theHitsGenerator->run(convRegion, region, *myEvent, *myEsetup);
0334
0335   unsigned int nHitss = hitss.size();
0336
0337   if (nHitss == 0)
0338     return false;
0339
0340 #ifdef debugTSPFSLA
0341   ss << "\n nHitss " << nHitss << "\n";
0342 #endif
0343
0344   if (seedCollection->empty())
0345     seedCollection->reserve(
0346         nHitss);  // don't do multiple reserves in the case of multiple regions: it would make things even worse
0347                   // as it will cause N re-allocations instead of the normal log(N)/log(2)
0348   for (unsigned int iHits = 0; iHits < nHitss; ++iHits) {
0349 #ifdef debugTSPFSLA
0350     ss << "\n iHits " << iHits << "\n";
0351 #endif
0352     const SeedingHitSet& hits = hitss[iHits];
0353     theSeedCreator->trajectorySeed(
0354         *seedCollection, hits, originPos, originBounds, ptmin, *myEsetup, convRegion.cotTheta(), ss);
0355   }
0356   return true;
0357 }