TrackHistory/plugins/QualityCutsAnalyzer.cc

0001 #include <algorithm>
0002 #include <cctype>
0003 #include <iomanip>
0004 #include <set>
0005 #include <sstream>
0006 #include <vector>
0007
0008 #include "TFile.h"
0009 #include "TH1F.h"
0010
0011 #include "HepPDT/ParticleID.hh"
0012
0013 // user include files
0014 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0015 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0016 #include "DataFormats/VertexReco/interface/Vertex.h"
0017
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/Frameworkfwd.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023
0024 #include "TrackingTools/IPTools/interface/IPTools.h"
0025
0026 #include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexSorter.h"
0027
0028 #include "SimTracker/TrackHistory/interface/TrackClassifier.h"
0029
0030 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0031
0032 //
0033 // class decleration
0034 //
0035
0036 class QualityCutsAnalyzer : public edm::one::EDAnalyzer<> {
0037 public:
0038   explicit QualityCutsAnalyzer(const edm::ParameterSet &);
0039   ~QualityCutsAnalyzer() override = default;
0040
0041 private:
0042   void beginJob() override;
0043   void analyze(const edm::Event &, const edm::EventSetup &) override;
0044   void endJob() override;
0045
0046   // Member data
0047   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> ttrkToken_;
0048   const edm::EDGetTokenT<edm::View<reco::Track>> trkToken_;
0049   const edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0050   const edm::EDGetTokenT<reco::JetTracksAssociationCollection> jtaToken_;
0051
0052   typedef std::vector<int> vint;
0053   typedef std::vector<std::string> vstring;
0054
0055   std::string rootFile_;
0056
0057   int minimumNumberOfHits_, minimumNumberOfPixelHits_;
0058   double minimumTransverseMomentum_, maximumChiSquared_;
0059
0060   bool useAllQualities_;
0061   reco::TrackBase::TrackQuality trackQuality_;
0062
0063   void LoopOverJetTracksAssociation(const edm::ESHandle<TransientTrackBuilder> &,
0064                                     const edm::Handle<reco::VertexCollection> &,
0065                                     const edm::Handle<reco::JetTracksAssociationCollection> &);
0066
0067   // Histograms for optimization
0068
0069   struct histogram_element_t {
0070     double sdl;             // Signed decay length
0071     double dta;             // Distance to jet axis
0072     double tip;             // Transverse impact parameter
0073     double lip;             // Longitudinal impact parameter
0074     double ips;             // Impact parameter significance.
0075     double pt;              // Transverse momentum
0076     double chi2;            // Chi^2
0077     std::size_t hits;       // Number of hits
0078     std::size_t pixelhits;  // Number of hits
0079
0080     histogram_element_t(
0081         double d, double a, double t, double l, double i, double p, double c, std::size_t h, std::size_t x) {
0082       sdl = d;
0083       dta = a;
0084       tip = t;
0085       lip = l;
0086       ips = i;
0087       pt = p;
0088       chi2 = c;
0089       hits = h;
0090       pixelhits = x;
0091     }
0092
0093     histogram_element_t(const histogram_element_t &orig) {
0094       sdl = orig.sdl;
0095       dta = orig.dta;
0096       tip = orig.tip;
0097       lip = orig.lip;
0098       ips = orig.ips;
0099       pt = orig.pt;
0100       chi2 = orig.chi2;
0101       hits = orig.hits;
0102       pixelhits = orig.pixelhits;
0103     }
0104   };
0105
0106   typedef std::vector<std::vector<histogram_element_t>> histogram_data_t;
0107   histogram_data_t histogram_data_;
0108
0109   class histogram_t {
0110     TH1F *sdl;
0111     TH1F *dta;
0112     TH1F *tip;
0113     TH1F *lip;
0114     TH1F *ips;
0115     TH1F *pixelhits;
0116     TH1F *pt_1gev;
0117     TH1F *chi2;
0118     TH1F *hits;
0119
0120   public:
0121     histogram_t(const std::string &particleType) {
0122       std::string name, title;
0123       name = std::string("hits_") + particleType;
0124       title = std::string("Hit distribution for ") + particleType;
0125       hits = new TH1F(name.c_str(), title.c_str(), 19, -0.5, 18.5);
0126
0127       name = std::string("chi2_") + particleType;
0128       title = std::string("Chi2 distribution for ") + particleType;
0129       chi2 = new TH1F(name.c_str(), title.c_str(), 100, 0., 30.);
0130
0131       name = std::string("pixelhits_") + particleType;
0132       title = std::string("Pixel hits distribution for ") + particleType;
0133       pixelhits = new TH1F(name.c_str(), title.c_str(), 21, -0.5, 20.5);
0134
0135       name = std::string("pt_1Gev_") + particleType;
0136       title = std::string("Pt distribution close 1Gev for ") + particleType;
0137       pt_1gev = new TH1F(name.c_str(), title.c_str(), 100, 0., 2.);
0138
0139       name = std::string("tip_") + particleType;
0140       title = std::string("Transverse impact parameter distribution for ") + particleType;
0141       tip = new TH1F(name.c_str(), title.c_str(), 100, -0.3, 0.3);
0142
0143       name = std::string("lip_") + particleType;
0144       title = std::string("Longitudinal impact parameter distribution for ") + particleType;
0145       lip = new TH1F(name.c_str(), title.c_str(), 100, -1., 1.);
0146
0147       name = std::string("ips_") + particleType;
0148       title = std::string("IPS distribution for ") + particleType;
0149       ips = new TH1F(name.c_str(), title.c_str(), 100, -25.0, 25.0);
0150
0151       name = std::string("sdl_") + particleType;
0152       title = std::string("Decay length distribution for ") + particleType;
0153       sdl = new TH1F(name.c_str(), title.c_str(), 100, -5., 5.);
0154
0155       name = std::string("dta_") + particleType;
0156       title = std::string("Distance to jet distribution for ") + particleType;
0157       dta = new TH1F(name.c_str(), title.c_str(), 100, 0.0, 0.2);
0158     }
0159
0160     ~histogram_t() {
0161       delete hits;
0162       delete chi2;
0163       delete pixelhits;
0164       delete pt_1gev;
0165       delete tip;
0166       delete lip;
0167       delete ips;
0168       delete sdl;
0169       delete dta;
0170     }
0171
0172     void Fill(const histogram_element_t &data) {
0173       hits->Fill(data.hits);
0174       chi2->Fill(data.chi2);
0175       pixelhits->Fill(data.pt);
0176       pt_1gev->Fill(data.pt);
0177       ips->Fill(data.ips);
0178       tip->Fill(data.tip);
0179       lip->Fill(data.lip);
0180       sdl->Fill(data.sdl);
0181       dta->Fill(data.dta);
0182     }
0183
0184     void Write() {
0185       hits->Write();
0186       chi2->Write();
0187       pixelhits->Write();
0188       pt_1gev->Write();
0189       ips->Write();
0190       tip->Write();
0191       lip->Write();
0192       sdl->Write();
0193       dta->Write();
0194     }
0195   };
0196
0197   // Track classification.
0198   TrackClassifier classifier_;
0199 };
0200
0201 //
0202 // constructors and destructor
0203 //
0204 QualityCutsAnalyzer::QualityCutsAnalyzer(const edm::ParameterSet &config)
0205     : ttrkToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0206       trkToken_(consumes<edm::View<reco::Track>>(config.getUntrackedParameter<edm::InputTag>("trackProducer"))),
0207       vtxToken_(consumes<reco::VertexCollection>(config.getUntrackedParameter<edm::InputTag>("primaryVertexProducer"))),
0208       jtaToken_(consumes<reco::JetTracksAssociationCollection>(
0209           config.getUntrackedParameter<edm::InputTag>("jetTracksAssociation"))),
0210       classifier_(config, consumesCollector()) {
0211   rootFile_ = config.getUntrackedParameter<std::string>("rootFile");
0212
0213   minimumNumberOfHits_ = config.getUntrackedParameter<int>("minimumNumberOfHits");
0214   minimumNumberOfPixelHits_ = config.getUntrackedParameter<int>("minimumNumberOfPixelHits");
0215   minimumTransverseMomentum_ = config.getUntrackedParameter<double>("minimumTransverseMomentum");
0216   maximumChiSquared_ = config.getUntrackedParameter<double>("maximumChiSquared");
0217
0218   std::string trackQualityType = config.getUntrackedParameter<std::string>("trackQualityClass");  // used
0219   trackQuality_ = reco::TrackBase::qualityByName(trackQualityType);
0220   useAllQualities_ = false;
0221
0222   std::transform(
0223       trackQualityType.begin(), trackQualityType.end(), trackQualityType.begin(), (int (*)(int))std::tolower);
0224   if (trackQualityType == "any") {
0225     std::cout << "Using any" << std::endl;
0226     useAllQualities_ = true;
0227   }
0228 }
0229
0230 //
0231 // member functions
0232 //
0233
0234 // ------------ method called to for each event  ------------
0235 void QualityCutsAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &setup) {
0236   // Track collection
0237   edm::Handle<edm::View<reco::Track>> trackCollection;
0238   event.getByToken(trkToken_, trackCollection);
0239   // Primary vertex
0240   edm::Handle<reco::VertexCollection> primaryVertexCollection;
0241   event.getByToken(vtxToken_, primaryVertexCollection);
0242   // Jet to tracks associator
0243   edm::Handle<reco::JetTracksAssociationCollection> jetTracks;
0244   event.getByToken(jtaToken_, jetTracks);
0245   // Trasient track builder
0246   const edm::ESHandle<TransientTrackBuilder> TTbuilder = setup.getHandle(ttrkToken_);
0247
0248   // Setting up event information for the track categories.
0249   classifier_.newEvent(event, setup);
0250
0251   LoopOverJetTracksAssociation(TTbuilder, primaryVertexCollection, jetTracks);
0252 }
0253
0254 // ------------ method called once each job just before starting event loop
0255 // ------------
0256 void QualityCutsAnalyzer::beginJob() { histogram_data_.resize(6); }
0257
0258 // ------------ method called once each job just after ending the event loop
0259 // ------------
0260 void QualityCutsAnalyzer::endJob() {
0261   TFile file(rootFile_.c_str(), "RECREATE");
0262   file.cd();
0263
0264   // saving the histograms
0265   for (std::size_t i = 0; i < 6; i++) {
0266     std::string particle;
0267     if (i == 0)
0268       particle = std::string("B_tracks");
0269     else if (i == 1)
0270       particle = std::string("C_tracks");
0271     else if (i == 2)
0272       particle = std::string("nonB_tracks");
0273     else if (i == 3)
0274       particle = std::string("displaced_tracks");
0275     else if (i == 4)
0276       particle = std::string("bad_tracks");
0277     else
0278       particle = std::string("fake_tracks");
0279
0280     histogram_t histogram(particle);
0281
0282     for (std::size_t j = 0; j < histogram_data_[i].size(); j++)
0283       histogram.Fill(histogram_data_[i][j]);
0284
0285     histogram.Write();
0286   }
0287
0288   file.Flush();
0289 }
0290
0291 void QualityCutsAnalyzer::LoopOverJetTracksAssociation(
0292     const edm::ESHandle<TransientTrackBuilder> &TTbuilder,
0293     const edm::Handle<reco::VertexCollection> &primaryVertexProducer_,
0294     const edm::Handle<reco::JetTracksAssociationCollection> &jetTracksAssociation) {
0295   const TransientTrackBuilder *bproduct = TTbuilder.product();
0296
0297   // getting the primary vertex
0298   // use first pv of the collection
0299   reco::Vertex pv;
0300
0301   if (!primaryVertexProducer_->empty()) {
0302     PrimaryVertexSorter pvs;
0303     std::vector<reco::Vertex> sortedList = pvs.sortedList(*(primaryVertexProducer_.product()));
0304     pv = (sortedList.front());
0305   } else {  // create a dummy PV
0306     // cout << "NO PV FOUND" << endl;
0307     reco::Vertex::Error e;
0308     e(0, 0) = 0.0015 * 0.0015;
0309     e(1, 1) = 0.0015 * 0.0015;
0310     e(2, 2) = 15. * 15.;
0311     reco::Vertex::Point p(0, 0, 0);
0312     pv = reco::Vertex(p, e, 1, 1, 1);
0313   }
0314
0315   reco::JetTracksAssociationCollection::const_iterator it = jetTracksAssociation->begin();
0316
0317   int i = 0;
0318
0319   for (; it != jetTracksAssociation->end(); it++, i++) {
0320     // get jetTracks object
0321     reco::JetTracksAssociationRef jetTracks(jetTracksAssociation, i);
0322
0323     double pvZ = pv.z();
0324     GlobalVector direction(jetTracks->first->px(), jetTracks->first->py(), jetTracks->first->pz());
0325
0326     // get the tracks associated to the jet
0327     reco::TrackRefVector tracks = jetTracks->second;
0328     for (std::size_t index = 0; index < tracks.size(); index++) {
0329       edm::RefToBase<reco::Track> track(tracks[index]);
0330
0331       double pt = tracks[index]->pt();
0332       double chi2 = tracks[index]->normalizedChi2();
0333       int hits = tracks[index]->hitPattern().numberOfValidHits();
0334       int pixelHits = tracks[index]->hitPattern().numberOfValidPixelHits();
0335
0336       if (hits < minimumNumberOfHits_ || pixelHits < minimumNumberOfPixelHits_ || pt < minimumTransverseMomentum_ ||
0337           chi2 > maximumChiSquared_ || (!useAllQualities_ && !tracks[index]->quality(trackQuality_)))
0338         continue;
0339
0340       const reco::TransientTrack transientTrack = bproduct->build(&(*tracks[index]));
0341       double dta = -IPTools::jetTrackDistance(transientTrack, direction, pv).second.value();
0342       double sdl = IPTools::signedDecayLength3D(transientTrack, direction, pv).second.value();
0343       double ips = IPTools::signedImpactParameter3D(transientTrack, direction, pv).second.value();
0344       double d0 = IPTools::signedTransverseImpactParameter(transientTrack, direction, pv).second.value();
0345       double dz = tracks[index]->dz() - pvZ;
0346
0347       // Classify the reco track;
0348       classifier_.evaluate(edm::RefToBase<reco::Track>(tracks[index]));
0349
0350       // Check for the different categories
0351       if (classifier_.is(TrackClassifier::Fake))
0352         histogram_data_[5].push_back(histogram_element_t(sdl, dta, d0, dz, ips, pt, chi2, hits, pixelHits));
0353       else if (classifier_.is(TrackClassifier::BWeakDecay))
0354         histogram_data_[0].push_back(histogram_element_t(sdl, dta, d0, dz, ips, pt, chi2, hits, pixelHits));
0355       else if (classifier_.is(TrackClassifier::Bad))
0356         histogram_data_[4].push_back(histogram_element_t(sdl, dta, d0, dz, ips, pt, chi2, hits, pixelHits));
0357       else if (!classifier_.is(TrackClassifier::CWeakDecay) && !classifier_.is(TrackClassifier::PrimaryVertex))
0358         histogram_data_[3].push_back(histogram_element_t(sdl, dta, d0, dz, ips, pt, chi2, hits, pixelHits));
0359       else if (classifier_.is(TrackClassifier::CWeakDecay))
0360         histogram_data_[1].push_back(histogram_element_t(sdl, dta, d0, dz, ips, pt, chi2, hits, pixelHits));
0361       else
0362         histogram_data_[2].push_back(histogram_element_t(sdl, dta, d0, dz, ips, pt, chi2, hits, pixelHits));
0363     }
0364   }
0365 }
0366
0367 DEFINE_FWK_MODULE(QualityCutsAnalyzer);