Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:     CommonAlignmentProducer
00004 // Class  :     AlignmentMonitorMuonVsCurvature
00005 // 
00006 // Implementation:
00007 //     <Notes on implementation>
00008 //
00009 // Original Author:  Jim Pivarski
00010 //         Created:  Fri Feb 19 21:45:02 CET 2010
00011 //
00012 
00013 // system include files
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 // user include files
00038 
00039 // 
00040 // class definition
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 // constants, enums and typedefs
00118 //
00119 
00120 //
00121 // static data member definitions
00122 //
00123 
00124 //
00125 // member functions
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                         // derivatives are in local coordinates, so these should be, too
00379                         resid_x = dt13->hitresid(m_layer);
00380                         resid_dxdz = dt13->resslope();
00381 
00382                         // calculate derivative
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.) {  // if in tracker (cheap, I know...)
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                      } // if it's a good segment
00428                   } // if on our chamber
00429                } // if DT
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                         // derivatives are in local coordinates, so these should be, too
00483                         resid_x = csc->hitresid(m_layer);
00484                         resid_dxdz = csc->resslope();
00485 
00486                         // calculate derivative
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                      } // if it's a good segment
00552                   } // if on our chamber
00553                } // if CSC
00554 
00555             } // end loop over chamberIds
00556          } // end if refit is okay
00557       } // end if track pT is within range
00558    } // end loop over tracks
00559 }
00560 
00561 void AlignmentMonitorMuonVsCurvature::afterAlignment(const edm::EventSetup &iSetup) {}
00562 
00563 //
00564 // constructors and destructor
00565 //
00566 
00567 // AlignmentMonitorMuonVsCurvature::AlignmentMonitorMuonVsCurvature(const AlignmentMonitorMuonVsCurvature& rhs)
00568 // {
00569 //    // do actual copying here;
00570 // }
00571 
00572 //
00573 // assignment operators
00574 //
00575 // const AlignmentMonitorMuonVsCurvature& AlignmentMonitorMuonVsCurvature::operator=(const AlignmentMonitorMuonVsCurvature& rhs)
00576 // {
00577 //   //An exception safe implementation is
00578 //   AlignmentMonitorMuonVsCurvature temp(rhs);
00579 //   swap(rhs);
00580 //
00581 //   return *this;
00582 // }
00583 
00584 //
00585 // const member functions
00586 //
00587 
00588 //
00589 // static member functions
00590 //
00591 
00592 //
00593 // SEAL definitions
00594 //
00595 
00596 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorMuonVsCurvature, "AlignmentMonitorMuonVsCurvature");
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/Alignment/CommonAlignmentMonitor/plugins/AlignmentMonitorMuonVsCurvature.cc
