CMS 3D CMS Logo

TrackBuildingAnalyzer.cc
Go to the documentation of this file.
11 
13 
18 #include <string>
19 #include "TMath.h"
20 
21 #include <iostream>
22 
24  doAllPlots( iConfig.getParameter<bool>("doAllPlots")),
25  doAllSeedPlots(iConfig.getParameter<bool>("doSeedParameterHistos")),
26  doTCPlots( iConfig.getParameter<bool>("doTrackCandHistos")),
27  doAllTCPlots( iConfig.getParameter<bool>("doAllTrackCandHistos")),
28  doPT( iConfig.getParameter<bool>("doSeedPTHisto")),
29  doETA( iConfig.getParameter<bool>("doSeedETAHisto")),
30  doPHI( iConfig.getParameter<bool>("doSeedPHIHisto")),
31  doPHIVsETA( iConfig.getParameter<bool>("doSeedPHIVsETAHisto")),
32  doTheta( iConfig.getParameter<bool>("doSeedThetaHisto")),
33  doQ( iConfig.getParameter<bool>("doSeedQHisto")),
34  doDxy( iConfig.getParameter<bool>("doSeedDxyHisto")),
35  doDz( iConfig.getParameter<bool>("doSeedDzHisto")),
36  doNRecHits( iConfig.getParameter<bool>("doSeedNRecHitsHisto")),
37  doProfPHI( iConfig.getParameter<bool>("doSeedNVsPhiProf")),
38  doProfETA( iConfig.getParameter<bool>("doSeedNVsEtaProf")),
39  doStopSource( iConfig.getParameter<bool>("doStopSource")),
40  doMVAPlots( iConfig.getParameter<bool>("doMVAPlots")),
41  doRegionPlots( iConfig.getParameter<bool>("doRegionPlots"))
42 {
43 }
44 
46 {
47 
48  // parameters from the configuration
49  std::string AlgoName = iConfig.getParameter<std::string>("AlgoName");
50  std::string MEFolderName = iConfig.getParameter<std::string>("FolderName");
51 
52  // std::cout << "[TrackBuildingAnalyzer::beginRun] AlgoName: " << AlgoName << std::endl;
53 
54  // use the AlgoName and Quality Name
55  const std::string& CatagoryName = AlgoName;
56 
57  // get binning from the configuration
58  int TrackPtBin = iConfig.getParameter<int>( "TrackPtBin");
59  double TrackPtMin = iConfig.getParameter<double>("TrackPtMin");
60  double TrackPtMax = iConfig.getParameter<double>("TrackPtMax");
61 
62  int PhiBin = iConfig.getParameter<int>( "PhiBin");
63  double PhiMin = iConfig.getParameter<double>("PhiMin");
64  double PhiMax = iConfig.getParameter<double>("PhiMax");
65 
66  int EtaBin = iConfig.getParameter<int>( "EtaBin");
67  double EtaMin = iConfig.getParameter<double>("EtaMin");
68  double EtaMax = iConfig.getParameter<double>("EtaMax");
69 
70  int ThetaBin = iConfig.getParameter<int>( "ThetaBin");
71  double ThetaMin = iConfig.getParameter<double>("ThetaMin");
72  double ThetaMax = iConfig.getParameter<double>("ThetaMax");
73 
74  int TrackQBin = iConfig.getParameter<int>( "TrackQBin");
75  double TrackQMin = iConfig.getParameter<double>("TrackQMin");
76  double TrackQMax = iConfig.getParameter<double>("TrackQMax");
77 
78  int SeedDxyBin = iConfig.getParameter<int>( "SeedDxyBin");
79  double SeedDxyMin = iConfig.getParameter<double>("SeedDxyMin");
80  double SeedDxyMax = iConfig.getParameter<double>("SeedDxyMax");
81 
82  int SeedDzBin = iConfig.getParameter<int>( "SeedDzBin");
83  double SeedDzMin = iConfig.getParameter<double>("SeedDzMin");
84  double SeedDzMax = iConfig.getParameter<double>("SeedDzMax");
85 
86  int SeedHitBin = iConfig.getParameter<int>( "SeedHitBin");
87  double SeedHitMin = iConfig.getParameter<double>("SeedHitMin");
88  double SeedHitMax = iConfig.getParameter<double>("SeedHitMax");
89 
90  int TCDxyBin = iConfig.getParameter<int>( "TCDxyBin");
91  double TCDxyMin = iConfig.getParameter<double>("TCDxyMin");
92  double TCDxyMax = iConfig.getParameter<double>("TCDxyMax");
93 
94  int TCDzBin = iConfig.getParameter<int>( "TCDzBin");
95  double TCDzMin = iConfig.getParameter<double>("TCDzMin");
96  double TCDzMax = iConfig.getParameter<double>("TCDzMax");
97 
98  int TCHitBin = iConfig.getParameter<int>( "TCHitBin");
99  double TCHitMin = iConfig.getParameter<double>("TCHitMin");
100  double TCHitMax = iConfig.getParameter<double>("TCHitMax");
101 
102  int MVABin = iConfig.getParameter<int>( "MVABin");
103  double MVAMin = iConfig.getParameter<double>("MVAMin");
104  double MVAMax = iConfig.getParameter<double>("MVAMax");
105 
106 
107  edm::InputTag seedProducer = iConfig.getParameter<edm::InputTag>("SeedProducer");
108  edm::InputTag tcProducer = iConfig.getParameter<edm::InputTag>("TCProducer");
109  std::vector<std::string> mvaProducers = iConfig.getParameter<std::vector<std::string> >("MVAProducers");
110  edm::InputTag regionProducer = iConfig.getParameter<edm::InputTag>("RegionProducer");
111 
112  // if (doAllPlots){doAllSeedPlots=true; doTCPlots=true;}
113 
114  ibooker.setCurrentFolder(MEFolderName);
115 
116  // book the Seed histograms
117  // ---------------------------------------------------------------------------------//
118  // std::cout << "[TrackBuildingAnalyzer::beginRun] MEFolderName: " << MEFolderName << std::endl;
119  ibooker.setCurrentFolder(MEFolderName+"/TrackBuilding");
120 
121  if (doAllSeedPlots || doPT) {
122  histname = "SeedPt_"+seedProducer.label() + "_";
123  SeedPt = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackPtBin, TrackPtMin, TrackPtMax);
124  SeedPt->setAxisTitle("Seed p_{T} (GeV/c)", 1);
125  SeedPt->setAxisTitle("Number of Seeds", 2);
126  }
127 
128  if (doAllSeedPlots || doETA) {
129  histname = "SeedEta_"+seedProducer.label() + "_";
130  SeedEta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax);
131  SeedEta->setAxisTitle("Seed #eta", 1);
132  SeedEta->setAxisTitle("Number of Seeds", 2);
133  }
134 
135  if (doAllSeedPlots || doPHI) {
136  histname = "SeedPhi_"+seedProducer.label() + "_";
137  SeedPhi = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax);
138  SeedPhi->setAxisTitle("Seed #phi", 1);
139  SeedPhi->setAxisTitle("Number of Seed", 2);
140  }
141 
142  if (doAllSeedPlots || doPHIVsETA) {
143  histname = "SeedPhiVsEta_"+seedProducer.label() + "_";
144  SeedPhiVsEta = ibooker.book2D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, PhiBin, PhiMin, PhiMax);
145  SeedPhiVsEta->setAxisTitle("Seed #eta", 1);
146  SeedPhiVsEta->setAxisTitle("Seed #phi", 2);
147  }
148 
149  if (doAllSeedPlots || doTheta){
150  histname = "SeedTheta_"+seedProducer.label() + "_";
151  SeedTheta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, ThetaBin, ThetaMin, ThetaMax);
152  SeedTheta->setAxisTitle("Seed #theta", 1);
153  SeedTheta->setAxisTitle("Number of Seeds", 2);
154  }
155 
156  if (doAllSeedPlots || doQ){
157  histname = "SeedQ_"+seedProducer.label() + "_";
158  SeedQ = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackQBin, TrackQMin, TrackQMax);
159  SeedQ->setAxisTitle("Seed Charge", 1);
160  SeedQ->setAxisTitle("Number of Seeds",2);
161  }
162 
163  if (doAllSeedPlots || doDxy){
164  histname = "SeedDxy_"+seedProducer.label() + "_";
165  SeedDxy = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, SeedDxyBin, SeedDxyMin, SeedDxyMax);
166  SeedDxy->setAxisTitle("Seed d_{xy} (cm)", 1);
167  SeedDxy->setAxisTitle("Number of Seeds",2);
168  }
169 
170  if (doAllSeedPlots || doDz){
171  histname = "SeedDz_"+seedProducer.label() + "_";
172  SeedDz = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, SeedDzBin, SeedDzMin, SeedDzMax);
173  SeedDz->setAxisTitle("Seed d_{z} (cm)", 1);
174  SeedDz->setAxisTitle("Number of Seeds",2);
175  }
176 
177  if (doAllSeedPlots || doNRecHits){
178  histname = "NumberOfRecHitsPerSeed_"+seedProducer.label() + "_";
179  NumberOfRecHitsPerSeed = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, SeedHitBin, SeedHitMin, SeedHitMax);
180  NumberOfRecHitsPerSeed->setAxisTitle("Number of RecHits per Seed", 1);
181  NumberOfRecHitsPerSeed->setAxisTitle("Number of Seeds",2);
182  }
183 
184  if (doAllSeedPlots || doProfPHI){
185  histname = "NumberOfRecHitsPerSeedVsPhiProfile_"+seedProducer.label() + "_";
186  NumberOfRecHitsPerSeedVsPhiProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax, SeedHitBin, SeedHitMin, SeedHitMax,"s");
188  NumberOfRecHitsPerSeedVsPhiProfile->setAxisTitle("Number of RecHits of each Seed",2);
189  }
190 
191  if (doAllSeedPlots || doProfETA){
192  histname = "NumberOfRecHitsPerSeedVsEtaProfile_"+seedProducer.label() + "_";
193  NumberOfRecHitsPerSeedVsEtaProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, SeedHitBin, SeedHitMin, SeedHitMax,"s");
195  NumberOfRecHitsPerSeedVsEtaProfile->setAxisTitle("Number of RecHits of each Seed",2);
196  }
197 
198  if (doRegionPlots) {
199  if (doAllSeedPlots || doPT) {
200  auto ptBin = iConfig.getParameter<int>( "RegionCandidatePtBin");
201  auto ptMin = iConfig.getParameter<double>("RegionCandidatePtMin");
202  auto ptMax = iConfig.getParameter<double>("RegionCandidatePtMax");
203 
204  histname = "TrackingRegionCandidatePt_"+seedProducer.label() + "_";
205  TrackingRegionCandidatePt = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, ptBin, ptMin, ptMax);
206  TrackingRegionCandidatePt->setAxisTitle("TrackingRegion Candidate p_{T} (GeV/c)", 1);
207  TrackingRegionCandidatePt->setAxisTitle("Number of TrackingRegion Candidates", 2);
208  }
209 
210  if (doAllSeedPlots || doETA) {
211  histname = "TrackingRegionCandidateEta_"+seedProducer.label() + "_";
212  TrackingRegionCandidateEta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax);
213  TrackingRegionCandidateEta->setAxisTitle("TrackingRegion Candidate #eta", 1);
214  TrackingRegionCandidateEta->setAxisTitle("Number of TrackingRegion Candidates", 2);
215  }
216 
217  if (doAllSeedPlots || doPHI) {
218  histname = "TrackingRegionCandidatePhi_"+seedProducer.label() + "_";
219  TrackingRegionCandidatePhi = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax);
220  TrackingRegionCandidatePhi->setAxisTitle("TrackingRegion Candidate #phi", 1);
221  TrackingRegionCandidatePhi->setAxisTitle("Number of TrackingRegion Candidates", 2);
222  }
223 
224  if (doAllSeedPlots || doPHIVsETA) {
225  histname = "TrackingRegionCandidatePhiVsEta_"+seedProducer.label() + "_";
226  TrackingRegionCandidatePhiVsEta = ibooker.book2D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, PhiBin, PhiMin, PhiMax);
227  TrackingRegionCandidatePhiVsEta->setAxisTitle("TrackingRegion Candidate #eta", 1);
228  TrackingRegionCandidatePhiVsEta->setAxisTitle("TrackingRegion Candidate #phi", 2);
229  }
230  }
231 
232  if (doAllSeedPlots || doStopSource) {
233  const auto stopReasonSize = static_cast<unsigned int>(SeedStopReason::SIZE);
234 
235  const auto candsBin = iConfig.getParameter<int>( "SeedCandBin");
236  const auto candsMin = iConfig.getParameter<double>("SeedCandMin");
237  const auto candsMax = iConfig.getParameter<double>("SeedCandMax");
238 
239  histname = "SeedStoppingSource_"+seedProducer.label() + "_";
240  seedStoppingSource = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, stopReasonSize, 0., stopReasonSize);
241  seedStoppingSource->setAxisTitle("Stopping reason",1);
242  seedStoppingSource->setAxisTitle("Number of seeds",2);
243 
244  histname = "StoppingSourceVsPhi_"+seedProducer.label() + "_";
245  seedStoppingSourceVsPhi = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax, 2, 0., 2.);
246  seedStoppingSourceVsPhi->setAxisTitle("seed #phi",1);
247  seedStoppingSourceVsPhi->setAxisTitle("fraction stopped",2);
248 
249  histname = "StoppingSourceVsEta_"+seedProducer.label() + "_";
250  seedStoppingSourceVsEta = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, 2, 0., 2.);
251  seedStoppingSourceVsEta->setAxisTitle("seed #eta",1);
252  seedStoppingSourceVsEta->setAxisTitle("fraction stopped",2);
253 
254  histname = "NumberOfTrajCandsPerSeed_"+seedProducer.label() + "_";
255  numberOfTrajCandsPerSeed = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, candsBin, candsMin, candsMax);
256  numberOfTrajCandsPerSeed->setAxisTitle("Number of Trajectory Candidate for each Seed", 1);
257  numberOfTrajCandsPerSeed->setAxisTitle("Number of Seeds", 2);
258 
259  histname = "NumberOfTrajCandsPerSeedVsPhi_"+seedProducer.label() + "_";
260  numberOfTrajCandsPerSeedVsPhi = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax, candsBin, candsMin, candsMax);
262  numberOfTrajCandsPerSeedVsPhi->setAxisTitle("Number of Trajectory Candidates for each Seed", 2);
263 
264  histname = "NumberOfTrajCandsPerSeedVsEta_"+seedProducer.label() + "_";
265  numberOfTrajCandsPerSeedVsEta = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, candsBin, candsMin, candsMax);
267  numberOfTrajCandsPerSeedVsEta->setAxisTitle("Number of Trajectory Candidates for each Seed", 2);
268 
269  histname = "SeedStoppingSourceVsNumberOfTrajCandsPerSeed_"+seedProducer.label() + "_";
270  seedStoppingSourceVsNumberOfTrajCandsPerSeed = ibooker.book2D(histname+CatagoryName, histname+CatagoryName, candsBin, candsMin, candsMax, stopReasonSize, 0, stopReasonSize);
271  seedStoppingSourceVsNumberOfTrajCandsPerSeed->setAxisTitle("Number of Trajectory Candidates for each Seed", 1);
272 
273  for(unsigned int i=0; i<stopReasonSize; ++i) {
276  }
277  }
278 
279  if (doAllTCPlots || doStopSource) {
280  // DataFormats/TrackReco/interface/TrajectoryStopReasons.h
281  size_t StopReasonNameSize = sizeof(StopReasonName::StopReasonName)/sizeof(std::string);
282 
283  histname = "StoppingSource_"+seedProducer.label() + "_";
284  stoppingSource = ibooker.book1D(histname+CatagoryName,
285  histname+CatagoryName,
286  StopReasonNameSize,
287  0., double(StopReasonNameSize));
288  stoppingSource->setAxisTitle("stopping reason",1);
289  stoppingSource->setAxisTitle("Number of Tracks",2);
290 
291  histname = "StoppingSourceVSeta_"+seedProducer.label() + "_";
292  stoppingSourceVSeta = ibooker.bookProfile(histname+CatagoryName,
293  histname+CatagoryName,
294  EtaBin,
295  EtaMin,
296  EtaMax,
297  2,
298  0., 2.);
299  stoppingSourceVSeta->setAxisTitle("track #eta",1);
300  stoppingSourceVSeta->setAxisTitle("fraction stopped",2);
301 
302  histname = "StoppingSourceVSphi_"+seedProducer.label() + "_";
303  stoppingSourceVSphi = ibooker.bookProfile(histname+CatagoryName,
304  histname+CatagoryName,
305  PhiBin,
306  PhiMin,
307  PhiMax,
308  2,
309  0., 2.);
310  stoppingSourceVSphi->setAxisTitle("track #phi",1);
311  stoppingSourceVSphi->setAxisTitle("fraction stopped",2);
312 
313  for (size_t ibin=0; ibin<StopReasonNameSize; ibin++) {
315  }
316  }
317 
318 
319 
320  // book the TrackCandidate histograms
321  // ---------------------------------------------------------------------------------//
322 
323  if (doTCPlots){
324 
325  ibooker.setCurrentFolder(MEFolderName+"/TrackBuilding");
326 
327  histname = "TrackCandPt_"+tcProducer.label() + "_";
328  TrackCandPt = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackPtBin, TrackPtMin, TrackPtMax);
329  TrackCandPt->setAxisTitle("Track Candidate p_{T} (GeV/c)", 1);
330  TrackCandPt->setAxisTitle("Number of Track Candidates", 2);
331 
332  histname = "TrackCandEta_"+tcProducer.label() + "_";
333  TrackCandEta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax);
334  TrackCandEta->setAxisTitle("Track Candidate #eta", 1);
335  TrackCandEta->setAxisTitle("Number of Track Candidates", 2);
336 
337  histname = "TrackCandPhi_"+tcProducer.label() + "_";
338  TrackCandPhi = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax);
339  TrackCandPhi->setAxisTitle("Track Candidate #phi", 1);
340  TrackCandPhi->setAxisTitle("Number of Track Candidates", 2);
341 
342  if (doTheta) {
343  histname = "TrackCandTheta_"+tcProducer.label() + "_";
344  TrackCandTheta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, ThetaBin, ThetaMin, ThetaMax);
345  TrackCandTheta->setAxisTitle("Track Candidate #theta", 1);
346  TrackCandTheta->setAxisTitle("Number of Track Candidates", 2);
347  }
348 
349  if (doAllTCPlots) {
350  histname = "TrackCandQ_"+tcProducer.label() + "_";
351  TrackCandQ = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackQBin, TrackQMin, TrackQMax);
352  TrackCandQ->setAxisTitle("Track Candidate Charge", 1);
353  TrackCandQ->setAxisTitle("Number of Track Candidates",2);
354 
355  histname = "TrackCandDxy_"+tcProducer.label() + "_";
356  TrackCandDxy = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TCDxyBin, TCDxyMin, TCDxyMax);
357  TrackCandDxy->setAxisTitle("Track Candidate d_{xy} (cm)", 1);
358  TrackCandDxy->setAxisTitle("Number of Track Candidates",2);
359 
360  histname = "TrackCandDz_"+tcProducer.label() + "_";
361  TrackCandDz = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TCDzBin, TCDzMin, TCDzMax);
362  TrackCandDz->setAxisTitle("Track Candidate d_{z} (cm)", 1);
363  TrackCandDz->setAxisTitle("Number of Track Candidates",2);
364 
365  histname = "NumberOfRecHitsPerTrackCand_"+tcProducer.label() + "_";
366  NumberOfRecHitsPerTrackCand = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TCHitBin, TCHitMin, TCHitMax);
367  NumberOfRecHitsPerTrackCand->setAxisTitle("Number of RecHits per Track Candidate", 1);
368  NumberOfRecHitsPerTrackCand->setAxisTitle("Number of Track Candidates",2);
369 
370  histname = "NumberOfRecHitsPerTrackCandVsPhiProfile_"+tcProducer.label() + "_";
371  NumberOfRecHitsPerTrackCandVsPhiProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax, TCHitBin, TCHitMin, TCHitMax,"s");
372  NumberOfRecHitsPerTrackCandVsPhiProfile->setAxisTitle("Track Candidate #phi",1);
373  NumberOfRecHitsPerTrackCandVsPhiProfile->setAxisTitle("Number of RecHits of each Track Candidate",2);
374 
375  histname = "NumberOfRecHitsPerTrackCandVsEtaProfile_"+tcProducer.label() + "_";
376  NumberOfRecHitsPerTrackCandVsEtaProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, TCHitBin, TCHitMin, TCHitMax,"s");
377  NumberOfRecHitsPerTrackCandVsEtaProfile->setAxisTitle("Track Candidate #eta",1);
378  NumberOfRecHitsPerTrackCandVsEtaProfile->setAxisTitle("Number of RecHits of each Track Candidate",2);
379  }
380 
381  histname = "TrackCandPhiVsEta_"+tcProducer.label() + "_";
382  TrackCandPhiVsEta = ibooker.book2D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, PhiBin, PhiMin, PhiMax);
383  TrackCandPhiVsEta->setAxisTitle("Track Candidate #eta", 1);
384  TrackCandPhiVsEta->setAxisTitle("Track Candidate #phi", 2);
385 
386  if(doAllTCPlots || doMVAPlots) {
387  for(size_t i=1, end=mvaProducers.size(); i<=end; ++i) {
388  auto num = std::to_string(i);
389  std::string pfix;
390 
391  if(i == 1) {
392  trackMVAsHP.push_back(nullptr);
393  trackMVAsHPVsPtProfile.push_back(nullptr);
394  trackMVAsHPVsEtaProfile.push_back(nullptr);
395  }
396  else {
397  pfix = " (not loose-selected)";
398  std::string pfix2 = " (not HP-selected)";
399  histname = "TrackMVA"+num+"HP_"+tcProducer.label() + "_";
400  trackMVAsHP.push_back(ibooker.book1D(histname+CatagoryName, histname+CatagoryName+pfix2, MVABin, MVAMin, MVAMax));
401  trackMVAsHP.back()->setAxisTitle("Track selection MVA"+num, 1);
402  trackMVAsHP.back()->setAxisTitle("Number of tracks", 2);
403 
404  histname = "TrackMVA"+num+"HPVsPtProfile_"+tcProducer.label() + "_";
405  trackMVAsHPVsPtProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix2, TrackPtBin, TrackPtMin, TrackPtMax, MVABin, MVAMin, MVAMax));
406  trackMVAsHPVsPtProfile.back()->setAxisTitle("Track p_{T} (GeV/c)", 1);
407  trackMVAsHPVsPtProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
408 
409  histname = "TrackMVA"+num+"HPVsEtaProfile_"+tcProducer.label() + "_";
410  trackMVAsHPVsEtaProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix2, EtaBin, EtaMin, EtaMax, MVABin, MVAMin, MVAMax));
411  trackMVAsHPVsEtaProfile.back()->setAxisTitle("Track #eta", 1);
412  trackMVAsHPVsEtaProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
413  }
414 
415  histname = "TrackMVA"+num+"_"+tcProducer.label() + "_";
416  trackMVAs.push_back(ibooker.book1D(histname+CatagoryName, histname+CatagoryName+pfix, MVABin, MVAMin, MVAMax));
417  trackMVAs.back()->setAxisTitle("Track selection MVA"+num, 1);
418  trackMVAs.back()->setAxisTitle("Number of tracks", 2);
419 
420  histname = "TrackMVA"+num+"VsPtProfile_"+tcProducer.label() + "_";
421  trackMVAsVsPtProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix, TrackPtBin, TrackPtMin, TrackPtMax, MVABin, MVAMin, MVAMax));
422  trackMVAsVsPtProfile.back()->setAxisTitle("Track p_{T} (GeV/c)", 1);
423  trackMVAsVsPtProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
424 
425  histname = "TrackMVA"+num+"VsEtaProfile_"+tcProducer.label() + "_";
426  trackMVAsVsEtaProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix, EtaBin, EtaMin, EtaMax, MVABin, MVAMin, MVAMax));
427  trackMVAsVsEtaProfile.back()->setAxisTitle("Track #eta", 1);
428  trackMVAsVsEtaProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
429  }
430  }
431  }
432 
433 }
434 
435 // -- Analyse
436 // ---------------------------------------------------------------------------------//
438 (
439  const edm::Event& iEvent,
440  const edm::EventSetup& iSetup,
441  const TrajectorySeed& candidate,
442  const SeedStopInfo& stopInfo,
443  const reco::BeamSpot& bs,
444  const edm::ESHandle<MagneticField>& theMF,
446 )
447 {
448  TSCBLBuilderNoMaterial tscblBuilder;
449 
450  //get parameters and errors from the candidate state
451  auto const & theG = ((TkTransientTrackingRecHitBuilder const *)(theTTRHBuilder.product()))->geometry();
452  auto const & candSS = candidate.startingState();
453  TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( candSS, &(theG->idToDet(candSS.detId())->surface()), theMF.product());
454  TrajectoryStateClosestToBeamLine tsAtClosestApproachSeed = tscblBuilder(*state.freeState(),bs);//as in TrackProducerAlgorithm
455  if(!(tsAtClosestApproachSeed.isValid())) {
456  edm::LogVerbatim("TrackBuilding") << "TrajectoryStateClosestToBeamLine not valid";
457  return;
458  }
459  GlobalPoint v0 = tsAtClosestApproachSeed.trackStateAtPCA().position();
460  GlobalVector p = tsAtClosestApproachSeed.trackStateAtPCA().momentum();
461  GlobalPoint v(v0.x()-bs.x0(),v0.y()-bs.y0(),v0.z()-bs.z0());
462 
463  double pt = sqrt(state.globalMomentum().perp2());
464  double eta = state.globalPosition().eta();
465  double phi = state.globalPosition().phi();
466  double theta = state.globalPosition().theta();
467  //double pm = sqrt(state.globalMomentum().mag2());
468  //double pz = state.globalMomentum().z();
469  //double qoverp = tsAtClosestApproachSeed.trackStateAtPCA().charge()/p.mag();
470  //double theta = p.theta();
471  //double lambda = M_PI/2-p.theta();
472  double numberOfHits = candidate.recHits().second-candidate.recHits().first;
473  double dxy = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
474  double dz = v.z() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.perp();
475 
476  // fill the ME's
477  if (doAllSeedPlots || doQ)SeedQ->Fill( state.charge() );
478  if (doAllSeedPlots || doPT) SeedPt->Fill( pt );
479  if (doAllSeedPlots || doETA) SeedEta->Fill( eta );
480  if (doAllSeedPlots || doPHI) SeedPhi->Fill( phi );
481  if (doAllSeedPlots || doPHIVsETA) SeedPhiVsEta->Fill( eta, phi);
482  if (doAllSeedPlots || doTheta) SeedTheta->Fill( theta );
483  if (doAllSeedPlots || doDxy) SeedDxy->Fill( dxy );
484  if (doAllSeedPlots || doDz) SeedDz->Fill( dz );
485  if (doAllSeedPlots || doNRecHits) NumberOfRecHitsPerSeed->Fill( numberOfHits );
486  if (doAllSeedPlots || doProfETA) NumberOfRecHitsPerSeedVsEtaProfile->Fill( eta, numberOfHits );
487  if (doAllSeedPlots || doProfPHI) NumberOfRecHitsPerSeedVsPhiProfile->Fill( phi, numberOfHits );
488  if (doAllSeedPlots || doStopSource) {
489  const double stopped = stopInfo.stopReason() == SeedStopReason::NOT_STOPPED ? 0. : 1.;
490  seedStoppingSource->Fill(stopInfo.stopReasonUC());
491  seedStoppingSourceVsPhi->Fill(phi, stopped);
492  seedStoppingSourceVsEta->Fill(eta, stopped);
493 
494  const auto ncands = stopInfo.candidatesPerSeed();
496  numberOfTrajCandsPerSeedVsPhi->Fill(phi, ncands);
497  numberOfTrajCandsPerSeedVsEta->Fill(eta, ncands);
498 
500  }
501 }
502 
503 // -- Analyse
504 // ---------------------------------------------------------------------------------//
506 (
507  const edm::Event& iEvent,
508  const edm::EventSetup& iSetup,
509  const TrackCandidate& candidate,
510  const reco::BeamSpot& bs,
511  const edm::ESHandle<MagneticField>& theMF,
513 )
514 {
515  TSCBLBuilderNoMaterial tscblBuilder;
516 
517  //get parameters and errors from the candidate state
518  auto const & theG = ((TkTransientTrackingRecHitBuilder const *)(theTTRHBuilder.product()))->geometry();
519  auto const & candSS = candidate.trajectoryStateOnDet();
520  TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( candSS, &(theG->idToDet(candSS.detId())->surface()), theMF.product());
521  TrajectoryStateClosestToBeamLine tsAtClosestApproachTrackCand = tscblBuilder(*state.freeState(),bs);//as in TrackProducerAlgorithm
522  if(!(tsAtClosestApproachTrackCand.isValid())) {
523  edm::LogVerbatim("TrackBuilding") << "TrajectoryStateClosestToBeamLine not valid";
524  return;
525  }
526  GlobalPoint v0 = tsAtClosestApproachTrackCand.trackStateAtPCA().position();
527  GlobalVector p = tsAtClosestApproachTrackCand.trackStateAtPCA().momentum();
528  GlobalPoint v(v0.x()-bs.x0(),v0.y()-bs.y0(),v0.z()-bs.z0());
529 
530  double pt = sqrt(state.globalMomentum().perp2());
531  double eta = state.globalPosition().eta();
532  double phi = state.globalPosition().phi();
533  double theta = state.globalPosition().theta();
534  //double pm = sqrt(state.globalMomentum().mag2());
535  //double pz = state.globalMomentum().z();
536  //double qoverp = tsAtClosestApproachTrackCand.trackStateAtPCA().charge()/p.mag();
537  //double theta = p.theta();
538  //double lambda = M_PI/2-p.theta();
539  double numberOfHits = candidate.recHits().second-candidate.recHits().first;
540  double dxy = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
541 
542  double dz = v.z() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.perp();
543 
544  if (doAllTCPlots || doStopSource) {
545  // stopping source
546  int max = stoppingSource->getNbinsX();
547  double stop = candidate.stopReason() > max ? double(max-1) : static_cast<double>(candidate.stopReason());
548  double stopped = int(StopReason::NOT_STOPPED)==candidate.stopReason() ? 0. : 1.;
549  stoppingSource ->Fill(stop);
550  stoppingSourceVSeta ->Fill(eta,stopped);
551  stoppingSourceVSphi ->Fill(phi,stopped);
552  }
553 
554  if (doTCPlots){
555  // fill the ME's
556  if (doAllTCPlots) TrackCandQ->Fill( state.charge() );
557  TrackCandPt->Fill( pt );
558  TrackCandEta->Fill( eta );
559  TrackCandPhi->Fill( phi );
560  TrackCandPhiVsEta->Fill( eta, phi );
561  if (doTheta) TrackCandTheta->Fill( theta );
562  if (doAllTCPlots) TrackCandDxy->Fill( dxy );
563  if (doAllTCPlots) TrackCandDz->Fill( dz );
564  if (doAllTCPlots) NumberOfRecHitsPerTrackCand->Fill( numberOfHits );
565  if (doAllTCPlots) NumberOfRecHitsPerTrackCandVsEtaProfile->Fill( eta, numberOfHits );
566  if (doAllTCPlots) NumberOfRecHitsPerTrackCandVsPhiProfile->Fill( phi, numberOfHits );
567  }
568 }
569 
570 namespace {
571  bool trackSelected(unsigned char mask, unsigned char qual) {
572  return mask & 1<<qual;
573  }
574 }
576  const std::vector<const MVACollection *>& mvaCollections,
577  const std::vector<const QualityMaskCollection *>& qualityMaskCollections) {
578  if(!(doAllTCPlots || doMVAPlots))
579  return;
580  if(trackCollection.empty())
581  return;
582 
583  const auto ntracks = trackCollection.size();
584  const auto nmva = mvaCollections.size();
585  for(const auto mva: mvaCollections) {
586  if(mva->size() != ntracks) {
587  edm::LogError("LogicError") << "TrackBuildingAnalyzer: Incompatible size of MVACollection, " << mva->size() << " differs from the size of the track collection " << ntracks;
588  return;
589  }
590  }
591  for(const auto qual: qualityMaskCollections) {
592  if(qual->size() != ntracks) {
593  edm::LogError("LogicError") << "TrackBuildingAnalyzer: Incompatible size of QualityMaskCollection, " << qual->size() << " differs from the size of the track collection " << ntracks;
594  return;
595  }
596  }
597 
598 
599  for(size_t iTrack=0; iTrack<ntracks; ++iTrack) {
600  // Fill MVA1 histos with all tracks, MVA2 histos only with tracks
601  // not selected by MVA1 etc
602  bool selectedLoose = false;
603  bool selectedHP = false;
604 
605  const auto pt = trackCollection[iTrack].pt();
606  const auto eta = trackCollection[iTrack].eta();
607 
608  for(size_t iMVA=0; iMVA<nmva; ++iMVA) {
609  const auto mva = (*(mvaCollections[iMVA]))[iTrack];
610  if(!selectedLoose) {
611  trackMVAs[iMVA]->Fill(mva);
612  trackMVAsVsPtProfile[iMVA]->Fill(pt, mva);
613  trackMVAsVsEtaProfile[iMVA]->Fill(eta, mva);
614  }
615  if(iMVA >= 1 && !selectedHP) {
616  trackMVAsHP[iMVA]->Fill(mva);
617  trackMVAsHPVsPtProfile[iMVA]->Fill(pt, mva);
618  trackMVAsHPVsEtaProfile[iMVA]->Fill(eta, mva);
619  }
620 
621  const auto qual = (*(qualityMaskCollections)[iMVA])[iTrack];
622  selectedLoose |= trackSelected(qual, reco::TrackBase::loose);
623  selectedHP |= trackSelected(qual, reco::TrackBase::highPurity);
624 
625  if(selectedLoose && selectedHP)
626  break;
627  }
628  }
629 }
630 
632  if(!doRegionPlots)
633  return;
634 
635  for(const auto& candidate: regionCandidates) {
636  const auto eta = candidate.eta();
637  const auto phi = candidate.phi();
638  if (doAllSeedPlots || doPT) TrackingRegionCandidatePt->Fill( candidate.pt() );
642  }
643 }
MonitorElement * NumberOfRecHitsPerTrackCandVsPhiProfile
T getParameter(std::string const &) const
MonitorElement * seedStoppingSourceVsPhi
double z0() const
z coordinate
Definition: BeamSpot.h:68
MonitorElement * numberOfTrajCandsPerSeed
MonitorElement * TrackingRegionCandidatePhiVsEta
void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup, const TrajectorySeed &seed, const SeedStopInfo &stopInfo, const reco::BeamSpot &bs, const edm::ESHandle< MagneticField > &theMF, const edm::ESHandle< TransientTrackingRecHitBuilder > &theTTRHBuilder)
T perp() const
Definition: PV3DBase.h:72
SeedStopReason stopReason() const
Definition: SeedStopInfo.h:15
range recHits() const
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
std::vector< MonitorElement * > trackMVAsHPVsEtaProfile
std::vector< MonitorElement * > trackMVAsVsPtProfile
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Geom::Theta< T > theta() const
T y() const
Definition: PV3DBase.h:63
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
MonitorElement * stoppingSourceVSphi
GlobalPoint globalPosition() const
size_type size() const
std::vector< MonitorElement * > trackMVAsVsEtaProfile
T perp2() const
Definition: PV3DBase.h:71
const double EtaMax[kNumberCalorimeter]
static const std::string StopReasonName[]
void initHisto(DQMStore::IBooker &ibooker, const edm::ParameterSet &)
MonitorElement * TrackCandPhiVsEta
MonitorElement * TrackingRegionCandidatePt
unsigned short candidatesPerSeed() const
Definition: SeedStopInfo.h:12
static const std::string SeedStopReasonName[]
void Fill(long long x)
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
MonitorElement * NumberOfRecHitsPerTrackCandVsEtaProfile
MonitorElement * stoppingSource
PTrajectoryStateOnDet const & trajectoryStateOnDet() const
MonitorElement * NumberOfRecHitsPerSeed
int iEvent
Definition: GenABIO.cc:230
MonitorElement * NumberOfRecHitsPerSeedVsPhiProfile
T sqrt(T t)
Definition: SSEVec.h:18
FreeTrajectoryState const * freeState(bool withErrors=true) const
T z() const
Definition: PV3DBase.h:64
std::vector< MonitorElement * > trackMVAsHP
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
bool empty() const
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
const double EtaMin[kNumberCalorimeter]
#define end
Definition: vmac.h:37
MonitorElement * numberOfTrajCandsPerSeedVsPhi
unsigned char stopReasonUC() const
Definition: SeedStopInfo.h:16
TrackBuildingAnalyzer(const edm::ParameterSet &)
MonitorElement * numberOfTrajCandsPerSeedVsEta
std::vector< MonitorElement * > trackMVAs
PTrajectoryStateOnDet const & startingState() const
MonitorElement * stoppingSourceVSeta
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:277
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field)
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
range recHits() const
MonitorElement * TrackCandTheta
MonitorElement * SeedPhiVsEta
T eta() const
Definition: PV3DBase.h:76
std::vector< MonitorElement * > trackMVAsHPVsPtProfile
MonitorElement * TrackingRegionCandidateEta
GlobalVector globalMomentum() const
MonitorElement * seedStoppingSourceVsEta
MonitorElement * TrackingRegionCandidatePhi
int getNbinsX(void) const
get # of bins in X-axis
double y0() const
y coordinate
Definition: BeamSpot.h:66
MonitorElement * seedStoppingSourceVsNumberOfTrajCandsPerSeed
MonitorElement * NumberOfRecHitsPerSeedVsEtaProfile
T x() const
Definition: PV3DBase.h:62
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
T const * product() const
Definition: ESHandle.h:86
MonitorElement * NumberOfRecHitsPerTrackCand
uint8_t stopReason() const
MonitorElement * seedStoppingSource
double x0() const
x coordinate
Definition: BeamSpot.h:64