Reference Manual

Go to the documentation of this file.
00001 #include "FWCore/Framework/interface/EDProducer.h"
00002 #include "FWCore/Framework/interface/Event.h"
00003 #include "FWCore/Framework/interface/EventSetup.h"
00004 #include "FWCore/Framework/interface/ESHandle.h"
00005 
00006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00007 #include "FWCore/Utilities/interface/InputTag.h"
00008 
00009 #include "DataFormats/Common/interface/View.h"
00010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00011 #include "DataFormats/TrackReco/interface/Track.h"
00012 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00013 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
00014 
00015 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00016 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00017 #include "MagneticField/Engine/interface/MagneticField.h"
00018 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00019 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00021 
00022 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00023 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00024 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00025 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00026 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00027 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00028 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00029 
00030 // added by me
00031 
00032 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00033 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00034 #include "DataFormats/GeometrySurface/interface/Surface.h"
00035 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
00036 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00037 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00038 
00039 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" 
00040 #include "DataFormats/GeometrySurface/interface/Surface.h" 
00041 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h" 
00042 #include "DataFormats/Math/interface/Vector.h" 
00043 #include "DataFormats/Math/interface/Error.h" 
00044 #include "TrackingTools/TrajectoryState/interface/CopyUsingClone.h" 
00045 #include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h" 
00046 
00047 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00048 
00049 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
00050 #include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h"
00051 
00052 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" 
00053 
00054 #include "FWCore/Framework/interface/ESHandle.h"
00055 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00056 
00057 
00058 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00059 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00060 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00061 
00062 
00063 #include <boost/regex.hpp>
00064 
00086 namespace reco { namespace modules {
00087         class CosmicTrackSplitter : public edm::EDProducer {
00088     public:
00089                 CosmicTrackSplitter(const edm::ParameterSet &iConfig) ; 
00090                 virtual void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) ;
00091                 
00092     private:
00093                 edm::InputTag tracks_;
00094                 edm::InputTag tjTag_;
00095                 int totalTracks_;
00096                 size_t minimumHits_;
00097                 
00098                 bool replaceWithInactiveHits_;
00099                 bool stripFrontInvalidHits_;
00100                 bool stripBackInvalidHits_;
00101                 bool stripAllInvalidHits_;
00102                 bool excludePixelHits_;
00103                 
00104                 double dZcut_;
00105                 double dXYcut_;
00106                 
00107                 std::vector<uint32_t> detsToIgnore_;
00108                 
00109                 edm::ESHandle<TrackerGeometry> theGeometry;
00110                 edm::ESHandle<MagneticField>   theMagField;
00111                 
00112                 TrackCandidate makeCandidate(const reco::Track &tk, std::vector<TrackingRecHit *>::iterator hitsBegin, std::vector<TrackingRecHit *>::iterator hitsEnd) ;
00113                 
00114         }; // class
00115         
00116         
00117         CosmicTrackSplitter::CosmicTrackSplitter(const edm::ParameterSet &iConfig) :
00118     tracks_(iConfig.getParameter<edm::InputTag>("tracks")),
00119     tjTag_(iConfig.getParameter<edm::InputTag>("tjTkAssociationMapTag") ),
00120     minimumHits_(iConfig.getParameter<uint32_t>("minimumHits")),
00121     replaceWithInactiveHits_(iConfig.getParameter<bool>("replaceWithInactiveHits")),
00122     stripFrontInvalidHits_(iConfig.getParameter<bool>("stripFrontInvalidHits")),
00123     stripBackInvalidHits_( iConfig.getParameter<bool>("stripBackInvalidHits") ),
00124     stripAllInvalidHits_(  iConfig.getParameter<bool>("stripAllInvalidHits")  ),
00125         excludePixelHits_(  iConfig.getParameter<bool>("excludePixelHits")  ),
00126     dZcut_(iConfig.getParameter<double>("dzCut") ),
00127     dXYcut_(iConfig.getParameter<double>("dxyCut") ),
00128     detsToIgnore_( iConfig.getParameter<std::vector<uint32_t> >("detsToIgnore") )
00129         {
00130                 // sanity check 
00131                 if (stripAllInvalidHits_ && replaceWithInactiveHits_) {
00132                         throw cms::Exception("Configuration") << "Inconsistent Configuration: you can't set both 'stripAllInvalidHits' and 'replaceWithInactiveHits' to true\n";
00133                 }
00134                 
00135                 LogDebug("CosmicTrackSplitter") << "sanity check";
00136                 
00137                 // sort detids to ignore
00138                 std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
00139                 
00140                 totalTracks_ = 0;
00141                 
00142                 // issue the produce<>
00143                 produces<TrackCandidateCollection>();
00144         }
00145         
00146         void 
00147         CosmicTrackSplitter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) 
00148         {
00149                 LogDebug("CosmicTrackSplitter") << "IN THE SPLITTER!!!!!";
00150                 
00151                 // read with View, so we can read also a TrackRefVector
00152                 edm::Handle<std::vector<reco::Track> > tracks;
00153                 iEvent.getByLabel(tracks_, tracks);
00154                 
00155                 // also need trajectories ...
00156                 // Retrieve trajectories and tracks from the event
00157                 // -> merely skip if collection is empty
00158                 edm::Handle<TrajTrackAssociationCollection> m_TrajTracksMap;
00159                 iEvent.getByLabel( tjTag_, m_TrajTracksMap ); 
00160                 
00161                 // read from EventSetup
00162                 iSetup.get<TrackerDigiGeometryRecord>().get(theGeometry);
00163                 iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
00164                 
00165                 // prepare output collection
00166                 std::auto_ptr<TrackCandidateCollection> output(new TrackCandidateCollection());
00167                 output->reserve(tracks->size());
00168                 
00169                 // working area and tools
00170                 std::vector<TrackingRecHit *> hits;
00171                 
00172                 // Form pairs of trajectories and tracks
00173                 //ConstTrajTrackPairCollection trajTracks;
00174                 LogDebug("CosmicTrackSplitter") << "size of map: " << m_TrajTracksMap->size();
00175                 int HITTOSPLITFROM = 0;
00176                 for ( TrajTrackAssociationCollection::const_iterator iPair = m_TrajTracksMap->begin(); iPair != m_TrajTracksMap->end(); iPair++ ){
00177                         const Trajectory* trajFromMap = &(*(*iPair).key);
00178                         const reco::Track* trackFromMap = &(*(*iPair).val);
00179                         
00180                         // loop to find the hit to split from (by taking dot product of pT and transverse position
00181                         std::vector<TrajectoryMeasurement> measurements = trajFromMap->measurements();
00182                         int totalNumberOfHits = measurements.size();
00183                         int numberOfHits = 0;
00184                         double previousDotProduct = 0;
00185                         for (trackingRecHit_iterator ith = trackFromMap->recHitsBegin(), edh = trackFromMap->recHitsEnd(); ith != edh; ++ith) {
00186                                 
00187                                 GlobalVector stateMomentum = measurements[numberOfHits].forwardPredictedState().globalMomentum();
00188                                 GlobalPoint statePosition = measurements[numberOfHits].forwardPredictedState().globalPosition();
00189                                 double dotProduct = stateMomentum.x()*statePosition.x() + stateMomentum.y()*statePosition.y();
00190                                 if ( dotProduct*previousDotProduct < 0 ){
00191                                         //found hit to split from...
00192                                         HITTOSPLITFROM = numberOfHits;
00193                                 }
00194 
00195                                 previousDotProduct = dotProduct;
00196                                 numberOfHits++;
00197 
00198                         }
00199                         LogDebug("CosmicTrackSplitter") << "number of rechits: " << numberOfHits;
00200                         
00201                         // check if the trajectories and rechits are in reverse order...
00202                         trackingRecHit_iterator bIt = trackFromMap->recHitsBegin();
00203                         trackingRecHit_iterator fIt = trackFromMap->recHitsEnd() - 1;
00204                         const TrackingRecHit* bHit = bIt->get();
00205                         const TrackingRecHit* fHit = fIt->get();
00206                         // hit type valid = 0, missing = 1, inactive = 2, bad = 3
00207                         if( bHit->type() != 0 || bHit->isValid() != 1){
00208                                 //loop over hits forwards until first Valid hit is found
00209                                 trackingRecHit_iterator ihit;
00210                                 for( ihit =  trackFromMap->recHitsBegin(); 
00211                                         ihit != trackFromMap->recHitsEnd(); ++ihit){
00212                                         const TrackingRecHit* iHit = ihit->get();
00213                                         if( iHit->type() == 0 && iHit->isValid() == 1){
00214                                                 bHit = iHit;
00215                                                 break;
00216                                         }
00217                                 }
00218                         }
00219                         DetId bdetid = bHit->geographicalId();
00220                         GlobalPoint bPosHit = theGeometry->idToDetUnit( bdetid)->surface().
00221                         toGlobal(bHit->localPosition());
00222                         if( fHit->type() != 0 || fHit->isValid() != 1){
00223                                 //loop over hits backwards until first Valid hit is found
00224                                 trackingRecHit_iterator ihitf;
00225                                 for( ihitf =  trackFromMap->recHitsEnd()-1; 
00226                                         ihitf != trackFromMap->recHitsBegin(); --ihitf){
00227                                         const TrackingRecHit* iHit = ihitf->get();
00228                                         if( iHit->type() == 0 && iHit->isValid() == 1){
00229                                                 fHit = iHit;
00230                                                 break;
00231                                         }
00232                                 }
00233                         }
00234                         DetId fdetid = fHit->geographicalId();
00235                         GlobalPoint fPosHit =  theGeometry->
00236                         idToDetUnit( fdetid )->surface().toGlobal(fHit->localPosition());
00237                         GlobalPoint bPosState = measurements[0].updatedState().globalPosition();
00238                         GlobalPoint fPosState = measurements[measurements.size()-1].
00239                         updatedState().globalPosition();
00240                         bool trajReversedFlag = false;
00241                         /*
00242                         DetId bdetid = bHit->geographicalId();
00243                         DetId fdetid = fHit->geographicalId();
00244                         GlobalPoint bPosHit =  theGeometry->idToDetUnit( bdetid )->surface().toGlobal(bHit->localPosition());
00245                         GlobalPoint fPosHit =  theGeometry->idToDetUnit( fdetid )->surface().toGlobal(fHit->localPosition());
00246                         GlobalPoint bPosState = measurements[0].updatedState().globalPosition();
00247                         GlobalPoint fPosState = measurements[measurements.size() - 1].updatedState().globalPosition();
00248                         bool trajReversedFlag = false;
00249                         */
00250                         if (( (bPosHit - bPosState).mag() > (bPosHit - fPosState).mag() ) && ( (fPosHit - fPosState).mag() > (fPosHit - bPosState).mag() ) ){
00251                                 trajReversedFlag = true;
00252                         }
00253                         if (trajReversedFlag){ int temp = HITTOSPLITFROM; HITTOSPLITFROM = totalNumberOfHits - temp; }
00254                 }
00255 
00256                 totalTracks_ = totalTracks_ + tracks->size();
00257                 // loop on tracks
00258                 for (std::vector<reco::Track>::const_iterator itt = tracks->begin(), edt = tracks->end(); itt != edt; ++itt) {
00259                         hits.clear(); // extra safety
00260                         
00261                         LogDebug("CosmicTrackSplitter") << "ntracks: " << tracks->size();
00262                         
00263                         // try to find distance of closest approach
00264                         math::XYZPoint refPoint = itt->referencePoint();
00265                         GlobalPoint v( itt->vx(), itt->vy(), itt->vz() );
00266                         
00267                         //checks on impact parameter
00268                         bool continueWithTrack = true;
00269                         if (fabs(v.z()) > dZcut_) continueWithTrack = false;
00270                         if (v.perp() > dXYcut_) continueWithTrack = false;
00271                         if (continueWithTrack == false) return;
00272                         
00273                         // LOOP TWICE, ONCE FOR TOP AND ONCE FOR BOTTOM
00274                         for (int i = 0; i < 2; ++i){
00275                                 hits.clear(); // extra safety
00276                                 LogDebug("CosmicTrackSplitter") << "   loop on hits of track #" << (itt - tracks->begin());
00277                                 int usedHitCtr = 0;
00278                                 int hitCtr = 0;
00279                                 for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
00280                                         //hitCtr++;
00281                                         const TrackingRecHit * hit = ith->get(); // ith is an iterator on edm::Ref to rechit
00282                                         LogDebug("CosmicTrackSplitter") << "         hit number " << (ith - itt->recHitsBegin());
00283                                         // let's look at valid hits
00284                                         if (hit->isValid()) { 
00285                                                 LogDebug("CosmicTrackSplitter") << "            valid, detid = " << hit->geographicalId().rawId();
00286                                                 DetId detid = hit->geographicalId();
00287                                                 
00288                                                 if (detid.det() == DetId::Tracker) {  // check for tracker hits
00289                                                         LogDebug("CosmicTrackSplitter") << "            valid, tracker ";
00290                                                         bool  verdict = false;
00291                                                         
00292                                                         //trying to get the global position of the hit
00293                                                         //const GeomDetUnit* geomDetUnit =  theGeometry->idToDetUnit( detid ).;
00294                                                         
00295                                                         const GlobalPoint pos =  theGeometry->idToDetUnit( detid )->surface().toGlobal(hit->localPosition());
00296                                                         LogDebug("CosmicTrackSplitter") << "hit pos: " << pos << ", dca pos: " << v;
00297                                                         
00298                                                         // top half
00299                                                         if ((i == 0)&&(hitCtr < HITTOSPLITFROM)){
00300                                                                 verdict = true;
00301                                                                 LogDebug("CosmicTrackSplitter") << "tophalf";
00302                                                         }
00303                                                         // bottom half
00304                                                         if ((i == 1)&&(hitCtr >= HITTOSPLITFROM)){
00305                                                                 verdict = true;
00306                                                                 LogDebug("CosmicTrackSplitter") << "bottomhalf";
00307                                                         }
00308                                                         
00309                                                         // if the hit is good, check again at module level
00310                                                         if ( verdict  && std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
00311                                                                 verdict = false;
00312                                                         }
00313                                                         
00314                                                         // if hit is good check to make sure that we are keeping pixel hits
00315                                                         if ( excludePixelHits_){
00316                                                                 if ((detid.det() == DetId::Tracker)&&((detid.subdetId() == 1)||(detid.subdetId() == 2))) {  // check for pixel hits
00317                                                                         verdict = false;
00318                                                                 }
00319                                                         }
00320                                                         
00321                                                         LogDebug("CosmicTrackSplitter") << "                   verdict after module list: " << (verdict ? "ok" : "no");
00322                                                         if (verdict == true) {
00323                                                                 // just copy the hit
00324                                                                 hits.push_back(hit->clone());
00325                                                                 usedHitCtr++;
00326                                                         } 
00327                                                         else {
00328                                                                 // still, if replaceWithInactiveHits is true we have to put a new hit
00329                                                                 if (replaceWithInactiveHits_) {
00330                                                                         hits.push_back(new InvalidTrackingRecHit(detid, TrackingRecHit::inactive));
00331                                                                 } 
00332                                                         }
00333                                                 } 
00334                                                 else { // just copy non tracker hits
00335                                                         hits.push_back(hit->clone());
00336                                                 }
00337                                         } 
00338                                         else {
00339                                                 if (!stripAllInvalidHits_) {
00340                                                         hits.push_back(hit->clone());
00341                                                 } 
00342                                         } // is valid hit
00343                                         LogDebug("CosmicTrackSplitter") << "         end of hit " << (ith - itt->recHitsBegin());
00344                                         hitCtr++;
00345                                 } // loop on hits
00346                                 LogDebug("CosmicTrackSplitter") << "   end of loop on hits of track #" << (itt - tracks->begin());
00347                                 
00348                                 std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
00349                                 
00350                                 LogDebug("CosmicTrackSplitter") << "   selected " << hits.size() << " hits ";
00351                                 
00352                                 // strip invalid hits at the beginning
00353                                 if (stripFrontInvalidHits_) {
00354                                         while ( (begin != end) && ( (*begin)->isValid() == false ) ) ++begin;
00355                                 }
00356                                 
00357                                 LogDebug("CosmicTrackSplitter") << "   after front stripping we have " << (end - begin) << " hits ";
00358                                 
00359                                 // strip invalid hits at the end
00360                                 if (stripBackInvalidHits_ && (begin != end)) {
00361                                         --end;
00362                                         while ( (begin != end) && ( (*end)->isValid() == false ) ) --end;
00363                                         ++end;
00364                                 }
00365                                 
00366                                 LogDebug("CosmicTrackSplitter") << "   after back stripping we have " << (end - begin) << " hits ";
00367                                 
00368                                 // if we still have some hits
00369                                 //if ((end - begin) >= int(minimumHits_)) {
00370                                 if ( usedHitCtr >= int(minimumHits_)) {
00371                                         output->push_back( makeCandidate ( *itt, begin, end ) );
00372                                         LogDebug("CosmicTrackSplitter") << "we made a candidate of " << hits.size() << " hits!";
00373                                 } 
00374                                 // now delete the hits not used by the candidate
00375                                 for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
00376                                         if (*begin) delete *begin;
00377                                 } 
00378                                 LogDebug("CosmicTrackSplitter") << "loop: " << i << " has " << usedHitCtr << " active hits and " << hits.size() << " total hits...";
00379                                 hits.clear();
00380                         } // loop twice for top and bottom
00381                 } // loop on tracks
00382                 LogDebug("CosmicTrackSplitter") << "totalTracks_ = " << totalTracks_;
00383                 iEvent.put(output);
00384         }
00385         
00386         TrackCandidate
00387         CosmicTrackSplitter::makeCandidate(const reco::Track &tk, std::vector<TrackingRecHit *>::iterator hitsBegin, std::vector<TrackingRecHit *>::iterator hitsEnd) {
00388                 
00389                 LogDebug("CosmicTrackSplitter") << "Making a candidate!";
00390                 
00391                 TrajectoryStateTransform transform;
00392                 PropagationDirection   pdir = tk.seedDirection();
00393                 PTrajectoryStateOnDet *state;
00394                 if ( pdir == anyDirection ) throw cms::Exception("UnimplementedFeature") << "Cannot work with tracks that have 'anyDirecton' \n";
00395                 //if ( (pdir == alongMomentum) == ( tk.p() >= tk.outerP() ) ) {
00396                 if ( (pdir == alongMomentum) == (  (tk.outerPosition()-tk.innerPosition()).Dot(tk.momentum()) >= 0    ) ) {
00397                         // use inner state
00398                         TrajectoryStateOnSurface originalTsosIn(transform.innerStateOnSurface(tk, *theGeometry, &*theMagField));
00399                         state = transform.persistentState( originalTsosIn, DetId(tk.innerDetId()) );
00400                 } else { 
00401                         // use outer state
00402                         TrajectoryStateOnSurface originalTsosOut(transform.outerStateOnSurface(tk, *theGeometry, &*theMagField));
00403                         state = transform.persistentState( originalTsosOut, DetId(tk.outerDetId()) );
00404                 }
00405                 
00406                 TrajectorySeed seed(*state, TrackCandidate::RecHitContainer(), pdir);
00407                 
00408                 TrackCandidate::RecHitContainer ownHits;
00409                 ownHits.reserve(hitsEnd - hitsBegin);
00410                 for ( ; hitsBegin != hitsEnd; ++hitsBegin) { ownHits.push_back( *hitsBegin ); }
00411                 
00412                 TrackCandidate cand(ownHits, seed, *state, tk.seedRef());
00413                 delete state;
00414                 
00415                 LogDebug("CosmicTrackSplitter") << "   dumping the hits now: ";
00416                 for (TrackCandidate::range hitR = cand.recHits(); hitR.first != hitR.second; ++hitR.first) {
00417                       LogTrace("CosmicTrackSplitter") << "     hit detid = " << hitR.first->geographicalId().rawId() <<
00418                         ", type  = " << typeid(*hitR.first).name();
00419                 }
00420                 
00421                 return cand;
00422         }
00423         
00424 }} //namespaces
00425 
00426 
00427 // ========= MODULE DEF ==============
00428 #include "FWCore/PluginManager/interface/ModuleDef.h"
00429 #include "FWCore/Framework/interface/MakerMacros.h"
00430 
00431 using reco::modules::CosmicTrackSplitter;
00432 DEFINE_FWK_MODULE(CosmicTrackSplitter);
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoTracker/FinalTrackSelectors/src/CosmicTrackSplitter.cc
