Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-08 23:10:00 

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                 LogDebug("CosmicTrackSplitter")
0292                     << "hit pos: " << theGeometry->idToDetUnit(detid)->surface().toGlobal(hit->localPosition())
0293                     << ", dca pos: " << v;
0294 
0295                 // top half
0296                 if ((i == 0) && (hitCtr < HITTOSPLITFROM)) {
0297                   verdict = true;
0298                   LogDebug("CosmicTrackSplitter") << "tophalf";
0299                 }
0300                 // bottom half
0301                 if ((i == 1) && (hitCtr >= HITTOSPLITFROM)) {
0302                   verdict = true;
0303                   LogDebug("CosmicTrackSplitter") << "bottomhalf";
0304                 }
0305 
0306                 // if the hit is good, check again at module level
0307                 if (verdict && std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
0308                   verdict = false;
0309                 }
0310 
0311                 // if hit is good check to make sure that we are keeping pixel hits
0312                 if (excludePixelHits_) {
0313                   // check for pixel hits
0314                   if ((detid.det() == DetId::Tracker) && ((detid.subdetId() == 1) || (detid.subdetId() == 2))) {
0315                     verdict = false;
0316                   }
0317                 }
0318 
0319                 LogDebug("CosmicTrackSplitter")
0320                     << "                   verdict after module list: " << (verdict ? "ok" : "no");
0321                 if (verdict == true) {
0322                   // just copy the hit
0323                   hits.push_back(hit->clone());
0324                   usedHitCtr++;
0325                 } else {
0326                   // still, if replaceWithInactiveHits is true we have to put a new hit
0327                   if (replaceWithInactiveHits_) {
0328                     hits.push_back(new InvalidTrackingRecHit(*hit->det(), TrackingRecHit::inactive));
0329                   }
0330                 }
0331               } else {  // just copy non tracker hits
0332                 hits.push_back(hit->clone());
0333               }
0334             } else {
0335               if (!stripAllInvalidHits_) {
0336                 hits.push_back(hit->clone());
0337               }
0338             }  // is valid hit
0339             LogDebug("CosmicTrackSplitter") << "         end of hit " << (ith - itt->recHitsBegin());
0340             hitCtr++;
0341           }  // loop on hits
0342           LogDebug("CosmicTrackSplitter") << "   end of loop on hits of track #" << (itt - tracks->begin());
0343 
0344           std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
0345 
0346           LogDebug("CosmicTrackSplitter") << "   selected " << hits.size() << " hits ";
0347 
0348           // strip invalid hits at the beginning
0349           if (stripFrontInvalidHits_) {
0350             while ((begin != end) && ((*begin)->isValid() == false))
0351               ++begin;
0352           }
0353 
0354           LogDebug("CosmicTrackSplitter") << "   after front stripping we have " << (end - begin) << " hits ";
0355 
0356           // strip invalid hits at the end
0357           if (stripBackInvalidHits_ && (begin != end)) {
0358             --end;
0359             while ((begin != end) && ((*end)->isValid() == false))
0360               --end;
0361             ++end;
0362           }
0363 
0364           LogDebug("CosmicTrackSplitter") << "   after back stripping we have " << (end - begin) << " hits ";
0365 
0366           // if we still have some hits
0367           //if ((end - begin) >= int(minimumHits_)) {
0368           if (usedHitCtr >= int(minimumHits_)) {
0369             output->push_back(makeCandidate(*itt, begin, end));
0370             LogDebug("CosmicTrackSplitter") << "we made a candidate of " << hits.size() << " hits!";
0371           }
0372           // now delete the hits not used by the candidate
0373           for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0374             if (*begin)
0375               delete *begin;
0376           }
0377           LogDebug("CosmicTrackSplitter")
0378               << "loop: " << i << " has " << usedHitCtr << " active hits and " << hits.size() << " total hits...";
0379           hits.clear();
0380         }  // loop twice for top and bottom
0381       }  // loop on tracks
0382       LogDebug("CosmicTrackSplitter") << "totalTracks_ = " << totalTracks_;
0383       iEvent.put(std::move(output));
0384     }
0385 
0386     TrackCandidate CosmicTrackSplitter::makeCandidate(const reco::Track &tk,
0387                                                       std::vector<TrackingRecHit *>::iterator hitsBegin,
0388                                                       std::vector<TrackingRecHit *>::iterator hitsEnd) {
0389       LogDebug("CosmicTrackSplitter") << "Making a candidate!";
0390 
0391       PropagationDirection pdir = tk.seedDirection();
0392       PTrajectoryStateOnDet state;
0393       if (pdir == anyDirection)
0394         throw cms::Exception("UnimplementedFeature") << "Cannot work with tracks that have 'anyDirecton' \n";
0395       //if ( (pdir == alongMomentum) == ( tk.p() >= tk.outerP() ) ) {
0396       if ((pdir == alongMomentum) == ((tk.outerPosition() - tk.innerPosition()).Dot(tk.momentum()) >= 0)) {
0397         // use inner state
0398         TrajectoryStateOnSurface originalTsosIn(
0399             trajectoryStateTransform::innerStateOnSurface(tk, *theGeometry, &*theMagField));
0400         state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(tk.innerDetId()));
0401       } else {
0402         // use outer state
0403         TrajectoryStateOnSurface originalTsosOut(
0404             trajectoryStateTransform::outerStateOnSurface(tk, *theGeometry, &*theMagField));
0405         state = trajectoryStateTransform::persistentState(originalTsosOut, DetId(tk.outerDetId()));
0406       }
0407 
0408       TrajectorySeed seed(state, TrackCandidate::RecHitContainer(), pdir);
0409 
0410       TrackCandidate::RecHitContainer ownHits;
0411       ownHits.reserve(hitsEnd - hitsBegin);
0412       for (; hitsBegin != hitsEnd; ++hitsBegin) {
0413         ownHits.push_back(*hitsBegin);
0414       }
0415 
0416       TrackCandidate cand(ownHits, seed, state, tk.seedRef());
0417 
0418 #ifdef EDM_ML_DEBUG
0419       LogDebug("CosmicTrackSplitter") << "   dumping the hits now: ";
0420       for (auto const &hit : cand.recHits()) {
0421         LogTrace("CosmicTrackSplitter") << "     hit detid = " << hit.geographicalId().rawId()
0422                                         << ", type  = " << typeid(hit).name();
0423       }
0424 #endif
0425 
0426       return cand;
0427     }
0428 
0429   }  // namespace modules
0430 }  // namespace reco
0431 
0432 // ========= MODULE DEF ==============
0433 #include "FWCore/PluginManager/interface/ModuleDef.h"
0434 #include "FWCore/Framework/interface/MakerMacros.h"
0435 
0436 using reco::modules::CosmicTrackSplitter;
0437 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!