SpecialSeedGenerators/src/CtfSpecialSeedGenerator.cc

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