CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackBuildingAnalyzer.cc
Go to the documentation of this file.
10 
15 #include <string>
16 #include "TMath.h"
17 
18 #include <iostream>
19 
21  : conf_( 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 {
34 }
35 
37 {
38 }
39 
41 {
42 
43  // parameters from the configuration
45  std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
46 
47  // std::cout << "[TrackBuildingAnalyzer::beginRun] AlgoName: " << AlgoName << std::endl;
48 
49  // use the AlgoName and Quality Name
50  std::string CatagoryName = AlgoName;
51 
52  // get binning from the configuration
53  int TrackPtBin = conf_.getParameter<int>( "TrackPtBin");
54  double TrackPtMin = conf_.getParameter<double>("TrackPtMin");
55  double TrackPtMax = conf_.getParameter<double>("TrackPtMax");
56 
57  int PhiBin = conf_.getParameter<int>( "PhiBin");
58  double PhiMin = conf_.getParameter<double>("PhiMin");
59  double PhiMax = conf_.getParameter<double>("PhiMax");
60 
61  int EtaBin = conf_.getParameter<int>( "EtaBin");
62  double EtaMin = conf_.getParameter<double>("EtaMin");
63  double EtaMax = conf_.getParameter<double>("EtaMax");
64 
65  int ThetaBin = conf_.getParameter<int>( "ThetaBin");
66  double ThetaMin = conf_.getParameter<double>("ThetaMin");
67  double ThetaMax = conf_.getParameter<double>("ThetaMax");
68 
69  int TrackQBin = conf_.getParameter<int>( "TrackQBin");
70  double TrackQMin = conf_.getParameter<double>("TrackQMin");
71  double TrackQMax = conf_.getParameter<double>("TrackQMax");
72 
73  int SeedDxyBin = conf_.getParameter<int>( "SeedDxyBin");
74  double SeedDxyMin = conf_.getParameter<double>("SeedDxyMin");
75  double SeedDxyMax = conf_.getParameter<double>("SeedDxyMax");
76 
77  int SeedDzBin = conf_.getParameter<int>( "SeedDzBin");
78  double SeedDzMin = conf_.getParameter<double>("SeedDzMin");
79  double SeedDzMax = conf_.getParameter<double>("SeedDzMax");
80 
81  int SeedHitBin = conf_.getParameter<int>( "SeedHitBin");
82  double SeedHitMin = conf_.getParameter<double>("SeedHitMin");
83  double SeedHitMax = conf_.getParameter<double>("SeedHitMax");
84 
85  int TCDxyBin = conf_.getParameter<int>( "TCDxyBin");
86  double TCDxyMin = conf_.getParameter<double>("TCDxyMin");
87  double TCDxyMax = conf_.getParameter<double>("TCDxyMax");
88 
89  int TCDzBin = conf_.getParameter<int>( "TCDzBin");
90  double TCDzMin = conf_.getParameter<double>("TCDzMin");
91  double TCDzMax = conf_.getParameter<double>("TCDzMax");
92 
93  int TCHitBin = conf_.getParameter<int>( "TCHitBin");
94  double TCHitMin = conf_.getParameter<double>("TCHitMin");
95  double TCHitMax = conf_.getParameter<double>("TCHitMax");
96 
97 
98  edm::InputTag seedProducer = conf_.getParameter<edm::InputTag>("SeedProducer");
99  edm::InputTag tcProducer = conf_.getParameter<edm::InputTag>("TCProducer");
100 
101  doAllPlots = conf_.getParameter<bool>("doAllPlots");
102  doAllSeedPlots = conf_.getParameter<bool>("doSeedParameterHistos");
103  doTCPlots = conf_.getParameter<bool>("doTrackCandHistos");
104  doPT = conf_.getParameter<bool>("doSeedPTHisto");
105  doETA = conf_.getParameter<bool>("doSeedETAHisto");
106  doPHI = conf_.getParameter<bool>("doSeedPHIHisto");
107  doPHIVsETA = conf_.getParameter<bool>("doSeedPHIVsETAHisto");
108  doTheta = conf_.getParameter<bool>("doSeedThetaHisto");
109  doQ = conf_.getParameter<bool>("doSeedQHisto");
110  doDxy = conf_.getParameter<bool>("doSeedDxyHisto");
111  doDz = conf_.getParameter<bool>("doSeedDzHisto");
112  doNRecHits = conf_.getParameter<bool>("doSeedNRecHitsHisto");
113  doProfPHI = conf_.getParameter<bool>("doSeedNVsPhiProf");
114  doProfETA = conf_.getParameter<bool>("doSeedNVsEtaProf");
115 
116  // if (doAllPlots){doAllSeedPlots=true; doTCPlots=true;}
117 
118  ibooker.setCurrentFolder(MEFolderName);
119 
120  // book the Seed histograms
121  // ---------------------------------------------------------------------------------//
122  // std::cout << "[TrackBuildingAnalyzer::beginRun] MEFolderName: " << MEFolderName << std::endl;
123  ibooker.setCurrentFolder(MEFolderName+"/TrackBuilding");
124 
125  if (doAllSeedPlots || doPT) {
126  histname = "SeedPt_"+seedProducer.label() + "_";
127  SeedPt = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackPtBin, TrackPtMin, TrackPtMax);
128  SeedPt->setAxisTitle("Seed p_{T} (GeV/c)", 1);
129  SeedPt->setAxisTitle("Number of Seeds", 2);
130  }
131 
132  if (doAllSeedPlots || doETA) {
133  histname = "SeedEta_"+seedProducer.label() + "_";
134  SeedEta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax);
135  SeedEta->setAxisTitle("Seed #eta", 1);
136  SeedEta->setAxisTitle("Number of Seeds", 2);
137  }
138 
139  if (doAllSeedPlots || doPHI) {
140  histname = "SeedPhi_"+seedProducer.label() + "_";
141  SeedPhi = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax);
142  SeedPhi->setAxisTitle("Seed #phi", 1);
143  SeedPhi->setAxisTitle("Number of Seed", 2);
144  }
145 
146  if (doAllSeedPlots || doPHIVsETA) {
147  histname = "SeedPhiVsEta_"+seedProducer.label() + "_";
148  SeedPhiVsEta = ibooker.book2D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, PhiBin, PhiMin, PhiMax);
149  SeedPhiVsEta->setAxisTitle("Seed #eta", 1);
150  SeedPhiVsEta->setAxisTitle("Seed #phi", 2);
151  }
152 
153  if (doAllSeedPlots || doTheta){
154  histname = "SeedTheta_"+seedProducer.label() + "_";
155  SeedTheta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, ThetaBin, ThetaMin, ThetaMax);
156  SeedTheta->setAxisTitle("Seed #theta", 1);
157  SeedTheta->setAxisTitle("Number of Seeds", 2);
158  }
159 
160  if (doAllSeedPlots || doQ){
161  histname = "SeedQ_"+seedProducer.label() + "_";
162  SeedQ = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackQBin, TrackQMin, TrackQMax);
163  SeedQ->setAxisTitle("Seed Charge", 1);
164  SeedQ->setAxisTitle("Number of Seeds",2);
165  }
166 
167  if (doAllSeedPlots || doDxy){
168  histname = "SeedDxy_"+seedProducer.label() + "_";
169  SeedDxy = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, SeedDxyBin, SeedDxyMin, SeedDxyMax);
170  SeedDxy->setAxisTitle("Seed d_{xy} (cm)", 1);
171  SeedDxy->setAxisTitle("Number of Seeds",2);
172  }
173 
174  if (doAllSeedPlots || doDz){
175  histname = "SeedDz_"+seedProducer.label() + "_";
176  SeedDz = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, SeedDzBin, SeedDzMin, SeedDzMax);
177  SeedDz->setAxisTitle("Seed d_{z} (cm)", 1);
178  SeedDz->setAxisTitle("Number of Seeds",2);
179  }
180 
181  if (doAllSeedPlots || doNRecHits){
182  histname = "NumberOfRecHitsPerSeed_"+seedProducer.label() + "_";
183  NumberOfRecHitsPerSeed = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, SeedHitBin, SeedHitMin, SeedHitMax);
184  NumberOfRecHitsPerSeed->setAxisTitle("Number of RecHits per Seed", 1);
185  NumberOfRecHitsPerSeed->setAxisTitle("Number of Seeds",2);
186  }
187 
188  if (doAllSeedPlots || doProfPHI){
189  histname = "NumberOfRecHitsPerSeedVsPhiProfile_"+seedProducer.label() + "_";
190  NumberOfRecHitsPerSeedVsPhiProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax, SeedHitBin, SeedHitMin, SeedHitMax,"s");
192  NumberOfRecHitsPerSeedVsPhiProfile->setAxisTitle("Number of RecHits of each Seed",2);
193  }
194 
195  if (doAllSeedPlots || doProfETA){
196  histname = "NumberOfRecHitsPerSeedVsEtaProfile_"+seedProducer.label() + "_";
197  NumberOfRecHitsPerSeedVsEtaProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, SeedHitBin, SeedHitMin, SeedHitMax,"s");
199  NumberOfRecHitsPerSeedVsEtaProfile->setAxisTitle("Number of RecHits of each Seed",2);
200  }
201 
202  // book the TrackCandidate histograms
203  // ---------------------------------------------------------------------------------//
204 
205  if (doTCPlots){
206 
207  ibooker.setCurrentFolder(MEFolderName+"/TrackBuilding");
208 
209  histname = "TrackCandPt_"+tcProducer.label() + "_";
210  TrackCandPt = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackPtBin, TrackPtMin, TrackPtMax);
211  TrackCandPt->setAxisTitle("Track Candidate p_{T} (GeV/c)", 1);
212  TrackCandPt->setAxisTitle("Number of Track Candidates", 2);
213 
214  histname = "TrackCandEta_"+tcProducer.label() + "_";
215  TrackCandEta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax);
216  TrackCandEta->setAxisTitle("Track Candidate #eta", 1);
217  TrackCandEta->setAxisTitle("Number of Track Candidates", 2);
218 
219  histname = "TrackCandPhi_"+tcProducer.label() + "_";
220  TrackCandPhi = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax);
221  TrackCandPhi->setAxisTitle("Track Candidate #phi", 1);
222  TrackCandPhi->setAxisTitle("Number of Track Candidates", 2);
223 
224  histname = "TrackCandTheta_"+tcProducer.label() + "_";
225  TrackCandTheta = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, ThetaBin, ThetaMin, ThetaMax);
226  TrackCandTheta->setAxisTitle("Track Candidate #theta", 1);
227  TrackCandTheta->setAxisTitle("Number of Track Candidates", 2);
228 
229  histname = "TrackCandQ_"+tcProducer.label() + "_";
230  TrackCandQ = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TrackQBin, TrackQMin, TrackQMax);
231  TrackCandQ->setAxisTitle("Track Candidate Charge", 1);
232  TrackCandQ->setAxisTitle("Number of Track Candidates",2);
233 
234  histname = "TrackCandDxy_"+tcProducer.label() + "_";
235  TrackCandDxy = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TCDxyBin, TCDxyMin, TCDxyMax);
236  TrackCandDxy->setAxisTitle("Track Candidate d_{xy} (cm)", 1);
237  TrackCandDxy->setAxisTitle("Number of Track Candidates",2);
238 
239  histname = "TrackCandDz_"+tcProducer.label() + "_";
240  TrackCandDz = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TCDzBin, TCDzMin, TCDzMax);
241  TrackCandDz->setAxisTitle("Track Candidate d_{z} (cm)", 1);
242  TrackCandDz->setAxisTitle("Number of Track Candidates",2);
243 
244  histname = "NumberOfRecHitsPerTrackCand_"+tcProducer.label() + "_";
245  NumberOfRecHitsPerTrackCand = ibooker.book1D(histname+CatagoryName, histname+CatagoryName, TCHitBin, TCHitMin, TCHitMax);
246  NumberOfRecHitsPerTrackCand->setAxisTitle("Number of RecHits per Track Candidate", 1);
247  NumberOfRecHitsPerTrackCand->setAxisTitle("Number of Track Candidates",2);
248 
249  histname = "NumberOfRecHitsPerTrackCandVsPhiProfile_"+tcProducer.label() + "_";
250  NumberOfRecHitsPerTrackCandVsPhiProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, PhiBin, PhiMin, PhiMax, TCHitBin, TCHitMin, TCHitMax,"s");
251  NumberOfRecHitsPerTrackCandVsPhiProfile->setAxisTitle("Track Candidate #phi",1);
252  NumberOfRecHitsPerTrackCandVsPhiProfile->setAxisTitle("Number of RecHits of each Track Candidate",2);
253 
254  histname = "NumberOfRecHitsPerTrackCandVsEtaProfile_"+tcProducer.label() + "_";
255  NumberOfRecHitsPerTrackCandVsEtaProfile = ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, TCHitBin, TCHitMin, TCHitMax,"s");
256  NumberOfRecHitsPerTrackCandVsEtaProfile->setAxisTitle("Track Candidate #eta",1);
257  NumberOfRecHitsPerTrackCandVsEtaProfile->setAxisTitle("Number of RecHits of each Track Candidate",2);
258  }
259 
260 }
261 
262 // -- Analyse
263 // ---------------------------------------------------------------------------------//
265 (
266  const edm::Event& iEvent,
267  const edm::EventSetup& iSetup,
268  const TrajectorySeed& candidate,
269  const reco::BeamSpot& bs,
270  const edm::ESHandle<MagneticField>& theMF,
272 )
273 {
274  TSCBLBuilderNoMaterial tscblBuilder;
275 
276  //get parameters and errors from the candidate state
277  TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(&*(candidate.recHits().second-1));
278  TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( candidate.startingState(), recHit->surface(), theMF.product());
279  TrajectoryStateClosestToBeamLine tsAtClosestApproachSeed = tscblBuilder(*state.freeState(),bs);//as in TrackProducerAlgorithm
280  if(!(tsAtClosestApproachSeed.isValid())) {
281  edm::LogVerbatim("TrackBuilding") << "TrajectoryStateClosestToBeamLine not valid";
282  return;
283  }
284  GlobalPoint v0 = tsAtClosestApproachSeed.trackStateAtPCA().position();
285  GlobalVector p = tsAtClosestApproachSeed.trackStateAtPCA().momentum();
286  GlobalPoint v(v0.x()-bs.x0(),v0.y()-bs.y0(),v0.z()-bs.z0());
287 
288  double pt = sqrt(state.globalMomentum().perp2());
289  double eta = state.globalMomentum().eta();
290  double phi = state.globalMomentum().phi();
291  double theta = state.globalMomentum().theta();
292  //double pm = sqrt(state.globalMomentum().mag2());
293  //double pz = state.globalMomentum().z();
294  //double qoverp = tsAtClosestApproachSeed.trackStateAtPCA().charge()/p.mag();
295  //double theta = p.theta();
296  //double lambda = M_PI/2-p.theta();
297  double numberOfHits = candidate.recHits().second-candidate.recHits().first;
298  double dxy = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
299  double dz = v.z() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.perp();
300 
301  // fill the ME's
302  if (doAllSeedPlots || doQ)SeedQ->Fill( state.charge() );
303  if (doAllSeedPlots || doPT) SeedPt->Fill( pt );
304  if (doAllSeedPlots || doETA) SeedEta->Fill( eta );
305  if (doAllSeedPlots || doPHI) SeedPhi->Fill( phi );
306  if (doAllSeedPlots || doPHIVsETA) SeedPhiVsEta->Fill( eta, phi);
307  if (doAllSeedPlots || doTheta) SeedTheta->Fill( theta );
308  if (doAllSeedPlots || doDxy) SeedDxy->Fill( dxy );
309  if (doAllSeedPlots || doDz) SeedDz->Fill( dz );
310  if (doAllSeedPlots || doNRecHits) NumberOfRecHitsPerSeed->Fill( numberOfHits );
311  if (doAllSeedPlots || doProfETA) NumberOfRecHitsPerSeedVsEtaProfile->Fill( eta, numberOfHits );
312  if (doAllSeedPlots || doProfPHI) NumberOfRecHitsPerSeedVsPhiProfile->Fill( phi, numberOfHits );
313 
314 }
315 
316 // -- Analyse
317 // ---------------------------------------------------------------------------------//
319 (
320  const edm::Event& iEvent,
321  const edm::EventSetup& iSetup,
322  const TrackCandidate& candidate,
323  const reco::BeamSpot& bs,
324  const edm::ESHandle<MagneticField>& theMF,
326 )
327 {
328  TSCBLBuilderNoMaterial tscblBuilder;
329 
330  //get parameters and errors from the candidate state
331  TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(&*(candidate.recHits().second-1));
332  TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( candidate.trajectoryStateOnDet(), recHit->surface(), theMF.product());
333  TrajectoryStateClosestToBeamLine tsAtClosestApproachTrackCand = tscblBuilder(*state.freeState(),bs);//as in TrackProducerAlgorithm
334  if(!(tsAtClosestApproachTrackCand.isValid())) {
335  edm::LogVerbatim("TrackBuilding") << "TrajectoryStateClosestToBeamLine not valid";
336  return;
337  }
338  GlobalPoint v0 = tsAtClosestApproachTrackCand.trackStateAtPCA().position();
339  GlobalVector p = tsAtClosestApproachTrackCand.trackStateAtPCA().momentum();
340  GlobalPoint v(v0.x()-bs.x0(),v0.y()-bs.y0(),v0.z()-bs.z0());
341 
342  double pt = sqrt(state.globalMomentum().perp2());
343  double eta = state.globalMomentum().eta();
344  double phi = state.globalMomentum().phi();
345  double theta = state.globalMomentum().theta();
346  //double pm = sqrt(state.globalMomentum().mag2());
347  //double pz = state.globalMomentum().z();
348  //double qoverp = tsAtClosestApproachTrackCand.trackStateAtPCA().charge()/p.mag();
349  //double theta = p.theta();
350  //double lambda = M_PI/2-p.theta();
351  double numberOfHits = candidate.recHits().second-candidate.recHits().first;
352  double dxy = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
353 
354  double dz = v.z() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.perp();
355 
356  if (doTCPlots){
357  // fill the ME's
358  TrackCandQ->Fill( state.charge() );
359  TrackCandPt->Fill( pt );
360  TrackCandEta->Fill( eta );
361  TrackCandPhi->Fill( phi );
362  TrackCandTheta->Fill( theta );
363  TrackCandDxy->Fill( dxy );
364  TrackCandDz->Fill( dz );
365  NumberOfRecHitsPerTrackCand->Fill( numberOfHits );
366  NumberOfRecHitsPerTrackCandVsEtaProfile->Fill( eta, numberOfHits );
367  NumberOfRecHitsPerTrackCandVsPhiProfile->Fill( phi, numberOfHits );
368  }
369 }
MonitorElement * NumberOfRecHitsPerTrackCandVsPhiProfile
MonitorElement * TrackCandPt
T getParameter(std::string const &) const
double z0() const
z coordinate
Definition: BeamSpot.h:68
T perp() const
Definition: PV3DBase.h:72
range recHits() const
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
MonitorElement * TrackCandDz
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
#define NULL
Definition: scimark2.h:8
T perp2() const
Definition: PV3DBase.h:71
MonitorElement * TrackCandPhi
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
MonitorElement * NumberOfRecHitsPerTrackCandVsEtaProfile
PTrajectoryStateOnDet const & trajectoryStateOnDet() const
MonitorElement * NumberOfRecHitsPerSeed
virtual 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)
int iEvent
Definition: GenABIO.cc:230
MonitorElement * NumberOfRecHitsPerSeedVsPhiProfile
MonitorElement * TrackCandDxy
T sqrt(T t)
Definition: SSEVec.h:48
FreeTrajectoryState const * freeState(bool withErrors=true) const
T z() const
Definition: PV3DBase.h:64
virtual void initHisto(DQMStore::IBooker &ibooker)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
MonitorElement * TrackCandEta
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
const double EtaMin[kNumberCalorimeter]
TrackBuildingAnalyzer(const edm::ParameterSet &)
std::shared_ptr< TrackingRecHit const > RecHitPointer
PTrajectoryStateOnDet const & startingState() const
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field)
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
T const * product() const
Definition: ESHandle.h:86
range recHits() const
MonitorElement * TrackCandTheta
MonitorElement * SeedPhiVsEta
Geom::Phi< T > phi() const
T eta() const
Definition: PV3DBase.h:76
GlobalVector globalMomentum() const
double y0() const
y coordinate
Definition: BeamSpot.h:66
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)
MonitorElement * NumberOfRecHitsPerTrackCand
double x0() const
x coordinate
Definition: BeamSpot.h:64