MuonSeedGenerator/plugins/SETMuonSeedProducer.cc

0001 /** \class SETMuonSeedProducer
0002     I. Bloch, E. James, S. Stoynev
0003  */
0004
0005 #include "RecoMuon/MuonSeedGenerator/plugins/SETMuonSeedProducer.h"
0006
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0011 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0012 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0013 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0014
0015 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0016 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0017
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019
0020 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0021 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0023 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0024 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0025 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0026 #include "TMath.h"
0027
0028 using namespace edm;
0029 using namespace std;
0030
0031 SETMuonSeedProducer::SETMuonSeedProducer(const ParameterSet& parameterSet)
0032     : theSeedFinder(parameterSet), theBeamSpotTag(parameterSet.getParameter<edm::InputTag>("beamSpotTag")) {
0033   edm::ConsumesCollector iC = consumesCollector();
0034   thePatternRecognition = new SETPatternRecognition(parameterSet, iC);
0035
0036   const string metname = "Muon|RecoMuon|SETMuonSeedSeed";
0037   //std::cout<<" The SET SEED producer started."<<std::endl;
0038
0039   ParameterSet serviceParameters = parameterSet.getParameter<ParameterSet>("ServiceParameters");
0040   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0041   thePatternRecognition->setServiceProxy(theService);
0042   theSeedFinder.setServiceProxy(theService);
0043   // Parameter set for the Builder
0044   ParameterSet trajectoryBuilderParameters = parameterSet.getParameter<ParameterSet>("SETTrajBuilderParameters");
0045
0046   LogTrace(metname) << "constructor called" << endl;
0047
0048   apply_prePruning = trajectoryBuilderParameters.getParameter<bool>("Apply_prePruning");
0049
0050   useSegmentsInTrajectory = trajectoryBuilderParameters.getParameter<bool>("UseSegmentsInTrajectory");
0051
0052   // The inward-outward fitter (starts from seed state)
0053   ParameterSet filterPSet = trajectoryBuilderParameters.getParameter<ParameterSet>("FilterParameters");
0054   filterPSet.addUntrackedParameter("UseSegmentsInTrajectory", useSegmentsInTrajectory);
0055   theFilter = new SETFilter(filterPSet, theService);
0056
0057   //----
0058
0059   beamspotToken = consumes<reco::BeamSpot>(theBeamSpotTag);
0060   produces<TrajectorySeedCollection>();
0061 }
0062
0063 SETMuonSeedProducer::~SETMuonSeedProducer() {
0064   LogTrace("Muon|RecoMuon|SETMuonSeedProducer") << "SETMuonSeedProducer destructor called" << endl;
0065
0066   if (theFilter)
0067     delete theFilter;
0068   if (theService)
0069     delete theService;
0070   if (thePatternRecognition)
0071     delete thePatternRecognition;
0072 }
0073
0074 void SETMuonSeedProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0075   //std::cout<<" start producing..."<<std::endl;
0076   const string metname = "Muon|RecoMuon|SETMuonSeedSeed";
0077
0078   MuonPatternRecoDumper debug;
0079
0080   //Get the CSC Geometry :
0081   theService->update(eventSetup);
0082
0083   auto output = std::make_unique<TrajectorySeedCollection>();
0084
0085   Handle<View<TrajectorySeed> > seeds;
0086
0087   setEvent(event);
0088
0089   reco::BeamSpot beamSpot;
0090   edm::Handle<reco::BeamSpot> beamSpotHandle;
0091   event.getByToken(beamspotToken, beamSpotHandle);
0092   if (beamSpotHandle.isValid()) {
0093     beamSpot = *beamSpotHandle;
0094
0095   } else {
0096     edm::LogInfo("MuonSeedGenerator") << "No beam spot available from EventSetup \n";
0097   }
0098
0099   // make it a vector so we can subtract it from position vectors
0100   GlobalVector gv(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0101   theSeedFinder.setBeamSpot(gv);
0102
0103   bool fwFitFailed = true;
0104
0105   std::vector<SeedCandidate> seedCandidates_AllChosen;
0106   std::vector<MuonRecHitContainer> MuonRecHitContainer_clusters;
0107   //---- this is "clustering"; later a trajectory can not use hits from different clusters
0108   thePatternRecognition->produce(event, eventSetup, MuonRecHitContainer_clusters);
0109
0110   //std::cout<<"We have formed "<<MuonRecHitContainer_clusters.size()<<" clusters"<<std::endl;
0111   //---- for each cluster,
0112   for (unsigned int cluster = 0; cluster < MuonRecHitContainer_clusters.size(); ++cluster) {
0113     //std::cout<<" This is cluster number : "<<cluster<<std::endl;
0114     std::vector<SeedCandidate> seedCandidates_inCluster;
0115     //---- group hits in detector layers (if in same layer); the idea is that
0116     //---- some hits could not belong to a track simultaneously - these will be in a
0117     //---- group; two hits from one and the same group will not go to the same track
0118     std::vector<MuonRecHitContainer> MuonRecHitContainer_perLayer =
0119         theSeedFinder.sortByLayer(MuonRecHitContainer_clusters[cluster]);
0120     //---- Add protection against huge memory consumption
0121     //---- Delete busy layers if needed (due to combinatorics)
0122     theSeedFinder.limitCombinatorics(MuonRecHitContainer_perLayer);
0123     //---- build all possible combinations (valid sets)
0124     std::vector<MuonRecHitContainer> allValidSets = theSeedFinder.findAllValidSets(MuonRecHitContainer_perLayer);
0125     if (apply_prePruning) {
0126       //---- remove "wild" segments from the combination
0127       theSeedFinder.validSetsPrePruning(allValidSets);
0128     }
0129
0130     //---- build the appropriate output: seedCandidates_inCluster
0131     //---- if too many (?) valid sets in a cluster - skip it
0132     if (allValidSets.size() < 500) {  // hardcoded - remove it
0133       seedCandidates_inCluster = theSeedFinder.fillSeedCandidates(allValidSets);
0134     }
0135     //---- find the best valid combinations using simple (no propagation errors) chi2-fit
0136     //std::cout<<"  Found "<<seedCandidates_inCluster.size()<<" valid sets in the current cluster."<<std::endl;
0137     if (!seedCandidates_inCluster.empty()) {
0138       //---- this is the forward fitter (segments); choose which of the SETs in a cluster to be considered further
0139       std::vector<SeedCandidate> bestSets_inCluster;
0140       fwFitFailed = !(filter()->fwfit_SET(seedCandidates_inCluster, bestSets_inCluster));
0141
0142       //---- has the fit failed? continue to the next cluster instead of returning the empty trajectoryContainer and stop the loop IBL 080903
0143       if (fwFitFailed) {
0144         //std::cout<<"  fwfit_SET failed!"<<std::endl;
0145         continue;
0146       }
0147       for (unsigned int iSet = 0; iSet < bestSets_inCluster.size(); ++iSet) {
0148         seedCandidates_AllChosen.push_back(bestSets_inCluster[iSet]);
0149       }
0150     }
0151   }
0152   //---- loop over all the SETs candidates
0153   for (unsigned int iMuon = 0; iMuon < seedCandidates_AllChosen.size(); ++iMuon) {
0154     //std::cout<<" chosen iMuon = "<<iMuon<<std::endl;
0155     Trajectory::DataContainer finalCandidate;
0156     SeedCandidate* aFinalSet = &(seedCandidates_AllChosen[iMuon]);
0157     fwFitFailed = !(filter()->buildTrajectoryMeasurements(aFinalSet, finalCandidate));
0158     if (fwFitFailed) {
0159       //std::cout<<"  buildTrajectoryMeasurements failed!"<<std::endl;
0160       continue;
0161     }
0162     //---- are there measurements (or detLayers) used at all?
0163     if (!filter()->layers().empty())
0164       LogTrace(metname) << debug.dumpLayer(filter()->lastDetLayer());
0165     else {
0166       continue;
0167     }
0168     //std::cout<<"  chambers used - all : "<<filter()->getTotalChamberUsed()<<", DT : "<<filter()->getDTChamberUsed()<<
0169     //", CSC : "<<filter()->getCSCChamberUsed()<<", RPC : "<<filter()->getRPCChamberUsed()<<std::endl;
0170     //---- ask for some "reasonable" conditions to build a STA muon;
0171     //---- (totalChambers >= 2, dtChambers + cscChambers >0)
0172     if (filter()->goodState()) {
0173       TransientTrackingRecHit::ConstRecHitContainer hitContainer;
0174       TrajectoryStateOnSurface firstTSOS;
0175       bool conversionPassed = false;
0176       if (!useSegmentsInTrajectory) {
0177         //---- transforms set of segment measurements to a set of rechit measurements
0178         conversionPassed = filter()->transform(finalCandidate, hitContainer, firstTSOS);
0179       } else {
0180         //---- transforms set of segment measurements to a set of segment measurements
0181         conversionPassed = filter()->transformLight(finalCandidate, hitContainer, firstTSOS);
0182       }
0183       if (conversionPassed && !finalCandidate.empty() && !hitContainer.empty()) {
0184         //---- doesn't work...
0185         //output->push_back( theSeedFinder.makeSeed(firstTSOS, hitContainer) );
0186
0187         edm::OwnVector<TrackingRecHit> recHitsContainer;
0188         for (unsigned int iHit = 0; iHit < hitContainer.size(); ++iHit) {
0189           recHitsContainer.push_back(hitContainer.at(iHit)->hit()->clone());
0190         }
0191         PropagationDirection dir = oppositeToMomentum;
0192         if (useSegmentsInTrajectory) {
0193           dir = alongMomentum;  // why forward (for rechits) later?
0194         }
0195
0196         PTrajectoryStateOnDet seedTSOS =
0197             trajectoryStateTransform::persistentState(firstTSOS, hitContainer.at(0)->geographicalId().rawId());
0198         TrajectorySeed seed(seedTSOS, recHitsContainer, dir);
0199
0200         //MuonPatternRecoDumper debug;
0201         //std::cout<<"  firstTSOS (not IP) = "<<debug.dumpTSOS(firstTSOS)<<std::endl;
0202         //std::cout<<"  hits(from range) = "<<range.second-range.first<<" hits (from size) = "<<hitContainer.size()<<std::endl;
0203         //for(unsigned int iRH=0;iRH<hitContainer.size();++iRH){
0204         //std::cout<<"  RH = "<<iRH+1<<" globPos = "<<hitContainer.at(iRH)->globalPosition()<<std::endl;
0205         //}
0206         output->push_back(seed);
0207       } else {
0208         //std::cout<<" Transformation from TrajectoryMeasurements to RecHitContainer faild - skip "<<std::endl;
0209         continue;
0210       }
0211     } else {
0212       //std::cout<<" Not enough (as defined) measurements to build trajectory - skip"<<std::endl;
0213       continue;
0214     }
0215   }
0216   event.put(std::move(output));
0217   theFilter->reset();
0218 }
0219
0220 //
0221 void SETMuonSeedProducer::setEvent(const edm::Event& event) { theFilter->setEvent(event); }