Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:27:14

0001 /** \class StandAloneTrajectoryBuilder
0002  *  Concrete class for the STA Muon reco 
0003  *
0004  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0005  *  \author Stefano Lacaprara - INFN Legnaro
0006  *  \author D. Trocino - INFN Torino <daniele.trocino@to.infn.it>
0007  *
0008  *  Modified by C. Calabria
0009  *  Modified by D. Nash
0010  */
0011 
0012 #include "RecoMuon/StandAloneTrackFinder/interface/StandAloneTrajectoryBuilder.h"
0013 
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0016 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0018 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0019 
0020 #include "RecoMuon/StandAloneTrackFinder/interface/StandAloneMuonFilter.h"
0021 #include "RecoMuon/StandAloneTrackFinder/interface/StandAloneMuonBackwardFilter.h"
0022 #include "RecoMuon/StandAloneTrackFinder/interface/StandAloneMuonRefitter.h"
0023 
0024 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0025 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0026 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0027 
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 
0030 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0032 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0033 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0034 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0035 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0036 #include "TrackingTools/TrackRefitter/interface/SeedTransformer.h"
0037 
0038 using namespace edm;
0039 using namespace std;
0040 
0041 StandAloneMuonTrajectoryBuilder::StandAloneMuonTrajectoryBuilder(const ParameterSet& par,
0042                                                                  const MuonServiceProxy* service,
0043                                                                  edm::ConsumesCollector& iC)
0044     : theService(service) {
0045   const std::string metname = "Muon|RecoMuon|StandAloneMuonTrajectoryBuilder";
0046 
0047   LogTrace(metname) << "constructor called" << endl;
0048 
0049   // The navigation type:
0050   // "Direct","Standard"
0051   theNavigationType = par.getParameter<string>("NavigationType");
0052 
0053   // The inward-outward fitter (starts from seed state)
0054   ParameterSet filterPSet = par.getParameter<ParameterSet>("FilterParameters");
0055   filterPSet.addParameter<string>("NavigationType", theNavigationType);
0056   theFilter = std::make_unique<StandAloneMuonFilter>(filterPSet, theService, iC);
0057 
0058   // Fit direction
0059   string seedPosition = par.getParameter<string>("SeedPosition");
0060 
0061   if (seedPosition == "in")
0062     theSeedPosition = recoMuon::in;
0063   else if (seedPosition == "out")
0064     theSeedPosition = recoMuon::out;
0065   else
0066     throw cms::Exception("StandAloneMuonFilter constructor")
0067         << "Wrong seed position chosen in StandAloneMuonFilter::StandAloneMuonFilter ParameterSet"
0068         << "\n"
0069         << "Possible choices are:"
0070         << "\n"
0071         << "SeedPosition = in or SeedPosition = out";
0072 
0073   // Propagator for the seed extrapolation
0074   theSeedPropagatorName = par.getParameter<string>("SeedPropagator");
0075 
0076   // Disable/Enable the backward filter
0077   doBackwardFilter = par.getParameter<bool>("DoBackwardFilter");
0078 
0079   // Disable/Enable the refit of the trajectory
0080   doRefit = par.getParameter<bool>("DoRefit");
0081 
0082   // Disable/Enable the refit of the trajectory seed
0083   doSeedRefit = par.getParameter<bool>("DoSeedRefit");
0084 
0085   if (doBackwardFilter) {
0086     // The outward-inward fitter (starts from theFilter outermost state)
0087     ParameterSet bwFilterPSet = par.getParameter<ParameterSet>("BWFilterParameters");
0088     //  theBWFilter = new StandAloneMuonBackwardFilter(bwFilterPSet,theService); // FIXME
0089     bwFilterPSet.addParameter<string>("NavigationType", theNavigationType);
0090     theBWFilter = std::make_unique<StandAloneMuonFilter>(bwFilterPSet, theService, iC);
0091 
0092     theBWSeedType = bwFilterPSet.getParameter<string>("BWSeedType");
0093   }
0094 
0095   if (doRefit) {
0096     // The outward-inward fitter (starts from theBWFilter innermost state)
0097     ParameterSet refitterPSet = par.getParameter<ParameterSet>("RefitterParameters");
0098     theRefitter = std::make_unique<StandAloneMuonRefitter>(refitterPSet, iC, theService);
0099   }
0100 
0101   // The seed transformer (used to refit the seed and get the seed transient state)
0102   //  ParameterSet seedTransformerPSet = par.getParameter<ParameterSet>("SeedTransformerParameters");
0103   ParameterSet seedTransformerParameters = par.getParameter<ParameterSet>("SeedTransformerParameters");
0104   theSeedTransformer = std::make_unique<SeedTransformer>(seedTransformerParameters, iC);
0105 }
0106 
0107 void StandAloneMuonTrajectoryBuilder::setEvent(const edm::Event& event) {
0108   theFilter->setEvent(event);
0109   if (doBackwardFilter)
0110     theBWFilter->setEvent(event);
0111 }
0112 
0113 StandAloneMuonTrajectoryBuilder::~StandAloneMuonTrajectoryBuilder() {
0114   LogTrace("Muon|RecoMuon|StandAloneMuonTrajectoryBuilder")
0115       << "StandAloneMuonTrajectoryBuilder destructor called" << endl;
0116 }
0117 
0118 namespace {
0119   struct Resetter {
0120     StandAloneMuonFilter* mf;
0121     explicit Resetter(StandAloneMuonFilter* imf) : mf(imf) {}
0122     ~Resetter() {
0123       if (mf)
0124         mf->reset();
0125     }
0126   };
0127 }  // namespace
0128 
0129 MuonTrajectoryBuilder::TrajectoryContainer StandAloneMuonTrajectoryBuilder::trajectories(const TrajectorySeed& seed) {
0130   Resetter fwReset(filter());
0131   Resetter bwReset(bwfilter());
0132 
0133   const std::string metname = "Muon|RecoMuon|StandAloneMuonTrajectoryBuilder";
0134   MuonPatternRecoDumper debug;
0135 
0136   // Set the services for the seed transformer
0137   theSeedTransformer->setServices(theService->eventSetup());
0138 
0139   // the trajectory container. In principle starting from one seed we can
0140   // obtain more than one trajectory. TODO: this feature is not yet implemented!
0141   TrajectoryContainer trajectoryContainer;
0142 
0143   PropagationDirection fwDirection = (theSeedPosition == recoMuon::in) ? alongMomentum : oppositeToMomentum;
0144   Trajectory trajectoryFW(seed, fwDirection);
0145 
0146   TrajectoryStateOnSurface lastTSOS;
0147   DetId lastDetId;
0148 
0149   vector<Trajectory> seedTrajectories;
0150 
0151   if (doSeedRefit) {
0152     seedTrajectories = theSeedTransformer->seedTransform(seed);
0153     if (!seedTrajectories.empty()) {
0154       TrajectoryMeasurement lastTM(seedTrajectories.front().lastMeasurement());
0155       lastTSOS = lastTM.updatedState();
0156       lastDetId = lastTM.recHit()->geographicalId();
0157     }
0158   }
0159 
0160   if (!doSeedRefit || seedTrajectories.empty()) {
0161     lastTSOS = theSeedTransformer->seedTransientState(seed);
0162     lastDetId = seed.startingState().detId();
0163   }
0164 
0165   LogTrace(metname) << "Trajectory State on Surface before the extrapolation" << endl;
0166   LogTrace(metname) << debug.dumpTSOS(lastTSOS);
0167 
0168   // Segment layer
0169   LogTrace(metname) << "The RecSegment relies on: " << endl;
0170   LogTrace(metname) << debug.dumpMuonId(lastDetId);
0171 
0172   DetLayerWithState inputFromSeed = propagateTheSeedTSOS(lastTSOS, lastDetId);
0173 
0174   // refine the FTS given by the seed
0175 
0176   // the trajectory is filled in the refitter::refit
0177   filter()->refit(inputFromSeed.second, inputFromSeed.first, trajectoryFW);
0178 
0179   // "0th order" check...
0180   if (trajectoryFW.empty()) {
0181     LogTrace(metname) << "Forward trajectory EMPTY! No trajectory will be loaded!" << endl;
0182     return trajectoryContainer;
0183   }
0184 
0185   // Get the last TSOS
0186   //  TrajectoryStateOnSurface tsosAfterRefit = filter()->lastUpdatedTSOS();     // this is the last UPDATED TSOS
0187   TrajectoryStateOnSurface tsosAfterRefit = filter()->lastCompatibleTSOS();  // this is the last COMPATIBLE TSOS
0188 
0189   LogTrace(metname) << "StandAloneMuonTrajectoryBuilder filter output " << endl;
0190   LogTrace(metname) << debug.dumpTSOS(tsosAfterRefit);
0191 
0192   /*
0193   // -- 1st attempt
0194   if( filter()->isCompatibilitySatisfied() ) {
0195     if( filter()->layers().size() )   //  OR   if( filter()->goodState() ) ???  Maybe when only RPC hits are used...
0196       LogTrace(metname) << debug.dumpLayer( filter()->lastDetLayer() );
0197     else {
0198       LogTrace(metname) << "Compatibility satisfied, but all RecHits are invalid! A trajectory with only invalid hits will be loaded!" << endl;
0199       trajectoryContainer.push_back(new Trajectory(trajectoryFW));
0200       return trajectoryContainer;
0201     }
0202   }
0203   else {
0204     LogTrace(metname) << "Compatibility NOT satisfied after Forward filter! No trajectory will be loaded!" << endl;
0205     LogTrace(metname) << "Total chambers: " << filter()->getTotalCompatibleChambers() << "; DT: " 
0206               << filter()->getDTCompatibleChambers() << "; CSC: " << filter()->getCSCCompatibleChambers() << endl;
0207     return trajectoryContainer; 
0208   }
0209   // -- end 1st attempt
0210   */
0211 
0212   // -- 2nd attempt
0213   if (filter()->goodState()) {
0214     LogTrace(metname) << debug.dumpLayer(filter()->lastDetLayer());
0215   } else if (filter()->isCompatibilitySatisfied()) {
0216     int foundValidRh = trajectoryFW.foundHits();  // number of valid hits
0217     LogTrace(metname) << "Compatibility satisfied after Forward filter, but too few valid RecHits (" << foundValidRh
0218                       << ")." << endl
0219                       << "A trajectory with only invalid RecHits will be loaded!" << endl;
0220     if (!foundValidRh) {
0221       trajectoryContainer.push_back(std::make_unique<Trajectory>(trajectoryFW));
0222       return trajectoryContainer;
0223     }
0224     Trajectory defaultTraj(seed, fwDirection);
0225     filter()->createDefaultTrajectory(trajectoryFW, defaultTraj);
0226     trajectoryContainer.push_back(std::make_unique<Trajectory>(defaultTraj));
0227     return trajectoryContainer;
0228   } else {
0229     LogTrace(metname) << "Compatibility NOT satisfied after Forward filter! No trajectory will be loaded!" << endl;
0230     LogTrace(metname) << "Total compatible chambers: " << filter()->getTotalCompatibleChambers()
0231                       << ";  DT: " << filter()->getDTCompatibleChambers()
0232                       << ";  CSC: " << filter()->getCSCCompatibleChambers()
0233                       << ";  RPC: " << filter()->getRPCCompatibleChambers()
0234                       << ";  GEM: " << filter()->getGEMCompatibleChambers()
0235                       << ";  ME0: " << filter()->getME0CompatibleChambers() << endl;
0236     return trajectoryContainer;
0237   }
0238   // -- end 2nd attempt
0239 
0240   LogTrace(metname) << "Number of DT/CSC/RPC/GEM/ME0 chamber used (fw): " << filter()->getDTChamberUsed() << "/"
0241                     << filter()->getCSCChamberUsed() << "/" << filter()->getRPCChamberUsed() << "/"
0242                     << filter()->getGEMChamberUsed() << "/" << filter()->getME0ChamberUsed() << endl;
0243   LogTrace(metname) << "Momentum: " << tsosAfterRefit.freeState()->momentum();
0244 
0245   if (!doBackwardFilter) {
0246     LogTrace(metname) << "Only forward refit requested. No backward refit will be performed!" << endl;
0247 
0248     // ***** Refit of fwd step *****
0249     //    if (doRefit && !trajectoryFW.empty() && filter()->goodState()){    // CHECK!!! Can trajectoryFW really be empty at this point??? And goodState...?
0250     if (doRefit) {
0251       pair<bool, Trajectory> refitResult = refitter()->refit(trajectoryFW);
0252       if (refitResult.first) {
0253         trajectoryContainer.push_back(std::make_unique<Trajectory>(refitResult.second));
0254         LogTrace(metname) << "StandAloneMuonTrajectoryBuilder refit output " << endl;
0255         LogTrace(metname) << debug.dumpTSOS(refitResult.second.lastMeasurement().updatedState());
0256       } else
0257         trajectoryContainer.push_back(std::make_unique<Trajectory>(trajectoryFW));
0258     } else
0259       trajectoryContainer.push_back(std::make_unique<Trajectory>(trajectoryFW));
0260 
0261     LogTrace(metname) << "Trajectory saved" << endl;
0262     return trajectoryContainer;
0263   }
0264 
0265   // ***** Backward filtering *****
0266 
0267   TrajectorySeed seedForBW;
0268 
0269   if (theBWSeedType == "noSeed") {
0270   } else if (theBWSeedType == "fromFWFit") {
0271     PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(
0272         tsosAfterRefit, trajectoryFW.lastMeasurement().recHit()->geographicalId().rawId());
0273 
0274     edm::OwnVector<TrackingRecHit> recHitsContainer;
0275     PropagationDirection seedDirection = (theSeedPosition == recoMuon::in) ? oppositeToMomentum : alongMomentum;
0276     TrajectorySeed fwFit(seedTSOS, recHitsContainer, seedDirection);
0277 
0278     seedForBW = fwFit;
0279   } else if (theBWSeedType == "fromGenerator") {
0280     seedForBW = seed;
0281   } else
0282     LogWarning(metname) << "Wrong seed type for the backward filter!";
0283 
0284   PropagationDirection bwDirection = (theSeedPosition == recoMuon::in) ? oppositeToMomentum : alongMomentum;
0285   Trajectory trajectoryBW(seedForBW, bwDirection);
0286 
0287   // FIXME! under check!
0288   bwfilter()->refit(tsosAfterRefit, filter()->lastDetLayer(), trajectoryBW);
0289 
0290   // Get the last TSOS
0291   TrajectoryStateOnSurface tsosAfterBWRefit = bwfilter()->lastUpdatedTSOS();
0292 
0293   LogTrace(metname) << "StandAloneMuonTrajectoryBuilder BW filter output " << endl;
0294   LogTrace(metname) << debug.dumpTSOS(tsosAfterBWRefit);
0295 
0296   LogTrace(metname) << "Number of RecHits: " << trajectoryBW.foundHits() << "\n"
0297                     << "Number of DT/CSC/RPC/GEM/ME0 chamber used (bw): " << bwfilter()->getDTChamberUsed() << "/"
0298                     << bwfilter()->getCSCChamberUsed() << "/" << bwfilter()->getRPCChamberUsed() << "/"
0299                     << bwfilter()->getGEMChamberUsed() << "/" << bwfilter()->getME0ChamberUsed();
0300 
0301   // -- The trajectory is "good" if there are at least 2 chambers used in total and at
0302   //    least 1 is "tracking" (DT or CSC)
0303   // -- The trajectory satisfies the "compatibility" requirements if there are at least
0304   //    2 compatible chambers (not necessarily used!) in total and at
0305   //    least 1 is "tracking" (DT or CSC)
0306   // 1st attempt
0307   /*
0308   if (bwfilter()->isCompatibilitySatisfied()) {
0309     
0310     if (doRefit && !trajectoryBW.empty() && bwfilter()->goodState()){
0311       pair<bool,Trajectory> refitResult = refitter()->refit(trajectoryBW);
0312       if (refitResult.first){
0313         trajectoryContainer.push_back(new Trajectory(refitResult.second));
0314     LogTrace(metname) << "StandAloneMuonTrajectoryBuilder Refit output " << endl;
0315     LogTrace(metname) << debug.dumpTSOS(refitResult.second.lastMeasurement().updatedState());
0316       }
0317       else
0318     trajectoryContainer.push_back(new Trajectory(trajectoryBW));
0319     }
0320     else
0321       trajectoryContainer.push_back(new Trajectory(trajectoryBW));
0322     
0323     LogTrace(metname)<< "Trajectory saved" << endl;
0324     
0325   }
0326   //if the trajectory is not saved, but at least two chamber are used in the
0327   //forward filtering, try to build a new trajectory starting from the old
0328   //trajectory w/o the latest measurement and a looser chi2 cut
0329   else if ( filter()->getTotalChamberUsed() >= 2 ) {
0330     LogTrace(metname)<< "Trajectory NOT saved. Second Attempt." << endl;
0331     // FIXME:
0332     // a better choice could be: identify the poorest one, exclude it, redo
0333     // the fw and bw filtering. Or maybe redo only the bw without the excluded
0334     // measure. As first step I will port the ORCA algo, then I will move to the
0335     // above pattern.
0336     
0337   }
0338 
0339   else {
0340     LogTrace(metname) << "Compatibility NOT satisfied after Backward filter!" << endl;
0341     LogTrace(metname) << "The result of Forward filter will be loaded!" << endl;
0342 
0343     trajectoryContainer.push_back(new Trajectory(trajectoryFW));
0344   }
0345   */
0346   // end 1st attempt
0347 
0348   // 2nd attempt
0349   if (bwfilter()->goodState()) {
0350     LogTrace(metname) << debug.dumpLayer(bwfilter()->lastDetLayer());
0351   } else if (bwfilter()->isCompatibilitySatisfied()) {
0352     LogTrace(metname) << "Compatibility satisfied after Backward filter, but too few valid RecHits ("
0353                       << trajectoryBW.foundHits() << ")." << endl
0354                       << "The (good) result of FW filter will be loaded!" << endl;
0355     trajectoryContainer.push_back(std::make_unique<Trajectory>(trajectoryFW));
0356     return trajectoryContainer;
0357   } else {
0358     LogTrace(metname) << "Compatibility NOT satisfied after Backward filter!" << endl
0359                       << "The Forward trajectory will be invalidated and then loaded!" << endl;
0360     Trajectory defaultTraj(seed, fwDirection);
0361     filter()->createDefaultTrajectory(trajectoryFW, defaultTraj);
0362     trajectoryContainer.push_back(std::make_unique<Trajectory>(defaultTraj));
0363     return trajectoryContainer;
0364   }
0365   // end 2nd attempt
0366 
0367   if (doRefit) {
0368     pair<bool, Trajectory> refitResult = refitter()->refit(trajectoryBW);
0369     if (refitResult.first) {
0370       trajectoryContainer.push_back(std::make_unique<Trajectory>(refitResult.second));
0371       LogTrace(metname) << "StandAloneMuonTrajectoryBuilder Refit output " << endl;
0372       LogTrace(metname) << debug.dumpTSOS(refitResult.second.lastMeasurement().updatedState());
0373     } else
0374       trajectoryContainer.push_back(std::make_unique<Trajectory>(trajectoryBW));
0375   } else
0376     trajectoryContainer.push_back(std::make_unique<Trajectory>(trajectoryBW));
0377 
0378   LogTrace(metname) << "Trajectory saved" << endl;
0379 
0380   return trajectoryContainer;
0381 }
0382 
0383 StandAloneMuonTrajectoryBuilder::DetLayerWithState StandAloneMuonTrajectoryBuilder::propagateTheSeedTSOS(
0384     TrajectoryStateOnSurface& aTSOS, DetId& aDetId) {
0385   const std::string metname = "Muon|RecoMuon|StandAloneMuonTrajectoryBuilder";
0386   MuonPatternRecoDumper debug;
0387 
0388   DetId seedDetId(aDetId);
0389   //  const GeomDet* gdet = theService->trackingGeometry()->idToDet( seedDetId );
0390 
0391   TrajectoryStateOnSurface initialState(aTSOS);
0392 
0393   LogTrace(metname) << "Seed's Pt: " << initialState.freeState()->momentum().perp() << endl;
0394 
0395   LogTrace(metname) << "Seed's id: " << endl;
0396   LogTrace(metname) << debug.dumpMuonId(seedDetId);
0397 
0398   // Get the layer on which the seed relies
0399   const DetLayer* initialLayer = theService->detLayerGeometry()->idToLayer(seedDetId);
0400 
0401   LogTrace(metname) << "Seed's detLayer: " << endl;
0402   LogTrace(metname) << debug.dumpLayer(initialLayer);
0403 
0404   LogTrace(metname) << "TrajectoryStateOnSurface before propagation:" << endl;
0405   LogTrace(metname) << debug.dumpTSOS(initialState);
0406 
0407   PropagationDirection detLayerOrder = (theSeedPosition == recoMuon::in) ? oppositeToMomentum : alongMomentum;
0408 
0409   // ask for compatible layers
0410   vector<const DetLayer*> detLayers;
0411 
0412   if (theNavigationType == "Standard")
0413     detLayers =
0414         theService->muonNavigationSchool()->compatibleLayers(*initialLayer, *initialState.freeState(), detLayerOrder);
0415   else if (theNavigationType == "Direct") {
0416     DirectMuonNavigation navigation(&*theService->detLayerGeometry());
0417     detLayers = navigation.compatibleLayers(*initialState.freeState(), detLayerOrder);
0418   } else
0419     edm::LogError(metname) << "No Properly Navigation Selected!!" << endl;
0420 
0421   LogTrace(metname) << "There are " << detLayers.size() << " compatible layers" << endl;
0422 
0423   DetLayerWithState result = DetLayerWithState(initialLayer, initialState);
0424 
0425   if (!detLayers.empty()) {
0426     LogTrace(metname) << "Compatible layers:" << endl;
0427     for (vector<const DetLayer*>::const_iterator layer = detLayers.begin(); layer != detLayers.end(); layer++) {
0428       LogTrace(metname) << debug.dumpMuonId((*layer)->basicComponents().front()->geographicalId());
0429       LogTrace(metname) << debug.dumpLayer(*layer);
0430     }
0431 
0432     const DetLayer* finalLayer = detLayers.back();
0433 
0434     if (theSeedPosition == recoMuon::in)
0435       LogTrace(metname) << "Most internal one:" << endl;
0436     else
0437       LogTrace(metname) << "Most external one:" << endl;
0438 
0439     LogTrace(metname) << debug.dumpLayer(finalLayer);
0440 
0441     const TrajectoryStateOnSurface propagatedState =
0442         theService->propagator(theSeedPropagatorName)->propagate(initialState, finalLayer->surface());
0443 
0444     if (propagatedState.isValid()) {
0445       result = DetLayerWithState(finalLayer, propagatedState);
0446 
0447       LogTrace(metname) << "Trajectory State on Surface after the extrapolation" << endl;
0448       LogTrace(metname) << debug.dumpTSOS(propagatedState);
0449       LogTrace(metname) << debug.dumpLayer(finalLayer);
0450     } else
0451       LogTrace(metname) << "Extrapolation failed. Keep the original seed state" << endl;
0452   } else
0453     LogTrace(metname) << "No compatible layers. Keep the original seed state" << endl;
0454 
0455   return result;
0456 }