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