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