File indexing completed on 2022-05-10 22:24:54
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "RecoLocalMuon/DTSegment/test/STAnalyzer.h"
0015
0016
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023
0024 using namespace edm;
0025
0026 #include "TFile.h"
0027 #include "TH1F.h"
0028 #include "TH2F.h"
0029
0030 #include "DataFormats/LTCDigi/interface/LTCDigi.h"
0031
0032 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0033 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0034 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0035 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0036 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0037
0038 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0039 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0040 #include "MagneticField/Engine/interface/MagneticField.h"
0041 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0042 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0043 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0044
0045 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0046 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0047 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0048 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0049
0050 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052
0053
0054 #include <iostream>
0055 #include <cmath>
0056 using namespace std;
0057
0058
0059
0060
0061 STAnalyzer::STAnalyzer(const ParameterSet& pset) : _ev(0) {
0062 if (debug)
0063 cout << "STAnalyzer::STAnalyzer" << endl;
0064
0065
0066 debug = pset.getUntrackedParameter<bool>("debug");
0067 theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0068
0069
0070 theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");
0071
0072
0073 theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");
0074
0075
0076 theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
0077
0078
0079 theSTAMuonLabel = pset.getParameter<string>("SALabel");
0080
0081 thePropagatorName = pset.getParameter<std::string>("PropagatorName");
0082 thePropagator = 0;
0083 thePropagatorToken = esConsumes(edm::ESInputTag("", thePropagatorName));
0084
0085 doSA = pset.getParameter<bool>("doSA");
0086
0087
0088 theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0089 bool dirStat = TH1::AddDirectoryStatus();
0090 TH1::AddDirectory(kTRUE);
0091
0092
0093 new TH1F("hNSA", "Num SA tracks in events", 6, -0.5, 5.5);
0094 new TH1F("hNSADifferentSectors", "Num SA tracks with hits in different sectors", 6, -0.5, 5.5);
0095 new TH1F("hNhitsSA", "Num hits in SA tracks", 50, .5, 50.5);
0096 new TH1F("hChi2SA", "#chi^{2}/NDoF for SA tracks", 100, 0, 100.);
0097 new TH1F("hPIPSA", "p for SA tracks @ IP", 100, 0., 100);
0098 new TH1F("hPtIPSA", "pt for SA tracks @ IP", 100, 0., 100);
0099 new TH1F("hPhiIPSA", "#phi for SA tracks @ IP", 100, -M_PI_2, M_PI_2);
0100 new TH1F("hEtaIPSA", "#eta for SA tracks @ IP", 100, -2.5, 2.5);
0101 new TH1F("hPInnerTSOSSA", "p for SA tracks @ innermost TSOS", 100, 0., 100);
0102 new TH1F("hPtInnerTSOSSA", "pt for SA tracks @ innermost TSOS", 100, 0., 100);
0103 new TH1F("hPhiInnerTSOSSA", "#phi for SA tracks @ innermost TSOS", 100, -M_PI_2, M_PI_2);
0104 new TH1F("hEtaInnerTSOSSA", "#eta for SA tracks @ innermost TSOS", 100, -2.5, 2.5);
0105 new TH1F("hInnerRSA", "Radius of innermost TSOS for SA tracks", 100, 400, 1000.);
0106 new TH1F("hOuterRSA", "Radius of outermost TSOS for SA tracks", 100, 400, 1000.);
0107 new TH2F("hInnerOuterRSA", "Radius of outermost TSOS vs innermost one for SA tracks", 40, 300, 1000., 40, 300, 1000.);
0108
0109 new TH2F("hNSAVsNHits", "Num 1d Hits vs N SA", 100, -0.5, 39.5, 5, -0.5, 4.5);
0110 new TH2F("hNSAVsNSegs2D", "Num 2d Segs vs N SA", 20, -0.5, 19.5, 5, -0.5, 4.5);
0111 new TH2F("hNSAVsNSegs4D", "Num 4d Segs vs N SA", 10, -0.5, 9.5, 5, -0.5, 4.5);
0112
0113 new TH2F("hNHitsSAVsNHits", "Num 1d Hits vs N Hits SA", 40, -0.5, 39.5, 100, 0, 100);
0114 new TH2F("hNHitsSAVsNSegs2D", "Num 2d Segs vs N Hits SA", 20, -0.5, 19.5, 100, 0, 100);
0115 new TH2F("hNHitsSAVsNSegs4D", "Num 4d Segs vs N Hits SA", 10, -0.5, 9.5, 100, 0, 100);
0116
0117 new TH1F("hSANLayers", "Num Layers used SA", 50, -0.5, 49.5);
0118 new TH1F("hSANSLs", "Num SuperLayers used SA", 15, -0.5, 14.5);
0119 new TH1F("hSANLayersPerSL", "Num Layers per SL used SA", 15, -0.5, 14.5);
0120 new TH1F("hSANChambers", "Num Chambers used SA", 6, -0.5, 5.5);
0121 new TH1F("hSANLayersPerChamber", "Num Layers per SL used SA", 15, -0.5, 14.5);
0122 new TH1F("hSANSLsPerChamber", "Num SLs per Chamber used SA", 5, -0.5, 4.5);
0123
0124 new TH2F("hSANLayersVsNChambers", "Num Layers Vs Num Chambers used SA", 6, -0.5, 5.5, 50, -0.5, 49.5);
0125 new TH2F("hSANSLsVsNChambers", "Num SuperLayers Vs Num Chambers used SA", 6, -0.5, 5.5, 16, -0.5, 15.5);
0126
0127 new TH2F("hSANHitsVsNLayers", "Num Hits vs Num Layers used SA", 50, -0.5, 49.5, 50, 0., 50.);
0128 new TH2F("hSANHitsVsNSLs", "Num Hits vs Num SuperLayers used SA", 16, -0.5, 15.5, 50, 0., 50.);
0129 new TH2F("hSANHitsVsNChambers", "Num Hits vs Num Chambers used SA", 6, -0.5, 5.5, 50, 0., 50.);
0130
0131 new TH2F("hHitsPosXYSA", "Hits position (x,y) SA", 100, -800, 800, 100, -800, 800);
0132 new TH2F("hHitsPosXZSA", "Hits position (x,z) SA", 100, -800, 800, 100, -670, 670);
0133 new TH2F("hHitsPosYZSA", "Hits position (y,z) SA", 100, -670, 670, 100, -800, 800);
0134
0135 new TH1F("nHitsSAInChamber", "Num STA recHits in each chamber", 4, 0.5, 4.5);
0136 new TH1F("nHitsSAInSL", "Num STA recHits in each sl", 50, 0.5, 50.5);
0137 new TH1F("nHitsSAInLayer", "Num STA recHits in each layer", 75, 25.5, 100.5);
0138
0139 new TH1F("hHitsLostChamber", "Num hits in det w/o associated hits chamber", 30, 0., 30.);
0140 new TH1F("hHitsLostSL", "Num hits in det w/o associated hits SL", 30, 0., 30.);
0141 new TH1F("hHitsLostLayer", "Num hits in det w/o associated hits Layer", 30, 0., 30.);
0142
0143 new TH1F("hMinDistChamberFound",
0144 "Min distance between extrap post and closest hit in det w associated hits chamber",
0145 100,
0146 0.,
0147 300.);
0148 new TH1F("hMinDistChamberNotFound",
0149 "Min distance between extrap post and closest hit in det w/o associated hits chamber",
0150 100,
0151 0.,
0152 300.);
0153 new TH1F("hMinPhiChamberFound",
0154 "Min phi between extrap post and closest hit in det w associated hits chamber",
0155 100,
0156 -1.,
0157 1.);
0158 new TH1F("hMinPhiChamberNotFound",
0159 "Min phi between extrap post and closest hit in det w/o associated hits chamber",
0160 100,
0161 -1.,
0162 1.);
0163 new TH1F("hMinThetaChamberFound",
0164 "Min theta between extrap post and closest hit in det w associated hits chamber",
0165 100,
0166 -1.,
0167 1.);
0168 new TH1F("hMinThetaChamberNotFound",
0169 "Min theta between extrap post and closest hit in det w/o associated hits chamber",
0170 100,
0171 -1.,
0172 1.);
0173
0174 new TH1F(
0175 "hMinDistSLFound", "Min distance between extrap post and closest hit in det w associated hits SL", 100, 0., 300.);
0176 new TH1F("hMinDistSLNotFound",
0177 "Min distance between extrap post and closest hit in det w/o associated hits SL",
0178 100,
0179 0.,
0180 300.);
0181 new TH1F(
0182 "hMinThetaSLFound", "Min theta between extrap post and closest hit in det w associated hits SL", 100, -1., 1.);
0183 new TH1F("hMinThetaSLNotFound",
0184 "Min theta between extrap post and closest hit in det w/o associated hits SL",
0185 100,
0186 -1.,
0187 1.);
0188
0189 new TH1F("hMinDistLayerFound",
0190 "Min distance between extrap post and closest hit in det w associated hits Layer",
0191 100,
0192 0.,
0193 300.);
0194 new TH1F("hMinDistLayerNotFound",
0195 "Min distance between extrap post and closest hit in det w/o associated hits Layer",
0196 100,
0197 0.,
0198 300.);
0199
0200 TH1::AddDirectory(dirStat);
0201
0202 theDtGeomTokenBR = esConsumes<edm::Transition::BeginRun>();
0203 theDtGeomToken = esConsumes();
0204 theMGFieldToken = esConsumes();
0205 theTrackingGeometryToken = esConsumes();
0206 }
0207
0208
0209 STAnalyzer::~STAnalyzer() {
0210 theFile->cd();
0211 theFile->Write();
0212 theFile->Close();
0213 }
0214
0215 void STAnalyzer::beginRun(const edm::Run& run, const EventSetup& setup) {
0216
0217 ESHandle<DTGeometry> dtGeom = setup.getHandle(theDtGeomTokenBR);
0218
0219 if (firstPass) {
0220 const std::vector<const DTChamber*>& chs = dtGeom->chambers();
0221 for (auto ch = chs.begin(); ch != chs.end(); ++ch)
0222 hitsPerChamber[(*ch)->id()] = 0;
0223 }
0224
0225 firstPass = false;
0226 }
0227
0228
0229 void STAnalyzer::beginJob() {
0230 if (debug)
0231 cout << "STAnalyzer::beginJob" << endl;
0232 }
0233
0234 void STAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0235 if (debug)
0236 cout << "STAnalyzer::analyze" << endl;
0237
0238 if (debug)
0239 cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev++ << endl;
0240
0241 if (doSA)
0242 analyzeSATrack(event, eventSetup);
0243 }
0244
0245 void STAnalyzer::analyzeSATrack(const Event& event, const EventSetup& eventSetup) {
0246 if (debug)
0247 cout << "STAnalyzer::analyzeSATrack" << endl;
0248 if (!thePropagator) {
0249 ESHandle<Propagator> prop = eventSetup.getHandle(thePropagatorToken);
0250 thePropagator = prop->clone();
0251 thePropagator->setPropagationDirection(anyDirection);
0252 }
0253
0254 MuonPatternRecoDumper muonDumper;
0255
0256
0257 edm::Handle<DTRecSegment4DCollection> segs;
0258 event.getByLabel(theRecHits4DLabel, segs);
0259
0260
0261 edm::Handle<DTRecSegment2DCollection> segs2d;
0262 event.getByLabel(theRecHits2DLabel, segs2d);
0263
0264
0265 Handle<DTRecHitCollection> dtRecHits;
0266 event.getByLabel(theRecHits1DLabel, dtRecHits);
0267
0268
0269 Handle<reco::TrackCollection> staTracks;
0270 event.getByLabel(theSTAMuonLabel, staTracks);
0271 if (debug)
0272 cout << "SA " << staTracks->size() << endl;
0273
0274 std::vector<Handle<reco::TrackCollection> > tracks;
0275 event.getManyByType(tracks);
0276 std::vector<Handle<reco::TrackCollection> >::iterator i(tracks.begin()), end(tracks.end());
0277 for (; i != end; ++i) {
0278 if (debug)
0279 cout << "ALL Tracks id " << (*i).id() << endl;
0280 if (debug)
0281 cout << "ALL Tracks size " << (*i).product()->size() << endl;
0282 }
0283
0284 ESHandle<MagneticField> theMGField = eventSetup.getHandle(theMGFieldToken);
0285
0286 ESHandle<GlobalTrackingGeometry> theTrackingGeometry = eventSetup.getHandle(theTrackingGeometryToken);
0287
0288 reco::TrackCollection::const_iterator staTrack;
0289
0290 double recPt = 0.;
0291 histo("hNSA")->Fill(staTracks->size());
0292 histo2d("hNSAVsNHits")->Fill(dtRecHits->size(), staTracks->size());
0293 histo2d("hNSAVsNSegs2D")->Fill(segs2d->size(), staTracks->size());
0294 histo2d("hNSAVsNSegs4D")->Fill(segs->size(), staTracks->size());
0295
0296 if (debug && staTracks->size())
0297 cout << endl << "R:E " << event.id().run() << ":" << event.id().event() << " SA " << staTracks->size() << endl;
0298
0299 for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack) {
0300 reco::TransientTrack track(*staTrack, &*theMGField, theTrackingGeometry);
0301
0302 if (debug) {
0303 cout << muonDumper.dumpFTS(track.impactPointTSCP().theState());
0304
0305 recPt = track.impactPointTSCP().momentum().perp();
0306 cout << " p: " << track.impactPointTSCP().momentum().mag() << " pT: " << recPt << endl;
0307 cout << " normalized chi2: " << track.normalizedChi2() << endl;
0308 }
0309
0310 histo("hChi2SA")->Fill(track.normalizedChi2());
0311
0312 histo("hPIPSA")->Fill(track.impactPointTSCP().momentum().mag());
0313 histo("hPtIPSA")->Fill(recPt);
0314 histo("hPhiIPSA")->Fill(track.impactPointTSCP().momentum().phi());
0315 histo("hEtaIPSA")->Fill(track.impactPointTSCP().momentum().eta());
0316
0317 TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0318 if (debug) {
0319 cout << "Inner TSOS:" << endl;
0320 cout << muonDumper.dumpTSOS(innerTSOS);
0321 cout << " p: " << innerTSOS.globalMomentum().mag() << " pT: " << innerTSOS.globalMomentum().perp() << endl;
0322 }
0323
0324
0325
0326
0327 float radiusECAL = 129.0;
0328 Cylinder::PositionType pos0;
0329 Cylinder::RotationType rot0;
0330 const Cylinder::CylinderPointer ecal = Cylinder::build(radiusECAL, pos0, rot0);
0331
0332
0333 TrajectoryStateOnSurface tsosAtEcal = thePropagator->propagate(*innerTSOS.freeState(), *ecal);
0334
0335
0336
0337
0338
0339
0340
0341 float radiusHCAL = 181.1;
0342 const Cylinder::CylinderPointer hcal = Cylinder::build(radiusHCAL, pos0, rot0);
0343
0344
0345 TrajectoryStateOnSurface tsosAtHcal = thePropagator->propagate(*innerTSOS.freeState(), *hcal);
0346
0347
0348
0349
0350
0351
0352 histo("hPInnerTSOSSA")->Fill(innerTSOS.globalMomentum().mag());
0353 histo("hPtInnerTSOSSA")->Fill(innerTSOS.globalMomentum().perp());
0354 histo("hPhiInnerTSOSSA")->Fill(innerTSOS.globalMomentum().phi());
0355 histo("hEtaInnerTSOSSA")->Fill(innerTSOS.globalMomentum().eta());
0356
0357 histo("hNhitsSA")->Fill(staTrack->recHitsSize());
0358
0359 histo2d("hNHitsSAVsNHits")->Fill(dtRecHits->size(), staTrack->recHitsSize());
0360 histo2d("hNHitsSAVsNSegs2D")->Fill(segs2d->size(), staTrack->recHitsSize());
0361 histo2d("hNHitsSAVsNSegs4D")->Fill(segs->size(), staTrack->recHitsSize());
0362
0363 histo("hInnerRSA")->Fill(sqrt(staTrack->innerPosition().perp2()));
0364 histo("hOuterRSA")->Fill(sqrt(staTrack->outerPosition().perp2()));
0365 histo2d("hInnerOuterRSA")->Fill(sqrt(staTrack->innerPosition().perp2()), sqrt(staTrack->outerPosition().perp2()));
0366
0367 if (debug)
0368 cout << "RecHits:" << endl;
0369
0370 trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
0371 trackingRecHit_iterator rhend = staTrack->recHitsEnd();
0372
0373
0374 for (std::map<DTChamberId, int>::iterator h = hitsPerChamber.begin(); h != hitsPerChamber.end(); ++h)
0375 (*h).second = 0;
0376 for (std::map<DTSuperLayerId, int>::iterator h = hitsPerSL.begin(); h != hitsPerSL.end(); ++h)
0377 (*h).second = 0;
0378 for (std::map<DTLayerId, int>::iterator h = hitsPerLayer.begin(); h != hitsPerLayer.end(); ++h)
0379 (*h).second = 0;
0380
0381 int firstHitWheel = 0;
0382 int lastHitWheel = 0;
0383 int firstHitSector = 0;
0384 int lastHitSector = 0;
0385 for (trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit) {
0386 const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
0387 GlobalPoint gpos = geomDet->toGlobal((*recHit)->localPosition());
0388 if (debug)
0389 cout << "r: " << gpos.perp() << " z: " << gpos.z() << endl;
0390 histo2d("hHitsPosXYSA")->Fill(gpos.x(), gpos.y());
0391 histo2d("hHitsPosXZSA")->Fill(gpos.z(), gpos.x());
0392 histo2d("hHitsPosYZSA")->Fill(gpos.z(), gpos.y());
0393
0394 if (const DTLayer* lay = dynamic_cast<const DTLayer*>(geomDet)) {
0395 if (firstHitSector == 0)
0396 firstHitSector = lay->id().sector();
0397 lastHitSector = lay->id().sector();
0398 if (firstHitWheel == 0)
0399 firstHitWheel = lay->id().wheel();
0400 lastHitWheel = lay->id().wheel();
0401
0402 if (debug)
0403 cout << "Layer " << lay->id() << endl;
0404 histo("nHitsSAInChamber")->Fill(lay->id().station());
0405 hitsPerChamber[lay->id()]++;
0406 histo("nHitsSAInSL")->Fill(lay->id().station() * 10 + lay->id().superlayer());
0407 hitsPerSL[lay->id()]++;
0408 histo("nHitsSAInLayer")->Fill(lay->id().station() * 20 + lay->id().superlayer() * 5 + lay->id().layer());
0409 hitsPerLayer[lay->id()]++;
0410 }
0411 if (const DTSuperLayer* lay = dynamic_cast<const DTSuperLayer*>(geomDet))
0412 cout << "SuperLayer " << lay->id() << endl;
0413 if (const DTChamber* lay = dynamic_cast<const DTChamber*>(geomDet))
0414 cout << "Chamber " << lay->id() << endl;
0415
0416 TrajectoryStateOnSurface extraptsos = thePropagator->propagate(innerTSOS, geomDet->specificSurface());
0417 }
0418 if (debug) {
0419 cout << "PerChamber " << muonDumper.dumpTSOS(innerTSOS) << endl;
0420 for (std::map<DTChamberId, int>::iterator h = hitsPerChamber.begin(); h != hitsPerChamber.end(); ++h)
0421 if ((*h).second)
0422 cout << (*h).first << ":" << (*h).second << endl;
0423 cout << "=====" << endl;
0424
0425 cout << "PerSL " << endl;
0426 for (std::map<DTSuperLayerId, int>::iterator h = hitsPerSL.begin(); h != hitsPerSL.end(); ++h)
0427 if ((*h).second)
0428 cout << (*h).first << ":" << (*h).second << endl;
0429 cout << "=====" << endl;
0430 cout << "PerLayer " << endl;
0431 for (std::map<DTLayerId, int>::iterator h = hitsPerLayer.begin(); h != hitsPerLayer.end(); ++h)
0432 if ((*h).second)
0433 cout << (*h).first << ":" << (*h).second << endl;
0434 cout << "=====" << endl;
0435 }
0436
0437
0438 ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDtGeomToken);
0439
0440
0441 Handle<DTRecHitCollection> hits1d;
0442 event.getByLabel(theRecHits1DLabel, hits1d);
0443
0444
0445 edm::Handle<DTRecSegment2DCollection> segs2d;
0446 event.getByLabel(theRecHits2DLabel, segs2d);
0447
0448
0449 edm::Handle<DTRecSegment4DCollection> segs;
0450 event.getByLabel(theRecHits4DLabel, segs);
0451
0452
0453
0454 if (firstHitWheel == lastHitWheel &&
0455 ((firstHitSector == lastHitSector) || (firstHitSector == 10 && lastHitSector == 14) ||
0456 (firstHitSector == 14 && lastHitSector == 10) || (firstHitSector == 13 && lastHitSector == 4) ||
0457 (firstHitSector == 4 && lastHitSector == 13))) {
0458 int nLay = 0, nSL = 0, nCh = 0;
0459 for (int c = 1; c <= 5; ++c) {
0460
0461
0462 if (c == 5 && firstHitSector == 14)
0463 firstHitSector = 10;
0464 else if (c == 5 && lastHitSector == 14)
0465 lastHitSector = 10;
0466 else if (c == 5 && firstHitSector == 13)
0467 firstHitSector = 4;
0468 else if (c == 5 && lastHitSector == 13)
0469 lastHitSector = 4;
0470 else if (c == 5)
0471 continue;
0472 int cc = c;
0473 if (c == 5)
0474 cc = 4;
0475
0476 int sect = min(firstHitSector, lastHitSector);
0477 DTChamberId chid(firstHitWheel, cc, sect);
0478 const DTChamber* ch = dtGeom->chamber(chid);
0479 if (hitsPerChamber[chid] == 0) {
0480 missingHit(dtGeom, segs, ch, innerTSOS, false);
0481 } else {
0482 missingHit(dtGeom, segs, ch, innerTSOS, true);
0483 nCh++;
0484
0485 }
0486 int nLayPerCh = 0, nSLPerCh = 0;
0487 for (int sl = 1; sl <= 3; ++sl) {
0488 if (sl == 2 && cc == 4)
0489 continue;
0490 DTSuperLayerId slid(chid, sl);
0491 const DTSuperLayer* dtsl = dtGeom->superLayer(slid);
0492 if (hitsPerSL[slid] == 0) {
0493
0494 missingHit(dtGeom, segs2d, dtsl, innerTSOS, false);
0495 } else {
0496 missingHit(dtGeom, segs2d, dtsl, innerTSOS, true);
0497 nSL++;
0498 nSLPerCh++;
0499
0500 }
0501 int nLayPerSL = 0;
0502 for (int l = 1; l <= 4; ++l) {
0503 DTLayerId lid(slid, l);
0504 const DTLayer* lay = dtGeom->layer(lid);
0505 if (hitsPerLayer[lid] == 0) {
0506
0507 missingHit(dtGeom, hits1d, lay, innerTSOS, false);
0508 } else {
0509 missingHit(dtGeom, hits1d, lay, innerTSOS, true);
0510 nLay++;
0511 nLayPerCh++;
0512 nLayPerSL++;
0513
0514 }
0515 }
0516 histo("hSANLayersPerSL")->Fill(nLayPerSL);
0517 }
0518 histo("hSANLayersPerChamber")->Fill(nLayPerCh);
0519 histo("hSANSLsPerChamber")->Fill(nSLPerCh);
0520 }
0521
0522
0523 histo("hSANLayers")->Fill(nLay);
0524 histo("hSANSLs")->Fill(nSL);
0525 histo("hSANChambers")->Fill(nCh);
0526 histo2d("hSANLayersVsNChambers")->Fill(nCh, nLay);
0527 histo2d("hSANSLsVsNChambers")->Fill(nCh, nSL);
0528
0529 histo2d("hSANHitsVsNLayers")->Fill(nLay, staTrack->recHitsSize());
0530 histo2d("hSANHitsVsNSLs")->Fill(nSL, staTrack->recHitsSize());
0531 histo2d("hSANHitsVsNChambers")->Fill(nCh, staTrack->recHitsSize());
0532 if (nCh > 4)
0533 cout << endl << "R:E " << event.id().run() << ":" << event.id().event() << " nCh " << nCh << endl;
0534 } else {
0535
0536
0537 histo("hNSADifferentSectors")->Fill(1);
0538 }
0539 }
0540 if (debug)
0541 cout << "---" << endl;
0542 }
0543
0544 template <typename T, typename C>
0545 void STAnalyzer::missingHit(const ESHandle<DTGeometry>& dtGeom,
0546 const Handle<C>& segs,
0547 const T* ch,
0548 const TrajectoryStateOnSurface& startTsos,
0549 bool found) {
0550 if (!ch)
0551 return;
0552 TrajectoryStateOnSurface extraptsos = thePropagator->propagate(startTsos, ch->specificSurface());
0553
0554 if (extraptsos.isValid()) {
0555
0556 bool inside = ch->specificSurface().bounds().inside(extraptsos.localPosition());
0557
0558 if (inside) {
0559
0560
0561 typename C::range segsch = segs->get(ch->id());
0562
0563
0564 fillMinDist(segsch, ch, extraptsos, found);
0565 }
0566 }
0567 }
0568
0569 void STAnalyzer::fillMinDist(const DTRecSegment4DCollection::range segsch,
0570 const DTChamber* ch,
0571 const TrajectoryStateOnSurface& extraptsos,
0572 bool found) {
0573 int nSegsCh = segsch.second - segsch.first;
0574 histo("hHitsLostChamber")->Fill(nSegsCh);
0575 if (nSegsCh) {
0576 LocalPoint extrapPos = extraptsos.localPosition();
0577 LocalVector extrapDir = extraptsos.localDirection();
0578
0579 double minDist = 99999.;
0580 double minPhi = 0., minTheta = 0.;
0581 for (DTRecSegment4DCollection::const_iterator hit = segsch.first; hit != segsch.second; ++hit) {
0582
0583 LocalVector dist = hit->localPosition() - extrapPos;
0584
0585 if (dist.perp() < minDist) {
0586 minDist = dist.perp();
0587 minPhi = hit->localDirection().phi() - extrapDir.phi();
0588 minTheta = hit->localDirection().theta() - extrapDir.theta();
0589 }
0590 }
0591 string inOut = (found) ? "Found" : "NotFound";
0592 histo(string("hMinDistChamber") + inOut)->Fill(minDist);
0593 histo(string("hMinPhiChamber") + inOut)->Fill(minPhi);
0594 histo(string("hMinThetaChamber") + inOut)->Fill(minTheta);
0595 }
0596 }
0597
0598 void STAnalyzer::fillMinDist(const DTRecSegment2DCollection::range segsch,
0599 const DTSuperLayer* ch,
0600 const TrajectoryStateOnSurface& extraptsos,
0601 bool found) {
0602 int nSegsCh = segsch.second - segsch.first;
0603 histo("hHitsLostSL")->Fill(nSegsCh);
0604 if (nSegsCh) {
0605 LocalPoint extrapPos = extraptsos.localPosition();
0606 LocalVector extrapDir = extraptsos.localDirection();
0607
0608 double minDist = 99999.;
0609 double minTheta = 0.;
0610 for (DTRecSegment2DCollection::const_iterator hit = segsch.first; hit != segsch.second; ++hit) {
0611
0612 LocalVector dist = hit->localPosition() - extrapPos;
0613
0614 if (dist.perp() < minDist) {
0615 minDist = dist.x();
0616 minTheta = hit->localDirection().theta() - extrapDir.theta();
0617 }
0618 }
0619 string inOut = (found) ? "Found" : "NotFound";
0620 histo(string("hMinDistSL") + inOut)->Fill(minDist);
0621 histo(string("hMinThetaSL") + inOut)->Fill(minTheta);
0622 }
0623 }
0624
0625 void STAnalyzer::fillMinDist(const DTRecHitCollection::range segsch,
0626 const DTLayer* ch,
0627 const TrajectoryStateOnSurface& extraptsos,
0628 bool found) {
0629 int nSegsCh = segsch.second - segsch.first;
0630 histo("hHitsLostLayer")->Fill(nSegsCh);
0631 if (nSegsCh) {
0632 LocalPoint extrapPos = extraptsos.localPosition();
0633
0634 double minDist = 99999.;
0635 for (DTRecHitCollection::const_iterator hit = segsch.first; hit != segsch.second; ++hit) {
0636
0637 LocalVector dist = hit->localPosition() - extrapPos;
0638
0639 if (dist.perp() < minDist)
0640 minDist = dist.x();
0641 }
0642 string inOut = (found) ? "Found" : "NotFound";
0643 histo(string("hMinDistLayer") + inOut)->Fill(minDist);
0644 }
0645 }
0646
0647 TH1F* STAnalyzer::histo(const string& name) const {
0648 if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0649 return h;
0650 else
0651 throw cms::Exception("STAnalyzer") << " Not a TH1F " << name;
0652 }
0653
0654 TH2F* STAnalyzer::histo2d(const string& name) const {
0655 if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0656 return h;
0657 else
0658 throw cms::Exception("STAnalyzer") << " Not a TH2F " << name;
0659 }
0660
0661 string STAnalyzer::toString(const DTLayerId& id) const {
0662 stringstream result;
0663 result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer() << "_Lay"
0664 << id.layer();
0665 return result.str();
0666 }
0667
0668 string STAnalyzer::toString(const DTSuperLayerId& id) const {
0669 stringstream result;
0670 result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer();
0671 return result.str();
0672 }
0673
0674 string STAnalyzer::toString(const DTChamberId& id) const {
0675 stringstream result;
0676 result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station();
0677 return result.str();
0678 }
0679
0680 template <class T>
0681 string STAnalyzer::hName(const string& s, const T& id) const {
0682 string name(toString(id));
0683 stringstream hName;
0684 hName << s << name;
0685 return hName.str();
0686 }
0687
0688 void STAnalyzer::createTH1F(const std::string& name,
0689 const std::string& title,
0690 const std::string& suffix,
0691 int nbin,
0692 const double& binMin,
0693 const double& binMax) const {
0694 stringstream hName;
0695 stringstream hTitle;
0696 hName << name << suffix;
0697 hTitle << title << suffix;
0698 new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin, binMin, binMax);
0699 }
0700
0701 void STAnalyzer::createTH2F(const std::string& name,
0702 const std::string& title,
0703 const std::string& suffix,
0704 int nBinX,
0705 const double& binXMin,
0706 const double& binXMax,
0707 int nBinY,
0708 const double& binYMin,
0709 const double& binYMax) const {
0710 stringstream hName;
0711 stringstream hTitle;
0712 hName << name << suffix;
0713 hTitle << title << suffix;
0714 new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX, binXMin, binXMax, nBinY, binYMin, binYMax);
0715 }
0716
0717 DEFINE_FWK_MODULE(STAnalyzer);