Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:22:16 

0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 
0009 #include "DataFormats/Common/interface/View.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0013 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
0014 
0015 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0016 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0017 #include "MagneticField/Engine/interface/MagneticField.h"
0018 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0019 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0021 
0022 // added by me
0023 
0024 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0025 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0026 #include "DataFormats/GeometrySurface/interface/Surface.h"
0027 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
0028 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0029 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0030 
0031 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0032 #include "DataFormats/GeometrySurface/interface/Surface.h"
0033 #include "DataFormats/Math/interface/Vector.h"
0034 #include "DataFormats/Math/interface/Error.h"
0035 #include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h"
0036 
0037 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0038 
0039 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0040 #include "DataFormats/GeometryCommonDetAlgo/interface/AlignmentPositionError.h"
0041 
0042 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0043 
0044 #include "FWCore/Framework/interface/ESHandle.h"
0045 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0046 #include "FWCore/Utilities/interface/ESGetToken.h"
0047 
0048 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0049 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0050 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0051 
0052 #include <boost/regex.hpp>
0053 
0054 /**
0055  * Configurables:
0056  *
0057  *   Generic:
0058  *     tracks                  = InputTag of a collection of tracks to read from
0059  *     minimumHits             = Minimum hits that the output TrackCandidate must have to be saved
0060  *     replaceWithInactiveHits = instead of discarding hits, replace them with a invalid "inactive" hits,
0061  *                               so multiple scattering is accounted for correctly.
0062  *     stripFrontInvalidHits   = strip invalid hits at the beginning of the track
0063  *     stripBackInvalidHits    = strip invalid hits at the end of the track
0064  *     stripAllInvalidHits     = remove ALL invald hits (might be a problem for multiple scattering, use with care!)
0065  *
0066  *   Per structure:
0067  *      commands = list of commands, to be applied in order as they are written
0068  *      commands can be:
0069  *          "keep XYZ"  , "drop XYZ"    (XYZ = PXB, PXE, TIB, TID, TOB, TEC)
0070  *          "keep XYZ l", "drop XYZ n"  (XYZ as above, n = layer, wheel or disk = 1 .. 6 ; positive and negative are the same )
0071  *
0072  *   Individual modules:
0073  *     detsToIgnore        = individual list of detids on which hits must be discarded
0074  */
0075 namespace reco {
0076   namespace modules {
0077     class CosmicTrackSplitter : public edm::stream::EDProducer<> {
0078     public:
0079       CosmicTrackSplitter(const edm::ParameterSet &iConfig);
0080       void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
0081 
0082     private:
0083       edm::EDGetTokenT<reco::TrackCollection> tokenTracks;
0084       edm::EDGetTokenT<TrajTrackAssociationCollection> tokenTrajTrack;
0085       int totalTracks_;
0086       size_t minimumHits_;
0087 
0088       bool replaceWithInactiveHits_;
0089       bool stripFrontInvalidHits_;
0090       bool stripBackInvalidHits_;
0091       bool stripAllInvalidHits_;
0092       bool excludePixelHits_;
0093 
0094       double dZcut_;
0095       double dXYcut_;
0096 
0097       std::vector<uint32_t> detsToIgnore_;
0098 
0099       edm::ESHandle<TrackerGeometry> theGeometry;
0100       edm::ESHandle<MagneticField> theMagField;
0101       edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tokenGeometry;
0102       edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tokenMagField;
0103 
0104       TrackCandidate makeCandidate(const reco::Track &tk,
0105                                    std::vector<TrackingRecHit *>::iterator hitsBegin,
0106                                    std::vector<TrackingRecHit *>::iterator hitsEnd);
0107 
0108     };  // class
0109 
0110     CosmicTrackSplitter::CosmicTrackSplitter(const edm::ParameterSet &iConfig)
0111         : minimumHits_(iConfig.getParameter<uint32_t>("minimumHits")),
0112           replaceWithInactiveHits_(iConfig.getParameter<bool>("replaceWithInactiveHits")),
0113           stripFrontInvalidHits_(iConfig.getParameter<bool>("stripFrontInvalidHits")),
0114           stripBackInvalidHits_(iConfig.getParameter<bool>("stripBackInvalidHits")),
0115           stripAllInvalidHits_(iConfig.getParameter<bool>("stripAllInvalidHits")),
0116           excludePixelHits_(iConfig.getParameter<bool>("excludePixelHits")),
0117           dZcut_(iConfig.getParameter<double>("dzCut")),
0118           dXYcut_(iConfig.getParameter<double>("dxyCut")),
0119           detsToIgnore_(iConfig.getParameter<std::vector<uint32_t> >("detsToIgnore")) {
0120       // sanity check
0121       if (stripAllInvalidHits_ && replaceWithInactiveHits_) {
0122         throw cms::Exception("Configuration") << "Inconsistent Configuration: you can't set both 'stripAllInvalidHits' "
0123                                                  "and 'replaceWithInactiveHits' to true\n";
0124       }
0125       tokenTracks = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0126       tokenTrajTrack =
0127           consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("tjTkAssociationMapTag"));
0128       tokenGeometry = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0129       tokenMagField = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0130 
0131       LogDebug("CosmicTrackSplitter") << "sanity check";
0132 
0133       // sort detids to ignore
0134       std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
0135 
0136       totalTracks_ = 0;
0137 
0138       // issue the produce<>
0139       produces<TrackCandidateCollection>();
0140     }
0141 
0142     void CosmicTrackSplitter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0143       LogDebug("CosmicTrackSplitter") << "IN THE SPLITTER!!!!!";
0144 
0145       // read with View, so we can read also a TrackRefVector
0146       edm::Handle<std::vector<reco::Track> > tracks;
0147       iEvent.getByToken(tokenTracks, tracks);
0148 
0149       // also need trajectories ...
0150       // Retrieve trajectories and tracks from the event
0151       // -> merely skip if collection is empty
0152       edm::Handle<TrajTrackAssociationCollection> m_TrajTracksMap;
0153       iEvent.getByToken(tokenTrajTrack, m_TrajTracksMap);
0154 
0155       // read from EventSetup
0156       theGeometry = iSetup.getHandle(tokenGeometry);
0157       theMagField = iSetup.getHandle(tokenMagField);
0158 
0159       // prepare output collection
0160       auto output = std::make_unique<TrackCandidateCollection>();
0161       output->reserve(tracks->size());
0162 
0163       // working area and tools
0164       std::vector<TrackingRecHit *> hits;
0165 
0166       // Form pairs of trajectories and tracks
0167       //ConstTrajTrackPairCollection trajTracks;
0168       LogDebug("CosmicTrackSplitter") << "size of map: " << m_TrajTracksMap->size();
0169       int HITTOSPLITFROM = 0;
0170       for (TrajTrackAssociationCollection::const_iterator iPair = m_TrajTracksMap->begin();
0171            iPair != m_TrajTracksMap->end();
0172            iPair++) {
0173         const Trajectory *trajFromMap = &(*(*iPair).key);
0174         const reco::Track *trackFromMap = &(*(*iPair).val);
0175 
0176         // loop to find the hit to split from (by taking dot product of pT and transverse position
0177         std::vector<TrajectoryMeasurement> measurements = trajFromMap->measurements();
0178         int totalNumberOfHits = measurements.size();
0179         int numberOfHits = 0;
0180         double previousDotProduct = 0;
0181         for (trackingRecHit_iterator ith = trackFromMap->recHitsBegin(), edh = trackFromMap->recHitsEnd(); ith != edh;
0182              ++ith) {
0183           GlobalVector stateMomentum = measurements[numberOfHits].forwardPredictedState().globalMomentum();
0184           GlobalPoint statePosition = measurements[numberOfHits].forwardPredictedState().globalPosition();
0185           double dotProduct = stateMomentum.x() * statePosition.x() + stateMomentum.y() * statePosition.y();
0186           if (dotProduct * previousDotProduct < 0) {
0187             //found hit to split from...
0188             HITTOSPLITFROM = numberOfHits;
0189           }
0190 
0191           previousDotProduct = dotProduct;
0192           numberOfHits++;
0193         }
0194         LogDebug("CosmicTrackSplitter") << "number of rechits: " << numberOfHits;
0195 
0196         // check if the trajectories and rechits are in reverse order...
0197         trackingRecHit_iterator bIt = trackFromMap->recHitsBegin();
0198         trackingRecHit_iterator fIt = trackFromMap->recHitsEnd() - 1;
0199         const TrackingRecHit *bHit = (*bIt);
0200         const TrackingRecHit *fHit = (*fIt);
0201         // hit type valid = 0, missing = 1, inactive = 2, bad = 3
0202         if (bHit->type() != 0 || bHit->isValid() != 1) {
0203           //loop over hits forwards until first Valid hit is found
0204           trackingRecHit_iterator ihit;
0205           for (ihit = trackFromMap->recHitsBegin(); ihit != trackFromMap->recHitsEnd(); ++ihit) {
0206             const TrackingRecHit *iHit = (*ihit);
0207             if (iHit->type() == 0 && iHit->isValid() == 1) {
0208               bHit = iHit;
0209               break;
0210             }
0211           }
0212         }
0213         DetId bdetid = bHit->geographicalId();
0214         GlobalPoint bPosHit = theGeometry->idToDetUnit(bdetid)->surface().toGlobal(bHit->localPosition());
0215         if (fHit->type() != 0 || fHit->isValid() != 1) {
0216           //loop over hits backwards until first Valid hit is found
0217           trackingRecHit_iterator ihitf;
0218           for (ihitf = trackFromMap->recHitsEnd() - 1; ihitf != trackFromMap->recHitsBegin(); --ihitf) {
0219             const TrackingRecHit *iHit = (*ihitf);
0220             if (iHit->type() == 0 && iHit->isValid() == 1) {
0221               fHit = iHit;
0222               break;
0223             }
0224           }
0225         }
0226         DetId fdetid = fHit->geographicalId();
0227         GlobalPoint fPosHit = theGeometry->idToDetUnit(fdetid)->surface().toGlobal(fHit->localPosition());
0228         GlobalPoint bPosState = measurements[0].updatedState().globalPosition();
0229         GlobalPoint fPosState = measurements[measurements.size() - 1].updatedState().globalPosition();
0230         bool trajReversedFlag = false;
0231         /*
0232             DetId bdetid = bHit->geographicalId();
0233             DetId fdetid = fHit->geographicalId();
0234             GlobalPoint bPosHit =  theGeometry->idToDetUnit( bdetid )->surface().toGlobal(bHit->localPosition());
0235             GlobalPoint fPosHit =  theGeometry->idToDetUnit( fdetid )->surface().toGlobal(fHit->localPosition());
0236             GlobalPoint bPosState = measurements[0].updatedState().globalPosition();
0237             GlobalPoint fPosState = measurements[measurements.size() - 1].updatedState().globalPosition();
0238             bool trajReversedFlag = false;
0239             */
0240         if (((bPosHit - bPosState).mag() > (bPosHit - fPosState).mag()) &&
0241             ((fPosHit - fPosState).mag() > (fPosHit - bPosState).mag())) {
0242           trajReversedFlag = true;
0243         }
0244         if (trajReversedFlag) {
0245           int temp = HITTOSPLITFROM;
0246           HITTOSPLITFROM = totalNumberOfHits - temp;
0247         }
0248       }
0249 
0250       totalTracks_ = totalTracks_ + tracks->size();
0251       // loop on tracks
0252       for (std::vector<reco::Track>::const_iterator itt = tracks->begin(), edt = tracks->end(); itt != edt; ++itt) {
0253         hits.clear();  // extra safety
0254 
0255         LogDebug("CosmicTrackSplitter") << "ntracks: " << tracks->size();
0256 
0257         // try to find distance of closest approach
0258         GlobalPoint v(itt->vx(), itt->vy(), itt->vz());
0259 
0260         //checks on impact parameter
0261         bool continueWithTrack = true;
0262         if (fabs(v.z()) > dZcut_)
0263           continueWithTrack = false;
0264         if (v.perp() > dXYcut_)
0265           continueWithTrack = false;
0266         if (continueWithTrack == false)
0267           return;
0268 
0269         // LOOP TWICE, ONCE FOR TOP AND ONCE FOR BOTTOM
0270         for (int i = 0; i < 2; ++i) {
0271           hits.clear();  // extra safety
0272           LogDebug("CosmicTrackSplitter") << "   loop on hits of track #" << (itt - tracks->begin());
0273           int usedHitCtr = 0;
0274           int hitCtr = 0;
0275           for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
0276             //hitCtr++;
0277             const TrackingRecHit *hit = (*ith);  // ith is an iterator on edm::Ref to rechit
0278             LogDebug("CosmicTrackSplitter") << "         hit number " << (ith - itt->recHitsBegin());
0279             // let's look at valid hits
0280             if (hit->isValid()) {
0281               LogDebug("CosmicTrackSplitter") << "            valid, detid = " << hit->geographicalId().rawId();
0282               DetId detid = hit->geographicalId();
0283 
0284               if (detid.det() == DetId::Tracker) {  // check for tracker hits
0285                 LogDebug("CosmicTrackSplitter") << "            valid, tracker ";
0286                 bool verdict = false;
0287 
0288                 //trying to get the global position of the hit
0289                 //const GeomDetUnit* geomDetUnit =  theGeometry->idToDetUnit( detid ).;
0290 
0291                 const GlobalPoint pos = theGeometry->idToDetUnit(detid)->surface().toGlobal(hit->localPosition());
0292                 LogDebug("CosmicTrackSplitter") << "hit pos: " << pos << ", dca pos: " << v;
0293 
0294                 // top half
0295                 if ((i == 0) && (hitCtr < HITTOSPLITFROM)) {
0296                   verdict = true;
0297                   LogDebug("CosmicTrackSplitter") << "tophalf";
0298                 }
0299                 // bottom half
0300                 if ((i == 1) && (hitCtr >= HITTOSPLITFROM)) {
0301                   verdict = true;
0302                   LogDebug("CosmicTrackSplitter") << "bottomhalf";
0303                 }
0304 
0305                 // if the hit is good, check again at module level
0306                 if (verdict && std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
0307                   verdict = false;
0308                 }
0309 
0310                 // if hit is good check to make sure that we are keeping pixel hits
0311                 if (excludePixelHits_) {
0312                   // check for pixel hits
0313                   if ((detid.det() == DetId::Tracker) && ((detid.subdetId() == 1) || (detid.subdetId() == 2))) {
0314                     verdict = false;
0315                   }
0316                 }
0317 
0318                 LogDebug("CosmicTrackSplitter")
0319                     << "                   verdict after module list: " << (verdict ? "ok" : "no");
0320                 if (verdict == true) {
0321                   // just copy the hit
0322                   hits.push_back(hit->clone());
0323                   usedHitCtr++;
0324                 } else {
0325                   // still, if replaceWithInactiveHits is true we have to put a new hit
0326                   if (replaceWithInactiveHits_) {
0327                     hits.push_back(new InvalidTrackingRecHit(*hit->det(), TrackingRecHit::inactive));
0328                   }
0329                 }
0330               } else {  // just copy non tracker hits
0331                 hits.push_back(hit->clone());
0332               }
0333             } else {
0334               if (!stripAllInvalidHits_) {
0335                 hits.push_back(hit->clone());
0336               }
0337             }  // is valid hit
0338             LogDebug("CosmicTrackSplitter") << "         end of hit " << (ith - itt->recHitsBegin());
0339             hitCtr++;
0340           }  // loop on hits
0341           LogDebug("CosmicTrackSplitter") << "   end of loop on hits of track #" << (itt - tracks->begin());
0342 
0343           std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
0344 
0345           LogDebug("CosmicTrackSplitter") << "   selected " << hits.size() << " hits ";
0346 
0347           // strip invalid hits at the beginning
0348           if (stripFrontInvalidHits_) {
0349             while ((begin != end) && ((*begin)->isValid() == false))
0350               ++begin;
0351           }
0352 
0353           LogDebug("CosmicTrackSplitter") << "   after front stripping we have " << (end - begin) << " hits ";
0354 
0355           // strip invalid hits at the end
0356           if (stripBackInvalidHits_ && (begin != end)) {
0357             --end;
0358             while ((begin != end) && ((*end)->isValid() == false))
0359               --end;
0360             ++end;
0361           }
0362 
0363           LogDebug("CosmicTrackSplitter") << "   after back stripping we have " << (end - begin) << " hits ";
0364 
0365           // if we still have some hits
0366           //if ((end - begin) >= int(minimumHits_)) {
0367           if (usedHitCtr >= int(minimumHits_)) {
0368             output->push_back(makeCandidate(*itt, begin, end));
0369             LogDebug("CosmicTrackSplitter") << "we made a candidate of " << hits.size() << " hits!";
0370           }
0371           // now delete the hits not used by the candidate
0372           for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0373             if (*begin)
0374               delete *begin;
0375           }
0376           LogDebug("CosmicTrackSplitter")
0377               << "loop: " << i << " has " << usedHitCtr << " active hits and " << hits.size() << " total hits...";
0378           hits.clear();
0379         }  // loop twice for top and bottom
0380       }    // loop on tracks
0381       LogDebug("CosmicTrackSplitter") << "totalTracks_ = " << totalTracks_;
0382       iEvent.put(std::move(output));
0383     }
0384 
0385     TrackCandidate CosmicTrackSplitter::makeCandidate(const reco::Track &tk,
0386                                                       std::vector<TrackingRecHit *>::iterator hitsBegin,
0387                                                       std::vector<TrackingRecHit *>::iterator hitsEnd) {
0388       LogDebug("CosmicTrackSplitter") << "Making a candidate!";
0389 
0390       PropagationDirection pdir = tk.seedDirection();
0391       PTrajectoryStateOnDet state;
0392       if (pdir == anyDirection)
0393         throw cms::Exception("UnimplementedFeature") << "Cannot work with tracks that have 'anyDirecton' \n";
0394       //if ( (pdir == alongMomentum) == ( tk.p() >= tk.outerP() ) ) {
0395       if ((pdir == alongMomentum) == ((tk.outerPosition() - tk.innerPosition()).Dot(tk.momentum()) >= 0)) {
0396         // use inner state
0397         TrajectoryStateOnSurface originalTsosIn(
0398             trajectoryStateTransform::innerStateOnSurface(tk, *theGeometry, &*theMagField));
0399         state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(tk.innerDetId()));
0400       } else {
0401         // use outer state
0402         TrajectoryStateOnSurface originalTsosOut(
0403             trajectoryStateTransform::outerStateOnSurface(tk, *theGeometry, &*theMagField));
0404         state = trajectoryStateTransform::persistentState(originalTsosOut, DetId(tk.outerDetId()));
0405       }
0406 
0407       TrajectorySeed seed(state, TrackCandidate::RecHitContainer(), pdir);
0408 
0409       TrackCandidate::RecHitContainer ownHits;
0410       ownHits.reserve(hitsEnd - hitsBegin);
0411       for (; hitsBegin != hitsEnd; ++hitsBegin) {
0412         ownHits.push_back(*hitsBegin);
0413       }
0414 
0415       TrackCandidate cand(ownHits, seed, state, tk.seedRef());
0416 
0417 #ifdef EDM_ML_DEBUG
0418       LogDebug("CosmicTrackSplitter") << "   dumping the hits now: ";
0419       for (auto const &hit : cand.recHits()) {
0420         LogTrace("CosmicTrackSplitter") << "     hit detid = " << hit.geographicalId().rawId()
0421                                         << ", type  = " << typeid(hit).name();
0422       }
0423 #endif
0424 
0425       return cand;
0426     }
0427 
0428   }  // namespace modules
0429 }  // namespace reco
0430 
0431 // ========= MODULE DEF ==============
0432 #include "FWCore/PluginManager/interface/ModuleDef.h"
0433 #include "FWCore/Framework/interface/MakerMacros.h"
0434 
0435 using reco::modules::CosmicTrackSplitter;
0436 DEFINE_FWK_MODULE(CosmicTrackSplitter);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!