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, sel = 0;
0110   _countSeedTracks = 0;
0111
0112 #ifdef debugTSPFSLA
0113   ss.str("");
0114 #endif
0115
0116   for (reco::TrackCollection::const_iterator tr = trackCollectionH->begin(); tr != trackCollectionH->end();
0117        tr++, idx++) {
0118     if (rejectTrack(*tr))
0119       continue;
0120     std::vector<reco::Vertex> selectedPriVtxCompatibleWithTrack;
0121     if (!_applyTkVtxConstraint) {
0122       selectedPriVtxCompatibleWithTrack.push_back(*(vertexHandle->begin()));  //Same approach as before
0123     } else {
0124       if (!selectPriVtxCompatibleWithTrack(*tr, selectedPriVtxCompatibleWithTrack))
0125         continue;
0126     }
0127
0128     sel++;
0129     loopOnPriVtx(*tr, selectedPriVtxCompatibleWithTrack);
0130   }
0131 #ifdef debugTSPFSLA
0132   edm::LogInfo("debugTrajSeedFromSingleLeg") << ss.str();
0133   edm::LogInfo("debugTrajSeedFromSingleLeg") << "Inspected " << sel << " tracks over " << idx
0134                                              << " tracks. \t # tracks providing at least one seed " << _countSeedTracks;
0135 #endif
0136 }
0137
0138 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::selectPriVtxCompatibleWithTrack(
0139     const reco::Track& tk, std::vector<reco::Vertex>& selectedPriVtxCompatibleWithTrack) {
0140   std::vector<std::pair<double, short> > idx;
0141   short count = -1;
0142
0143   double cosPhi = tk.px() / tk.pt();
0144   double sinPhi = tk.py() / tk.pt();
0145   double sphi2 = tk.covariance(2, 2);
0146   double stheta2 = tk.covariance(1, 1);
0147
0148   for (const reco::Vertex& vtx : *vertexHandle) {
0149     count++;
0150     if (vtx.ndof() <= _vtxMinDoF)
0151       continue;
0152
0153     double _dz = tk.dz(vtx.position());
0154     double _dzError = tk.dzError();
0155
0156     double cotTheta = tk.pz() / tk.pt();
0157     double dx = vtx.position().x();
0158     double dy = vtx.position().y();
0159     double sx2 = vtx.covariance(0, 0);
0160     double sy2 = vtx.covariance(1, 1);
0161
0162     double sxy2 = sqr(cosPhi * cotTheta) * sx2 + sqr(sinPhi * cotTheta) * sy2 +
0163                   sqr(cotTheta * (-dx * sinPhi + dy * cosPhi)) * sphi2 +
0164                   sqr((1 + cotTheta * cotTheta) * (dx * cosPhi + dy * sinPhi)) * stheta2;
0165
0166     _dzError = sqrt(
0167         _dzError * _dzError + vtx.covariance(2, 2) +
0168         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.
0169
0170 #ifdef debugTSPFSLA
0171     ss << " primary vtx " << vtx.position() << " \tk vz " << tk.vz() << " vx " << tk.vx() << " vy " << tk.vy()
0172        << " pz/pt " << tk.pz() / tk.pt() << " \t dz " << _dz << " \t " << _dzError << " sxy2 " << sxy2
0173        << " \t dz/dzErr " << _dz / _dzError << std::endl;
0174 #endif
0175
0176     if (fabs(_dz) / _dzError > _maxDZSigmas)
0177       continue;
0178
0179     idx.push_back(std::pair<double, short>(fabs(_dz), count));
0180   }
0181   if (idx.empty()) {
0182 #ifdef debugTSPFSLA
0183     ss << "no vertex selected " << std::endl;
0184 #endif
0185     return false;
0186   }
0187
0188   std::stable_sort(
0189       idx.begin(), idx.end(), [](std::pair<double, short> a, std::pair<double, short> b) { return a.first < b.first; });
0190
0191   for (size_t i = 0; i < _maxNumSelVtx && i < idx.size(); ++i) {
0192     selectedPriVtxCompatibleWithTrack.push_back((*vertexHandle)[idx[i].second]);
0193 #ifdef debugTSPFSLA
0194     ss << "selected vtx dz " << idx[0].first << "  position" << idx[0].second << std::endl;
0195 #endif
0196   }
0197
0198   return true;
0199 }
0200
0201 void PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::loopOnPriVtx(
0202     const reco::Track& tk, const std::vector<reco::Vertex>& selectedPriVtxCompatibleWithTrack) {
0203   bool foundAtLeastASeedCand = false;
0204   for (auto const& vtx : selectedPriVtxCompatibleWithTrack) {
0205     math::XYZPoint primaryVertexPoint = math::XYZPoint(vtx.position());
0206
0207     for (IR ir = regions.begin(), irEnd = regions.end(); ir < irEnd; ++ir) {
0208       const TrackingRegion& region = **ir;
0209
0210 #ifdef debugTSPFSLA
0211       ss << "[PrintRegion] " << region.print() << std::endl;
0212 #endif
0213
0214       //This if is just for the _countSeedTracks. otherwise
0215       //inspectTrack(&tk,region, primaryVertexPoint);
0216       //would be enough
0217
0218       if (inspectTrack(&tk, region, primaryVertexPoint) and !foundAtLeastASeedCand) {
0219         foundAtLeastASeedCand = true;
0220         _countSeedTracks++;
0221       }
0222     }
0223   }
0224 }
0225
0226 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::rejectTrack(const reco::Track& track) {
0227   math::XYZVector beamSpot;
0228   if (recoBeamSpotHandle.isValid()) {
0229     beamSpot = math::XYZVector(recoBeamSpotHandle->position());
0230
0231     _IdealHelixParameters.setData(&track, beamSpot);
0232     if (_IdealHelixParameters.GetTangentPoint().r() == 0) {
0233       //this case means a null results on the _IdealHelixParameters side
0234       return true;
0235     }
0236
0237     float rMin = 2.;  //cm
0238     if (_IdealHelixParameters.GetTangentPoint().rho() < rMin) {
0239       //this case means a track that has the tangent point nearby the primary vertex
0240       // if the track is primary, this number tends to be the primary vertex itself
0241       //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
0242       //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
0243       return true;
0244     }
0245   }
0246
0247   //-------------------------------------------------------
0248   /*
0249   float maxPt2=64.; //Cut on pt^2 Indeed doesn't do almost nothing
0250   if(track.momentum().Perp2() > maxPt2)
0251     return true;
0252   */
0253   //-------------------------------------------------------
0254   //Cut in the barrel eta region FIXME: to be extended to endcaps
0255   /*
0256   float maxEta=1.3;
0257   if(fabs(track.eta()) > maxEta)
0258     return true;
0259   */
0260   //-------------------------------------------------------
0261   //Reject tracks that have a first valid hit in the pixel barrel/endcap layer/disk 1
0262   //assume that the hits are aligned along momentum
0263   /*
0264   const reco::HitPattern& p=track.hitPattern();
0265   for (int i=0; i<p.numberOfHits(); i++) {
0266     uint32_t hit = p.getHitPattern(i);
0267     // if the hit is valid and in pixel barrel, print out the layer
0268     if (! p.validHitFilter(hit) ) continue;
0269     if( (p.pixelBarrelHitFilter(hit) || p.pixelEndcapHitFilter(hit))
0270     &&
0271     p.getLayer(hit) == 1
0272     )
0273       return true;
0274     else
0275       break; //because the first valid hit is in a different layer
0276   }
0277   */
0278   //-------------------------------------------------------
0279
0280   return false;
0281 }
0282
0283 bool PhotonConversionTrajectorySeedProducerFromSingleLegAlgo::inspectTrack(const reco::Track* track,
0284                                                                            const TrackingRegion& region,
0285                                                                            math::XYZPoint& primaryVertexPoint) {
0286   _IdealHelixParameters.setData(track, primaryVertexPoint);
0287
0288   if (edm::isNotFinite(_IdealHelixParameters.GetTangentPoint().r()) ||
0289       (_IdealHelixParameters.GetTangentPoint().r() == 0)) {
0290     //this case means a null results on the _IdealHelixParameters side
0291     return false;
0292   }
0293
0294   float rMin = 3.;  //cm
0295   if (_IdealHelixParameters.GetTangentPoint().rho() < rMin) {
0296     //this case means a track that has the tangent point nearby the primary vertex
0297     // if the track is primary, this number tends to be the primary vertex itself
0298     //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
0299     //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
0300     return false;
0301   }
0302
0303   float ptmin = 0.5;
0304   float originRBound = 3;
0305   float originZBound = 3.;
0306
0307   GlobalPoint originPos;
0308   originPos = GlobalPoint(_IdealHelixParameters.GetTangentPoint().x(),
0309                           _IdealHelixParameters.GetTangentPoint().y(),
0310                           _IdealHelixParameters.GetTangentPoint().z());
0311   float cotTheta;
0312   if (std::abs(_IdealHelixParameters.GetMomentumAtTangentPoint().rho()) > 1.e-4f) {
0313     cotTheta =
0314         _IdealHelixParameters.GetMomentumAtTangentPoint().z() / _IdealHelixParameters.GetMomentumAtTangentPoint().rho();
0315   } else {
0316     if (_IdealHelixParameters.GetMomentumAtTangentPoint().z() > 0)
0317       cotTheta = 99999.f;
0318     else
0319       cotTheta = -99999.f;
0320   }
0321   GlobalVector originBounds(originRBound, originRBound, originZBound);
0322
0323   GlobalPoint pvtxPoint(primaryVertexPoint.x(), primaryVertexPoint.y(), primaryVertexPoint.z());
0324
0325   ConversionRegion convRegion(originPos, pvtxPoint, cotTheta, track->thetaError(), -1 * track->charge());
0326
0327 #ifdef debugTSPFSLA
0328   ss << "\nConversion Point " << originPos << " " << originPos.perp() << "\n";
0329 #endif
0330
0331   const OrderedSeedingHits& hitss = theHitsGenerator->run(convRegion, region, *myEvent, *myEsetup);
0332
0333   unsigned int nHitss = hitss.size();
0334
0335   if (nHitss == 0)
0336     return false;
0337
0338 #ifdef debugTSPFSLA
0339   ss << "\n nHitss " << nHitss << "\n";
0340 #endif
0341
0342   if (seedCollection->empty())
0343     seedCollection->reserve(
0344         nHitss);  // don't do multiple reserves in the case of multiple regions: it would make things even worse
0345                   // as it will cause N re-allocations instead of the normal log(N)/log(2)
0346   for (unsigned int iHits = 0; iHits < nHitss; ++iHits) {
0347 #ifdef debugTSPFSLA
0348     ss << "\n iHits " << iHits << "\n";
0349 #endif
0350     const SeedingHitSet& hits = hitss[iHits];
0351     theSeedCreator->trajectorySeed(
0352         *seedCollection, hits, originPos, originBounds, ptmin, *myEsetup, convRegion.cotTheta(), ss);
0353   }
0354   return true;
0355 }