CMS 3D CMS Logo

TrackCount.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TrackCount
4 // Class: TrackCount
5 //
13 //
14 // Original Author: Andrea Venturi
15 // Created: Thu Sep 25 16:32:56 CEST 2008
16 // $Id: TrackCount.cc,v 1.15 2011/11/15 10:09:24 venturia Exp $
17 //
18 //
19 
20 // system include files
21 #include <memory>
22 
23 // user include files
26 
31 
33 
34 // my includes
35 
36 #include <iostream>
37 #include <map>
38 #include <string>
39 #include <numeric>
40 
41 #include "TH1F.h"
42 #include "TH1D.h"
43 #include "TH2D.h"
44 #include "TProfile.h"
45 #include "TProfile2D.h"
46 
48 
51 
53 
56 
58 
60 //
61 // class decleration
62 //
63 
64 class TrackCount : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
65 public:
66  explicit TrackCount(const edm::ParameterSet&);
67  ~TrackCount() override;
68 
69 private:
70  void beginRun(const edm::Run&, const edm::EventSetup&) override;
71  void endRun(const edm::Run&, const edm::EventSetup&) override {}
72  void analyze(const edm::Event&, const edm::EventSetup&) override;
73 
74  // ----------member data ---------------------------
75 
76  TH1F* m_ntrk;
77  TProfile* m_ntrkvslumi;
79  TH1F* m_nhptrk;
80  TH1F* m_hhpfrac;
81  TH1F* m_hsqsumptsq;
82  TH1F* m_hphi;
83  TH1F* m_heta;
84  TH1F* m_hcos;
85  TH1F* m_hpt;
86  TH1F* m_chisq;
87  TH1F* m_chisqnorm;
88  TH1F* m_ndof;
89  TH2F* m_chisqvseta;
91  TH2F* m_ndofvseta;
92  TProfile2D* m_hptphieta;
93  TH1D* m_hnlosthits;
94  TH1D* m_hnrhits;
98  TH1D* m_hnlayers;
101  TH1D* m_halgo;
102  TH2D* m_hphieta;
103  TProfile2D* m_hnhitphieta;
104  TProfile2D* m_hnlayerphieta;
105 
106  TProfile** m_ntrkvsorbrun;
107 
108  // std::map<int,int> m_multiplicity;
110  const unsigned int m_maxLS;
111  const unsigned int m_LSfrac;
112  const bool m_2dhistos;
113  const bool m_runHisto;
114  const bool m_dump;
117  const unsigned int m_nptbin;
118  const double m_ptmin;
119  const double m_ptmax;
120 };
121 
122 //
123 // constants, enums and typedefs
124 //
125 
126 //
127 // static data member definitions
128 //
129 
130 //
131 // constructors and destructor
132 //
134  : m_rhm(consumesCollector()),
135  m_maxLS(100),
136  m_LSfrac(16),
137  m_2dhistos(iConfig.getUntrackedParameter<bool>("wanted2DHistos", false)),
138  m_runHisto(iConfig.getUntrackedParameter<bool>("runHisto", false)),
139  m_dump(iConfig.getUntrackedParameter<bool>("dumpTracks", false)),
140  m_trkcollToken(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trackCollection"))),
141  m_lumiProducerToken(consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"))),
142  m_nptbin(iConfig.getUntrackedParameter<unsigned int>("numberPtBins", 200)),
143  m_ptmin(iConfig.getUntrackedParameter<double>("ptMin", 0.)),
144  m_ptmax(iConfig.getUntrackedParameter<double>("ptMax", 20.))
145 
146 {
147  //now do what ever initialization is needed
148 
149  usesResource(TFileService::kSharedResource);
150  // histogram parameters
151 
152  const unsigned int netabin1d = iConfig.getUntrackedParameter<unsigned int>("netabin1D", 120);
153  const unsigned int netabin2d = iConfig.getUntrackedParameter<unsigned int>("netabin2D", 40);
154  const float etamin = iConfig.getUntrackedParameter<double>("etaMin", -3.);
155  const float etamax = iConfig.getUntrackedParameter<double>("etaMax", 3.);
156  const unsigned int nchisqbin1d = iConfig.getUntrackedParameter<unsigned int>("nchi2bin1D", 50);
157  const unsigned int nndofbin1d = iConfig.getUntrackedParameter<unsigned int>("nndofbin1D", 50);
158  const unsigned int nchisqbin2d = iConfig.getUntrackedParameter<unsigned int>("nchi2bin2D", 50);
159  const unsigned int nndofbin2d = iConfig.getUntrackedParameter<unsigned int>("nndofbin2D", 50);
160 
161  edm::LogInfo("TrackCollection") << "Using collection "
162  << iConfig.getParameter<edm::InputTag>("trackCollection").label().c_str();
163 
165 
166  m_ntrk = tfserv->make<TH1F>("ntrk", "Number of Tracks", 2500, -0.5, 2499.5);
167  m_ntrk->GetXaxis()->SetTitle("Tracks");
168  m_ntrk->GetYaxis()->SetTitle("Events");
169  m_ntrkvslumi = tfserv->make<TProfile>("ntrkvslumi", "Number of Tracks vs Luminosity", 250, 0., 10.);
170  m_ntrkvslumi->GetXaxis()->SetTitle("BX lumi [10^{30}cm^{-2}s^{-1}]");
171  m_ntrkvslumi->GetYaxis()->SetTitle("Ntracks");
172  m_ntrkvslumi2D = tfserv->make<TH2D>("ntrkvslumi2D", "Number of Tracks vs Luminosity", 80, 0., 10., 125, -0.5, 2499.5);
173  m_ntrkvslumi2D->GetXaxis()->SetTitle("BX lumi [10^{30}cm^{-2}s^{-1}]");
174  m_ntrkvslumi2D->GetYaxis()->SetTitle("Ntracks");
175 
176  m_nhptrk = tfserv->make<TH1F>("nhptrk", "Number of High Purity Tracks", 2500, -0.5, 2499.5);
177  m_nhptrk->GetXaxis()->SetTitle("Tracks");
178  m_nhptrk->GetYaxis()->SetTitle("Events");
179  m_hhpfrac = tfserv->make<TH1F>("hhpfrac", "Fraction of High Purity Tracks", 51, 0., 1.02);
180  m_hhpfrac->GetXaxis()->SetTitle("hp/all");
181  m_hhpfrac->GetYaxis()->SetTitle("Events");
182  m_hsqsumptsq = tfserv->make<TH1F>("hsqsumptsq", "Sqrt(Sum pt**2)", 1000, 0., 200.);
183  m_hsqsumptsq->GetXaxis()->SetTitle("#sqrt(#Sigma pt^2) (GeV)");
184  m_hsqsumptsq->GetYaxis()->SetTitle("Events");
185 
186  m_hphi = tfserv->make<TH1F>("phi", "Track azimuth", 40, -M_PI, M_PI);
187  m_hphi->GetXaxis()->SetTitle("#phi (rad)");
188  m_hphi->GetYaxis()->SetTitle("Tracks");
189  m_heta = tfserv->make<TH1F>("eta", "Track pseudorapidity", netabin1d, etamin, etamax);
190  m_heta->GetXaxis()->SetTitle("#eta");
191  m_heta->GetYaxis()->SetTitle("Tracks");
192  m_hcos = tfserv->make<TH1F>("cos", "Track polar angle", 50, -1., 1.);
193  m_hcos->GetXaxis()->SetTitle("cos(#theta)");
194  m_hcos->GetYaxis()->SetTitle("Tracks");
195  m_hpt = tfserv->make<TH1F>("pt", "Track pt", m_nptbin, m_ptmin, m_ptmax);
196  m_hpt->GetXaxis()->SetTitle("pt (GeV)");
197  m_hpt->GetYaxis()->SetTitle("Tracks");
198 
199  m_chisq = tfserv->make<TH1F>("chisq", "Track Chi2", nchisqbin1d, 0., 100.);
200  m_chisq->GetXaxis()->SetTitle("chi2");
201  m_chisq->GetYaxis()->SetTitle("Tracks");
202  m_chisqnorm = tfserv->make<TH1F>("chisqnorm", "Track normalized Chi2", nchisqbin1d, 0., 10.);
203  m_chisqnorm->GetXaxis()->SetTitle("normalized chi2");
204  m_chisqnorm->GetYaxis()->SetTitle("Tracks");
205  m_ndof = tfserv->make<TH1F>("ndof", "Track ndof", nndofbin1d, 0., 100.);
206  m_ndof->GetXaxis()->SetTitle("ndof");
207  m_ndof->GetYaxis()->SetTitle("Tracks");
208  if (m_2dhistos) {
209  m_chisqvseta =
210  tfserv->make<TH2F>("chisqvseta", "Track Chi2 vs #eta", netabin2d, etamin, etamax, nchisqbin2d, 0., 100.);
211  m_chisqvseta->GetXaxis()->SetTitle("#eta");
212  m_chisqvseta->GetYaxis()->SetTitle("#chi2");
213  m_chisqnormvseta = tfserv->make<TH2F>(
214  "chisqnormvseta", "Track normalized Chi2 vs #eta", netabin2d, etamin, etamax, nchisqbin2d, 0., 10.);
215  m_chisqnormvseta->GetXaxis()->SetTitle("#eta");
216  m_chisqnormvseta->GetYaxis()->SetTitle("normalized #chi2");
217  m_ndofvseta =
218  tfserv->make<TH2F>("ndofvseta", "Track ndof vs #eta", netabin2d, etamin, etamax, nndofbin2d, 0., 100.);
219  m_ndofvseta->GetXaxis()->SetTitle("#eta");
220  m_ndofvseta->GetYaxis()->SetTitle("ndof");
221  }
222 
223  m_hptphieta =
224  tfserv->make<TProfile2D>("ptphivseta", "Average pt vs #phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
225  m_hptphieta->GetXaxis()->SetTitle("#eta");
226  m_hptphieta->GetYaxis()->SetTitle("#phi");
227 
228  m_hnlosthits = tfserv->make<TH1D>("nlosthits", "Number of Lost Hits", 10, -0.5, 9.5);
229  m_hnlosthits->GetXaxis()->SetTitle("Nlost");
230  m_hnlosthits->GetYaxis()->SetTitle("Tracks");
231 
232  m_hnrhits = tfserv->make<TH1D>("nrhits", "Number of Valid Hits", 55, -0.5, 54.5);
233  m_hnrhits->GetXaxis()->SetTitle("Nvalid");
234  m_hnrhits->GetYaxis()->SetTitle("Tracks");
235  m_hnpixelrhits = tfserv->make<TH1D>("npixelrhits", "Number of Valid Pixel Hits", 20, -0.5, 19.5);
236  m_hnpixelrhits->GetXaxis()->SetTitle("Nvalid");
237  m_hnpixelrhits->GetYaxis()->SetTitle("Tracks");
238  m_hnstriprhits = tfserv->make<TH1D>("nstriprhits", "Number of Valid Strip Hits", 45, -0.5, 44.5);
239  m_hnstriprhits->GetXaxis()->SetTitle("Nvalid");
240  m_hnstriprhits->GetYaxis()->SetTitle("Tracks");
241 
242  m_hnlostlayers = tfserv->make<TH1D>("nlostlayers", "Number of Layers w/o measurement", 10, -0.5, 9.5);
243  m_hnlostlayers->GetXaxis()->SetTitle("Nlostlay");
244  m_hnlostlayers->GetYaxis()->SetTitle("Tracks");
245 
246  m_hnlayers = tfserv->make<TH1D>("nlayers", "Number of Layers", 20, -0.5, 19.5);
247  m_hnlayers->GetXaxis()->SetTitle("Nlayers");
248  m_hnlayers->GetYaxis()->SetTitle("Tracks");
249  m_hnpixellayers = tfserv->make<TH1D>("npixellayers", "Number of Pixel Layers", 10, -0.5, 9.5);
250  m_hnpixellayers->GetXaxis()->SetTitle("Nlayers");
251  m_hnpixellayers->GetYaxis()->SetTitle("Tracks");
252  m_hnstriplayers = tfserv->make<TH1D>("nstriplayers", "Number of Strip Layers", 20, -0.5, 19.5);
253  m_hnstriplayers->GetXaxis()->SetTitle("Nlayers");
254  m_hnstriplayers->GetYaxis()->SetTitle("Tracks");
255 
256  m_hnhitphieta = tfserv->make<TProfile2D>(
257  "nhitphivseta", "N valid hits vs #phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
258  m_hnhitphieta->GetXaxis()->SetTitle("#eta");
259  m_hnhitphieta->GetYaxis()->SetTitle("#phi");
260  m_hnlayerphieta = tfserv->make<TProfile2D>(
261  "nlayerphivseta", "N layers vs #phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
262  m_hnlayerphieta->GetXaxis()->SetTitle("#eta");
263  m_hnlayerphieta->GetYaxis()->SetTitle("#phi");
264 
265  m_halgo =
266  tfserv->make<TH1D>("algo", "Algorithm number", reco::TrackBase::algoSize, -0.5, reco::TrackBase::algoSize - 0.5);
267  m_halgo->GetXaxis()->SetTitle("algorithm");
268  m_halgo->GetYaxis()->SetTitle("Tracks");
269  m_hphieta = tfserv->make<TH2D>("phivseta", "#phi vs #eta", netabin2d, etamin, etamax, 40, -M_PI, M_PI);
270  m_hphieta->GetXaxis()->SetTitle("#eta");
271  m_hphieta->GetYaxis()->SetTitle("#phi");
272 
273  if (m_runHisto) {
275  m_rhm.makeTProfile("ntrkvsorbrun", "Number of tracks vs orbit number", m_maxLS * m_LSfrac, 0, m_maxLS * 262144);
276  }
277 }
278 
280  // do anything here that needs to be done at desctruction time
281  // (e.g. close files, deallocate resources etc.)
282 }
283 
284 //
285 // member functions
286 //
287 
288 // ------------ method called to for each event ------------
290  using namespace edm;
292 
294  iEvent.getByToken(m_trkcollToken, tracks);
295 
296  //
297 
298  m_ntrk->Fill(tracks->size());
299 
300  // get luminosity
301 
303  iEvent.getLuminosityBlock().getByToken(m_lumiProducerToken, ld);
304 
305  if (ld.isValid()) {
306  if (ld->isValid()) {
307  float lumi = ld->lumiValue(LumiDetails::kOCC1, iEvent.bunchCrossing()) * 6.37;
308  m_ntrkvslumi->Fill(lumi, tracks->size());
309  m_ntrkvslumi2D->Fill(lumi, tracks->size());
310  }
311  }
312 
313  if (m_runHisto) {
315  (*m_ntrkvsorbrun)->Fill(iEvent.orbitNumber(), tracks->size());
316  }
317 
318  unsigned int nhptrk = 0;
319  double sumptsq = 0.;
320 
322 
323  if (m_dump)
324  edm::LogInfo("TrackDump") << " isHP algo pt eta phi chi2N chi2 ndof nlay npxl n3dl nlost ";
325 
326  for (reco::TrackCollection::const_iterator it = tracks->begin(); it != tracks->end(); it++) {
327  if (m_dump) {
328  edm::LogVerbatim("TrackDump") << it->quality(quality) << " " << it->algo() << " " << it->pt() << " " << it->eta()
329  << " " << it->phi() << " " << it->normalizedChi2() << " " << it->chi2() << " "
330  << it->ndof() << " " << it->hitPattern().trackerLayersWithMeasurement() << " "
331  << it->hitPattern().pixelLayersWithMeasurement() << " "
332  << it->hitPattern().numberOfValidStripLayersWithMonoAndStereo() << " "
333  << it->hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS)
334  << " ";
335  }
336 
337  m_hnlosthits->Fill(it->hitPattern().numberOfLostTrackerHits(reco::HitPattern::TRACK_HITS));
338 
339  m_hnrhits->Fill(it->hitPattern().numberOfValidTrackerHits());
340  m_hnpixelrhits->Fill(it->hitPattern().numberOfValidPixelHits());
341  m_hnstriprhits->Fill(it->hitPattern().numberOfValidStripHits());
342  m_hnhitphieta->Fill(it->eta(), it->phi(), it->hitPattern().numberOfValidTrackerHits());
343 
344  m_hnlostlayers->Fill(it->hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS));
345 
346  m_hnlayers->Fill(it->hitPattern().trackerLayersWithMeasurement());
347  m_hnpixellayers->Fill(it->hitPattern().pixelLayersWithMeasurement());
348  m_hnstriplayers->Fill(it->hitPattern().stripLayersWithMeasurement());
349  m_hnlayerphieta->Fill(it->eta(), it->phi(), it->hitPattern().trackerLayersWithMeasurement());
350 
351  m_halgo->Fill(it->algo());
352 
353  m_hphi->Fill(it->phi());
354  m_heta->Fill(it->eta());
355  m_hphieta->Fill(it->eta(), it->phi());
356 
357  double pt = it->pt();
358  m_hpt->Fill(pt);
359  m_chisq->Fill(it->chi2());
360  m_chisqnorm->Fill(it->normalizedChi2());
361  m_ndof->Fill(it->ndof());
362  if (m_2dhistos) {
363  m_chisqvseta->Fill(it->eta(), it->chi2());
364  m_chisqnormvseta->Fill(it->eta(), it->normalizedChi2());
365  m_ndofvseta->Fill(it->eta(), it->ndof());
366  }
367  m_hptphieta->Fill(it->eta(), it->phi(), pt);
368  sumptsq += pt * pt;
369  if (it->p())
370  m_hcos->Fill(it->pz() / it->p());
371  if (it->quality(quality))
372  nhptrk++;
373  }
374 
375  m_nhptrk->Fill(nhptrk);
376 
377  const double hpfrac = !tracks->empty() ? double(nhptrk) / double(tracks->size()) : 0.;
378  m_hhpfrac->Fill(hpfrac);
379  m_hsqsumptsq->Fill(sqrt(sumptsq));
380 }
381 
382 void TrackCount::beginRun(const edm::Run& iRun, const edm::EventSetup&) {
383  m_rhm.beginRun(iRun);
384 
385  if (m_runHisto) {
386  (*m_ntrkvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
387  (*m_ntrkvsorbrun)->GetYaxis()->SetTitle("Ntracks");
388  (*m_ntrkvsorbrun)->SetCanExtend(TH1::kXaxis);
389  }
390 }
391 
392 //define this as a plug-in
static const std::string kSharedResource
Definition: TFileService.h:76
Log< level::Info, true > LogVerbatim
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
TH1F * m_ntrk
Definition: TrackCount.cc:76
TH2F * m_chisqvseta
Definition: TrackCount.cc:89
RunHistogramManager m_rhm
Definition: TrackCount.cc:109
TProfile2D * m_hnlayerphieta
Definition: TrackCount.cc:104
TrackQuality
track quality
Definition: TrackBase.h:150
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
TProfile2D * m_hnhitphieta
Definition: TrackCount.cc:103
TProfile2D * m_hptphieta
Definition: TrackCount.cc:92
TH1D * m_halgo
Definition: TrackCount.cc:101
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
~TrackCount() override
Definition: TrackCount.cc:279
const double m_ptmax
Definition: TrackCount.cc:119
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: TrackCount.cc:289
edm::EDGetTokenT< reco::TrackCollection > m_trkcollToken
Definition: TrackCount.cc:115
T getUntrackedParameter(std::string const &, T const &) const
char const * label
TH1D * m_hnstriplayers
Definition: TrackCount.cc:100
TH1F * m_chisqnorm
Definition: TrackCount.cc:87
int iEvent
Definition: GenABIO.cc:224
TH1F * m_chisq
Definition: TrackCount.cc:86
TH2F * m_chisqnormvseta
Definition: TrackCount.cc:90
void beginRun(const edm::Run &, const edm::EventSetup &) override
Definition: TrackCount.cc:382
TH1D * m_hnlostlayers
Definition: TrackCount.cc:97
void endRun(const edm::Run &, const edm::EventSetup &) override
Definition: TrackCount.cc:71
T sqrt(T t)
Definition: SSEVec.h:19
TH1F * m_hpt
Definition: TrackCount.cc:85
TProfile ** makeTProfile(const char *name, const char *title, const unsigned int nbinx, const double xmin, const double xmax)
TProfile ** m_ntrkvsorbrun
Definition: TrackCount.cc:106
bool isValid() const
Definition: LumiDetails.cc:44
TH1F * m_nhptrk
Definition: TrackCount.cc:79
TH1F * m_hsqsumptsq
Definition: TrackCount.cc:81
float lumiValue(AlgoType algo, unsigned int bx) const
Definition: LumiDetails.cc:72
#define M_PI
const unsigned int m_maxLS
Definition: TrackCount.cc:110
Log< level::Info, false > LogInfo
TH1F * m_heta
Definition: TrackCount.cc:83
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:126
TH1D * m_hnlayers
Definition: TrackCount.cc:98
auto const & tracks
cannot be loose
TrackCount(const edm::ParameterSet &)
Definition: TrackCount.cc:133
TH1D * m_hnpixelrhits
Definition: TrackCount.cc:95
void beginRun(const edm::Run &iRun)
TH1D * m_hnrhits
Definition: TrackCount.cc:94
TH1D * m_hnpixellayers
Definition: TrackCount.cc:99
bool isValid() const
Definition: HandleBase.h:70
fixed size matrix
HLT enums.
TH2F * m_ndofvseta
Definition: TrackCount.cc:91
TH1D * m_hnstriprhits
Definition: TrackCount.cc:96
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
const unsigned int m_LSfrac
Definition: TrackCount.cc:111
const double m_ptmin
Definition: TrackCount.cc:118
TH1F * m_hcos
Definition: TrackCount.cc:84
TProfile * m_ntrkvslumi
Definition: TrackCount.cc:77
string quality
const bool m_dump
Definition: TrackCount.cc:114
TH2D * m_ntrkvslumi2D
Definition: TrackCount.cc:78
const bool m_2dhistos
Definition: TrackCount.cc:112
TH1D * m_hnlosthits
Definition: TrackCount.cc:93
TH1F * m_ndof
Definition: TrackCount.cc:88
TH1F * m_hphi
Definition: TrackCount.cc:82
TH2D * m_hphieta
Definition: TrackCount.cc:102
Definition: Run.h:45
const bool m_runHisto
Definition: TrackCount.cc:113
edm::EDGetTokenT< LumiDetails > m_lumiProducerToken
Definition: TrackCount.cc:116
const unsigned int m_nptbin
Definition: TrackCount.cc:117
TH1F * m_hhpfrac
Definition: TrackCount.cc:80