ConversionSeedGenerators/plugins/PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo.cc

0001 #include "PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo.h"
0002 #include "Quad.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 // ClusterShapeIncludes:::
0005 #include "RecoTracker/TkTrackingRegions/interface/OrderedHitsGenerator.h"
0006 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0007 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
0008 #include "RecoTracker/TkSeedGenerator/interface/SeedCreator.h"
0009 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitorFactory.h"
0010 #include "RecoTracker/TkSeedGenerator/interface/SeedCreatorFactory.h"
0011 /*
0012 To Do:
0013
0014 assign the parameters to some data member to avoid search at every event
0015
0016  */
0017
0018 //#define debugTSPFSLA
0019 //#define mydebug_knuenz
0020
0021 PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo::PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo(
0022     const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0023     : _conf(conf),
0024       seedCollection(nullptr),
0025       theClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet"), iC),
0026       QuadCutPSet(conf.getParameter<edm::ParameterSet>("QuadCutPSet")),
0027       theSilentOnClusterCheck(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet")
0028                                   .getUntrackedParameter<bool>("silentClusterCheck", false)),
0029       theHitsGenerator(new CombinedHitQuadrupletGeneratorForPhotonConversion(
0030           conf.getParameter<edm::ParameterSet>("OrderedHitsFactoryPSet"), iC)),
0031       theSeedCreator(
0032           new SeedForPhotonConversionFromQuadruplets(conf.getParameter<edm::ParameterSet>("SeedCreatorPSet"),
0033                                                      iC,
0034                                                      conf.getParameter<edm::ParameterSet>("SeedComparitorPSet"))),
0035       theRegionProducer(
0036           new GlobalTrackingRegionProducerFromBeamSpot(conf.getParameter<edm::ParameterSet>("RegionFactoryPSet"), iC)) {
0037   token_vertex = iC.consumes<reco::VertexCollection>(_conf.getParameter<edm::InputTag>("primaryVerticesTag"));
0038 }
0039
0040 PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo::~PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo() {
0041 }
0042
0043 void PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo::analyze(const edm::Event& event,
0044                                                                         const edm::EventSetup& setup) {
0045   myEsetup = &setup;
0046   myEvent = &event;
0047
0048   if (seedCollection != nullptr)
0049     delete seedCollection;
0050
0051   seedCollection = new TrajectorySeedCollection();
0052
0053   size_t clustsOrZero = theClusterCheck.tooManyClusters(event);
0054   if (clustsOrZero) {
0055     if (!theSilentOnClusterCheck)
0056       edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
0057     return;
0058   }
0059
0060   // Why is the regions variable (and theRegionProducer) needed at all?
0061   regions = theRegionProducer->regions(event, setup);
0062
0063   event.getByToken(token_vertex, vertexHandle);
0064   if (!vertexHandle.isValid()) {
0065     edm::LogError("PhotonConversionFinderFromTracks")
0066         << "Error! Can't get the product primary Vertex Collection "
0067         << _conf.getParameter<edm::InputTag>("primaryVerticesTag") << "\n";
0068     return;
0069   }
0070
0071   //Do the analysis
0072   loop();
0073
0074 #ifdef mydebug_knuenz
0075   std::cout << "Running PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo" << std::endl;
0076 #endif
0077
0078 #ifdef debugTSPFSLA
0079   std::stringstream ss;
0080   ss.str("");
0081   ss << "\n++++++++++++++++++\n";
0082   ss << "seed collection size " << seedCollection->size();
0083   for (auto const& tjS : *seedCollection) {
0084     po.print(ss, tjS);
0085   }
0086   edm::LogInfo("debugTrajSeedFromQuadruplets") << ss.str();
0087   //-------------------------------------------------
0088 #endif
0089 }
0090
0091 void PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo::loop() {
0092   ss.str("");
0093
0094   float ptMin = 0.1;
0095
0096   for (auto const& primaryVertex : *vertexHandle) {
0097     //FIXME currnetly using the first primaryVertex, should loop on the promaryVertexes
0098     GlobalTrackingRegion region(
0099         ptMin,
0100         GlobalPoint(primaryVertex.position().x(), primaryVertex.position().y(), primaryVertex.position().z()),
0101         primaryVertex.xError() * 10,
0102         primaryVertex.zError() * 10,
0103         true);
0104
0105 #ifdef debugTSPFSLA
0106     ss << "[PrintRegion] " << region.print() << std::endl;
0107 #endif
0108
0109     inspect(region);
0110   }
0111 #ifdef debugTSPFSLA
0112   edm::LogInfo("debugTrajSeedFromQuadruplets") << ss.str();
0113 #endif
0114 }
0115
0116 bool PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo::inspect(const TrackingRegion& region) {
0117   const OrderedSeedingHits& hitss = theHitsGenerator->run(region, *myEvent, *myEsetup);
0118
0119   unsigned int nHitss = hitss.size();
0120
0121 #ifdef debugTSPFSLA
0122   ss << "\n nHitss " << nHitss << "\n";
0123 #endif
0124
0125   // don't do multiple reserves in the case of multiple regions: it would make things even worse
0126   // as it will cause N re-allocations instead of the normal log(N)/log(2)
0127   if (seedCollection->empty())
0128     seedCollection->reserve(nHitss / 2);
0129
0130   unsigned int iHits = 0, jHits = 1;
0131
0132   //
0133   // Trivial arbitration
0134   //
0135   // Vector to store old quad values
0136   std::vector<Quad> quadVector;
0137
0138   for (; iHits < nHitss && jHits < nHitss; iHits += 2, jHits += 2) {
0139 #ifdef debugTSPFSLA
0140     //    ss << "\n iHits " << iHits << " " << jHits << "\n";
0141 #endif
0142     //phits is the pair of hits for the positron
0143     const SeedingHitSet& phits = hitss[iHits];
0144     //mhits is the pair of hits for the electron
0145     const SeedingHitSet& mhits = hitss[jHits];
0146
0147     try {
0148       //FIXME (modify the interface of the seed generator if needed)
0149       //passing the region, that is centered around the primary vertex
0150       theSeedCreator->trajectorySeed(
0151           *seedCollection, phits, mhits, region, *myEvent, *myEsetup, ss, quadVector, QuadCutPSet);
0152     } catch (cms::Exception& er) {
0153       edm::LogError("SeedingConversion") << " Problem in the Quad Seed creator " << er.what() << std::endl;
0154     } catch (std::exception& er) {
0155       edm::LogError("SeedingConversion") << " Problem in the Quad Seed creator " << er.what() << std::endl;
0156     }
0157   }
0158   quadVector.clear();
0159   return true;
0160 }