SETMuonSeedProducer.cc

Go to the documentation of this file.

#include "RecoMuon/MuonSeedGenerator/plugins/SETMuonSeedProducer.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"

#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
#include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "TMath.h"

using namespace edm;
using namespace std;

SETMuonSeedProducer::SETMuonSeedProducer(const ParameterSet& parameterSet)
: theSeedFinder(parameterSet),
  theBeamSpotTag(parameterSet.getParameter<edm::InputTag>("beamSpotTag"))
{
  
  edm::ConsumesCollector iC =  consumesCollector();
  thePatternRecognition = new SETPatternRecognition(parameterSet,iC);
  
  const string metname = "Muon|RecoMuon|SETMuonSeedSeed";  
  //std::cout<<" The SET SEED producer started."<<std::endl;

  ParameterSet serviceParameters = parameterSet.getParameter<ParameterSet>("ServiceParameters");
  theService        = new MuonServiceProxy(serviceParameters);
  thePatternRecognition->setServiceProxy(theService);
  theSeedFinder.setServiceProxy(theService);
  // Parameter set for the Builder
  ParameterSet trajectoryBuilderParameters = parameterSet.getParameter<ParameterSet>("SETTrajBuilderParameters");

  LogTrace(metname) << "constructor called" << endl;
  
  apply_prePruning = trajectoryBuilderParameters.getParameter<bool>("Apply_prePruning");

  useSegmentsInTrajectory = trajectoryBuilderParameters.getParameter<bool>("UseSegmentsInTrajectory");

  // The inward-outward fitter (starts from seed state)
  ParameterSet filterPSet = trajectoryBuilderParameters.getParameter<ParameterSet>("FilterParameters");
  filterPSet.addUntrackedParameter("UseSegmentsInTrajectory", useSegmentsInTrajectory);
  theFilter = new SETFilter(filterPSet,theService);

  //----


  beamspotToken = consumes<reco::BeamSpot>(theBeamSpotTag);
  produces<TrajectorySeedCollection>();

} 

SETMuonSeedProducer::~SETMuonSeedProducer(){

  LogTrace("Muon|RecoMuon|SETMuonSeedProducer") 
    << "SETMuonSeedProducer destructor called" << endl;
  
  if(theFilter) delete theFilter;
   if (theService) delete theService;
   if(thePatternRecognition) delete  thePatternRecognition;

}

void SETMuonSeedProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup){
  //std::cout<<" start producing..."<<std::endl;  
  const string metname = "Muon|RecoMuon|SETMuonSeedSeed";  

  MuonPatternRecoDumper debug;

  //Get the CSC Geometry :
  theService->update(eventSetup);

  std::auto_ptr<TrajectorySeedCollection> output(new TrajectorySeedCollection());
  
  Handle<View<TrajectorySeed> > seeds; 

  setEvent(event);

  reco::BeamSpot beamSpot;
  edm::Handle<reco::BeamSpot> beamSpotHandle;
  event.getByToken(beamspotToken, beamSpotHandle);
  if ( beamSpotHandle.isValid() )
  {
    beamSpot = *beamSpotHandle;

  } else
  {
    edm::LogInfo("MuonSeedGenerator")
      << "No beam spot available from EventSetup \n";
  }

  // make it a vector so we can subtract it from position vectors
  GlobalVector gv(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
  theSeedFinder.setBeamSpot(gv);

  bool fwFitFailed = true;

  std::vector <SeedCandidate> seedCandidates_AllChosen;
  std::vector< MuonRecHitContainer > MuonRecHitContainer_clusters;
  //---- this is "clustering"; later a trajectory can not use hits from different clusters
  thePatternRecognition->produce(event, eventSetup, MuonRecHitContainer_clusters);

  //std::cout<<"We have formed "<<MuonRecHitContainer_clusters.size()<<" clusters"<<std::endl;
  //---- for each cluster,
  for(unsigned int cluster = 0; cluster < MuonRecHitContainer_clusters.size(); ++cluster) {
    //std::cout<<" This is cluster number : "<<cluster<<std::endl;
    std::vector <SeedCandidate> seedCandidates_inCluster; 
    //---- group hits in detector layers (if in same layer); the idea is that
    //---- some hits could not belong to a track simultaneously - these will be in a 
    //---- group; two hits from one and the same group will not go to the same track 
    std::vector< MuonRecHitContainer > MuonRecHitContainer_perLayer = theSeedFinder.sortByLayer(MuonRecHitContainer_clusters[cluster]);
    //---- Add protection against huge memory consumption
    //---- Delete busy layers if needed (due to combinatorics)
    theSeedFinder.limitCombinatorics(MuonRecHitContainer_perLayer);
    //---- build all possible combinations (valid sets)
    std::vector <MuonRecHitContainer> allValidSets = theSeedFinder.findAllValidSets(MuonRecHitContainer_perLayer);
    if(apply_prePruning){
      //---- remove "wild" segments from the combination
      theSeedFinder.validSetsPrePruning(allValidSets);
    }
    
    //---- build the appropriate output: seedCandidates_inCluster
    //---- if too many (?) valid sets in a cluster - skip it 
    if(allValidSets.size()<500){// hardcoded - remove it
      seedCandidates_inCluster = theSeedFinder.fillSeedCandidates(allValidSets);
    }
    //---- find the best valid combinations using simple (no propagation errors) chi2-fit 
    //std::cout<<"  Found "<<seedCandidates_inCluster.size()<<" valid sets in the current cluster."<<std::endl;
    if(!seedCandidates_inCluster.empty()){
      //---- this is the forward fitter (segments); choose which of the SETs in a cluster to be considered further
      std::vector < SeedCandidate> bestSets_inCluster;
      fwFitFailed = !(filter()->fwfit_SET(seedCandidates_inCluster, bestSets_inCluster));
      
      //---- has the fit failed? continue to the next cluster instead of returning the empty trajectoryContainer and stop the loop IBL 080903
      if(fwFitFailed){
    //std::cout<<"  fwfit_SET failed!"<<std::endl;
    continue;
      }  
      for(unsigned int iSet = 0; iSet<bestSets_inCluster.size();++iSet){
        seedCandidates_AllChosen.push_back(bestSets_inCluster[iSet]);
      }
    }
  }
  //---- loop over all the SETs candidates
  for(unsigned int iMuon = 0;iMuon<seedCandidates_AllChosen.size();++iMuon){
    //std::cout<<" chosen iMuon = "<<iMuon<<std::endl;
    Trajectory::DataContainer finalCandidate;
    SeedCandidate * aFinalSet = &(seedCandidates_AllChosen[iMuon]);
    fwFitFailed = !(filter()->buildTrajectoryMeasurements(aFinalSet, finalCandidate));
    if(fwFitFailed ){
      //std::cout<<"  buildTrajectoryMeasurements failed!"<<std::endl;
      continue;
    }
    //---- are there measurements (or detLayers) used at all?
    if( filter()->layers().size() )
      LogTrace(metname) << debug.dumpLayer( filter()->lastDetLayer());
    else {
      continue;
    }
    //std::cout<<"  chambers used - all : "<<filter()->getTotalChamberUsed()<<", DT : "<<filter()->getDTChamberUsed()<<
    //", CSC : "<<filter()->getCSCChamberUsed()<<", RPC : "<<filter()->getRPCChamberUsed()<<std::endl;
    //---- ask for some "reasonable" conditions to build a STA muon; 
    //---- (totalChambers >= 2, dtChambers + cscChambers >0)
    if (filter()->goodState()) {
      TransientTrackingRecHit::ConstRecHitContainer hitContainer;
      TrajectoryStateOnSurface firstTSOS;
      bool conversionPassed = false;
      if(!useSegmentsInTrajectory){
    //---- transforms set of segment measurements to a set of rechit measurements
    conversionPassed = filter()->transform(finalCandidate, hitContainer, firstTSOS);
      }
      else{
    //---- transforms set of segment measurements to a set of segment measurements
    conversionPassed = filter()->transformLight(finalCandidate, hitContainer, firstTSOS);
      }
      if ( conversionPassed && !finalCandidate.empty() && !hitContainer.empty()) {
    //---- doesn't work... 
    //output->push_back( theSeedFinder.makeSeed(firstTSOS, hitContainer) );
    
    edm::OwnVector<TrackingRecHit> recHitsContainer;
    for(unsigned int iHit = 0;iHit < hitContainer.size();++iHit){
      recHitsContainer.push_back(hitContainer.at(iHit)->hit()->clone());
    }
    PropagationDirection dir = oppositeToMomentum;
    if(useSegmentsInTrajectory){
      dir = alongMomentum;// why forward (for rechits) later?
    }
    
    PTrajectoryStateOnDet seedTSOS =
     trajectoryStateTransform::persistentState( firstTSOS, hitContainer.at(0)->geographicalId().rawId());
    TrajectorySeed seed(seedTSOS,recHitsContainer,dir);
    
    //MuonPatternRecoDumper debug;
    //std::cout<<"  firstTSOS (not IP) = "<<debug.dumpTSOS(firstTSOS)<<std::endl;
    //std::cout<<"  hits(from range) = "<<range.second-range.first<<" hits (from size) = "<<hitContainer.size()<<std::endl;
    //for(unsigned int iRH=0;iRH<hitContainer.size();++iRH){
    //std::cout<<"  RH = "<<iRH+1<<" globPos = "<<hitContainer.at(iRH)->globalPosition()<<std::endl;
    //}
    output->push_back(seed);
      }
      else{
    //std::cout<<" Transformation from TrajectoryMeasurements to RecHitContainer faild - skip "<<std::endl;
    continue;
      }
    }else{
      //std::cout<<" Not enough (as defined) measurements to build trajectory - skip"<<std::endl;
      continue;
    }
  }
  event.put(output);
  theFilter->reset();
}

//
void SETMuonSeedProducer::setEvent(const edm::Event& event){
  theFilter->setEvent(event);
}