SpecialSeedGenerators/src/CtfSpecialSeedGenerator.cc

0001 #include "RecoTracker/SpecialSeedGenerators/interface/CtfSpecialSeedGenerator.h"
0002 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0003 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0004 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0005 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0006
0007 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducerFactory.h"
0008 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0009 #include "RecoTracker/TkTrackingRegions/interface/OrderedHitsGeneratorFactory.h"
0010 #include "RecoTracker/TkSeedingLayers/interface/OrderedSeedingHits.h"
0011 #include "FWCore/Framework/interface/ConsumesCollector.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015
0016 using namespace ctfseeding;
0017
0018 CtfSpecialSeedGenerator::CtfSpecialSeedGenerator(const edm::ParameterSet& conf)
0019     : conf_(conf),
0020       theMFToken(esConsumes<edm::Transition::BeginRun>()),
0021       theBuilderToken(
0022           esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", conf_.getParameter<std::string>("TTRHBuilder")))),
0023       theTrackerToken(esConsumes<edm::Transition::BeginRun>()),
0024       thePropAlongToken(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "PropagatorWithMaterial"))),
0025       thePropOppositeToken(
0026           esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "PropagatorWithMaterialOpposite"))),
0027       theTopoToken(esConsumes()),
0028       requireBOFF(conf.getParameter<bool>("requireBOFF")),
0029       theMaxSeeds(conf.getParameter<int32_t>("maxSeeds")),
0030       check(conf, consumesCollector())
0031
0032 {
0033   useScintillatorsConstraint = conf_.getParameter<bool>("UseScintillatorsConstraint");
0034   edm::LogVerbatim("CtfSpecialSeedGenerator") << "Constructing CtfSpecialSeedGenerator";
0035   produces<TrajectorySeedCollection>();
0036
0037   edm::ParameterSet regfactoryPSet = conf_.getParameter<edm::ParameterSet>("RegionFactoryPSet");
0038   std::string regfactoryName = regfactoryPSet.getParameter<std::string>("ComponentName");
0039   theRegionProducer = std::unique_ptr<TrackingRegionProducer>{
0040       TrackingRegionProducerFactory::get()->create(regfactoryName, regfactoryPSet, consumesCollector())};
0041
0042   std::vector<edm::ParameterSet> pSets = conf_.getParameter<std::vector<edm::ParameterSet>>("OrderedHitsFactoryPSets");
0043   std::vector<edm::ParameterSet>::const_iterator iPSet;
0044   edm::ConsumesCollector iC = consumesCollector();
0045   for (iPSet = pSets.begin(); iPSet != pSets.end(); iPSet++) {
0046     std::string hitsfactoryName = iPSet->getParameter<std::string>("ComponentName");
0047     theGenerators.emplace_back(OrderedHitsGeneratorFactory::get()->create(hitsfactoryName, *iPSet, iC));
0048   }
0049 }
0050
0051 CtfSpecialSeedGenerator::~CtfSpecialSeedGenerator() = default;
0052
0053 void CtfSpecialSeedGenerator::endRun(edm::Run const&, edm::EventSetup const&) { theSeedBuilder.reset(); }
0054
0055 void CtfSpecialSeedGenerator::beginRun(edm::Run const&, const edm::EventSetup& iSetup) {
0056   theMagfield = iSetup.getHandle(theMFToken);
0057   theBuilder = iSetup.getHandle(theBuilderToken);
0058   theTracker = iSetup.getHandle(theTrackerToken);
0059
0060   edm::LogVerbatim("CtfSpecialSeedGenerator") << "Initializing...";
0061   if (useScintillatorsConstraint) {
0062     edm::ParameterSet upperScintPar = conf_.getParameter<edm::ParameterSet>("UpperScintillatorParameters");
0063     edm::ParameterSet lowerScintPar = conf_.getParameter<edm::ParameterSet>("LowerScintillatorParameters");
0064     RectangularPlaneBounds upperBounds(
0065         upperScintPar.getParameter<double>("WidthInX"), upperScintPar.getParameter<double>("LenghtInZ"), 1);
0066     GlobalPoint upperPosition(upperScintPar.getParameter<double>("GlobalX"),
0067                               upperScintPar.getParameter<double>("GlobalY"),
0068                               upperScintPar.getParameter<double>("GlobalZ"));
0069     edm::LogVerbatim("CtfSpecialSeedGenerator") << "Upper Scintillator position x, y, z " << upperPosition.x() << ", "
0070                                                 << upperPosition.y() << ", " << upperPosition.z();
0071     RectangularPlaneBounds lowerBounds(
0072         lowerScintPar.getParameter<double>("WidthInX"), lowerScintPar.getParameter<double>("LenghtInZ"), 1);
0073     GlobalPoint lowerPosition(lowerScintPar.getParameter<double>("GlobalX"),
0074                               lowerScintPar.getParameter<double>("GlobalY"),
0075                               lowerScintPar.getParameter<double>("GlobalZ"));
0076     edm::LogVerbatim("CtfSpecialSeedGenerator") << "Lower Scintillator position x, y, z " << lowerPosition.x() << ", "
0077                                                 << lowerPosition.y() << ", " << lowerPosition.z();
0078     TkRotation<float> rot(1, 0, 0, 0, 0, 1, 0, 1, 0);
0079     upperScintillator = BoundPlane::build(upperPosition, rot, &upperBounds);
0080     lowerScintillator = BoundPlane::build(lowerPosition, rot, &lowerBounds);
0081   }
0082
0083   edm::ESHandle<Propagator> propagatorAlongHandle = iSetup.getHandle(thePropAlongToken);
0084   edm::ESHandle<Propagator> propagatorOppositeHandle = iSetup.getHandle(thePropOppositeToken);
0085
0086   std::vector<edm::ParameterSet> pSets = conf_.getParameter<std::vector<edm::ParameterSet>>("OrderedHitsFactoryPSets");
0087   std::vector<edm::ParameterSet>::const_iterator iPSet;
0088   for (iPSet = pSets.begin(); iPSet != pSets.end(); iPSet++) {
0089     std::string propagationDirection = iPSet->getParameter<std::string>("PropagationDirection");
0090     if (propagationDirection == "alongMomentum")
0091       thePropDirs.push_back(alongMomentum);
0092     else
0093       thePropDirs.push_back(oppositeToMomentum);
0094     std::string navigationDirection = iPSet->getParameter<std::string>("NavigationDirection");
0095     if (navigationDirection == "insideOut")
0096       theNavDirs.push_back(insideOut);
0097     else
0098       theNavDirs.push_back(outsideIn);
0099     edm::LogVerbatim("CtfSpecialSeedGenerator") << "hitsGenerator done";
0100   }
0101   bool setMomentum = conf_.getParameter<bool>("SetMomentum");
0102   std::vector<int> charges;
0103   if (setMomentum) {
0104     charges = conf_.getParameter<std::vector<int>>("Charges");
0105   }
0106   theSeedBuilder = std::make_unique<SeedFromGenericPairOrTriplet>(theMagfield.product(),
0107                                                                   theTracker.product(),
0108                                                                   theBuilder.product(),
0109                                                                   propagatorAlongHandle.product(),
0110                                                                   propagatorOppositeHandle.product(),
0111                                                                   charges,
0112                                                                   setMomentum,
0113                                                                   conf_.getParameter<double>("ErrorRescaling"));
0114   double p = 1;
0115   if (setMomentum) {
0116     p = conf_.getParameter<double>("SeedMomentum");
0117     theSeedBuilder->setMomentumTo(p);
0118   }
0119 }
0120
0121 void CtfSpecialSeedGenerator::produce(edm::Event& e, const edm::EventSetup& iSetup) {
0122   // get Inputs
0123   auto output = std::make_unique<TrajectorySeedCollection>();
0124
0125   //check on the number of clusters
0126   if (!requireBOFF || (theMagfield->inTesla(GlobalPoint(0, 0, 0)).mag() == 0.00)) {
0127     size_t clustsOrZero = check.tooManyClusters(e);
0128     if (!clustsOrZero) {
0129       bool ok = run(iSetup, e, *output);
0130       if (!ok) {
0131       }  // nothing to do
0132     } else
0133       edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
0134   }
0135
0136   edm::LogVerbatim("CtfSpecialSeedGenerator") << " number of seeds = " << output->size();
0137   e.put(std::move(output));
0138 }
0139
0140 bool CtfSpecialSeedGenerator::run(const edm::EventSetup& iSetup,
0141                                   const edm::Event& e,
0142                                   TrajectorySeedCollection& output) {
0143   std::vector<std::unique_ptr<TrackingRegion>> regions = theRegionProducer->regions(e, iSetup);
0144   bool ok = true;
0145   for (auto iReg = regions.begin(); iReg != regions.end(); iReg++) {
0146     if (!theSeedBuilder->momentumFromPSet())
0147       theSeedBuilder->setMomentumTo((*iReg)->ptMin());
0148     int i = 0;
0149     for (auto iGen = theGenerators.begin(); iGen != theGenerators.end(); iGen++) {
0150       ok = buildSeeds(iSetup, e, (*iGen)->run(**iReg, e, iSetup), theNavDirs[i], thePropDirs[i], output);
0151       i++;
0152       if (!ok)
0153         break;
0154     }
0155     if (!ok)
0156       break;
0157   }
0158   return ok;
0159 }
0160
0161 bool CtfSpecialSeedGenerator::buildSeeds(const edm::EventSetup& iSetup,
0162                                          const edm::Event& e,
0163                                          const OrderedSeedingHits& osh,
0164                                          const NavigationDirection& navdir,
0165                                          const PropagationDirection& dir,
0166                                          TrajectorySeedCollection& output) {
0167   //SeedFromGenericPairOrTriplet seedBuilder(conf_, magfield.product(), tracker.product(), theBuilder.product());
0168   edm::LogInfo("CtfSpecialSeedGenerator") << "osh.size() " << osh.size();
0169   for (unsigned int i = 0; i < osh.size(); i++) {
0170     const SeedingHitSet& shs = osh[i];
0171     if (preliminaryCheck(shs, iSetup)) {
0172       std::vector<TrajectorySeed*> seeds = theSeedBuilder->seed(shs, dir, navdir, iSetup);
0173       for (std::vector<TrajectorySeed*>::const_iterator iSeed = seeds.begin(); iSeed != seeds.end(); iSeed++) {
0174         if (!*iSeed) {
0175           edm::LogError("CtfSpecialSeedGenerator") << "a seed pointer is null. skipping.";
0176           continue;
0177         }
0178         if (postCheck(**iSeed)) {
0179           output.push_back(**iSeed);
0180         }
0181         delete *iSeed;
0182         edm::LogVerbatim("CtfSpecialSeedGenerator") << "Seed built";
0183       }
0184     }
0185   }
0186   if ((theMaxSeeds > 0) && (output.size() > size_t(theMaxSeeds))) {
0187     edm::LogWarning("TooManySeeds") << "Too many seeds (" << output.size() << "), bailing out.\n";
0188     output.clear();
0189     return false;
0190   }
0191   return true;
0192 }
0193 //checks the hits are on diffrent layers
0194 bool CtfSpecialSeedGenerator::preliminaryCheck(const SeedingHitSet& shs, const edm::EventSetup& es) {
0195   edm::ESHandle<TrackerTopology> tTopo = es.getHandle(theTopoToken);
0196
0197   std::vector<std::pair<unsigned int, unsigned int>> vSubdetLayer;
0198   //std::vector<std::string> vSeedLayerNames;
0199   bool checkHitsAtPositiveY = conf_.getParameter<bool>("SeedsFromPositiveY");
0200   //***top-bottom
0201   bool checkHitsAtNegativeY = conf_.getParameter<bool>("SeedsFromNegativeY");
0202   //***
0203   bool checkHitsOnDifferentLayers = conf_.getParameter<bool>("CheckHitsAreOnDifferentLayers");
0204   unsigned int nHits = shs.size();
0205   for (unsigned int iHit = 0; iHit < nHits; ++iHit) {
0206     //hits for the seeds must be at positive y
0207     auto trh = shs[iHit];
0208     auto recHit = trh;
0209     GlobalPoint hitPos = recHit->globalPosition();
0210     //GlobalPoint point =
0211     //  theTracker->idToDet(iHits->geographicalId() )->surface().toGlobal(iHits->localPosition());
0212     if (checkHitsAtPositiveY) {
0213       if (hitPos.y() < 0)
0214         return false;
0215     }
0216     //***top-bottom
0217     if (checkHitsAtNegativeY) {
0218       if (hitPos.y() > 0)
0219         return false;
0220     }
0221     //***
0222     //std::string name = iHits->seedinglayer().name();
0223     //hits for the seeds must be in different layers
0224     unsigned int subid = (*trh).geographicalId().subdetId();
0225     unsigned int layer = tTopo->layer((*trh).geographicalId());
0226     std::vector<std::pair<unsigned int, unsigned int>>::const_iterator iter;
0227     //std::vector<std::string>::const_iterator iNames;
0228     if (checkHitsOnDifferentLayers) {
0229       for (iter = vSubdetLayer.begin(); iter != vSubdetLayer.end(); iter++) {
0230         if (iter->first == subid && iter->second == layer)
0231           return false;
0232       }
0233       /*
0234         for (iNames = vSeedLayerNames.begin(); iNames != vSeedLayerNames.end(); iNames++){
0235         if (*iNames == name) return false;
0236         }
0237       */
0238     }
0239     //vSeedLayerNames.push_back(iHits->seedinglayer().name());
0240     vSubdetLayer.push_back(std::make_pair(subid, layer));
0241   }
0242   return true;
0243 }
0244
0245 bool CtfSpecialSeedGenerator::postCheck(const TrajectorySeed& seed) {
0246   if (!useScintillatorsConstraint)
0247     return true;
0248
0249   PTrajectoryStateOnDet pstate = seed.startingState();
0250   TrajectoryStateOnSurface theTSOS = trajectoryStateTransform::transientState(
0251       pstate, &(theTracker->idToDet(DetId(pstate.detId()))->surface()), &(*theMagfield));
0252   const FreeTrajectoryState* state = theTSOS.freeState();
0253   StraightLinePlaneCrossing planeCrossingLower(
0254       Basic3DVector<float>(state->position()), Basic3DVector<float>(state->momentum()), alongMomentum);
0255   StraightLinePlaneCrossing planeCrossingUpper(
0256       Basic3DVector<float>(state->position()), Basic3DVector<float>(state->momentum()), oppositeToMomentum);
0257   std::pair<bool, StraightLinePlaneCrossing::PositionType> positionUpper =
0258       planeCrossingUpper.position(*upperScintillator);
0259   std::pair<bool, StraightLinePlaneCrossing::PositionType> positionLower =
0260       planeCrossingLower.position(*lowerScintillator);
0261   if (!(positionUpper.first && positionLower.first)) {
0262     edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection") << "Scintillator plane not crossed";
0263     return false;
0264   }
0265   LocalPoint positionUpperLocal = upperScintillator->toLocal((GlobalPoint)(positionUpper.second));
0266   LocalPoint positionLowerLocal = lowerScintillator->toLocal((GlobalPoint)(positionLower.second));
0267   if (upperScintillator->bounds().inside(positionUpperLocal) &&
0268       lowerScintillator->bounds().inside(positionLowerLocal)) {
0269     edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection")
0270         << "position on Upper scintillator " << positionUpper.second;
0271     edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection")
0272         << "position on Lower scintillator " << positionLower.second;
0273
0274     return true;
0275   }
0276   edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection")
0277       << "scintillator not crossed in bounds: position on Upper scintillator " << positionUpper.second
0278       << " position on Lower scintillator " << positionLower.second;
0279   return false;
0280 }