TrackerHitAssociation/test/myTrackAnalyzer.cc

0001
0002 #include "SimTracker/TrackerHitAssociation/test/myTrackAnalyzer.h"
0003 #include "Math/GenVector/BitReproducible.h"
0004
0005 #include <memory>
0006 #include <iostream>
0007 #include <string>
0008
0009 class TrackerHitAssociator;
0010
0011 myTrackAnalyzer::myTrackAnalyzer(edm::ParameterSet const& conf)
0012     : trackerHitAssociatorConfig_(conf, consumesCollector()),
0013       doPixel_(conf.getParameter<bool>("associatePixel")),
0014       doStrip_(conf.getParameter<bool>("associateStrip")),
0015       trackCollectionTag_(conf.getParameter<edm::InputTag>("trackCollectionTag")),
0016       tokGeo_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0017       tokTrack_(consumes<reco::TrackCollection>(trackCollectionTag_)),
0018       tokSimTk_(consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"))),
0019       tokSimVtx_(consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"))) {}
0020
0021 void myTrackAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0022   //
0023   // extract tracker geometry
0024   //
0025   auto const& theG = setup.getHandle(tokGeo_);
0026
0027   std::cout << "\nEvent ID = " << event.id() << std::endl;
0028
0029   const edm::Handle<reco::TrackCollection>& trackCollection = event.getHandle(tokTrack_);
0030
0031   //get simtrack info
0032   std::vector<SimTrack> theSimTracks;
0033   std::vector<SimVertex> theSimVertexes;
0034
0035   const edm::Handle<edm::SimTrackContainer>& SimTk = event.getHandle(tokSimTk_);
0036   const edm::Handle<edm::SimVertexContainer>& SimVtx = event.getHandle(tokSimVtx_);
0037   theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
0038   theSimVertexes.insert(theSimVertexes.end(), SimVtx->begin(), SimVtx->end());
0039
0040   if (!doPixel_ && !doStrip_)
0041     throw edm::Exception(edm::errors::Configuration, "Strip and pixel association disabled");
0042   //NEW
0043   std::vector<PSimHit> matched;
0044   TrackerHitAssociator associate(event, trackerHitAssociatorConfig_);
0045   std::vector<unsigned int> SimTrackIds;
0046
0047   const reco::TrackCollection tC = *(trackCollection.product());
0048
0049   std::cout << "Reconstructed " << tC.size() << " tracks" << std::endl;
0050
0051   int i = 1;
0052   for (reco::TrackCollection::const_iterator track = tC.begin(); track != tC.end(); track++) {
0053     std::cout << "Track number " << i << std::endl;
0054     std::cout << "\tmomentum: " << track->momentum() << std::endl;
0055     std::cout << "\tPT: " << track->pt() << std::endl;
0056     std::cout << "\tvertex: " << track->vertex() << std::endl;
0057     std::cout << "\timpact parameter: " << track->d0() << std::endl;
0058     std::cout << "\tcharge: " << track->charge() << std::endl;
0059     std::cout << "\tnormalizedChi2: " << track->normalizedChi2() << std::endl;
0060     std::cout << "\tFrom EXTRA : " << std::endl;
0061     std::cout << "\t\touter PT " << track->outerPt() << std::endl;
0062     //
0063     // try and access Hits
0064     //
0065     SimTrackIds.clear();
0066     std::cout << "\t\tNumber of RecHits " << track->recHitsSize() << std::endl;
0067     for (trackingRecHit_iterator it = track->recHitsBegin(); it != track->recHitsEnd(); it++) {
0068       if ((*it)->isValid()) {
0069         std::cout << "\t\t\tRecHit on det " << (*it)->geographicalId().rawId() << std::endl;
0070         std::cout << "\t\t\tRecHit in LP " << (*it)->localPosition() << std::endl;
0071         std::cout << "\t\t\tRecHit in GP "
0072                   << theG->idToDet((*it)->geographicalId())->surface().toGlobal((*it)->localPosition()) << std::endl;
0073         //try SimHit matching
0074         float mindist = 999999;
0075         float dist;
0076         PSimHit closest;
0077         matched.clear();
0078         matched = associate.associateHit((**it));
0079         if (!matched.empty()) {
0080           std::cout << "\t\t\tmatched  " << matched.size() << std::endl;
0081           for (std::vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); m++) {
0082             //        std::cout << "\t\t\tSimhit  ID  " << (*m).trackId()
0083             //         << "\t\t\tSimhit  LP  " << (*m).localPosition()
0084             //         << "\t\t\tSimhit  GP  " << theG->idToDet((*it)->geographicalId())->surface().toGlobal((*m).localPosition()) << std::endl;
0085             //                std::cout << "Track parameters " << theSimTracks[(*m).trackId()].momentum() << std::endl;
0086             //            //do the majority of the simtrack here properly.
0087             //      }
0088             //      std::cout << "\t\t\tSimhit  ID  " << matched[0].trackId() << std::endl;
0089             // << "\t\t\tSimhit  LP  " << matched[0].localPosition()
0090             //  << "\t\t\tSimhit  GP  " << theG->idToDet((*it)->geographicalId())->surface().toGlobal(matched[0].localPosition()) << std::endl;
0091             //std::cout << "Track parameters " << theSimTracks[matched[0].trackId()].momentum() << std::endl;
0092             //now figure out which is the majority of the ids
0093             dist = (*it)->localPosition().x() - (*m).localPosition().x();
0094             if (dist < mindist) {
0095               mindist = dist;
0096               closest = (*m);
0097             }
0098           }
0099           SimTrackIds.push_back(closest.trackId());
0100         }
0101       } else {
0102         std::cout << "\t\t Invalid Hit On " << (*it)->geographicalId().rawId() << std::endl;
0103       }
0104     }
0105
0106     int nmax = 0;
0107     int idmax = -1;
0108     for (size_t j = 0; j < SimTrackIds.size(); j++) {
0109       int n = 0;
0110       n = std::count(SimTrackIds.begin(), SimTrackIds.end(), SimTrackIds[j]);
0111       //    std::cout << " Tracks # of rechits = " << track->recHitsSize() << " found match = " << SimTrackIds.size() << std::endl;
0112       //        std::cout << " rechit = " << i << " sim ID = " << SimTrackIds[i] << " Occurrence = " << n << std::endl;
0113       if (n > nmax) {
0114         nmax = n;
0115         idmax = SimTrackIds[i];
0116       }
0117     }
0118     float totsim = nmax;
0119     float tothits = track->recHitsSize();  //include pixel as well..
0120     float fraction = totsim / tothits;
0121
0122     std::cout << " Track id # " << i << "# of rechits = " << track->recHitsSize() << " matched simtrack id= " << idmax
0123               << " fraction = " << fraction << std::endl;
0124     std::cout << " sim track mom = " << theSimTracks[idmax].momentum() << std::endl;
0125     i++;
0126   }
0127 }