MuonCSCDigis/src/CSCStubResolutionValidation.cc

0001 #include <memory>
0002
0003 #include "Validation/MuonCSCDigis/interface/CSCStubResolutionValidation.h"
0004 #include "Validation/MuonCSCDigis/interface/CSCStubMatcher.h"
0005
0006 CSCStubResolutionValidation::CSCStubResolutionValidation(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0007     : CSCBaseValidation(pset) {
0008   const auto& simVertex = pset.getParameter<edm::ParameterSet>("simVertex");
0009   simVertexInput_ = iC.consumes<edm::SimVertexContainer>(simVertex.getParameter<edm::InputTag>("inputTag"));
0010   const auto& simTrack = pset.getParameter<edm::ParameterSet>("simTrack");
0011   simTrackInput_ = iC.consumes<edm::SimTrackContainer>(simTrack.getParameter<edm::InputTag>("inputTag"));
0012
0013   // Initialize stub matcher
0014   cscStubMatcher_ = std::make_unique<CSCStubMatcher>(pset, std::move(iC));
0015 }
0016
0017 CSCStubResolutionValidation::~CSCStubResolutionValidation() {}
0018
0019 //create folder for resolution histograms and book them
0020 void CSCStubResolutionValidation::bookHistograms(DQMStore::IBooker& iBooker) {
0021   for (int i = 1; i <= 10; ++i) {
0022     int j = i - 1;
0023     const std::string cn(CSCDetId::chamberName(i));
0024
0025     //Position resolution; CLCT
0026     std::string t1 = "CLCTPosRes_hs_" + cn;
0027     std::string t2 = "CLCTPosRes_qs_" + cn;
0028     std::string t3 = "CLCTPosRes_es_" + cn;
0029
0030     iBooker.setCurrentFolder("MuonCSCDigisV/CSCDigiTask/CLCT/Resolution/");
0031     posresCLCT_hs[j] = iBooker.book1D(
0032         t1, cn + " CLCT Position Resolution (1/2-strip prec.); Strip_{L1T} - Strip_{SIM}; Entries", 50, -1, 1);
0033     posresCLCT_qs[j] = iBooker.book1D(
0034         t2, cn + " CLCT Position Resolution (1/4-strip prec.); Strip_{L1T} - Strip_{SIM}; Entries", 50, -1, 1);
0035     posresCLCT_es[j] = iBooker.book1D(
0036         t3, cn + " CLCT Position Resolution (1/8-strip prec.); Strip_{L1T} - Strip_{SIM}; Entries", 50, -1, 1);
0037
0038     //Slope resolution; CLCT
0039     std::string t4 = "CLCTBendRes_" + cn;
0040
0041     bendresCLCT[j] =
0042         iBooker.book1D(t4, cn + " CLCT Bend Resolution; Slope_{L1T} - Slope_{SIM}; Entries", 50, -0.5, 0.5);
0043   }
0044 }
0045
0046 void CSCStubResolutionValidation::analyze(const edm::Event& e, const edm::EventSetup& eventSetup) {
0047   // Define handles
0048   edm::Handle<edm::SimTrackContainer> sim_tracks;
0049   edm::Handle<edm::SimVertexContainer> sim_vertices;
0050
0051   // Use token to retreive event information
0052   e.getByToken(simTrackInput_, sim_tracks);
0053   e.getByToken(simVertexInput_, sim_vertices);
0054
0055   // Initialize StubMatcher
0056   cscStubMatcher_->init(e, eventSetup);
0057
0058   const edm::SimTrackContainer& sim_track = *sim_tracks.product();
0059   const edm::SimVertexContainer& sim_vert = *sim_vertices.product();
0060
0061   // select simtracks for true muons
0062   edm::SimTrackContainer sim_track_selected;
0063   for (const auto& t : sim_track) {
0064     if (!isSimTrackGood(t))
0065       continue;
0066     sim_track_selected.push_back(t);
0067   }
0068
0069   // Skip events with no selected simtracks
0070   if (sim_track_selected.empty())
0071     return;
0072
0073   // Loop through good tracks, use corresponding vertex to match stubs, then fill hists of chambers where the stub appears.
0074   for (const auto& t : sim_track_selected) {
0075     std::vector<bool> hitCLCT(10);
0076
0077     std::vector<float> delta_fhs_clct(10);
0078     std::vector<float> delta_fqs_clct(10);
0079     std::vector<float> delta_fes_clct(10);
0080
0081     std::vector<float> dslope_clct(10);
0082
0083     // Match track to stubs with appropriate vertex
0084     cscStubMatcher_->match(t, sim_vert[t.vertIndex()]);
0085
0086     // Store matched stubs.
0087     // Key: ChamberID, Value : CSCStubDigiContainer
0088     const auto& clcts = cscStubMatcher_->clcts();
0089
0090     // CLCTs
0091     for (auto& [id, container] : clcts) {
0092       const CSCDetId cscId(id);
0093
0094       // get the best clct in chamber
0095       const auto& clct = cscStubMatcher_->bestClctInChamber(id);
0096       if (!clct.isValid())
0097         continue;
0098
0099       // ME1a CLCTs are saved in ME1b container. So the DetId need to be specified
0100       const bool isME11(cscId.station() == 1 and (cscId.ring() == 4 or cscId.ring() == 1));
0101       const bool isME1a(isME11 and clct.getKeyStrip() > CSCConstants::MAX_HALF_STRIP_ME1B);
0102       int ring = cscId.ring();
0103       if (isME1a)
0104         ring = 4;
0105       else if (isME11)
0106         ring = 1;
0107       CSCDetId cscId2(cscId.endcap(), cscId.station(), ring, cscId.chamber(), 0);
0108       auto id2 = cscId2.rawId();
0109
0110       // calculate deltastrip for ME1/a. Basically, we need to subtract 64 from the CLCT key strip to
0111       // compare with key strip as obtained through the fit to simhits positions.
0112       int deltaStrip = 0;
0113       if (isME1a)
0114         deltaStrip = CSCConstants::NUM_STRIPS_ME1B;
0115
0116       // fractional strip
0117       const float fhs_clct = clct.getFractionalStrip(2);
0118       const float fqs_clct = clct.getFractionalStrip(4);
0119       const float fes_clct = clct.getFractionalStrip(8);
0120
0121       // in half-strips per layer
0122       const float slopeHalfStrip(clct.getFractionalSlope());
0123       const float slopeStrip(slopeHalfStrip / 2.);
0124
0125       // get the fit hits in chamber for true value
0126       float stripIntercept, stripSlope;
0127       cscStubMatcher_->cscDigiMatcher()->muonSimHitMatcher()->fitHitsInChamber(id2, stripIntercept, stripSlope);
0128
0129       // add offset of +0.25 strips for non-ME1/1 chambers
0130       if (!isME11) {
0131         stripIntercept -= 0.25;
0132       }
0133
0134       const float strip_csc_sh = stripIntercept;
0135       const float bend_csc_sh = stripSlope;
0136
0137       const unsigned chamberType(cscId2.iChamberType());
0138       hitCLCT[chamberType - 1] = true;
0139
0140       delta_fhs_clct[chamberType - 1] = fhs_clct - deltaStrip - strip_csc_sh;
0141       delta_fqs_clct[chamberType - 1] = fqs_clct - deltaStrip - strip_csc_sh;
0142       delta_fes_clct[chamberType - 1] = fes_clct - deltaStrip - strip_csc_sh;
0143
0144       dslope_clct[chamberType - 1] = slopeStrip - bend_csc_sh;
0145     }
0146
0147     for (int i = 0; i < 10; i++) {
0148       if (hitCLCT[i]) {
0149         //fill histograms
0150         posresCLCT_hs[i]->Fill(delta_fhs_clct[i]);
0151         posresCLCT_qs[i]->Fill(delta_fqs_clct[i]);
0152         posresCLCT_es[i]->Fill(delta_fes_clct[i]);
0153
0154         bendresCLCT[i]->Fill(dslope_clct[i]);
0155       }
0156     }
0157   }
0158 }