File indexing completed on 2024-04-06 12:06:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include <memory>
0022
0023
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
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
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
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
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
0124
0125
0126
0127
0128
0129
0130
0131
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
0148
0149 usesResource(TFileService::kSharedResource);
0150
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 = tfserv->make<TH1D>(
0266 "algo", "Algorithm number", reco::TrackBase::algoSize, -0.5, static_cast<int>(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
0281
0282 }
0283
0284
0285
0286
0287
0288
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
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
0393 DEFINE_FWK_MODULE(TrackCount);
