CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonTrackAnalyzer.cc
Go to the documentation of this file.
1 
8 
9 // Collaborating Class Header
14 
16 
20 
22 
24 
28 
31 
32 #include "TFile.h"
33 #include "TH1F.h"
34 #include "TH2F.h"
35 
36 using namespace std;
37 using namespace edm;
38 
41  // service parameters
42  pset = ps;
43  ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
44  // the services
45  theService = new MuonServiceProxy(serviceParameters);
46 
47  theSimTracksLabel = edm::InputTag("g4SimHits");
48  theSimTracksToken = consumes<edm::SimTrackContainer>(theSimTracksLabel);
49 
50  theTracksLabel = pset.getParameter<InputTag>("Tracks");
51  theTracksToken = consumes<reco::TrackCollection>(theTracksLabel);
52  doTracksAnalysis = pset.getUntrackedParameter<bool>("DoTracksAnalysis", true);
53 
54  doSeedsAnalysis = pset.getUntrackedParameter<bool>("DoSeedsAnalysis", false);
55  if (doSeedsAnalysis) {
56  theSeedsLabel = pset.getParameter<InputTag>("MuonSeed");
57  theSeedsToken = consumes<TrajectorySeedCollection>(theSeedsLabel);
58  ParameterSet updatorPar = pset.getParameter<ParameterSet>("MuonUpdatorAtVertexParameters");
59  theSeedPropagatorName = updatorPar.getParameter<string>("Propagator");
60 
61  theUpdator = new MuonUpdatorAtVertex(updatorPar, theService);
62  }
63 
64  theCSCSimHitLabel = pset.getParameter<InputTag>("CSCSimHit");
65  theDTSimHitLabel = pset.getParameter<InputTag>("DTSimHit");
66  theRPCSimHitLabel = pset.getParameter<InputTag>("RPCSimHit");
67  theCSCSimHitToken = consumes<std::vector<PSimHit> >(theCSCSimHitLabel);
68  theDTSimHitToken = consumes<std::vector<PSimHit> >(theDTSimHitLabel);
69  theRPCSimHitToken = consumes<std::vector<PSimHit> >(theRPCSimHitLabel);
70 
71  theEtaRange = (EtaRange)pset.getParameter<int>("EtaRange");
72 
73  // number of sim tracks
74  numberOfSimTracks = 0;
75  // number of reco tracks
76  numberOfRecTracks = 0;
77 
79  out = pset.getUntrackedParameter<string>("rootFileName");
80  dirName_ = pset.getUntrackedParameter<std::string>("dirName");
81  subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
82 }
83 
86  if (theService)
87  delete theService;
88 }
89 
91  //theFile->cd();
92 }
94  edm::Run const &iRun,
95  edm::EventSetup const & /* iSetup */) {
96  ibooker.cd();
97 
98  InputTag algo = theTracksLabel;
99  string dirName = dirName_;
100  if (!algo.process().empty())
101  dirName += algo.process() + "_";
102  if (!algo.label().empty())
103  dirName += algo.label() + "_";
104  if (!algo.instance().empty())
105  dirName += algo.instance() + "";
106  if (dirName.find("Tracks") < dirName.length()) {
107  dirName.replace(dirName.find("Tracks"), 6, "");
108  }
109  std::replace(dirName.begin(), dirName.end(), ':', '_');
110  ibooker.setCurrentFolder(dirName);
111 
112  //ibooker.goUp();
113  std::string simName = dirName;
114  simName += "/SimTracks";
115  hSimTracks = new HTrackVariables(ibooker, simName, "SimTracks");
116 
117  ibooker.cd();
118  ibooker.setCurrentFolder(dirName);
119 
120  // Create the root file
121  //theFile = new TFile(theRootFileName.c_str(), "RECREATE");
122 
123  if (doSeedsAnalysis) {
124  ibooker.cd();
125  ibooker.setCurrentFolder(dirName);
126  hRecoSeedInner = new HTrack(ibooker, dirName, "RecoSeed", "Inner");
127  hRecoSeedPCA = new HTrack(ibooker, dirName, "RecoSeed", "PCA");
128  }
129 
130  if (doTracksAnalysis) {
131  ibooker.cd();
132  ibooker.setCurrentFolder(dirName);
133  hRecoTracksPCA = new HTrack(ibooker, dirName, "RecoTracks", "PCA");
134  hRecoTracksInner = new HTrack(ibooker, dirName, "RecoTracks", "Inner");
135  hRecoTracksOuter = new HTrack(ibooker, dirName, "RecoTracks", "Outer");
136 
137  ibooker.cd();
138  ibooker.setCurrentFolder(dirName);
139 
140  // General Histos
141 
142  hChi2 = ibooker.book1D("chi2", "#chi^2", 200, 0, 200);
143  hChi2VsEta = ibooker.book2D("chi2VsEta", "#chi^2 VS #eta", 120, -3., 3., 200, 0, 200);
144 
145  hChi2Norm = ibooker.book1D("chi2Norm", "Normalized #chi^2", 400, 0, 100);
146  hChi2NormVsEta = ibooker.book2D("chi2NormVsEta", "Normalized #chi^2 VS #eta", 120, -3., 3., 400, 0, 100);
147 
148  hHitsPerTrack = ibooker.book1D("HitsPerTrack", "Number of hits per track", 55, 0, 55);
149  hHitsPerTrackVsEta =
150  ibooker.book2D("HitsPerTrackVsEta", "Number of hits per track VS #eta", 120, -3., 3., 55, 0, 55);
151 
152  hDof = ibooker.book1D("dof", "Number of Degree of Freedom", 55, 0, 55);
153  hDofVsEta = ibooker.book2D("dofVsEta", "Number of Degree of Freedom VS #eta", 120, -3., 3., 55, 0, 55);
154 
155  hChi2Prob = ibooker.book1D("chi2Prob", "#chi^2 probability", 200, 0, 1);
156  hChi2ProbVsEta = ibooker.book2D("chi2ProbVsEta", "#chi^2 probability VS #eta", 120, -3., 3., 200, 0, 1);
157 
158  hNumberOfTracks = ibooker.book1D("NumberOfTracks", "Number of reconstructed tracks per event", 200, 0, 200);
159  hNumberOfTracksVsEta = ibooker.book2D(
160  "NumberOfTracksVsEta", "Number of reconstructed tracks per event VS #eta", 120, -3., 3., 10, 0, 10);
161 
162  hChargeVsEta = ibooker.book2D("ChargeVsEta", "Charge vs #eta gen", 120, -3., 3., 4, -2., 2.);
163  hChargeVsPt = ibooker.book2D("ChargeVsPt", "Charge vs P_{T} gen", 250, 0, 200, 4, -2., 2.);
164  hPtRecVsPtGen = ibooker.book2D("PtRecVsPtGen", "P_{T} rec vs P_{T} gen", 250, 0, 200, 250, 0, 200);
165 
166  hDeltaPtVsEta = ibooker.book2D("DeltaPtVsEta", "#Delta P_{t} vs #eta gen", 120, -3., 3., 500, -250., 250.);
167  hDeltaPt_In_Out_VsEta =
168  ibooker.book2D("DeltaPt_In_Out_VsEta_", "P^{in}_{t} - P^{out}_{t} vs #eta gen", 120, -3., 3., 500, -250., 250.);
169  }
170 }
171 
172 void MuonTrackAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
173  LogDebug("MuonTrackAnalyzer") << "Run: " << event.id().run() << " Event: " << event.id().event();
174 
175  // Update the services
176  theService->update(eventSetup);
177 
179  event.getByToken(theSimTracksToken, simTracks);
180  fillPlots(event, simTracks);
181 
182  if (doTracksAnalysis)
183  tracksAnalysis(event, eventSetup, simTracks);
184 
185  if (doSeedsAnalysis)
186  seedsAnalysis(event, eventSetup, simTracks);
187 }
188 
190  const EventSetup &eventSetup,
193 
194  // Get the RecTrack collection from the event
196  event.getByToken(theSeedsToken, seeds);
197 
198  LogTrace("MuonTrackAnalyzer") << "Number of reconstructed seeds: " << seeds->size() << endl;
199 
200  for (TrajectorySeedCollection::const_iterator seed = seeds->begin(); seed != seeds->end(); ++seed) {
201  TrajectoryStateOnSurface seedTSOS = getSeedTSOS(*seed);
202  pair<SimTrack, double> sim = getSimTrack(seedTSOS, simTracks);
203  fillPlots(seedTSOS, sim.first, hRecoSeedInner, debug);
204 
205  std::pair<bool, FreeTrajectoryState> propSeed = theUpdator->propagateToNominalLine(seedTSOS);
206  if (propSeed.first)
207  fillPlots(propSeed.second, sim.first, hRecoSeedPCA, debug);
208  else
209  LogTrace("MuonTrackAnalyzer") << "Error in seed propagation" << endl;
210  }
211 }
212 
214  const EventSetup &eventSetup,
217 
218  // Get the RecTrack collection from the event
220  event.getByToken(theTracksToken, tracks);
221 
222  LogTrace("MuonTrackAnalyzer") << "Reconstructed tracks: " << tracks->size() << endl;
223  hNumberOfTracks->Fill(tracks->size());
224 
225  if (!tracks->empty())
226  numberOfRecTracks++;
227 
228  // Loop over the Rec tracks
229  for (reco::TrackCollection::const_iterator t = tracks->begin(); t != tracks->end(); ++t) {
230  reco::TransientTrack track(*t, &*theService->magneticField(), theService->trackingGeometry());
231 
232  TrajectoryStateOnSurface outerTSOS = track.outermostMeasurementState();
233  TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
234  TrajectoryStateOnSurface pcaTSOS = track.impactPointState();
235 
236  pair<SimTrack, double> sim = getSimTrack(pcaTSOS, simTracks);
237  SimTrack simTrack = sim.first;
238  hNumberOfTracksVsEta->Fill(simTrack.momentum().eta(), tracks->size());
239  fillPlots(track, simTrack);
240 
241  LogTrace("MuonTrackAnalyzer") << "State at the outer surface: " << endl;
242  fillPlots(outerTSOS, simTrack, hRecoTracksOuter, debug);
243 
244  LogTrace("MuonTrackAnalyzer") << "State at the inner surface: " << endl;
245  fillPlots(innerTSOS, simTrack, hRecoTracksInner, debug);
246 
247  LogTrace("MuonTrackAnalyzer") << "State at PCA: " << endl;
248  fillPlots(pcaTSOS, simTrack, hRecoTracksPCA, debug);
249 
250  double deltaPt_in_out = innerTSOS.globalMomentum().perp() - outerTSOS.globalMomentum().perp();
251  hDeltaPt_In_Out_VsEta->Fill(simTrack.momentum().eta(), deltaPt_in_out);
252 
253  double deltaPt_pca_sim = pcaTSOS.globalMomentum().perp() - sqrt(simTrack.momentum().Perp2());
254  hDeltaPtVsEta->Fill(simTrack.momentum().eta(), deltaPt_pca_sim);
255 
256  hChargeVsEta->Fill(simTrack.momentum().eta(), pcaTSOS.charge());
257 
258  hChargeVsPt->Fill(sqrt(simTrack.momentum().perp2()), pcaTSOS.charge());
259 
260  hPtRecVsPtGen->Fill(sqrt(simTrack.momentum().perp2()), pcaTSOS.globalMomentum().perp());
261  }
262  LogTrace("MuonTrackAnalyzer") << "--------------------------------------------" << endl;
263 }
264 
266  if (!checkMuonSimHitPresence(event, simTracks))
267  return;
268 
269  // Loop over the Sim tracks
270  SimTrackContainer::const_iterator simTrack;
271  LogTrace("MuonTrackAnalyzer") << "Simulated tracks: " << simTracks->size() << endl;
272 
273  for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
274  if (abs((*simTrack).type()) == 13) {
275  if (!isInTheAcceptance((*simTrack).momentum().eta()))
276  continue; // FIXME!!
277 
278  numberOfSimTracks++;
279 
280  LogTrace("MuonTrackAnalyzer") << "Simualted muon:" << endl;
281  LogTrace("MuonTrackAnalyzer") << "Sim pT: " << sqrt((*simTrack).momentum().perp2()) << endl;
282  LogTrace("MuonTrackAnalyzer") << "Sim Eta: " << (*simTrack).momentum().eta() << endl; // FIXME
283 
284  hSimTracks->Fill((*simTrack).momentum().mag(),
285  sqrt((*simTrack).momentum().perp2()),
286  (*simTrack).momentum().eta(),
287  (*simTrack).momentum().phi(),
288  -(*simTrack).type() / abs((*simTrack).type())); // Double FIXME
289  LogTrace("MuonTrackAnalyzer") << "hSimTracks filled" << endl;
290  }
291 
292  LogTrace("MuonTrackAnalyzer") << endl;
293 }
294 
296  LogTrace("MuonTrackAnalyzer") << "Analizer: New track, chi2: " << track.chi2() << " dof: " << track.ndof() << endl;
297  hChi2->Fill(track.chi2());
298  hDof->Fill(track.ndof());
299  hChi2Norm->Fill(track.normalizedChi2());
300  hHitsPerTrack->Fill(track.recHitsSize());
301 
302  hChi2Prob->Fill(ChiSquaredProbability(track.chi2(), track.ndof()));
303 
304  hChi2VsEta->Fill(simTrack.momentum().eta(), track.chi2());
305  hChi2NormVsEta->Fill(simTrack.momentum().eta(), track.normalizedChi2());
306  hChi2ProbVsEta->Fill(simTrack.momentum().eta(), ChiSquaredProbability(track.chi2(), track.ndof()));
307  hHitsPerTrackVsEta->Fill(simTrack.momentum().eta(), track.recHitsSize());
308  hDofVsEta->Fill(simTrack.momentum().eta(), track.ndof());
309 }
310 
313  HTrack *histo,
315  LogTrace("MuonTrackAnalyzer") << debug.dumpTSOS(recoTSOS) << endl;
316  histo->Fill(recoTSOS);
317 
318  GlobalVector tsosVect = recoTSOS.globalMomentum();
319  math::XYZVectorD reco(tsosVect.x(), tsosVect.y(), tsosVect.z());
320  double deltaRVal = deltaR<double>(reco.eta(), reco.phi(), simTrack.momentum().eta(), simTrack.momentum().phi());
321  histo->FillDeltaR(deltaRVal);
322 
323  histo->computeResolutionAndPull(recoTSOS, simTrack);
324 }
325 
328  HTrack *histo,
330  LogTrace("MuonTrackAnalyzer") << debug.dumpFTS(recoFTS) << endl;
331  histo->Fill(recoFTS);
332 
333  GlobalVector ftsVect = recoFTS.momentum();
334  math::XYZVectorD reco(ftsVect.x(), ftsVect.y(), ftsVect.z());
335  double deltaRVal = deltaR<double>(reco.eta(), reco.phi(), simTrack.momentum().eta(), simTrack.momentum().phi());
336  histo->FillDeltaR(deltaRVal);
337 
338  histo->computeResolutionAndPull(recoFTS, simTrack);
339 }
340 
343  // // Loop over the Sim tracks
344  // SimTrackContainer::const_iterator simTrack;
345 
346  // SimTrack result;
347  // int mu=0;
348  // for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
349  // if (abs((*simTrack).type()) == 13) {
350  // result = *simTrack;
351  // ++mu;
352  // }
353 
354  // if(mu != 1) LogTrace("MuonTrackAnalyzer") << "WARNING!! more than 1 simulated muon!!" <<endl;
355  // return result;
356 
357  // Loop over the Sim tracks
358  SimTrackContainer::const_iterator simTrack;
359 
361 
362  double bestDeltaR = 10e5;
363  for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
364  if (abs((*simTrack).type()) != 13)
365  continue;
366 
367  // double newDeltaR = tsos.globalMomentum().basicVector().deltaR(simTrack->momentum().vect());
368  GlobalVector tsosVect = tsos.globalMomentum();
369  math::XYZVectorD vect(tsosVect.x(), tsosVect.y(), tsosVect.z());
370  double newDeltaR = deltaR<double>(vect.eta(), vect.phi(), simTrack->momentum().eta(), simTrack->momentum().phi());
371 
372  if (newDeltaR < bestDeltaR) {
373  LogTrace("MuonTrackAnalyzer") << "Matching Track with DeltaR = " << newDeltaR << endl;
374  bestDeltaR = newDeltaR;
375  result = *simTrack;
376  }
377  }
378  return pair<SimTrack, double>(result, bestDeltaR);
379 }
380 
382  switch (theEtaRange) {
383  case all:
384  return (abs(eta) <= 2.4) ? true : false;
385  case barrel:
386  return (abs(eta) < 1.1) ? true : false;
387  case endcap:
388  return (abs(eta) >= 1.1 && abs(eta) <= 2.4) ? true : false;
389  default: {
390  LogTrace("MuonTrackAnalyzer") << "No correct Eta range selected!! " << endl;
391  return false;
392  }
393  }
394 }
395 
397  // Get the SimHit collection from the event
398  Handle<PSimHitContainer> dtSimHits;
399  event.getByToken(theDTSimHitToken, dtSimHits);
400 
401  Handle<PSimHitContainer> cscSimHits;
402  event.getByToken(theCSCSimHitToken, cscSimHits);
403 
404  Handle<PSimHitContainer> rpcSimHits;
405  event.getByToken(theRPCSimHitToken, rpcSimHits);
406 
407  map<unsigned int, vector<const PSimHit *> > mapOfMuonSimHits;
408 
409  for (PSimHitContainer::const_iterator simhit = dtSimHits->begin(); simhit != dtSimHits->end(); ++simhit) {
410  if (abs(simhit->particleType()) != 13)
411  continue;
412  mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
413  }
414 
415  for (PSimHitContainer::const_iterator simhit = cscSimHits->begin(); simhit != cscSimHits->end(); ++simhit) {
416  if (abs(simhit->particleType()) != 13)
417  continue;
418  mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
419  }
420 
421  for (PSimHitContainer::const_iterator simhit = rpcSimHits->begin(); simhit != rpcSimHits->end(); ++simhit) {
422  if (abs(simhit->particleType()) != 13)
423  continue;
424  mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
425  }
426 
427  bool presence = false;
428 
429  for (SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
430  if (abs(simTrack->type()) != 13)
431  continue;
432 
433  map<unsigned int, vector<const PSimHit *> >::const_iterator mapIterator =
434  mapOfMuonSimHits.find(simTrack->trackId());
435 
436  if (mapIterator != mapOfMuonSimHits.end())
437  presence = true;
438  }
439 
440  return presence;
441 }
442 
444  // Get the Trajectory State on Det (persistent version of a TSOS) from the seed
445  PTrajectoryStateOnDet pTSOD = seed.startingState();
446 
447  // Transform it in a TrajectoryStateOnSurface
448 
449  DetId seedDetId(pTSOD.detId());
450 
451  const GeomDet *gdet = theService->trackingGeometry()->idToDet(seedDetId);
452 
453  TrajectoryStateOnSurface initialState =
454  trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*theService->magneticField());
455 
456  // Get the layer on which the seed relies
457  const DetLayer *initialLayer = theService->detLayerGeometry()->idToLayer(seedDetId);
458 
460 
461  // ask for compatible layers
462  vector<const DetLayer *> detLayers;
463  detLayers =
464  theService->muonNavigationSchool()->compatibleLayers(*initialLayer, *initialState.freeState(), detLayerOrder);
465 
466  TrajectoryStateOnSurface result = initialState;
467  if (!detLayers.empty()) {
468  const DetLayer *finalLayer = detLayers.back();
469  const TrajectoryStateOnSurface propagatedState =
470  theService->propagator(theSeedPropagatorName)->propagate(initialState, finalLayer->surface());
471  if (propagatedState.isValid())
472  result = propagatedState;
473  }
474 
475  return result;
476 }
#define LogDebug(id)
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX)
Definition: DQMStore.cc:239
T getParameter(std::string const &) const
void computeResolutionAndPull(TrajectoryStateOnSurface &vtx, SimTrack &simTrack)
Definition: HTrack.cc:117
T perp() const
Definition: PV3DBase.h:69
std::pair< SimTrack, double > getSimTrack(TrajectoryStateOnSurface &tsos, edm::Handle< edm::SimTrackContainer > simTracks)
bool checkMuonSimHitPresence(const edm::Event &event, edm::Handle< edm::SimTrackContainer > simTracks)
void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:418
dqm::legacy::DQMStore * dbe_
T y() const
Definition: PV3DBase.h:60
def replace(string, replacements)
size_t recHitsSize() const
number of RecHits
PropagationDirection
TrajectoryStateOnSurface getSeedTSOS(const TrajectorySeed &seed)
void seedsAnalysis(const edm::Event &event, const edm::EventSetup &eventSetup, edm::Handle< edm::SimTrackContainer > simTracks)
double chi2() const
chi-squared of the fit
MuonTrackAnalyzer(const edm::ParameterSet &pset)
Constructor.
std::string dumpFTS(const FreeTrajectoryState &fts) const
std::string dumpTSOS(const TrajectoryStateOnSurface &tsos) const
void beginJob() override
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
Definition: Vector3D.h:8
T sqrt(T t)
Definition: SSEVec.h:19
bool isInTheAcceptance(double eta)
FreeTrajectoryState const * freeState(bool withErrors=true) const
T z() const
Definition: PV3DBase.h:61
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
float ChiSquaredProbability(double chiSquared, double nrDOF)
unsigned int detId() const
TH1F * hChi2
void Fill(TrajectoryStateOnSurface &)
Definition: HTrack.cc:104
GlobalVector momentum() const
#define LogTrace(id)
Definition: DetId.h:17
PTrajectoryStateOnDet const & startingState() const
#define debug
Definition: HDRShower.cc:19
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field)
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
double normalizedChi2() const
chi-squared divided by n.d.o.f.
std::string const & label() const
Definition: InputTag.h:36
double ndof() const
number of degrees of freedom of the fit
std::string const & process() const
Definition: InputTag.h:40
void tracksAnalysis(const edm::Event &event, const edm::EventSetup &eventSetup, edm::Handle< edm::SimTrackContainer > simTracks)
fixed size matrix
HLT enums.
GlobalVector globalMomentum() const
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Definition: DQMStore.cc:266
const math::XYZTLorentzVectorD & momentum() const
Definition: CoreSimTrack.h:19
~MuonTrackAnalyzer() override
Destructor.
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
void FillDeltaR(double)
Definition: HTrack.cc:111
T x() const
Definition: PV3DBase.h:59
Definition: HTrack.h:12
std::string const & instance() const
Definition: InputTag.h:37
void fillPlots(const edm::Event &event, edm::Handle< edm::SimTrackContainer > &simTracks)
Definition: event.py:1
Definition: Run.h:45