48 theSimTracksToken = consumes<edm::SimTrackContainer>(theSimTracksLabel);
51 theTracksToken = consumes<reco::TrackCollection>(theTracksLabel);
52 doTracksAnalysis =
pset.getUntrackedParameter<
bool>(
"DoTracksAnalysis",
true);
54 doSeedsAnalysis =
pset.getUntrackedParameter<
bool>(
"DoSeedsAnalysis",
false);
55 if (doSeedsAnalysis) {
57 theSeedsToken = consumes<TrajectorySeedCollection>(theSeedsLabel);
59 theSeedPropagatorName = updatorPar.
getParameter<
string>(
"Propagator");
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);
71 theEtaRange = (
EtaRange)
pset.getParameter<
int>(
"EtaRange");
74 numberOfSimTracks = 0;
76 numberOfRecTracks = 0;
79 out =
pset.getUntrackedParameter<
string>(
"rootFileName");
81 subsystemname_ =
pset.getUntrackedParameter<
std::string>(
"subSystemFolder",
"YourSubsystem");
101 dirName += algo.
process() +
"_";
102 if (!algo.
label().empty())
103 dirName += algo.
label() +
"_";
106 if (dirName.find(
"Tracks") < dirName.length()) {
107 dirName.replace(dirName.find(
"Tracks"), 6,
"");
114 simName +=
"/SimTracks";
123 if (doSeedsAnalysis) {
126 hRecoSeedInner =
new HTrack(ibooker, dirName,
"RecoSeed",
"Inner");
127 hRecoSeedPCA =
new HTrack(ibooker, dirName,
"RecoSeed",
"PCA");
130 if (doTracksAnalysis) {
133 hRecoTracksPCA =
new HTrack(ibooker, dirName,
"RecoTracks",
"PCA");
134 hRecoTracksInner =
new HTrack(ibooker, dirName,
"RecoTracks",
"Inner");
135 hRecoTracksOuter =
new HTrack(ibooker, dirName,
"RecoTracks",
"Outer");
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);
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);
148 hHitsPerTrack = ibooker.
book1D(
"HitsPerTrack",
"Number of hits per track", 55, 0, 55);
150 ibooker.
book2D(
"HitsPerTrackVsEta",
"Number of hits per track VS #eta", 120, -3., 3., 55, 0, 55);
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);
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);
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);
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);
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.);
173 LogDebug(
"MuonTrackAnalyzer") <<
"Run: " <<
event.id().run() <<
" Event: " <<
event.id().event();
176 theService->update(eventSetup);
179 event.getByToken(theSimTracksToken, simTracks);
180 fillPlots(event, simTracks);
182 if (doTracksAnalysis)
183 tracksAnalysis(event, eventSetup, simTracks);
186 seedsAnalysis(event, eventSetup, simTracks);
196 event.getByToken(theSeedsToken, seeds);
198 LogTrace(
"MuonTrackAnalyzer") <<
"Number of reconstructed seeds: " << seeds->size() << endl;
200 for (TrajectorySeedCollection::const_iterator
seed = seeds->begin();
seed != seeds->end(); ++
seed) {
202 pair<SimTrack, double>
sim = getSimTrack(seedTSOS, simTracks);
203 fillPlots(seedTSOS, sim.first, hRecoSeedInner, debug);
205 std::pair<bool, FreeTrajectoryState> propSeed = theUpdator->propagateToNominalLine(seedTSOS);
207 fillPlots(propSeed.second, sim.first, hRecoSeedPCA, debug);
209 LogTrace(
"MuonTrackAnalyzer") <<
"Error in seed propagation" << endl;
220 event.getByToken(theTracksToken, tracks);
222 LogTrace(
"MuonTrackAnalyzer") <<
"Reconstructed tracks: " << tracks->size() << endl;
223 hNumberOfTracks->Fill(tracks->size());
225 if (!tracks->empty())
229 for (reco::TrackCollection::const_iterator
t = tracks->begin();
t != tracks->end(); ++
t) {
236 pair<SimTrack, double>
sim = getSimTrack(pcaTSOS, simTracks);
238 hNumberOfTracksVsEta->Fill(simTrack.
momentum().eta(), tracks->size());
239 fillPlots(
track, simTrack);
241 LogTrace(
"MuonTrackAnalyzer") <<
"State at the outer surface: " << endl;
242 fillPlots(outerTSOS, simTrack, hRecoTracksOuter, debug);
244 LogTrace(
"MuonTrackAnalyzer") <<
"State at the inner surface: " << endl;
245 fillPlots(innerTSOS, simTrack, hRecoTracksInner, debug);
247 LogTrace(
"MuonTrackAnalyzer") <<
"State at PCA: " << endl;
248 fillPlots(pcaTSOS, simTrack, hRecoTracksPCA, debug);
250 double deltaPt_in_out = innerTSOS.
globalMomentum().
perp() - outerTSOS.globalMomentum().perp();
251 hDeltaPt_In_Out_VsEta->Fill(simTrack.
momentum().eta(), deltaPt_in_out);
254 hDeltaPtVsEta->Fill(simTrack.
momentum().eta(), deltaPt_pca_sim);
262 LogTrace(
"MuonTrackAnalyzer") <<
"--------------------------------------------" << endl;
266 if (!checkMuonSimHitPresence(event, simTracks))
270 SimTrackContainer::const_iterator
simTrack;
271 LogTrace(
"MuonTrackAnalyzer") <<
"Simulated tracks: " << simTracks->size() << endl;
273 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++
simTrack)
274 if (
abs((*simTrack).type()) == 13) {
275 if (!isInTheAcceptance((*simTrack).momentum().eta()))
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;
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()));
289 LogTrace(
"MuonTrackAnalyzer") <<
"hSimTracks filled" << endl;
292 LogTrace(
"MuonTrackAnalyzer") << endl;
296 LogTrace(
"MuonTrackAnalyzer") <<
"Analizer: New track, chi2: " << track.
chi2() <<
" dof: " << track.
ndof() << endl;
298 hDof->Fill(track.
ndof());
304 hChi2VsEta->Fill(simTrack.
momentum().eta(), track.
chi2());
308 hDofVsEta->Fill(simTrack.
momentum().eta(), track.
ndof());
316 histo->
Fill(recoTSOS);
331 histo->
Fill(recoFTS);
358 SimTrackContainer::const_iterator
simTrack;
362 double bestDeltaR = 10e5;
363 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++
simTrack) {
364 if (
abs((*simTrack).type()) != 13)
370 double newDeltaR = deltaR<double>(vect.eta(), vect.phi(), simTrack->momentum().eta(), simTrack->momentum().phi());
372 if (newDeltaR < bestDeltaR) {
373 LogTrace(
"MuonTrackAnalyzer") <<
"Matching Track with DeltaR = " << newDeltaR << endl;
374 bestDeltaR = newDeltaR;
378 return pair<SimTrack, double>(
result, bestDeltaR);
382 switch (theEtaRange) {
384 return (
abs(eta) <= 2.4) ?
true :
false;
386 return (
abs(eta) < 1.1) ?
true :
false;
388 return (
abs(eta) >= 1.1 &&
abs(eta) <= 2.4) ?
true :
false;
390 LogTrace(
"MuonTrackAnalyzer") <<
"No correct Eta range selected!! " << endl;
399 event.getByToken(theDTSimHitToken, dtSimHits);
402 event.getByToken(theCSCSimHitToken, cscSimHits);
405 event.getByToken(theRPCSimHitToken, rpcSimHits);
407 map<unsigned int, vector<const PSimHit *> > mapOfMuonSimHits;
409 for (PSimHitContainer::const_iterator simhit = dtSimHits->begin(); simhit != dtSimHits->end(); ++simhit) {
410 if (
abs(simhit->particleType()) != 13)
412 mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
415 for (PSimHitContainer::const_iterator simhit = cscSimHits->begin(); simhit != cscSimHits->end(); ++simhit) {
416 if (
abs(simhit->particleType()) != 13)
418 mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
421 for (PSimHitContainer::const_iterator simhit = rpcSimHits->begin(); simhit != rpcSimHits->end(); ++simhit) {
422 if (
abs(simhit->particleType()) != 13)
424 mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
427 bool presence =
false;
433 map<unsigned int, vector<const PSimHit *> >::const_iterator mapIterator =
434 mapOfMuonSimHits.find(
simTrack->trackId());
436 if (mapIterator != mapOfMuonSimHits.end())
451 const GeomDet *gdet = theService->trackingGeometry()->idToDet(seedDetId);
457 const DetLayer *initialLayer = theService->detLayerGeometry()->idToLayer(seedDetId);
462 vector<const DetLayer *> detLayers;
464 theService->muonNavigationSchool()->compatibleLayers(*initialLayer, *initialState.
freeState(), detLayerOrder);
467 if (!detLayers.empty()) {
468 const DetLayer *finalLayer = detLayers.back();
470 theService->propagator(theSeedPropagatorName)->propagate(initialState, finalLayer->
surface());
472 result = propagatedState;
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX)
T getParameter(std::string const &) const
void computeResolutionAndPull(TrajectoryStateOnSurface &vtx, SimTrack &simTrack)
TrackCharge charge() const
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)
dqm::legacy::DQMStore * dbe_
def replace(string, replacements)
size_t recHitsSize() const
number of RecHits
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
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
bool isInTheAcceptance(double eta)
FreeTrajectoryState const * freeState(bool withErrors=true) const
Abs< T >::type abs(const T &t)
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
float ChiSquaredProbability(double chiSquared, double nrDOF)
unsigned int detId() const
void Fill(TrajectoryStateOnSurface &)
GlobalVector momentum() const
PTrajectoryStateOnDet const & startingState() const
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
double normalizedChi2() const
chi-squared divided by n.d.o.f.
double ndof() const
number of degrees of freedom of the fit
void tracksAnalysis(const edm::Event &event, const edm::EventSetup &eventSetup, edm::Handle< edm::SimTrackContainer > simTracks)
GlobalVector globalMomentum() const
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
const math::XYZTLorentzVectorD & momentum() const
~MuonTrackAnalyzer() override
Destructor.
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
void fillPlots(const edm::Event &event, edm::Handle< edm::SimTrackContainer > &simTracks)