00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
00015 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00020 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00021
00022 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
00023 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsPositionFitter.h"
00024 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsAngleFitter.h"
00025 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"
00026
00027 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00028 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00029 #include "MagneticField/Engine/interface/MagneticField.h"
00030 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00031
00032 #include <sstream>
00033 #include "TProfile.h"
00034 #include "TH2F.h"
00035 #include "TH1F.h"
00036
00037
00038
00039
00040
00041
00042
00043 class AlignmentMonitorMuonVsCurvature: public AlignmentMonitorBase {
00044 public:
00045 AlignmentMonitorMuonVsCurvature(const edm::ParameterSet& cfg);
00046 ~AlignmentMonitorMuonVsCurvature() {};
00047
00048 void book();
00049 void event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& iTrajTracks);
00050 void afterAlignment(const edm::EventSetup &iSetup);
00051
00052 private:
00053 double m_minTrackPt;
00054 int m_minTrackerHits;
00055 double m_maxTrackerRedChi2;
00056 bool m_allowTIDTEC;
00057 double m_maxDxy;
00058 int m_minDT13Hits;
00059 int m_minDT2Hits;
00060 int m_minCSCHits;
00061 int m_layer;
00062 std::string m_propagator;
00063
00064 enum {
00065 kWheelMinus2 = 0,
00066 kWheelMinus1,
00067 kWheelZero,
00068 kWheelPlus1,
00069 kWheelPlus2,
00070 kWheelMEm11,
00071 kWheelMEm12,
00072 kWheelMEm13,
00073 kWheelMEm14,
00074 kWheelMEp11,
00075 kWheelMEp12,
00076 kWheelMEp13,
00077 kWheelMEp14,
00078 kNumWheels
00079 };
00080
00081 enum {
00082 kEndcapMEm11 = 0,
00083 kEndcapMEm12,
00084 kEndcapMEm13,
00085 kEndcapMEm14,
00086 kEndcapMEp11,
00087 kEndcapMEp12,
00088 kEndcapMEp13,
00089 kEndcapMEp14,
00090 kNumEndcap
00091 };
00092
00093 enum {
00094 kDeltaX = 0,
00095 kDeltaDxDz,
00096 kPtErr,
00097 kCurvErr,
00098 kNumComponents
00099 };
00100
00101 TH2F *th2f_wheelsector[kNumWheels][12][kNumComponents];
00102 TProfile *tprofile_wheelsector[kNumWheels][12][kNumComponents];
00103
00104 TH2F *th2f_evens[kNumEndcap][kNumComponents];
00105 TH2F *th2f_odds[kNumEndcap][kNumComponents];
00106 TProfile *tprofile_evens[kNumEndcap][kNumComponents];
00107 TProfile *tprofile_odds[kNumEndcap][kNumComponents];
00108
00109 TH2F *th2f_endcap[kNumEndcap][36][kNumComponents];
00110 TProfile *tprofile_endcap[kNumEndcap][36][kNumComponents];
00111
00112 TH1F *th1f_trackerRedChi2;
00113 TH1F *th1f_trackerRedChi2Diff;
00114 };
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 AlignmentMonitorMuonVsCurvature::AlignmentMonitorMuonVsCurvature(const edm::ParameterSet& cfg)
00129 : AlignmentMonitorBase(cfg, "AlignmentMonitorMuonVsCurvature")
00130 , m_minTrackPt(cfg.getParameter<double>("minTrackPt"))
00131 , m_minTrackerHits(cfg.getParameter<int>("minTrackerHits"))
00132 , m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2"))
00133 , m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC"))
00134 , m_maxDxy(cfg.getParameter<double>("maxDxy"))
00135 , m_minDT13Hits(cfg.getParameter<int>("minDT13Hits"))
00136 , m_minDT2Hits(cfg.getParameter<int>("minDT2Hits"))
00137 , m_minCSCHits(cfg.getParameter<int>("minCSCHits"))
00138 , m_layer(cfg.getParameter<int>("layer"))
00139 , m_propagator(cfg.getParameter<std::string>("propagator"))
00140 {}
00141
00142 void AlignmentMonitorMuonVsCurvature::book() {
00143 for (int wheel = 0; wheel < kNumWheels; wheel++) {
00144 std::stringstream wheelname;
00145 if (wheel == kWheelMinus2) wheelname << "wheelm2_";
00146 else if (wheel == kWheelMinus1) wheelname << "wheelm1_";
00147 else if (wheel == kWheelZero) wheelname << "wheelz_";
00148 else if (wheel == kWheelPlus1) wheelname << "wheelp1_";
00149 else if (wheel == kWheelPlus2) wheelname << "wheelp2_";
00150 else if (wheel == kWheelMEm11) wheelname << "wheelmem11_";
00151 else if (wheel == kWheelMEm12) wheelname << "wheelmem12_";
00152 else if (wheel == kWheelMEm13) wheelname << "wheelmem13_";
00153 else if (wheel == kWheelMEm14) wheelname << "wheelmem14_";
00154 else if (wheel == kWheelMEp11) wheelname << "wheelmep11_";
00155 else if (wheel == kWheelMEp12) wheelname << "wheelmep12_";
00156 else if (wheel == kWheelMEp13) wheelname << "wheelmep13_";
00157 else if (wheel == kWheelMEp14) wheelname << "wheelmep14_";
00158
00159 for (int sector = 0; sector < 12; sector++) {
00160 std::stringstream sectorname;
00161 if (sector == 0) sectorname << "sector01_";
00162 else if (sector == 1) sectorname << "sector02_";
00163 else if (sector == 2) sectorname << "sector03_";
00164 else if (sector == 3) sectorname << "sector04_";
00165 else if (sector == 4) sectorname << "sector05_";
00166 else if (sector == 5) sectorname << "sector06_";
00167 else if (sector == 6) sectorname << "sector07_";
00168 else if (sector == 7) sectorname << "sector08_";
00169 else if (sector == 8) sectorname << "sector09_";
00170 else if (sector == 9) sectorname << "sector10_";
00171 else if (sector == 10) sectorname << "sector11_";
00172 else if (sector == 11) sectorname << "sector12_";
00173
00174 for (int component = 0; component < kNumComponents; component++) {
00175 std::stringstream th2f_name, tprofile_name;
00176 th2f_name << "th2f_" << wheelname.str() << sectorname.str();
00177 tprofile_name << "tprofile_" << wheelname.str() << sectorname.str();
00178
00179 double minmax = 15.;
00180 if (component == kDeltaX) {
00181 th2f_name << "deltax";
00182 tprofile_name << "deltax";
00183 minmax = 15.;
00184 }
00185 else if (component == kDeltaDxDz) {
00186 th2f_name << "deltadxdz";
00187 tprofile_name << "deltadxdz";
00188 minmax = 15.;
00189 }
00190 else if (component == kPtErr) {
00191 th2f_name << "pterr";
00192 tprofile_name << "pterr";
00193 minmax = 15.;
00194 }
00195 else if (component == kCurvErr) {
00196 th2f_name << "curverr";
00197 tprofile_name << "curverr";
00198 minmax = 0.0015;
00199 }
00200
00201 if (component == kPtErr) {
00202 th2f_wheelsector[wheel][sector][component] = book2D("/iterN/", th2f_name.str().c_str(), "", 25, -200., 200., 30, -minmax, minmax);
00203 tprofile_wheelsector[wheel][sector][component] = bookProfile("/iterN/", tprofile_name.str().c_str(), "", 25, -200., 200.);
00204 }
00205 else {
00206 th2f_wheelsector[wheel][sector][component] = book2D("/iterN/", th2f_name.str().c_str(), "", 25, -0.05, 0.05, 30, -minmax, minmax);
00207 tprofile_wheelsector[wheel][sector][component] = bookProfile("/iterN/", tprofile_name.str().c_str(), "", 25, -0.05, 0.05);
00208 }
00209 }
00210 }
00211 }
00212
00213 for (int endcap = 0; endcap < kNumEndcap; endcap++) {
00214 std::stringstream endcapname;
00215 if (endcap == kEndcapMEm11) endcapname << "endcapmem11_";
00216 else if (endcap == kEndcapMEm12) endcapname << "endcapmem12_";
00217 else if (endcap == kEndcapMEm13) endcapname << "endcapmem13_";
00218 else if (endcap == kEndcapMEm14) endcapname << "endcapmem14_";
00219 else if (endcap == kEndcapMEp11) endcapname << "endcapmep11_";
00220 else if (endcap == kEndcapMEp12) endcapname << "endcapmep12_";
00221 else if (endcap == kEndcapMEp13) endcapname << "endcapmep13_";
00222 else if (endcap == kEndcapMEp14) endcapname << "endcapmep14_";
00223
00224 for (int component = 0; component < kNumComponents; component++) {
00225 std::stringstream componentname;
00226 double minmax = 15.;
00227 if (component == kDeltaX) {
00228 componentname << "deltax";
00229 minmax = 15.;
00230 }
00231 else if (component == kDeltaDxDz) {
00232 componentname << "deltadxdz";
00233 minmax = 15.;
00234 }
00235 else if (component == kPtErr) {
00236 componentname << "pterr";
00237 minmax = 15.;
00238 }
00239 else if (component == kCurvErr) {
00240 componentname << "curverr";
00241 minmax = 0.0015;
00242 }
00243
00244 std::stringstream th2f_evens_name, th2f_odds_name, tprofile_evens_name, tprofile_odds_name;
00245 th2f_evens_name << "th2f_" << endcapname.str() << "evens_" << componentname.str();
00246 th2f_odds_name << "th2f_" << endcapname.str() << "odds_" << componentname.str();
00247 tprofile_evens_name << "tprofile_" << endcapname.str() << "evens_" << componentname.str();
00248 tprofile_odds_name << "tprofile_" << endcapname.str() << "odds_" << componentname.str();
00249
00250 if (component == kPtErr) {
00251 th2f_evens[endcap][component] = book2D("/iterN/", th2f_evens_name.str().c_str(), "", 25, -200., 200., 30, -minmax, minmax);
00252 th2f_odds[endcap][component] = book2D("/iterN/", th2f_odds_name.str().c_str(), "", 25, -200., 200., 30, -minmax, minmax);
00253 tprofile_evens[endcap][component] = bookProfile("/iterN/", tprofile_evens_name.str().c_str(), "", 25, -200., 200.);
00254 tprofile_odds[endcap][component] = bookProfile("/iterN/", tprofile_odds_name.str().c_str(), "", 25, -200., 200.);
00255 }
00256 else {
00257 th2f_evens[endcap][component] = book2D("/iterN/", th2f_evens_name.str().c_str(), "", 25, -0.05, 0.05, 30, -minmax, minmax);
00258 th2f_odds[endcap][component] = book2D("/iterN/", th2f_odds_name.str().c_str(), "", 25, -0.05, 0.05, 30, -minmax, minmax);
00259 tprofile_evens[endcap][component] = bookProfile("/iterN/", tprofile_evens_name.str().c_str(), "", 25, -0.05, 0.05);
00260 tprofile_odds[endcap][component] = bookProfile("/iterN/", tprofile_odds_name.str().c_str(), "", 25, -0.05, 0.05);
00261 }
00262
00263 for (int chamber = 0; chamber < 36; chamber++) {
00264 std::stringstream chambername;
00265 if (chamber == 0) chambername << "chamber01_";
00266 else if (chamber == 1) chambername << "chamber02_";
00267 else if (chamber == 2) chambername << "chamber03_";
00268 else if (chamber == 3) chambername << "chamber04_";
00269 else if (chamber == 4) chambername << "chamber05_";
00270 else if (chamber == 5) chambername << "chamber06_";
00271 else if (chamber == 6) chambername << "chamber07_";
00272 else if (chamber == 7) chambername << "chamber08_";
00273 else if (chamber == 8) chambername << "chamber09_";
00274 else if (chamber == 9) chambername << "chamber10_";
00275 else if (chamber == 10) chambername << "chamber11_";
00276 else if (chamber == 11) chambername << "chamber12_";
00277 else if (chamber == 12) chambername << "chamber13_";
00278 else if (chamber == 13) chambername << "chamber14_";
00279 else if (chamber == 14) chambername << "chamber15_";
00280 else if (chamber == 15) chambername << "chamber16_";
00281 else if (chamber == 16) chambername << "chamber17_";
00282 else if (chamber == 17) chambername << "chamber18_";
00283 else if (chamber == 18) chambername << "chamber19_";
00284 else if (chamber == 19) chambername << "chamber20_";
00285 else if (chamber == 20) chambername << "chamber21_";
00286 else if (chamber == 21) chambername << "chamber22_";
00287 else if (chamber == 22) chambername << "chamber23_";
00288 else if (chamber == 23) chambername << "chamber24_";
00289 else if (chamber == 24) chambername << "chamber25_";
00290 else if (chamber == 25) chambername << "chamber26_";
00291 else if (chamber == 26) chambername << "chamber27_";
00292 else if (chamber == 27) chambername << "chamber28_";
00293 else if (chamber == 28) chambername << "chamber29_";
00294 else if (chamber == 29) chambername << "chamber30_";
00295 else if (chamber == 30) chambername << "chamber31_";
00296 else if (chamber == 31) chambername << "chamber32_";
00297 else if (chamber == 32) chambername << "chamber33_";
00298 else if (chamber == 33) chambername << "chamber34_";
00299 else if (chamber == 34) chambername << "chamber35_";
00300 else if (chamber == 35) chambername << "chamber36_";
00301
00302 std::stringstream th2f_name, tprofile_name;
00303 th2f_name << "th2f_" << endcapname.str() << chambername.str() << componentname.str();
00304 tprofile_name << "tprofile_" << endcapname.str() << chambername.str() << componentname.str();
00305
00306 if (component == kPtErr) {
00307 th2f_endcap[endcap][chamber][component] = book2D("/iterN/", th2f_name.str().c_str(), "", 25, -200., 200., 30, -minmax, minmax);
00308 tprofile_endcap[endcap][chamber][component] = bookProfile("/iterN/", tprofile_name.str().c_str(), "", 25, -200., 200.);
00309 }
00310 else {
00311 th2f_endcap[endcap][chamber][component] = book2D("/iterN/", th2f_name.str().c_str(), "", 25, -0.05, 0.05, 30, -minmax, minmax);
00312 tprofile_endcap[endcap][chamber][component] = bookProfile("/iterN/", tprofile_name.str().c_str(), "", 25, -0.05, 0.05);
00313 }
00314
00315 }
00316 }
00317 }
00318
00319 th1f_trackerRedChi2 = book1D("/iterN/", "trackerRedChi2", "Refit tracker reduced chi^2", 100, 0., 30.);
00320 th1f_trackerRedChi2Diff = book1D("/iterN/", "trackerRedChi2Diff", "Fit-minus-refit tracker reduced chi^2", 100, -5., 5.);
00321 }
00322
00323 void AlignmentMonitorMuonVsCurvature::event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& trajtracks) {
00324 edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
00325 iSetup.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
00326
00327 edm::ESHandle<Propagator> propagator;
00328 iSetup.get<TrackingComponentsRecord>().get(m_propagator, propagator);
00329
00330 edm::ESHandle<MagneticField> magneticField;
00331 iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
00332
00333 for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end(); ++trajtrack) {
00334 const Trajectory* traj = (*trajtrack).first;
00335 const reco::Track* track = (*trajtrack).second;
00336
00337 if (track->pt() > m_minTrackPt && fabs(track->dxy()) < m_maxDxy) {
00338 MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, traj, pNavigator(), 1000.);
00339
00340 if (muonResidualsFromTrack.trackerNumHits() >= m_minTrackerHits && muonResidualsFromTrack.trackerRedChi2() < m_maxTrackerRedChi2 && (m_allowTIDTEC || !muonResidualsFromTrack.contains_TIDTEC())) {
00341 std::vector<DetId> chamberIds = muonResidualsFromTrack.chamberIds();
00342
00343 double qoverpt = track->charge() / track->pt();
00344 double px = track->px();
00345 double py = track->py();
00346 double pz = track->pz();
00347
00348 th1f_trackerRedChi2->Fill(muonResidualsFromTrack.trackerRedChi2());
00349 th1f_trackerRedChi2Diff->Fill(track->normalizedChi2() - muonResidualsFromTrack.trackerRedChi2());
00350
00351 for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
00352 if (chamberId->det() == DetId::Muon && chamberId->subdetId() == MuonSubdetId::DT) {
00353 DTChamberId dtid(chamberId->rawId());
00354 if (dtid.station() == 1) {
00355
00356 MuonChamberResidual *dt13 = muonResidualsFromTrack.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
00357
00358 if (dt13 != NULL && dt13->numHits() >= m_minDT13Hits) {
00359 int wheel = -1;
00360 if (dtid.wheel() == -2) wheel = kWheelMinus2;
00361 else if (dtid.wheel() == -1) wheel = kWheelMinus1;
00362 else if (dtid.wheel() == 0) wheel = kWheelZero;
00363 else if (dtid.wheel() == 1) wheel = kWheelPlus1;
00364 else if (dtid.wheel() == 2) wheel = kWheelPlus2;
00365
00366 int sector = dtid.sector() - 1;
00367
00368 double resid_x = dt13->global_hitresid(m_layer);
00369 double resid_dxdz = dt13->global_resslope();
00370
00371 if (fabs(resid_x) < 10.) {
00372 th2f_wheelsector[wheel][sector][kDeltaX]->Fill(qoverpt, resid_x*10.);
00373 tprofile_wheelsector[wheel][sector][kDeltaX]->Fill(qoverpt, resid_x*10.);
00374 th2f_wheelsector[wheel][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00375 tprofile_wheelsector[wheel][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00376 }
00377
00378
00379 resid_x = dt13->hitresid(m_layer);
00380 resid_dxdz = dt13->resslope();
00381
00382
00383 TrajectoryStateOnSurface last_tracker_tsos;
00384 double last_tracker_R = 0.;
00385 std::vector<TrajectoryMeasurement> measurements = traj->measurements();
00386 for (std::vector<TrajectoryMeasurement>::const_iterator im = measurements.begin(); im != measurements.end(); ++im) {
00387 TrajectoryStateOnSurface tsos = im->forwardPredictedState();
00388 if (tsos.isValid()) {
00389 GlobalPoint pos = tsos.globalPosition();
00390 if (pos.perp() < 200. && fabs(pos.z()) < 400.) {
00391 if (pos.perp() > last_tracker_R) {
00392 last_tracker_tsos = tsos;
00393 last_tracker_R = pos.perp();
00394 }
00395 }
00396 }
00397 }
00398 if (last_tracker_R > 0.) {
00399 FreeTrajectoryState ts_rebuilt(last_tracker_tsos.globalPosition(),
00400 last_tracker_tsos.globalMomentum(),
00401 last_tracker_tsos.charge(),
00402 &*magneticField);
00403
00404 double factor = (last_tracker_tsos.globalMomentum().mag() + 1.) / last_tracker_tsos.globalMomentum().mag();
00405 FreeTrajectoryState ts_plus1GeV(last_tracker_tsos.globalPosition(),
00406 GlobalVector(factor*last_tracker_tsos.globalMomentum().x(), factor*last_tracker_tsos.globalMomentum().y(), factor*last_tracker_tsos.globalMomentum().z()),
00407 last_tracker_tsos.charge(),
00408 &*magneticField);
00409
00410 TrajectoryStateOnSurface extrapolation_rebuilt = propagator->propagate(ts_rebuilt, globalGeometry->idToDet(dt13->localid(m_layer))->surface());
00411 TrajectoryStateOnSurface extrapolation_plus1GeV = propagator->propagate(ts_plus1GeV, globalGeometry->idToDet(dt13->localid(m_layer))->surface());
00412
00413 if (extrapolation_rebuilt.isValid() && extrapolation_plus1GeV.isValid()) {
00414 double rebuiltx = extrapolation_rebuilt.localPosition().x();
00415 double plus1x = extrapolation_plus1GeV.localPosition().x();
00416
00417 if (fabs(resid_x) < 10.) {
00418 th2f_wheelsector[wheel][sector][kPtErr]->Fill(1./qoverpt, resid_x/(rebuiltx-plus1x)/sqrt(1. + pz*pz/(px*px + py*py))*fabs(qoverpt)*100.);
00419 tprofile_wheelsector[wheel][sector][kPtErr]->Fill(1./qoverpt, resid_x/(rebuiltx-plus1x)/sqrt(1. + pz*pz/(px*px + py*py))*fabs(qoverpt)*100.);
00420
00421 th2f_wheelsector[wheel][sector][kCurvErr]->Fill(qoverpt, resid_x/(rebuiltx-plus1x)*qoverpt/sqrt(px*px + py*py + pz*pz));
00422 tprofile_wheelsector[wheel][sector][kCurvErr]->Fill(qoverpt, resid_x/(rebuiltx-plus1x)*qoverpt/sqrt(px*px + py*py + pz*pz));
00423 }
00424 }
00425 }
00426
00427 }
00428 }
00429 }
00430
00431 if (chamberId->det() == DetId::Muon && chamberId->subdetId() == MuonSubdetId::CSC) {
00432 CSCDetId cscid(chamberId->rawId());
00433 if (cscid.station() == 1) {
00434
00435 MuonChamberResidual *csc = muonResidualsFromTrack.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
00436
00437 if (csc != NULL && csc->numHits() >= m_minCSCHits) {
00438 int wheel = -1;
00439 int endcap = -1;
00440 if (cscid.endcap() == 1 && cscid.ring() == 1) { wheel = kWheelMEp11; endcap = kEndcapMEp11; }
00441 else if (cscid.endcap() == 1 && cscid.ring() == 2) { wheel = kWheelMEp12; endcap = kEndcapMEp12; }
00442 else if (cscid.endcap() == 1 && cscid.ring() == 3) { wheel = kWheelMEp13; endcap = kEndcapMEp13; }
00443 else if (cscid.endcap() == 1 && cscid.ring() == 4) { wheel = kWheelMEp14; endcap = kEndcapMEp14; }
00444 else if (cscid.endcap() != 1 && cscid.ring() == 1) { wheel = kWheelMEm11; endcap = kEndcapMEm11; }
00445 else if (cscid.endcap() != 1 && cscid.ring() == 2) { wheel = kWheelMEm12; endcap = kEndcapMEm12; }
00446 else if (cscid.endcap() != 1 && cscid.ring() == 3) { wheel = kWheelMEm13; endcap = kEndcapMEm13; }
00447 else if (cscid.endcap() != 1 && cscid.ring() == 4) { wheel = kWheelMEm14; endcap = kEndcapMEm14; }
00448
00449 int chamber = cscid.chamber() - 1;
00450
00451 int sector = cscid.chamber() / 3;
00452 if (cscid.chamber() == 36) sector = 0;
00453
00454 double resid_x = csc->global_hitresid(m_layer);
00455 double resid_dxdz = csc->global_resslope();
00456
00457 if (fabs(resid_x) < 10.) {
00458 th2f_wheelsector[wheel][sector][kDeltaX]->Fill(qoverpt, resid_x*10.);
00459 tprofile_wheelsector[wheel][sector][kDeltaX]->Fill(qoverpt, resid_x*10.);
00460 th2f_wheelsector[wheel][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00461 tprofile_wheelsector[wheel][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00462
00463 th2f_endcap[endcap][chamber][kDeltaX]->Fill(qoverpt, resid_x*10.);
00464 tprofile_endcap[endcap][chamber][kDeltaX]->Fill(qoverpt, resid_x*10.);
00465 th2f_endcap[endcap][chamber][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00466 tprofile_endcap[endcap][chamber][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00467
00468 if (cscid.chamber() % 2 == 0) {
00469 th2f_evens[endcap][kDeltaX]->Fill(qoverpt, resid_x*10.);
00470 tprofile_evens[endcap][kDeltaX]->Fill(qoverpt, resid_x*10.);
00471 th2f_evens[endcap][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00472 tprofile_evens[endcap][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00473 }
00474 else {
00475 th2f_odds[endcap][kDeltaX]->Fill(qoverpt, resid_x*10.);
00476 tprofile_odds[endcap][kDeltaX]->Fill(qoverpt, resid_x*10.);
00477 th2f_odds[endcap][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00478 tprofile_odds[endcap][kDeltaDxDz]->Fill(qoverpt, resid_dxdz*1000.);
00479 }
00480 }
00481
00482
00483 resid_x = csc->hitresid(m_layer);
00484 resid_dxdz = csc->resslope();
00485
00486
00487 TrajectoryStateOnSurface last_tracker_tsos;
00488 double last_tracker_absZ = 0.;
00489 std::vector<TrajectoryMeasurement> measurements = traj->measurements();
00490 for (std::vector<TrajectoryMeasurement>::const_iterator im = measurements.begin(); im != measurements.end(); ++im) {
00491 TrajectoryStateOnSurface tsos = im->forwardPredictedState();
00492 if (tsos.isValid()) {
00493 GlobalPoint pos = tsos.globalPosition();
00494 if (pos.perp() < 200. && fabs(pos.z()) < 400.) {
00495 if (fabs(pos.z()) > last_tracker_absZ) {
00496 last_tracker_tsos = tsos;
00497 last_tracker_absZ = fabs(pos.z());
00498 }
00499 }
00500 }
00501 }
00502 if (last_tracker_absZ > 0.) {
00503 FreeTrajectoryState ts_rebuilt(last_tracker_tsos.globalPosition(),
00504 last_tracker_tsos.globalMomentum(),
00505 last_tracker_tsos.charge(),
00506 &*magneticField);
00507
00508 double factor = (last_tracker_tsos.globalMomentum().mag() + 1.) / last_tracker_tsos.globalMomentum().mag();
00509 FreeTrajectoryState ts_plus1GeV(last_tracker_tsos.globalPosition(),
00510 GlobalVector(factor*last_tracker_tsos.globalMomentum().x(), factor*last_tracker_tsos.globalMomentum().y(), factor*last_tracker_tsos.globalMomentum().z()),
00511 last_tracker_tsos.charge(),
00512 &*magneticField);
00513
00514 TrajectoryStateOnSurface extrapolation_rebuilt = propagator->propagate(ts_rebuilt, globalGeometry->idToDet(csc->localid(m_layer))->surface());
00515 TrajectoryStateOnSurface extrapolation_plus1GeV = propagator->propagate(ts_plus1GeV, globalGeometry->idToDet(csc->localid(m_layer))->surface());
00516
00517 if (extrapolation_rebuilt.isValid() && extrapolation_plus1GeV.isValid()) {
00518
00519 double rebuiltx = extrapolation_rebuilt.localPosition().x();
00520 double plus1x = extrapolation_plus1GeV.localPosition().x();
00521
00522 if (fabs(resid_x) < 10.) {
00523 double pterroverpt = resid_x/(rebuiltx-plus1x)/sqrt(1. + pz*pz/(px*px + py*py))*fabs(qoverpt)*100.;
00524 double curverror = resid_x/(rebuiltx-plus1x)*qoverpt/sqrt(px*px + py*py + pz*pz);
00525
00526 th2f_wheelsector[wheel][sector][kPtErr]->Fill(1./qoverpt, pterroverpt);
00527 tprofile_wheelsector[wheel][sector][kPtErr]->Fill(1./qoverpt, pterroverpt);
00528 th2f_wheelsector[wheel][sector][kCurvErr]->Fill(qoverpt, curverror);
00529 tprofile_wheelsector[wheel][sector][kCurvErr]->Fill(qoverpt, curverror);
00530
00531 th2f_endcap[endcap][chamber][kPtErr]->Fill(1./qoverpt, pterroverpt);
00532 tprofile_endcap[endcap][chamber][kPtErr]->Fill(1./qoverpt, pterroverpt);
00533 th2f_endcap[endcap][chamber][kCurvErr]->Fill(qoverpt, curverror);
00534 tprofile_endcap[endcap][chamber][kCurvErr]->Fill(qoverpt, curverror);
00535
00536 if (cscid.chamber() % 2 == 0) {
00537 th2f_evens[endcap][kPtErr]->Fill(1./qoverpt, pterroverpt);
00538 tprofile_evens[endcap][kPtErr]->Fill(1./qoverpt, pterroverpt);
00539 th2f_evens[endcap][kCurvErr]->Fill(qoverpt, curverror);
00540 tprofile_evens[endcap][kCurvErr]->Fill(qoverpt, curverror);
00541 }
00542 else {
00543 th2f_odds[endcap][kPtErr]->Fill(1./qoverpt, pterroverpt);
00544 tprofile_odds[endcap][kPtErr]->Fill(1./qoverpt, pterroverpt);
00545 th2f_odds[endcap][kCurvErr]->Fill(qoverpt, curverror);
00546 tprofile_odds[endcap][kCurvErr]->Fill(qoverpt, curverror);
00547 }
00548 }
00549 }
00550 }
00551 }
00552 }
00553 }
00554
00555 }
00556 }
00557 }
00558 }
00559 }
00560
00561 void AlignmentMonitorMuonVsCurvature::afterAlignment(const edm::EventSetup &iSetup) {}
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorMuonVsCurvature, "AlignmentMonitorMuonVsCurvature");
1.7.3