Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-08 08:16:04

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackCount
0004 // Class:      TrackCount
0005 //
0006 /**\class TrackCount TrackCount.cc trackCount/TrackCount/src/TrackCount.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Andrea Venturi
0015 //         Created:  Thu Sep 25 16:32:56 CEST 2008
0016 // $Id: TrackCount.cc,v 1.15 2011/11/15 10:09:24 venturia Exp $
0017 //
0018 //
0019 
0020 // system include files
0021 #include <memory>
0022 
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/LuminosityBlock.h"
0029 #include "FWCore/Framework/interface/Run.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 
0034 // my includes
0035 
0036 #include <iostream>
0037 #include <map>
0038 #include <string>
0039 #include <numeric>
0040 
0041 #include "TH1F.h"
0042 #include "TH1D.h"
0043 #include "TH2D.h"
0044 #include "TProfile.h"
0045 #include "TProfile2D.h"
0046 
0047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0048 
0049 #include "FWCore/ServiceRegistry/interface/Service.h"
0050 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0051 
0052 #include "FWCore/Utilities/interface/InputTag.h"
0053 
0054 #include "DataFormats/TrackReco/interface/Track.h"
0055 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0056 
0057 #include "DPGAnalysis/SiStripTools/interface/RunHistogramManager.h"
0058 
0059 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0060 //
0061 // class decleration
0062 //
0063 
0064 class TrackCount : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0065 public:
0066   explicit TrackCount(const edm::ParameterSet&);
0067   ~TrackCount() override;
0068 
0069 private:
0070   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0071   void endRun(const edm::Run&, const edm::EventSetup&) override {}
0072   void analyze(const edm::Event&, const edm::EventSetup&) override;
0073 
0074   // ----------member data ---------------------------
0075 
0076   TH1F* m_ntrk;
0077   TProfile* m_ntrkvslumi;
0078   TH2D* m_ntrkvslumi2D;
0079   TH1F* m_nhptrk;
0080   TH1F* m_hhpfrac;
0081   TH1F* m_hsqsumptsq;
0082   TH1F* m_hphi;
0083   TH1F* m_heta;
0084   TH1F* m_hcos;
0085   TH1F* m_hpt;
0086   TH1F* m_chisq;
0087   TH1F* m_chisqnorm;
0088   TH1F* m_ndof;
0089   TH2F* m_chisqvseta;
0090   TH2F* m_chisqnormvseta;
0091   TH2F* m_ndofvseta;
0092   TProfile2D* m_hptphieta;
0093   TH1D* m_hnlosthits;
0094   TH1D* m_hnrhits;
0095   TH1D* m_hnpixelrhits;
0096   TH1D* m_hnstriprhits;
0097   TH1D* m_hnlostlayers;
0098   TH1D* m_hnlayers;
0099   TH1D* m_hnpixellayers;
0100   TH1D* m_hnstriplayers;
0101   TH1D* m_halgo;
0102   TH2D* m_hphieta;
0103   TProfile2D* m_hnhitphieta;
0104   TProfile2D* m_hnlayerphieta;
0105 
0106   TProfile** m_ntrkvsorbrun;
0107 
0108   //  std::map<int,int> m_multiplicity;
0109   RunHistogramManager m_rhm;
0110   const unsigned int m_maxLS;
0111   const unsigned int m_LSfrac;
0112   const bool m_2dhistos;
0113   const bool m_runHisto;
0114   const bool m_dump;
0115   edm::EDGetTokenT<reco::TrackCollection> m_trkcollToken;
0116   edm::EDGetTokenT<LumiDetails> m_lumiProducerToken;
0117   const unsigned int m_nptbin;
0118   const double m_ptmin;
0119   const double m_ptmax;
0120 };
0121 
0122 //
0123 // constants, enums and typedefs
0124 //
0125 
0126 //
0127 // static data member definitions
0128 //
0129 
0130 //
0131 // constructors and destructor
0132 //
0133 TrackCount::TrackCount(const edm::ParameterSet& iConfig)
0134     : m_rhm(consumesCollector()),
0135       m_maxLS(100),
0136       m_LSfrac(16),
0137       m_2dhistos(iConfig.getUntrackedParameter<bool>("wanted2DHistos", false)),
0138       m_runHisto(iConfig.getUntrackedParameter<bool>("runHisto", false)),
0139       m_dump(iConfig.getUntrackedParameter<bool>("dumpTracks", false)),
0140       m_trkcollToken(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trackCollection"))),
0141       m_lumiProducerToken(consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"))),
0142       m_nptbin(iConfig.getUntrackedParameter<unsigned int>("numberPtBins", 200)),
0143       m_ptmin(iConfig.getUntrackedParameter<double>("ptMin", 0.)),
0144       m_ptmax(iConfig.getUntrackedParameter<double>("ptMax", 20.))
0145 
0146 {
0147   //now do what ever initialization is needed
0148 
0149   usesResource(TFileService::kSharedResource);
0150   // histogram parameters
0151 
0152   const unsigned int netabin1d = iConfig.getUntrackedParameter<unsigned int>("netabin1D", 120);
0153   const unsigned int netabin2d = iConfig.getUntrackedParameter<unsigned int>("netabin2D", 40);
0154   const float etamin = iConfig.getUntrackedParameter<double>("etaMin", -3.);
0155   const float etamax = iConfig.getUntrackedParameter<double>("etaMax", 3.);
0156   const unsigned int nchisqbin1d = iConfig.getUntrackedParameter<unsigned int>("nchi2bin1D", 50);
0157   const unsigned int nndofbin1d = iConfig.getUntrackedParameter<unsigned int>("nndofbin1D", 50);
0158   const unsigned int nchisqbin2d = iConfig.getUntrackedParameter<unsigned int>("nchi2bin2D", 50);
0159   const unsigned int nndofbin2d = iConfig.getUntrackedParameter<unsigned int>("nndofbin2D", 50);
0160 
0161   edm::LogInfo("TrackCollection") << "Using collection "
0162                                   << iConfig.getParameter<edm::InputTag>("trackCollection").label().c_str();
0163 
0164   edm::Service<TFileService> tfserv;
0165 
0166   m_ntrk = tfserv->make<TH1F>("ntrk", "Number of Tracks", 2500, -0.5, 2499.5);
0167   m_ntrk->GetXaxis()->SetTitle("Tracks");
0168   m_ntrk->GetYaxis()->SetTitle("Events");
0169   m_ntrkvslumi = tfserv->make<TProfile>("ntrkvslumi", "Number of Tracks vs Luminosity", 250, 0., 10.);
0170   m_ntrkvslumi->GetXaxis()->SetTitle("BX lumi [10^{30}cm^{-2}s^{-1}]");
0171   m_ntrkvslumi->GetYaxis()->SetTitle("Ntracks");
0172   m_ntrkvslumi2D = tfserv->make<TH2D>("ntrkvslumi2D", "Number of Tracks vs Luminosity", 80, 0., 10., 125, -0.5, 2499.5);
0173   m_ntrkvslumi2D->GetXaxis()->SetTitle("BX lumi [10^{30}cm^{-2}s^{-1}]");
0174   m_ntrkvslumi2D->GetYaxis()->SetTitle("Ntracks");
0175 
0176   m_nhptrk = tfserv->make<TH1F>("nhptrk", "Number of High Purity Tracks", 2500, -0.5, 2499.5);
0177   m_nhptrk->GetXaxis()->SetTitle("Tracks");
0178   m_nhptrk->GetYaxis()->SetTitle("Events");
0179   m_hhpfrac = tfserv->make<TH1F>("hhpfrac", "Fraction of High Purity Tracks", 51, 0., 1.02);
0180   m_hhpfrac->GetXaxis()->SetTitle("hp/all");
0181   m_hhpfrac->GetYaxis()->SetTitle("Events");
0182   m_hsqsumptsq = tfserv->make<TH1F>("hsqsumptsq", "Sqrt(Sum pt**2)", 1000, 0., 200.);
0183   m_hsqsumptsq->GetXaxis()->SetTitle("#sqrt(#Sigma pt^2) (GeV)");
0184   m_hsqsumptsq->GetYaxis()->SetTitle("Events");
0185 
0186   m_hphi = tfserv->make<TH1F>("phi", "Track azimuth", 40, -M_PI, M_PI);
0187   m_hphi->GetXaxis()->SetTitle("#phi (rad)");
0188   m_hphi->GetYaxis()->SetTitle("Tracks");
0189   m_heta = tfserv->make<TH1F>("eta", "Track pseudorapidity", netabin1d, etamin, etamax);
0190   m_heta->GetXaxis()->SetTitle("#eta");
0191   m_heta->GetYaxis()->SetTitle("Tracks");
0192   m_hcos = tfserv->make<TH1F>("cos", "Track polar angle", 50, -1., 1.);
0193   m_hcos->GetXaxis()->SetTitle("cos(#theta)");
0194   m_hcos->GetYaxis()->SetTitle("Tracks");
0195   m_hpt = tfserv->make<TH1F>("pt", "Track pt", m_nptbin, m_ptmin, m_ptmax);
0196   m_hpt->GetXaxis()->SetTitle("pt (GeV)");
0197   m_hpt->GetYaxis()->SetTitle("Tracks");
0198 
0199   m_chisq = tfserv->make<TH1F>("chisq", "Track Chi2", nchisqbin1d, 0., 100.);
0200   m_chisq->GetXaxis()->SetTitle("chi2");
0201   m_chisq->GetYaxis()->SetTitle("Tracks");
0202   m_chisqnorm = tfserv->make<TH1F>("chisqnorm", "Track normalized Chi2", nchisqbin1d, 0., 10.);
0203   m_chisqnorm->GetXaxis()->SetTitle("normalized chi2");
0204   m_chisqnorm->GetYaxis()->SetTitle("Tracks");
0205   m_ndof = tfserv->make<TH1F>("ndof", "Track ndof", nndofbin1d, 0., 100.);
0206   m_ndof->GetXaxis()->SetTitle("ndof");
0207   m_ndof->GetYaxis()->SetTitle("Tracks");
0208   if (m_2dhistos) {
0209     m_chisqvseta =
0210         tfserv->make<TH2F>("chisqvseta", "Track Chi2 vs #eta", netabin2d, etamin, etamax, nchisqbin2d, 0., 100.);
0211     m_chisqvseta->GetXaxis()->SetTitle("#eta");
0212     m_chisqvseta->GetYaxis()->SetTitle("#chi2");
0213     m_chisqnormvseta = tfserv->make<TH2F>(
0214         "chisqnormvseta", "Track normalized Chi2 vs #eta", netabin2d, etamin, etamax, nchisqbin2d, 0., 10.);
0215     m_chisqnormvseta->GetXaxis()->SetTitle("#eta");
0216     m_chisqnormvseta->GetYaxis()->SetTitle("normalized #chi2");
0217     m_ndofvseta =
0218         tfserv->make<TH2F>("ndofvseta", "Track ndof vs #eta", netabin2d, etamin, etamax, nndofbin2d, 0., 100.);
0219     m_ndofvseta->GetXaxis()->SetTitle("#eta");
0220     m_ndofvseta->GetYaxis()->SetTitle("ndof");
0221   }
0222 
0223   m_hptphieta =
0224       tfserv->make<TProfile2D>("ptphivseta", "Average pt vs #phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
0225   m_hptphieta->GetXaxis()->SetTitle("#eta");
0226   m_hptphieta->GetYaxis()->SetTitle("#phi");
0227 
0228   m_hnlosthits = tfserv->make<TH1D>("nlosthits", "Number of Lost Hits", 10, -0.5, 9.5);
0229   m_hnlosthits->GetXaxis()->SetTitle("Nlost");
0230   m_hnlosthits->GetYaxis()->SetTitle("Tracks");
0231 
0232   m_hnrhits = tfserv->make<TH1D>("nrhits", "Number of Valid Hits", 55, -0.5, 54.5);
0233   m_hnrhits->GetXaxis()->SetTitle("Nvalid");
0234   m_hnrhits->GetYaxis()->SetTitle("Tracks");
0235   m_hnpixelrhits = tfserv->make<TH1D>("npixelrhits", "Number of Valid Pixel Hits", 20, -0.5, 19.5);
0236   m_hnpixelrhits->GetXaxis()->SetTitle("Nvalid");
0237   m_hnpixelrhits->GetYaxis()->SetTitle("Tracks");
0238   m_hnstriprhits = tfserv->make<TH1D>("nstriprhits", "Number of Valid Strip Hits", 45, -0.5, 44.5);
0239   m_hnstriprhits->GetXaxis()->SetTitle("Nvalid");
0240   m_hnstriprhits->GetYaxis()->SetTitle("Tracks");
0241 
0242   m_hnlostlayers = tfserv->make<TH1D>("nlostlayers", "Number of Layers w/o measurement", 10, -0.5, 9.5);
0243   m_hnlostlayers->GetXaxis()->SetTitle("Nlostlay");
0244   m_hnlostlayers->GetYaxis()->SetTitle("Tracks");
0245 
0246   m_hnlayers = tfserv->make<TH1D>("nlayers", "Number of Layers", 20, -0.5, 19.5);
0247   m_hnlayers->GetXaxis()->SetTitle("Nlayers");
0248   m_hnlayers->GetYaxis()->SetTitle("Tracks");
0249   m_hnpixellayers = tfserv->make<TH1D>("npixellayers", "Number of Pixel Layers", 10, -0.5, 9.5);
0250   m_hnpixellayers->GetXaxis()->SetTitle("Nlayers");
0251   m_hnpixellayers->GetYaxis()->SetTitle("Tracks");
0252   m_hnstriplayers = tfserv->make<TH1D>("nstriplayers", "Number of Strip Layers", 20, -0.5, 19.5);
0253   m_hnstriplayers->GetXaxis()->SetTitle("Nlayers");
0254   m_hnstriplayers->GetYaxis()->SetTitle("Tracks");
0255 
0256   m_hnhitphieta = tfserv->make<TProfile2D>(
0257       "nhitphivseta", "N valid hits vs #phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
0258   m_hnhitphieta->GetXaxis()->SetTitle("#eta");
0259   m_hnhitphieta->GetYaxis()->SetTitle("#phi");
0260   m_hnlayerphieta = tfserv->make<TProfile2D>(
0261       "nlayerphivseta", "N layers vs #phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
0262   m_hnlayerphieta->GetXaxis()->SetTitle("#eta");
0263   m_hnlayerphieta->GetYaxis()->SetTitle("#phi");
0264 
0265   m_halgo =
0266       tfserv->make<TH1D>("algo", "Algorithm number", reco::TrackBase::algoSize, -0.5, reco::TrackBase::algoSize - 0.5);
0267   m_halgo->GetXaxis()->SetTitle("algorithm");
0268   m_halgo->GetYaxis()->SetTitle("Tracks");
0269   m_hphieta = tfserv->make<TH2D>("phivseta", "#phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
0270   m_hphieta->GetXaxis()->SetTitle("#eta");
0271   m_hphieta->GetYaxis()->SetTitle("#phi");
0272 
0273   if (m_runHisto) {
0274     m_ntrkvsorbrun =
0275         m_rhm.makeTProfile("ntrkvsorbrun", "Number of tracks vs orbit number", m_maxLS * m_LSfrac, 0, m_maxLS * 262144);
0276   }
0277 }
0278 
0279 TrackCount::~TrackCount() {
0280   // do anything here that needs to be done at desctruction time
0281   // (e.g. close files, deallocate resources etc.)
0282 }
0283 
0284 //
0285 // member functions
0286 //
0287 
0288 // ------------ method called to for each event  ------------
0289 void TrackCount::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0290   using namespace edm;
0291   edm::Service<TFileService> tfserv;
0292 
0293   Handle<reco::TrackCollection> tracks;
0294   iEvent.getByToken(m_trkcollToken, tracks);
0295 
0296   //
0297 
0298   m_ntrk->Fill(tracks->size());
0299 
0300   // get luminosity
0301 
0302   edm::Handle<LumiDetails> ld;
0303   iEvent.getLuminosityBlock().getByToken(m_lumiProducerToken, ld);
0304 
0305   if (ld.isValid()) {
0306     if (ld->isValid()) {
0307       float lumi = ld->lumiValue(LumiDetails::kOCC1, iEvent.bunchCrossing()) * 6.37;
0308       m_ntrkvslumi->Fill(lumi, tracks->size());
0309       m_ntrkvslumi2D->Fill(lumi, tracks->size());
0310     }
0311   }
0312 
0313   if (m_runHisto) {
0314     if (m_ntrkvsorbrun && *m_ntrkvsorbrun)
0315       (*m_ntrkvsorbrun)->Fill(iEvent.orbitNumber(), tracks->size());
0316   }
0317 
0318   unsigned int nhptrk = 0;
0319   double sumptsq = 0.;
0320 
0321   reco::TrackBase::TrackQuality quality = reco::TrackBase::qualityByName("highPurity");
0322 
0323   if (m_dump)
0324     edm::LogInfo("TrackDump") << " isHP algo pt eta phi chi2N chi2 ndof nlay npxl n3dl nlost ";
0325 
0326   for (reco::TrackCollection::const_iterator it = tracks->begin(); it != tracks->end(); it++) {
0327     if (m_dump) {
0328       edm::LogVerbatim("TrackDump") << it->quality(quality) << " " << it->algo() << " " << it->pt() << " " << it->eta()
0329                                     << " " << it->phi() << " " << it->normalizedChi2() << " " << it->chi2() << " "
0330                                     << it->ndof() << " " << it->hitPattern().trackerLayersWithMeasurement() << " "
0331                                     << it->hitPattern().pixelLayersWithMeasurement() << " "
0332                                     << it->hitPattern().numberOfValidStripLayersWithMonoAndStereo() << " "
0333                                     << it->hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS)
0334                                     << " ";
0335     }
0336 
0337     m_hnlosthits->Fill(it->hitPattern().numberOfLostTrackerHits(reco::HitPattern::TRACK_HITS));
0338 
0339     m_hnrhits->Fill(it->hitPattern().numberOfValidTrackerHits());
0340     m_hnpixelrhits->Fill(it->hitPattern().numberOfValidPixelHits());
0341     m_hnstriprhits->Fill(it->hitPattern().numberOfValidStripHits());
0342     m_hnhitphieta->Fill(it->eta(), it->phi(), it->hitPattern().numberOfValidTrackerHits());
0343 
0344     m_hnlostlayers->Fill(it->hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS));
0345 
0346     m_hnlayers->Fill(it->hitPattern().trackerLayersWithMeasurement());
0347     m_hnpixellayers->Fill(it->hitPattern().pixelLayersWithMeasurement());
0348     m_hnstriplayers->Fill(it->hitPattern().stripLayersWithMeasurement());
0349     m_hnlayerphieta->Fill(it->eta(), it->phi(), it->hitPattern().trackerLayersWithMeasurement());
0350 
0351     m_halgo->Fill(it->algo());
0352 
0353     m_hphi->Fill(it->phi());
0354     m_heta->Fill(it->eta());
0355     m_hphieta->Fill(it->eta(), it->phi());
0356 
0357     double pt = it->pt();
0358     m_hpt->Fill(pt);
0359     m_chisq->Fill(it->chi2());
0360     m_chisqnorm->Fill(it->normalizedChi2());
0361     m_ndof->Fill(it->ndof());
0362     if (m_2dhistos) {
0363       m_chisqvseta->Fill(it->eta(), it->chi2());
0364       m_chisqnormvseta->Fill(it->eta(), it->normalizedChi2());
0365       m_ndofvseta->Fill(it->eta(), it->ndof());
0366     }
0367     m_hptphieta->Fill(it->eta(), it->phi(), pt);
0368     sumptsq += pt * pt;
0369     if (it->p())
0370       m_hcos->Fill(it->pz() / it->p());
0371     if (it->quality(quality))
0372       nhptrk++;
0373   }
0374 
0375   m_nhptrk->Fill(nhptrk);
0376 
0377   const double hpfrac = !tracks->empty() ? double(nhptrk) / double(tracks->size()) : 0.;
0378   m_hhpfrac->Fill(hpfrac);
0379   m_hsqsumptsq->Fill(sqrt(sumptsq));
0380 }
0381 
0382 void TrackCount::beginRun(const edm::Run& iRun, const edm::EventSetup&) {
0383   m_rhm.beginRun(iRun);
0384 
0385   if (m_runHisto) {
0386     (*m_ntrkvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0387     (*m_ntrkvsorbrun)->GetYaxis()->SetTitle("Ntracks");
0388     (*m_ntrkvsorbrun)->SetCanExtend(TH1::kXaxis);
0389   }
0390 }
0391 
0392 //define this as a plug-in
0393 DEFINE_FWK_MODULE(TrackCount);