Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:31:53

0001 /*
0002  * Package:     CommonAlignmentProducer
0003  * Class  :     AlignmentMonitorMuonSystemMap1D
0004  *
0005  * Original Author:  Jim Pivarski
0006  *         Created:  Mon Nov 12 13:30:14 CST 2007
0007  *
0008  * $Id: AlignmentMonitorMuonSystemMap1D.cc,v 1.5 2011/04/15 23:09:37 khotilov Exp $
0009  */
0010 
0011 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
0012 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
0013 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
0014 
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/ServiceRegistry/interface/Service.h"
0021 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0022 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0023 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0024 
0025 #include "TH1F.h"
0026 #include "TH2F.h"
0027 
0028 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0029 #include "TrackingTools/Records/interface/DetIdAssociatorRecord.h"
0030 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0031 #include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
0032 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0033 #include "MagneticField/Engine/interface/MagneticField.h"
0034 
0035 class AlignmentMonitorMuonSystemMap1D : public AlignmentMonitorBase {
0036 public:
0037   AlignmentMonitorMuonSystemMap1D(const edm::ParameterSet &cfg, edm::ConsumesCollector iC);
0038   ~AlignmentMonitorMuonSystemMap1D() override = default;
0039   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0040 
0041   void book() override;
0042 
0043   void event(const edm::Event &iEvent,
0044              const edm::EventSetup &iSetup,
0045              const ConstTrajTrackPairCollection &iTrajTracks) override;
0046   void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft, const edm::Event &iEvent);
0047 
0048   void afterAlignment() override;
0049 
0050 private:
0051   // es token
0052   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_esTokenGBTGeom;
0053   const edm::ESGetToken<DetIdAssociator, DetIdAssociatorRecord> m_esTokenDetId;
0054   const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_esTokenProp;
0055   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_esTokenMF;
0056   const MuonResidualsFromTrack::BuilderToken m_esTokenBuilder;
0057 
0058   // parameters
0059   const edm::InputTag m_muonCollectionTag;
0060   const double m_minTrackPt;
0061   const double m_maxTrackPt;
0062   const double m_minTrackP;
0063   const double m_maxTrackP;
0064   const double m_maxDxy;
0065   const int m_minTrackerHits;
0066   const double m_maxTrackerRedChi2;
0067   const bool m_allowTIDTEC;
0068   const int m_minNCrossedChambers;
0069   const int m_minDT13Hits;
0070   const int m_minDT2Hits;
0071   const int m_minCSCHits;
0072   const bool m_doDT;
0073   const bool m_doCSC;
0074   const bool m_useStubPosition;
0075   const bool m_createNtuple;
0076   const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0077   const edm::EDGetTokenT<reco::MuonCollection> muonToken_;
0078 
0079   // counter
0080   long m_counter_event;
0081   long m_counter_track;
0082   long m_counter_trackmoment;
0083   long m_counter_trackdxy;
0084   long m_counter_trackokay;
0085   long m_counter_dt;
0086   long m_counter_13numhits;
0087   long m_counter_2numhits;
0088   long m_counter_csc;
0089   long m_counter_cscnumhits;
0090 
0091   // histogram helper
0092   class MuonSystemMapPlot1D {
0093   public:
0094     MuonSystemMapPlot1D(std::string name,
0095                         AlignmentMonitorMuonSystemMap1D *module,
0096                         int bins,
0097                         double low,
0098                         double high,
0099                         bool xy,
0100                         bool add_1d);
0101 
0102     void fill_x_1d(double residx, double chi2, int dof);
0103     void fill_x(char charge, double abscissa, double residx, double chi2, int dof);
0104     void fill_y(char charge, double abscissa, double residy, double chi2, int dof);
0105     void fill_dxdz(char charge, double abscissa, double slopex, double chi2, int dof);
0106     void fill_dydz(char charge, double abscissa, double slopey, double chi2, int dof);
0107 
0108   private:
0109     std::string m_name;
0110     int m_bins;
0111     bool m_xy;
0112     bool m_1d;
0113     TH1F *m_x_1d;
0114     TH2F *m_x_2d, *m_y_2d, *m_dxdz_2d, *m_dydz_2d;
0115   };
0116 
0117   MuonSystemMapPlot1D *m_DTvsz_station[4][14];   // [station][sector]
0118   MuonSystemMapPlot1D *m_CSCvsr_me[2][4][36];    // [endcap][station][chamber]
0119   MuonSystemMapPlot1D *m_DTvsphi_station[4][5];  // [station][wheel]
0120   MuonSystemMapPlot1D *m_CSCvsphi_me[2][4][3];   // [endcap][station][ring]
0121 
0122   std::vector<MuonSystemMapPlot1D *> m_plots;
0123 
0124   std::string num02d(int num);
0125 
0126   // optional debug ntuple
0127   TTree *m_cscnt;
0128 
0129   struct MyCSCDetId {
0130     void init(CSCDetId &id) {
0131       e = id.endcap();
0132       s = id.station();
0133       r = id.ring();
0134       c = id.chamber();
0135       t = id.iChamberType();
0136     }
0137     Short_t e, s, r, c;
0138     Short_t t;  // type 1-10: ME1/a,1/b,1/2,1/3,2/1...4/2
0139   };
0140   MyCSCDetId m_id;
0141 
0142   struct MyTrack {
0143     Int_t q;
0144     Float_t pt, pz;
0145   };
0146   MyTrack m_tr;
0147 
0148   struct MyResidual {
0149     Float_t res, slope, rho, phi, z;
0150   };
0151   MyResidual m_re;
0152 
0153   UInt_t m_run;
0154 };
0155 
0156 AlignmentMonitorMuonSystemMap1D::AlignmentMonitorMuonSystemMap1D(const edm::ParameterSet &cfg,
0157                                                                  edm::ConsumesCollector iC)
0158     : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorMuonSystemMap1D"),
0159       m_esTokenGBTGeom(iC.esConsumes()),
0160       m_esTokenDetId(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
0161       m_esTokenProp(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0162       m_esTokenMF(iC.esConsumes()),
0163       m_esTokenBuilder(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
0164       m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
0165       m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
0166       m_maxTrackPt(cfg.getParameter<double>("maxTrackPt")),
0167       m_minTrackP(cfg.getParameter<double>("minTrackP")),
0168       m_maxTrackP(cfg.getParameter<double>("maxTrackP")),
0169       m_maxDxy(cfg.getParameter<double>("maxDxy")),
0170       m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
0171       m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
0172       m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
0173       m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
0174       m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
0175       m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
0176       m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
0177       m_doDT(cfg.getParameter<bool>("doDT")),
0178       m_doCSC(cfg.getParameter<bool>("doCSC")),
0179       m_useStubPosition(cfg.getParameter<bool>("useStubPosition")),
0180       m_createNtuple(cfg.getParameter<bool>("createNtuple")),
0181       bsToken_(iC.consumes<reco::BeamSpot>(m_beamSpotTag)),
0182       muonToken_(iC.consumes<reco::MuonCollection>(m_muonCollectionTag)) {
0183   if (m_createNtuple) {
0184     edm::Service<TFileService> fs;
0185     m_cscnt = fs->make<TTree>("mualNtuple", "mualNtuple");
0186     m_cscnt->Branch("id", &m_id.e, "e/S:s:r:c:t");
0187     m_cscnt->Branch("tr", &m_tr.q, "q/I:pt/F:pz");
0188     m_cscnt->Branch("re", &m_re.res, "res/F:slope:rho:phi:z");
0189     m_cscnt->Branch("run", &m_run, "run/i");
0190   }
0191 }
0192 
0193 void AlignmentMonitorMuonSystemMap1D::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0194   edm::ParameterSetDescription desc;
0195   desc.add<edm::InputTag>("muonCollectionTag", edm::InputTag(""));
0196   desc.addUntracked<edm::InputTag>("beamSpotTag", edm::InputTag("offlineBeamSpot"));
0197   desc.add<double>("minTrackPt", 100.);
0198   desc.add<double>("maxTrackPt", 200.);
0199   desc.add<double>("minTrackP", 0.);
0200   desc.add<double>("maxTrackP", 99999.);
0201   desc.add<double>("maxDxy", 100.);
0202   desc.add<int>("minTrackerHits", 15);
0203   desc.add<double>("maxTrackerRedChi2", 10.);
0204   desc.add<bool>("allowTIDTEC", true);
0205   desc.add<int>("minNCrossedChambers", 3);
0206   desc.add<int>("minDT13Hits", 8);
0207   desc.add<int>("minDT2Hits", 4);
0208   desc.add<int>("minCSCHits", 6);
0209   desc.add<bool>("doDT", true);
0210   desc.add<bool>("doCSC", true);
0211   desc.add<bool>("useStubPosition", false);
0212   desc.add<bool>("createNtuple", false);
0213   descriptions.add("alignmentMonitorMuonSystemMap1D", desc);
0214 }
0215 
0216 std::string AlignmentMonitorMuonSystemMap1D::num02d(int num) {
0217   assert(num >= 0 && num < 100);
0218   char tmp[4];
0219   sprintf(tmp, "%02d", num);
0220   return std::string(tmp);
0221 }
0222 
0223 void AlignmentMonitorMuonSystemMap1D::book() {
0224   std::string wheel_label[5] = {"A", "B", "C", "D", "E"};
0225 
0226   for (unsigned char station = 1; station <= 4; station++) {
0227     std::string s_station = std::to_string(station);
0228 
0229     bool do_y = true;
0230     if (station == 4)
0231       do_y = false;
0232 
0233     // *** DT ***
0234     if (m_doDT)
0235       for (int sector = 1; sector <= 14; sector++) {
0236         if ((station < 4 && sector <= 12) || station == 4) {
0237           m_DTvsz_station[station - 1][sector - 1] = new MuonSystemMapPlot1D(
0238               "DTvsz_st" + s_station + "sec" + num02d(sector), this, 60, -660., 660., do_y, false);
0239           m_plots.push_back(m_DTvsz_station[station - 1][sector - 1]);
0240         }
0241       }
0242 
0243     if (m_doDT)
0244       for (int wheel = -2; wheel <= 2; wheel++) {
0245         m_DTvsphi_station[station - 1][wheel + 2] = new MuonSystemMapPlot1D(
0246             "DTvsphi_st" + s_station + "wh" + wheel_label[wheel + 2], this, 180, -M_PI, M_PI, do_y, false);
0247         m_plots.push_back(m_DTvsphi_station[station - 1][wheel + 2]);
0248       }
0249 
0250     // *** CSC ***
0251     if (m_doCSC)
0252       for (int endcap = 1; endcap <= 2; endcap++) {
0253         std::string s_endcap("m");
0254         if (endcap == 1)
0255           s_endcap = "p";
0256 
0257         for (int chamber = 1; chamber <= 36; chamber++) {
0258           m_CSCvsr_me[endcap - 1][station - 1][chamber - 1] = new MuonSystemMapPlot1D(
0259               "CSCvsr_me" + s_endcap + s_station + "ch" + num02d(chamber), this, 60, 100., 700., false, false);
0260           m_plots.push_back(m_CSCvsr_me[endcap - 1][station - 1][chamber - 1]);
0261         }
0262 
0263         for (int ring = 1; ring <= 3; ring++)  // the ME1/a (ring4) is not independent from ME1/b (ring1)
0264         {
0265           std::string s_ring = std::to_string(ring);
0266           if ((station > 1 && ring <= 2) || station == 1) {
0267             m_CSCvsphi_me[endcap - 1][station - 1][ring - 1] =
0268                 new MuonSystemMapPlot1D("CSCvsphi_me" + s_endcap + s_station + s_ring,
0269                                         this,
0270                                         180,
0271                                         -M_PI / 180. * 5.,
0272                                         M_PI * (2. - 5. / 180.),
0273                                         false,
0274                                         true);
0275             m_plots.push_back(m_CSCvsphi_me[endcap - 1][station - 1][ring - 1]);
0276           }
0277         }
0278       }  // endcaps
0279   }      // stations
0280 
0281   m_counter_event = 0;
0282   m_counter_track = 0;
0283   m_counter_trackmoment = 0;
0284   m_counter_trackdxy = 0;
0285   m_counter_trackokay = 0;
0286   m_counter_dt = 0;
0287   m_counter_13numhits = 0;
0288   m_counter_2numhits = 0;
0289   m_counter_csc = 0;
0290   m_counter_cscnumhits = 0;
0291 }
0292 
0293 void AlignmentMonitorMuonSystemMap1D::event(const edm::Event &iEvent,
0294                                             const edm::EventSetup &iSetup,
0295                                             const ConstTrajTrackPairCollection &trajtracks) {
0296   m_counter_event++;
0297 
0298   const edm::Handle<reco::BeamSpot> &beamSpot = iEvent.getHandle(bsToken_);
0299 
0300   const GlobalTrackingGeometry *globalGeometry = &iSetup.getData(m_esTokenGBTGeom);
0301   const DetIdAssociator *muonDetIdAssociator_ = &iSetup.getData(m_esTokenDetId);
0302   const Propagator *prop = &iSetup.getData(m_esTokenProp);
0303   const MagneticField *magneticField = &iSetup.getData(m_esTokenMF);
0304   auto builder = iSetup.getHandle(m_esTokenBuilder);
0305 
0306   if (m_muonCollectionTag.label().empty())  // use trajectories
0307   {
0308     for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
0309          ++trajtrack) {
0310       m_counter_track++;
0311       const Trajectory *traj = (*trajtrack).first;
0312       const reco::Track *track = (*trajtrack).second;
0313 
0314       if (m_minTrackPt < track->pt() && track->pt() < m_maxTrackPt && m_minTrackP < track->p() &&
0315           track->p() < m_maxTrackP) {
0316         m_counter_trackmoment++;
0317         if (fabs(track->dxy(beamSpot->position())) < m_maxDxy) {
0318           m_counter_trackdxy++;
0319 
0320           MuonResidualsFromTrack muonResidualsFromTrack(
0321               builder, magneticField, globalGeometry, muonDetIdAssociator_, prop, traj, track, pNavigator(), 1000.);
0322           processMuonResidualsFromTrack(muonResidualsFromTrack, iEvent);
0323         }
0324       }  // end if track has acceptable momentum
0325     }    // end loop over tracks
0326   } else {
0327     const edm::Handle<reco::MuonCollection> &muons = iEvent.getHandle(muonToken_);
0328 
0329     for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0330       if (!(muon->isTrackerMuon() && muon->innerTrack().isNonnull()))
0331         continue;
0332 
0333       m_counter_track++;
0334 
0335       if (m_minTrackPt < muon->pt() && muon->pt() < m_maxTrackPt && m_minTrackP < muon->p() &&
0336           muon->p() < m_maxTrackP) {
0337         m_counter_trackmoment++;
0338         if (fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy) {
0339           m_counter_trackdxy++;
0340 
0341           MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
0342           processMuonResidualsFromTrack(muonResidualsFromTrack, iEvent);
0343         }
0344       }
0345     }
0346   }
0347 }
0348 
0349 void AlignmentMonitorMuonSystemMap1D::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft,
0350                                                                     const edm::Event &iEvent) {
0351   if (mrft.trackerNumHits() < m_minTrackerHits)
0352     return;
0353   if (!m_allowTIDTEC && mrft.contains_TIDTEC())
0354     return;
0355   if (mrft.normalizedChi2() > m_maxTrackerRedChi2)
0356     return;
0357 
0358   int nMuChambers = 0;
0359   std::vector<DetId> chamberIds = mrft.chamberIds();
0360   for (unsigned ch = 0; ch < chamberIds.size(); ch++)
0361     if (chamberIds[ch].det() == DetId::Muon)
0362       nMuChambers++;
0363   if (nMuChambers < m_minNCrossedChambers)
0364     return;
0365 
0366   char charge = (mrft.getTrack()->charge() > 0 ? 1 : -1);
0367   // double qoverpt = track->charge() / track->pt();
0368   // double qoverpz = track->charge() / track->pz();
0369 
0370   m_counter_trackokay++;
0371 
0372   for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
0373     if (chamberId->det() != DetId::Muon)
0374       continue;
0375 
0376     if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT) {
0377       MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
0378       MuonChamberResidual *dt2 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT2);
0379       DTChamberId id(chamberId->rawId());
0380 
0381       m_counter_dt++;
0382 
0383       if (id.station() < 4 && dt13 != nullptr && dt13->numHits() >= m_minDT13Hits && dt2 != nullptr &&
0384           dt2->numHits() >= m_minDT2Hits && (dt2->chi2() / double(dt2->ndof())) < 2.0) {
0385         m_counter_13numhits++;
0386 
0387         double residual = dt13->global_residual();
0388         double resslope = dt13->global_resslope();
0389         double chi2 = dt13->chi2();
0390         int dof = dt13->ndof();
0391 
0392         align::GlobalPoint gpos;
0393         if (m_useStubPosition)
0394           gpos = dt13->global_stubpos();
0395         else
0396           gpos = dt13->global_trackpos();
0397         double phi = atan2(gpos.y(), gpos.x());
0398         double z = gpos.z();
0399 
0400         assert(1 <= id.sector() && id.sector() <= 14);
0401 
0402         m_DTvsz_station[id.station() - 1][id.sector() - 1]->fill_x(charge, z, residual, chi2, dof);
0403         m_DTvsz_station[id.station() - 1][id.sector() - 1]->fill_dxdz(charge, z, resslope, chi2, dof);
0404         m_DTvsphi_station[id.station() - 1][id.wheel() + 2]->fill_x(charge, phi, residual, chi2, dof);
0405         m_DTvsphi_station[id.station() - 1][id.wheel() + 2]->fill_dxdz(charge, phi, resslope, chi2, dof);
0406 
0407         m_counter_2numhits++;
0408 
0409         residual = dt2->global_residual();
0410         resslope = dt2->global_resslope();
0411         chi2 = dt2->chi2();
0412         dof = dt2->ndof();
0413 
0414         if (m_useStubPosition)
0415           gpos = dt2->global_stubpos();
0416         else
0417           gpos = dt2->global_trackpos();
0418         phi = atan2(gpos.y(), gpos.x());
0419         z = gpos.z();
0420 
0421         assert(1 <= id.sector() && id.sector() <= 14);
0422 
0423         m_DTvsz_station[id.station() - 1][id.sector() - 1]->fill_y(charge, z, residual, chi2, dof);
0424         m_DTvsz_station[id.station() - 1][id.sector() - 1]->fill_dydz(charge, z, resslope, chi2, dof);
0425         m_DTvsphi_station[id.station() - 1][id.wheel() + 2]->fill_y(charge, phi, residual, chi2, dof);
0426         m_DTvsphi_station[id.station() - 1][id.wheel() + 2]->fill_dydz(charge, phi, resslope, chi2, dof);
0427       }
0428 
0429       if (id.station() == 4 && dt13 != nullptr && dt13->numHits() >= m_minDT13Hits) {
0430         m_counter_13numhits++;
0431 
0432         double residual = dt13->global_residual();
0433         double resslope = dt13->global_resslope();
0434         double chi2 = dt13->chi2();
0435         int dof = dt13->ndof();
0436 
0437         align::GlobalPoint gpos;
0438         if (m_useStubPosition)
0439           gpos = dt13->global_stubpos();
0440         else
0441           gpos = dt13->global_trackpos();
0442         double phi = atan2(gpos.y(), gpos.x());
0443         double z = gpos.z();
0444 
0445         assert(1 <= id.sector() && id.sector() <= 14);
0446 
0447         m_DTvsz_station[id.station() - 1][id.sector() - 1]->fill_x(charge, z, residual, chi2, dof);
0448         m_DTvsz_station[id.station() - 1][id.sector() - 1]->fill_dxdz(charge, z, resslope, chi2, dof);
0449         m_DTvsphi_station[id.station() - 1][id.wheel() + 2]->fill_x(charge, phi, residual, chi2, dof);
0450         m_DTvsphi_station[id.station() - 1][id.wheel() + 2]->fill_dxdz(charge, phi, resslope, chi2, dof);
0451       }
0452     }
0453 
0454     else if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
0455       MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
0456       CSCDetId id(chamberId->rawId());
0457 
0458       int ring = id.ring();
0459       if (id.ring() == 4)
0460         ring = 1;  // combine ME1/a + ME1/b
0461 
0462       m_counter_csc++;
0463 
0464       if (csc != nullptr && csc->numHits() >= m_minCSCHits) {
0465         m_counter_cscnumhits++;
0466 
0467         double residual = csc->global_residual();
0468         double resslope = csc->global_resslope();
0469         double chi2 = csc->chi2();
0470         int dof = csc->ndof();
0471 
0472         align::GlobalPoint gpos;
0473         if (m_useStubPosition)
0474           gpos = csc->global_stubpos();
0475         else
0476           gpos = csc->global_trackpos();
0477         double phi = atan2(gpos.y(), gpos.x());
0478         // start phi from -5deg
0479         if (phi < -M_PI / 180. * 5.)
0480           phi += 2. * M_PI;
0481         double R = sqrt(pow(gpos.x(), 2) + pow(gpos.y(), 2));
0482 
0483         int chamber = id.chamber() - 1;
0484         if (id.station() > 1 && ring == 1)
0485           chamber *= 2;
0486 
0487         assert(1 <= id.endcap() && id.endcap() <= 2 && 0 <= chamber && chamber <= 35);
0488 
0489         if (R > 0.)
0490           m_CSCvsphi_me[id.endcap() - 1][id.station() - 1][ring - 1]->fill_x_1d(residual / R, chi2, dof);
0491 
0492         m_CSCvsr_me[id.endcap() - 1][id.station() - 1][chamber]->fill_x(charge, R, residual, chi2, dof);
0493         m_CSCvsr_me[id.endcap() - 1][id.station() - 1][chamber]->fill_dxdz(charge, R, resslope, chi2, dof);
0494         m_CSCvsphi_me[id.endcap() - 1][id.station() - 1][ring - 1]->fill_x(charge, phi, residual, chi2, dof);
0495         m_CSCvsphi_me[id.endcap() - 1][id.station() - 1][ring - 1]->fill_dxdz(charge, phi, resslope, chi2, dof);
0496 
0497         if (m_createNtuple && chi2 > 0.) {  //  &&  TMath::Prob(chi2, dof) < 0.95)
0498           m_id.init(id);
0499           m_tr.q = charge;
0500           m_tr.pt = mrft.getTrack()->pt();
0501           m_tr.pz = mrft.getTrack()->pz();
0502           m_re.res = residual;
0503           m_re.slope = resslope;
0504           m_re.rho = R;
0505           m_re.phi = phi;
0506           m_re.z = gpos.z();
0507           m_run = iEvent.id().run();
0508           m_cscnt->Fill();
0509         }
0510       }
0511     }
0512 
0513     //else { assert(false); }
0514   }  // end loop over chambers
0515 }
0516 
0517 void AlignmentMonitorMuonSystemMap1D::afterAlignment() {
0518   edm::LogVerbatim("AlignmentMuonSystemMap") << "AlignmentMonitorMuonSystemMap1D counters:";
0519   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_event      = " << m_counter_event;
0520   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_track      = " << m_counter_track;
0521   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_trackppt   = " << m_counter_trackmoment;
0522   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_trackdxy   = " << m_counter_trackdxy;
0523   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_trackokay  = " << m_counter_trackokay;
0524   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_dt         = " << m_counter_dt;
0525   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_13numhits  = " << m_counter_13numhits;
0526   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_2numhits   = " << m_counter_2numhits;
0527   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_csc        = " << m_counter_csc;
0528   edm::LogVerbatim("AlignmentMuonSystemMap") << " monitor m_counter_cscnumhits = " << m_counter_cscnumhits;
0529 }
0530 
0531 AlignmentMonitorMuonSystemMap1D::MuonSystemMapPlot1D::MuonSystemMapPlot1D(
0532     std::string name, AlignmentMonitorMuonSystemMap1D *module, int bins, double low, double high, bool xy, bool add_1d)
0533     : m_name(name), m_bins(bins), m_xy(xy), m_1d(add_1d) {
0534   m_x_2d = m_y_2d = m_dxdz_2d = m_dydz_2d = nullptr;
0535   std::stringstream name_x_2d, name_y_2d, name_dxdz_2d, name_dydz_2d;
0536   name_x_2d << m_name << "_x_2d";
0537   name_y_2d << m_name << "_y_2d";
0538   name_dxdz_2d << m_name << "_dxdz_2d";
0539   name_dydz_2d << m_name << "_dydz_2d";
0540 
0541   const int nbins = 200;
0542   const double window = 100.;
0543 
0544   m_x_2d = module->book2D("/iterN/", name_x_2d.str(), "", m_bins, low, high, nbins, -window, window);
0545   if (m_xy)
0546     m_y_2d = module->book2D("/iterN/", name_y_2d.str(), "", m_bins, low, high, nbins, -window, window);
0547   m_dxdz_2d = module->book2D("/iterN/", name_dxdz_2d.str(), "", m_bins, low, high, nbins, -window, window);
0548   if (m_xy)
0549     m_dydz_2d = module->book2D("/iterN/", name_dydz_2d.str(), "", m_bins, low, high, nbins, -window, window);
0550 
0551   m_x_1d = nullptr;
0552   if (m_1d) {
0553     std::stringstream name_x_1d;  //, name_y_1d, name_dxdz_1d, name_dydz_1d;
0554     name_x_1d << m_name << "_x_1d";
0555     m_x_1d = module->book1D("/iterN/", name_x_1d.str(), "", nbins, -window, window);
0556   }
0557 }
0558 
0559 void AlignmentMonitorMuonSystemMap1D::MuonSystemMapPlot1D::fill_x_1d(double residx, double chi2, int dof) {
0560   if (m_1d && chi2 > 0.) {
0561     // assume that residx was in radians
0562     double residual = residx * 1000.;
0563     m_x_1d->Fill(residual);
0564   }
0565 }
0566 
0567 void AlignmentMonitorMuonSystemMap1D::MuonSystemMapPlot1D::fill_x(
0568     char charge, double abscissa, double residx, double chi2, int dof) {
0569   if (chi2 > 0.) {
0570     double residual = residx * 10.;
0571     //double weight = dof / chi2;
0572     m_x_2d->Fill(abscissa, residual);
0573   }
0574 }
0575 
0576 void AlignmentMonitorMuonSystemMap1D::MuonSystemMapPlot1D::fill_y(
0577     char charge, double abscissa, double residy, double chi2, int dof) {
0578   if (m_xy && chi2 > 0.) {
0579     double residual = residy * 10.;
0580     //double weight = dof / chi2;
0581     m_y_2d->Fill(abscissa, residual);
0582   }
0583 }
0584 
0585 void AlignmentMonitorMuonSystemMap1D::MuonSystemMapPlot1D::fill_dxdz(
0586     char charge, double abscissa, double slopex, double chi2, int dof) {
0587   if (chi2 > 0.) {
0588     double residual = slopex * 1000.;
0589     //double weight = dof / chi2;
0590     m_dxdz_2d->Fill(abscissa, residual);
0591   }
0592 }
0593 
0594 void AlignmentMonitorMuonSystemMap1D::MuonSystemMapPlot1D::fill_dydz(
0595     char charge, double abscissa, double slopey, double chi2, int dof) {
0596   if (m_xy && chi2 > 0.) {
0597     double residual = slopey * 1000.;
0598     //double weight = dof / chi2;
0599     m_dydz_2d->Fill(abscissa, residual);
0600   }
0601 }
0602 
0603 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorMuonSystemMap1D, "AlignmentMonitorMuonSystemMap1D");