Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:27:16

0001 /**
0002   \class    TSGForOIFromL2
0003   \brief    Create L3MuonTrajectorySeeds from L2 Muons updated at vertex in an outside-in manner
0004   \author   Benjamin Radburn-Smith, Santiago Folgueras, Bibhuprasad Mahakud, Jan Frederik Schulte (Purdue University, West Lafayette, USA)
0005  */
0006 
0007 #include "RecoMuon/TrackerSeedGenerator/plugins/TSGForOIFromL2.h"
0008 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0011 
0012 #include <memory>
0013 
0014 TSGForOIFromL2::TSGForOIFromL2(const edm::ParameterSet& iConfig)
0015     : src_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0016       maxSeeds_(iConfig.getParameter<uint32_t>("maxSeeds")),
0017       maxHitlessSeeds_(iConfig.getParameter<uint32_t>("maxHitlessSeeds")),
0018       maxHitSeeds_(iConfig.getParameter<uint32_t>("maxHitSeeds")),
0019       numOfLayersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
0020       numOfHitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
0021       numL2ValidHitsCutAllEta_(iConfig.getParameter<uint32_t>("numL2ValidHitsCutAllEta")),
0022       numL2ValidHitsCutAllEndcap_(iConfig.getParameter<uint32_t>("numL2ValidHitsCutAllEndcap")),
0023       fixedErrorRescalingForHits_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHits")),
0024       fixedErrorRescalingForHitless_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHitless")),
0025       adjustErrorsDynamicallyForHits_(iConfig.getParameter<bool>("adjustErrorsDynamicallyForHits")),
0026       adjustErrorsDynamicallyForHitless_(iConfig.getParameter<bool>("adjustErrorsDynamicallyForHitless")),
0027       estimatorName_(iConfig.getParameter<std::string>("estimator")),
0028       minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
0029       maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
0030       useHitLessSeeds_(iConfig.getParameter<bool>("UseHitLessSeeds")),
0031       updator_(new KFUpdator()),
0032       measurementTrackerTag_(
0033           consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
0034       pT1_(iConfig.getParameter<double>("pT1")),
0035       pT2_(iConfig.getParameter<double>("pT2")),
0036       pT3_(iConfig.getParameter<double>("pT3")),
0037       eta1_(iConfig.getParameter<double>("eta1")),
0038       eta2_(iConfig.getParameter<double>("eta2")),
0039       eta3_(iConfig.getParameter<double>("eta3")),
0040       eta4_(iConfig.getParameter<double>("eta4")),
0041       eta5_(iConfig.getParameter<double>("eta5")),
0042       eta6_(iConfig.getParameter<double>("eta6")),
0043       eta7_(iConfig.getParameter<double>("eta7")),
0044       SF1_(iConfig.getParameter<double>("SF1")),
0045       SF2_(iConfig.getParameter<double>("SF2")),
0046       SF3_(iConfig.getParameter<double>("SF3")),
0047       SF4_(iConfig.getParameter<double>("SF4")),
0048       SF5_(iConfig.getParameter<double>("SF5")),
0049       SF6_(iConfig.getParameter<double>("SF6")),
0050       SFHld_(iConfig.getParameter<double>("SFHld")),
0051       SFHd_(iConfig.getParameter<double>("SFHd")),
0052       tsosDiff1_(iConfig.getParameter<double>("tsosDiff1")),
0053       tsosDiff2_(iConfig.getParameter<double>("tsosDiff2")),
0054       displacedReco_(iConfig.getParameter<bool>("displacedReco")),
0055       propagatorName_(iConfig.getParameter<std::string>("propagatorName")),
0056       theCategory_(std::string("Muon|RecoMuon|TSGForOIFromL2")),
0057       estimatorToken_(esConsumes(edm::ESInputTag("", estimatorName_))),
0058       magfieldToken_(esConsumes()),
0059       propagatorToken_(esConsumes(edm::ESInputTag("", propagatorName_))),
0060       tmpTkGeometryToken_(esConsumes()),
0061       geometryToken_(esConsumes()),
0062       sHPOppositeToken_(esConsumes(edm::ESInputTag("", "hltESPSteppingHelixPropagatorOpposite"))) {
0063   produces<std::vector<TrajectorySeed> >();
0064 }
0065 
0066 TSGForOIFromL2::~TSGForOIFromL2() {}
0067 
0068 //
0069 // Produce seeds
0070 //
0071 void TSGForOIFromL2::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0072   // Initialize variables
0073   unsigned int numSeedsMade = 0;
0074   unsigned int layerCount = 0;
0075   unsigned int hitlessSeedsMadeIP = 0;
0076   unsigned int hitlessSeedsMadeMuS = 0;
0077   unsigned int hitSeedsMade = 0;
0078   unsigned int hitSeedsMadeMuS = 0;
0079 
0080   // Surface used to make a TSOS at the PCA to the beamline
0081   Plane::PlanePointer dummyPlane = Plane::build(Plane::PositionType(), Plane::RotationType());
0082 
0083   // Read ESHandles
0084   edm::Handle<MeasurementTrackerEvent> measurementTrackerH;
0085   const edm::ESHandle<Chi2MeasurementEstimatorBase> estimatorH = iSetup.getHandle(estimatorToken_);
0086   const edm::ESHandle<MagneticField> magfieldH = iSetup.getHandle(magfieldToken_);
0087   const edm::ESHandle<Propagator> propagatorAlongH = iSetup.getHandle(propagatorToken_);
0088   const edm::ESHandle<Propagator>& propagatorOppositeH = propagatorAlongH;
0089   const edm::ESHandle<TrackerGeometry> tmpTkGeometryH = iSetup.getHandle(tmpTkGeometryToken_);
0090   const edm::ESHandle<GlobalTrackingGeometry> geometryH = iSetup.getHandle(geometryToken_);
0091 
0092   iEvent.getByToken(measurementTrackerTag_, measurementTrackerH);
0093 
0094   // Read L2 track collection
0095   edm::Handle<reco::TrackCollection> l2TrackCol;
0096   iEvent.getByToken(src_, l2TrackCol);
0097 
0098   // The product
0099   std::unique_ptr<std::vector<TrajectorySeed> > result(new std::vector<TrajectorySeed>());
0100 
0101   // Get vector of Detector layers
0102   std::vector<BarrelDetLayer const*> const& tob = measurementTrackerH->geometricSearchTracker()->tobLayers();
0103   std::vector<ForwardDetLayer const*> const& tecPositive =
0104       tmpTkGeometryH->isThere(GeomDetEnumerators::P2OTEC)
0105           ? measurementTrackerH->geometricSearchTracker()->posTidLayers()
0106           : measurementTrackerH->geometricSearchTracker()->posTecLayers();
0107   std::vector<ForwardDetLayer const*> const& tecNegative =
0108       tmpTkGeometryH->isThere(GeomDetEnumerators::P2OTEC)
0109           ? measurementTrackerH->geometricSearchTracker()->negTidLayers()
0110           : measurementTrackerH->geometricSearchTracker()->negTecLayers();
0111 
0112   // Get suitable propagators
0113   std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(*propagatorAlongH, alongMomentum);
0114   std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(*propagatorOppositeH, oppositeToMomentum);
0115 
0116   // Stepping Helix Propagator for propogation from muon system to tracker
0117   const edm::ESHandle<Propagator> SHPOpposite = iSetup.getHandle(sHPOppositeToken_);
0118 
0119   // Loop over the L2's and make seeds for all of them
0120   LogTrace(theCategory_) << "TSGForOIFromL2::produce: Number of L2's: " << l2TrackCol->size();
0121   for (unsigned int l2TrackColIndex(0); l2TrackColIndex != l2TrackCol->size(); ++l2TrackColIndex) {
0122     const reco::TrackRef l2(l2TrackCol, l2TrackColIndex);
0123 
0124     // Container of Seeds
0125     std::vector<TrajectorySeed> out;
0126     LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: L2 muon pT, eta, phi --> " << l2->pt() << " , " << l2->eta()
0127                                << " , " << l2->phi() << std::endl;
0128 
0129     FreeTrajectoryState fts = trajectoryStateTransform::initialFreeState(*l2, magfieldH.product());
0130 
0131     dummyPlane->move(fts.position() - dummyPlane->position());
0132     TrajectoryStateOnSurface tsosAtIP = TrajectoryStateOnSurface(fts, *dummyPlane);
0133     LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: Created TSOSatIP: " << tsosAtIP << std::endl;
0134 
0135     // Get the TSOS on the innermost layer of the L2
0136     TrajectoryStateOnSurface tsosAtMuonSystem =
0137         trajectoryStateTransform::innerStateOnSurface(*l2, *geometryH, magfieldH.product());
0138     LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: Created TSOSatMuonSystem: " << tsosAtMuonSystem
0139                                << std::endl;
0140 
0141     LogTrace("TSGForOIFromL2")
0142         << "TSGForOIFromL2::produce: Check the error of the L2 parameter and use hit seeds if big errors" << std::endl;
0143 
0144     StateOnTrackerBound fromInside(propagatorAlong.get());
0145     TrajectoryStateOnSurface outerTkStateInside = fromInside(fts);
0146 
0147     StateOnTrackerBound fromOutside(&*SHPOpposite);
0148     TrajectoryStateOnSurface outerTkStateOutside = fromOutside(tsosAtMuonSystem);
0149 
0150     // Check if the two positions (using updated and not-updated TSOS) agree withing certain extent.
0151     // If both TSOSs agree, use only the one at vertex, as it uses more information. If they do not agree, search for seeds based on both.
0152     double L2muonEta = l2->eta();
0153     double absL2muonEta = std::abs(L2muonEta);
0154     bool useBoth = false;
0155     if (outerTkStateInside.isValid() && outerTkStateOutside.isValid()) {
0156       //following commented out variables dist1 (5 par compatibility of tsos at outertracker surface)
0157       //dist2 (angle between two tsos) could further be explored in combination of L2 valid hits for seeding. So kept for
0158       //future developers
0159       //auto dist1 = match_Chi2(outerTkStateInside,outerTkStateOutside);//for future developers
0160       //auto dist2 = deltaR(outerTkStateInside.globalMomentum(),outerTkStateOutside.globalMomentum());//for future developers
0161       //if ((dist1 > tsosDiff1_ || dist2 > tsosDiff2_) && l2->numberOfValidHits() < 20) useBoth = true;//for future developers
0162       if (l2->numberOfValidHits() < numL2ValidHitsCutAllEta_)
0163         useBoth = true;
0164       if (l2->numberOfValidHits() < numL2ValidHitsCutAllEndcap_ && absL2muonEta > eta7_)
0165         useBoth = true;
0166       if (absL2muonEta > eta1_ && absL2muonEta < eta1_)
0167         useBoth = true;
0168     }
0169 
0170     numSeedsMade = 0;
0171     hitlessSeedsMadeIP = 0;
0172     hitlessSeedsMadeMuS = 0;
0173     hitSeedsMade = 0;
0174     hitSeedsMadeMuS = 0;
0175 
0176     // calculate scale factors
0177     double errorSFHits = (adjustErrorsDynamicallyForHits_ ? calculateSFFromL2(l2) : fixedErrorRescalingForHits_);
0178     double errorSFHitless =
0179         (adjustErrorsDynamicallyForHitless_ ? calculateSFFromL2(l2) : fixedErrorRescalingForHitless_);
0180 
0181     // BARREL
0182     if (absL2muonEta < maxEtaForTOB_) {
0183       layerCount = 0;
0184       for (auto it = tob.rbegin(); it != tob.rend(); ++it) {
0185         LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: looping in TOB layer " << layerCount << std::endl;
0186         if (useHitLessSeeds_ && hitlessSeedsMadeIP < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0187           makeSeedsWithoutHits(**it,
0188                                tsosAtIP,
0189                                *(propagatorAlong.get()),
0190                                estimatorH,
0191                                errorSFHitless,
0192                                hitlessSeedsMadeIP,
0193                                numSeedsMade,
0194                                out);
0195 
0196         // Do not create hitbased seeds in barrel region
0197         if (absL2muonEta > 1.0 && hitSeedsMade < maxHitSeeds_ && numSeedsMade < maxSeeds_)
0198           makeSeedsFromHits(**it,
0199                             tsosAtIP,
0200                             *(propagatorAlong.get()),
0201                             estimatorH,
0202                             measurementTrackerH,
0203                             errorSFHits,
0204                             hitSeedsMade,
0205                             numSeedsMade,
0206                             layerCount,
0207                             out);
0208 
0209         if (useBoth && !displacedReco_) {
0210           if (useHitLessSeeds_ && hitlessSeedsMadeMuS < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0211             makeSeedsWithoutHits(**it,
0212                                  outerTkStateOutside,
0213                                  *(propagatorOpposite.get()),
0214                                  estimatorH,
0215                                  errorSFHitless,
0216                                  hitlessSeedsMadeMuS,
0217                                  numSeedsMade,
0218                                  out);
0219         }
0220       }
0221       LogTrace("TSGForOIFromL2") << "TSGForOIFromL2:::produce: NumSeedsMade = " << numSeedsMade
0222                                  << " , layerCount = " << layerCount << std::endl;
0223     }
0224 
0225     // Reset number of seeds if in overlap region
0226     if (absL2muonEta > minEtaForTEC_ && absL2muonEta < maxEtaForTOB_ && !displacedReco_) {
0227       numSeedsMade = 0;
0228       hitlessSeedsMadeIP = 0;
0229       hitlessSeedsMadeMuS = 0;
0230       hitSeedsMade = 0;
0231     }
0232 
0233     // ENDCAP+
0234     if (L2muonEta > minEtaForTEC_) {
0235       layerCount = 0;
0236       for (auto it = tecPositive.rbegin(); it != tecPositive.rend(); ++it) {
0237         LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: looping in TEC+ layer " << layerCount << std::endl;
0238         if (useHitLessSeeds_ && hitlessSeedsMadeIP < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0239           makeSeedsWithoutHits(**it,
0240                                tsosAtIP,
0241                                *(propagatorAlong.get()),
0242                                estimatorH,
0243                                errorSFHitless,
0244                                hitlessSeedsMadeIP,
0245                                numSeedsMade,
0246                                out);
0247 
0248         if (absL2muonEta > 1.0 && hitSeedsMade < maxHitSeeds_ && numSeedsMade < maxSeeds_)
0249           makeSeedsFromHits(**it,
0250                             tsosAtIP,
0251                             *(propagatorAlong.get()),
0252                             estimatorH,
0253                             measurementTrackerH,
0254                             errorSFHits,
0255                             hitSeedsMade,
0256                             numSeedsMade,
0257                             layerCount,
0258                             out);
0259 
0260         if (useBoth && !displacedReco_) {
0261           if (useHitLessSeeds_ && hitlessSeedsMadeMuS < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0262             makeSeedsWithoutHits(**it,
0263                                  outerTkStateOutside,
0264                                  *(propagatorOpposite.get()),
0265                                  estimatorH,
0266                                  errorSFHitless,
0267                                  hitlessSeedsMadeMuS,
0268                                  numSeedsMade,
0269                                  out);
0270         }
0271       }
0272       LogTrace("TSGForOIFromL2") << "TSGForOIFromL2:::produce: NumSeedsMade = " << numSeedsMade
0273                                  << " , layerCount = " << layerCount << std::endl;
0274     }
0275 
0276     // ENDCAP-
0277     if (L2muonEta < -minEtaForTEC_) {
0278       layerCount = 0;
0279       for (auto it = tecNegative.rbegin(); it != tecNegative.rend(); ++it) {
0280         LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: looping in TEC- layer " << layerCount << std::endl;
0281         if (useHitLessSeeds_ && hitlessSeedsMadeIP < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0282           makeSeedsWithoutHits(**it,
0283                                tsosAtIP,
0284                                *(propagatorAlong.get()),
0285                                estimatorH,
0286                                errorSFHitless,
0287                                hitlessSeedsMadeIP,
0288                                numSeedsMade,
0289                                out);
0290 
0291         if (absL2muonEta > 1.0 && hitSeedsMade < maxHitSeeds_ && numSeedsMade < maxSeeds_)
0292           makeSeedsFromHits(**it,
0293                             tsosAtIP,
0294                             *(propagatorAlong.get()),
0295                             estimatorH,
0296                             measurementTrackerH,
0297                             errorSFHits,
0298                             hitSeedsMade,
0299                             numSeedsMade,
0300                             layerCount,
0301                             out);
0302 
0303         if (useBoth && !displacedReco_) {
0304           if (useHitLessSeeds_ && hitlessSeedsMadeMuS < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0305             makeSeedsWithoutHits(**it,
0306                                  outerTkStateOutside,
0307                                  *(propagatorOpposite.get()),
0308                                  estimatorH,
0309                                  errorSFHitless,
0310                                  hitlessSeedsMadeMuS,
0311                                  numSeedsMade,
0312                                  out);
0313         }
0314       }
0315       LogTrace("TSGForOIFromL2") << "TSGForOIFromL2:::produce: NumSeedsMade = " << numSeedsMade
0316                                  << " , layerCount = " << layerCount << std::endl;
0317     }
0318 
0319     // Displaced Reconstruction
0320     if (displacedReco_ && outerTkStateOutside.isValid()) {
0321       layerCount = 0;
0322       for (auto it = tob.rbegin(); it != tob.rend(); ++it) {
0323         LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: looping in TOB layer " << layerCount;
0324         if (useHitLessSeeds_ && hitlessSeedsMadeMuS < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0325           makeSeedsWithoutHits(**it,
0326                                outerTkStateOutside,
0327                                *(propagatorOpposite.get()),
0328                                estimatorH,
0329                                errorSFHitless * SFHld_,
0330                                hitlessSeedsMadeMuS,
0331                                numSeedsMade,
0332                                out);
0333         if (hitSeedsMadeMuS < maxHitSeeds_ && numSeedsMade < maxSeeds_)
0334           makeSeedsFromHits(**it,
0335                             outerTkStateOutside,
0336                             *(propagatorOpposite.get()),
0337                             estimatorH,
0338                             measurementTrackerH,
0339                             errorSFHits * SFHd_,
0340                             hitSeedsMadeMuS,
0341                             numSeedsMade,
0342                             layerCount,
0343                             out);
0344       }
0345       LogTrace("TSGForOIFromL2") << "TSGForOIFromL2:::produce: NumSeedsMade = " << numSeedsMade
0346                                  << " , layerCount = " << layerCount;
0347       if (L2muonEta >= 0.0) {
0348         layerCount = 0;
0349         for (auto it = tecPositive.rbegin(); it != tecPositive.rend(); ++it) {
0350           LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: looping in TEC+ layer " << layerCount << std::endl;
0351           if (useHitLessSeeds_ && hitlessSeedsMadeMuS < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0352             makeSeedsWithoutHits(**it,
0353                                  outerTkStateOutside,
0354                                  *(propagatorOpposite.get()),
0355                                  estimatorH,
0356                                  errorSFHitless * SFHld_,
0357                                  hitlessSeedsMadeMuS,
0358                                  numSeedsMade,
0359                                  out);
0360           if (hitSeedsMadeMuS < maxHitSeeds_ && numSeedsMade < maxSeeds_)
0361             makeSeedsFromHits(**it,
0362                               outerTkStateOutside,
0363                               *(propagatorOpposite.get()),
0364                               estimatorH,
0365                               measurementTrackerH,
0366                               errorSFHits * SFHd_,
0367                               hitSeedsMadeMuS,
0368                               numSeedsMade,
0369                               layerCount,
0370                               out);
0371         }
0372         LogTrace("TSGForOIFromL2") << "TSGForOIFromL2:::produce: NumSeedsMade = " << numSeedsMade
0373                                    << " , layerCount = " << layerCount;
0374       }
0375 
0376       else {
0377         layerCount = 0;
0378         for (auto it = tecNegative.rbegin(); it != tecNegative.rend(); ++it) {
0379           LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::produce: looping in TEC- layer " << layerCount;
0380           if (useHitLessSeeds_ && hitlessSeedsMadeMuS < maxHitlessSeeds_ && numSeedsMade < maxSeeds_)
0381             makeSeedsWithoutHits(**it,
0382                                  outerTkStateOutside,
0383                                  *(propagatorOpposite.get()),
0384                                  estimatorH,
0385                                  errorSFHitless * SFHld_,
0386                                  hitlessSeedsMadeMuS,
0387                                  numSeedsMade,
0388                                  out);
0389           if (hitSeedsMadeMuS < maxHitSeeds_ && numSeedsMade < maxSeeds_)
0390             makeSeedsFromHits(**it,
0391                               outerTkStateOutside,
0392                               *(propagatorOpposite.get()),
0393                               estimatorH,
0394                               measurementTrackerH,
0395                               errorSFHits * SFHd_,
0396                               hitSeedsMadeMuS,
0397                               numSeedsMade,
0398                               layerCount,
0399                               out);
0400         }
0401         LogTrace("TSGForOIFromL2") << "TSGForOIFromL2:::produce: NumSeedsMade = " << numSeedsMade
0402                                    << " , layerCount = " << layerCount;
0403       }
0404     }
0405 
0406     for (std::vector<TrajectorySeed>::iterator it = out.begin(); it != out.end(); ++it) {
0407       result->push_back(*it);
0408     }
0409 
0410   }  // L2Collection
0411 
0412   edm::LogInfo(theCategory_) << "TSGForOIFromL2::produce: number of seeds made: " << result->size();
0413 
0414   iEvent.put(std::move(result));
0415 }
0416 
0417 //
0418 // Create seeds without hits on a given layer (TOB or TEC)
0419 //
0420 void TSGForOIFromL2::makeSeedsWithoutHits(const GeometricSearchDet& layer,
0421                                           const TrajectoryStateOnSurface& tsos,
0422                                           const Propagator& propagatorAlong,
0423                                           const edm::ESHandle<Chi2MeasurementEstimatorBase>& estimator,
0424                                           double errorSF,
0425                                           unsigned int& hitlessSeedsMade,
0426                                           unsigned int& numSeedsMade,
0427                                           std::vector<TrajectorySeed>& out) const {
0428   // create hitless seeds
0429   LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsWithoutHits: Start hitless" << std::endl;
0430   std::vector<GeometricSearchDet::DetWithState> dets;
0431   layer.compatibleDetsV(tsos, propagatorAlong, *estimator, dets);
0432   if (!dets.empty()) {
0433     auto const& detOnLayer = dets.front().first;
0434     auto const& tsosOnLayer = dets.front().second;
0435     LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsWithoutHits: tsosOnLayer " << tsosOnLayer << std::endl;
0436     if (!tsosOnLayer.isValid()) {
0437       edm::LogInfo(theCategory_) << "ERROR!: Hitless TSOS is not valid!";
0438     } else {
0439       dets.front().second.rescaleError(errorSF);
0440       PTrajectoryStateOnDet const& ptsod =
0441           trajectoryStateTransform::persistentState(tsosOnLayer, detOnLayer->geographicalId().rawId());
0442       TrajectorySeed::RecHitContainer rHC;
0443       out.push_back(TrajectorySeed(ptsod, rHC, oppositeToMomentum));
0444       LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsWithoutHits: TSOS (Hitless) done " << std::endl;
0445       hitlessSeedsMade++;
0446       numSeedsMade++;
0447     }
0448   }
0449 }
0450 
0451 //
0452 // Find hits on a given layer (TOB or TEC) and create seeds from updated TSOS with hit
0453 //
0454 void TSGForOIFromL2::makeSeedsFromHits(const GeometricSearchDet& layer,
0455                                        const TrajectoryStateOnSurface& tsos,
0456                                        const Propagator& propagatorAlong,
0457                                        const edm::ESHandle<Chi2MeasurementEstimatorBase>& estimator,
0458                                        const edm::Handle<MeasurementTrackerEvent>& measurementTracker,
0459                                        double errorSF,
0460                                        unsigned int& hitSeedsMade,
0461                                        unsigned int& numSeedsMade,
0462                                        unsigned int& layerCount,
0463                                        std::vector<TrajectorySeed>& out) const {
0464   if (layerCount > numOfLayersToTry_)
0465     return;
0466 
0467   // Error Rescaling
0468   TrajectoryStateOnSurface onLayer(tsos);
0469   onLayer.rescaleError(errorSF);
0470 
0471   std::vector<GeometricSearchDet::DetWithState> dets;
0472   layer.compatibleDetsV(onLayer, propagatorAlong, *estimator, dets);
0473 
0474   // Find Measurements on each DetWithState
0475   LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsFromHits: Find measurements on each detWithState  "
0476                              << dets.size() << std::endl;
0477   std::vector<TrajectoryMeasurement> meas;
0478   for (std::vector<GeometricSearchDet::DetWithState>::iterator it = dets.begin(); it != dets.end(); ++it) {
0479     MeasurementDetWithData det = measurementTracker->idToDet(it->first->geographicalId());
0480     if (det.isNull())
0481       continue;
0482     if (!it->second.isValid())
0483       continue;  // Skip if TSOS is not valid
0484 
0485     std::vector<TrajectoryMeasurement> mymeas =
0486         det.fastMeasurements(it->second, onLayer, propagatorAlong, *estimator);  // Second TSOS is not used
0487     for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2;
0488          ++it2) {
0489       if (it2->recHit()->isValid())
0490         meas.push_back(*it2);  // Only save those which are valid
0491     }
0492   }
0493 
0494   // Update TSOS using TMs after sorting, then create Trajectory Seed and put into vector
0495   LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsFromHits: Update TSOS using TMs after sorting, then create "
0496                                 "Trajectory Seed, number of TM = "
0497                              << meas.size() << std::endl;
0498   std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0499 
0500   unsigned int found = 0;
0501   for (std::vector<TrajectoryMeasurement>::const_iterator it = meas.begin(); it != meas.end(); ++it) {
0502     TrajectoryStateOnSurface updatedTSOS = updator_->update(it->forwardPredictedState(), *it->recHit());
0503     LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsFromHits: TSOS for TM " << found << std::endl;
0504     if (not updatedTSOS.isValid())
0505       continue;
0506 
0507     edm::OwnVector<TrackingRecHit> seedHits;
0508     seedHits.push_back(*it->recHit()->hit());
0509     PTrajectoryStateOnDet const& pstate =
0510         trajectoryStateTransform::persistentState(updatedTSOS, it->recHit()->geographicalId().rawId());
0511     LogTrace("TSGForOIFromL2") << "TSGForOIFromL2::makeSeedsFromHits: Number of seedHits: " << seedHits.size()
0512                                << std::endl;
0513     TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0514     out.push_back(seed);
0515     found++;
0516     numSeedsMade++;
0517     hitSeedsMade++;
0518     if (found == numOfHitsToTry_)
0519       break;
0520     if (hitSeedsMade > maxHitSeeds_)
0521       return;
0522   }
0523 
0524   if (found)
0525     layerCount++;
0526 }
0527 
0528 //
0529 // Calculate the dynamic error SF by analysing the L2
0530 //
0531 double TSGForOIFromL2::calculateSFFromL2(const reco::TrackRef track) const {
0532   double theSF = 1.0;
0533   // L2 direction vs pT blowup - as was previously done:
0534   // Split into 4 pT ranges: <pT1_, pT1_<pT2_, pT2_<pT3_, <pT4_: 13,30,70
0535   // Split into different eta ranges depending in pT
0536   double abseta = std::abs(track->eta());
0537   if (track->pt() <= pT1_)
0538     theSF = SF1_;
0539   else if (track->pt() > pT1_ && track->pt() <= pT2_) {
0540     if (abseta <= eta3_)
0541       theSF = SF3_;
0542     else if (abseta > eta3_ && abseta <= eta6_)
0543       theSF = SF2_;
0544     else if (abseta > eta6_)
0545       theSF = SF3_;
0546   } else if (track->pt() > pT2_ && track->pt() <= pT3_) {
0547     if (abseta <= eta1_)
0548       theSF = SF6_;
0549     else if (abseta > eta1_ && abseta <= eta2_)
0550       theSF = SF4_;
0551     else if (abseta > eta2_ && abseta <= eta3_)
0552       theSF = SF6_;
0553     else if (abseta > eta3_ && abseta <= eta4_)
0554       theSF = SF1_;
0555     else if (abseta > eta4_ && abseta <= eta5_)
0556       theSF = SF1_;
0557     else if (abseta > eta5_)
0558       theSF = SF5_;
0559   } else if (track->pt() > pT3_) {
0560     if (abseta <= eta3_)
0561       theSF = SF5_;
0562     else if (abseta > eta3_ && abseta <= eta4_)
0563       theSF = SF4_;
0564     else if (abseta > eta4_ && abseta <= eta5_)
0565       theSF = SF4_;
0566     else if (abseta > eta5_)
0567       theSF = SF5_;
0568   }
0569 
0570   LogTrace(theCategory_) << "TSGForOIFromL2::calculateSFFromL2: SF has been calculated as: " << theSF;
0571 
0572   return theSF;
0573 }
0574 
0575 //
0576 // calculate Chi^2 of two trajectory states
0577 //
0578 double TSGForOIFromL2::match_Chi2(const TrajectoryStateOnSurface& tsos1, const TrajectoryStateOnSurface& tsos2) const {
0579   if (!tsos1.isValid() || !tsos2.isValid())
0580     return -1.;
0581 
0582   AlgebraicVector5 v(tsos1.localParameters().vector() - tsos2.localParameters().vector());
0583   AlgebraicSymMatrix55 m(tsos1.localError().matrix() + tsos2.localError().matrix());
0584 
0585   bool ierr = !m.Invert();
0586 
0587   if (ierr) {
0588     edm::LogInfo("TSGForOIFromL2") << "Error inverting covariance matrix";
0589     return -1;
0590   }
0591 
0592   double est = ROOT::Math::Similarity(v, m);
0593 
0594   return est;
0595 }
0596 
0597 //
0598 //
0599 //
0600 void TSGForOIFromL2::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0601   edm::ParameterSetDescription desc;
0602   desc.add<edm::InputTag>("src", edm::InputTag("hltL2Muons", "UpdatedAtVtx"));
0603   desc.add<int>("layersToTry", 2);
0604   desc.add<double>("fixedErrorRescaleFactorForHitless", 2.0);
0605   desc.add<int>("hitsToTry", 1);
0606   desc.add<bool>("adjustErrorsDynamicallyForHits", false);
0607   desc.add<bool>("adjustErrorsDynamicallyForHitless", true);
0608   desc.add<edm::InputTag>("MeasurementTrackerEvent", edm::InputTag("hltSiStripClusters"));
0609   desc.add<bool>("UseHitLessSeeds", true);
0610   desc.add<std::string>("estimator", "hltESPChi2MeasurementEstimator100");
0611   desc.add<double>("maxEtaForTOB", 1.8);
0612   desc.add<double>("minEtaForTEC", 0.7);
0613   desc.addUntracked<bool>("debug", false);
0614   desc.add<double>("fixedErrorRescaleFactorForHits", 1.0);
0615   desc.add<unsigned int>("maxSeeds", 20);
0616   desc.add<unsigned int>("maxHitlessSeeds", 5);
0617   desc.add<unsigned int>("maxHitSeeds", 1);
0618   desc.add<unsigned int>("numL2ValidHitsCutAllEta", 20);
0619   desc.add<unsigned int>("numL2ValidHitsCutAllEndcap", 30);
0620   desc.add<double>("pT1", 13.0);
0621   desc.add<double>("pT2", 30.0);
0622   desc.add<double>("pT3", 70.0);
0623   desc.add<double>("eta1", 0.2);
0624   desc.add<double>("eta2", 0.3);
0625   desc.add<double>("eta3", 1.0);
0626   desc.add<double>("eta4", 1.2);
0627   desc.add<double>("eta5", 1.6);
0628   desc.add<double>("eta6", 1.4);
0629   desc.add<double>("eta7", 2.1);
0630   desc.add<double>("SF1", 3.0);
0631   desc.add<double>("SF2", 4.0);
0632   desc.add<double>("SF3", 5.0);
0633   desc.add<double>("SF4", 7.0);
0634   desc.add<double>("SF5", 10.0);
0635   desc.add<double>("SF6", 2.0);
0636   desc.add<double>("SFHld", 2.0)->setComment("Scale Factor used to rescale the TSOS error of the hitless seeds");
0637   desc.add<double>("SFHd", 4.0)->setComment("Scale Factor used to rescale the TSOS error of the hit based seeds");
0638   desc.add<double>("tsosDiff1", 0.2);
0639   desc.add<double>("tsosDiff2", 0.02);
0640   desc.add<bool>("displacedReco", false)->setComment("Flag to turn on the displaced seeding");
0641   desc.add<std::string>("propagatorName", "PropagatorWithMaterialParabolicMf");
0642   descriptions.add("TSGForOIFromL2", desc);
0643 }
0644 
0645 DEFINE_FWK_MODULE(TSGForOIFromL2);