Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:33:22

0001 #include "DataFormats/Scalers/interface/DcsStatus.h"
0002 #include "Validation/MuonME0Validation/interface/ME0SegmentsValidation.h"
0003 #include <TMath.h>
0004 
0005 ME0SegmentsValidation::ME0SegmentsValidation(const edm::ParameterSet &cfg) : ME0BaseValidation(cfg) {
0006   InputTagToken_Segments = consumes<ME0SegmentCollection>(cfg.getParameter<edm::InputTag>("segmentInputLabel"));
0007   InputTagToken_Digis = consumes<ME0DigiPreRecoCollection>(cfg.getParameter<edm::InputTag>("digiInputLabel"));
0008   InputTagToken_ = consumes<edm::PSimHitContainer>(cfg.getParameter<edm::InputTag>("simInputLabel"));
0009   InputTagTokenST_ = consumes<edm::SimTrackContainer>(cfg.getParameter<edm::InputTag>("simInputLabelST"));
0010   sigma_x_ = cfg.getParameter<double>("sigma_x");
0011   sigma_y_ = cfg.getParameter<double>("sigma_y");
0012   eta_max_ = cfg.getParameter<double>("eta_max");
0013   eta_min_ = cfg.getParameter<double>("eta_min");
0014   pt_min_ = cfg.getParameter<double>("pt_min");
0015   isMuonGun_ = cfg.getParameter<bool>("isMuonGun");
0016 }
0017 
0018 void ME0SegmentsValidation::bookHistograms(DQMStore::IBooker &ibooker,
0019                                            edm::Run const &Run,
0020                                            edm::EventSetup const &iSetup) {
0021   LogDebug("MuonME0SegmentsValidation") << "Info : Loading Geometry information\n";
0022   ibooker.setCurrentFolder("MuonME0RecHitsV/ME0SegmentsTask");
0023 
0024   unsigned int nregion = 2;
0025 
0026   edm::LogInfo("MuonME0SegmentsValidation") << "+++ Info : # of region : " << nregion << std::endl;
0027 
0028   LogDebug("MuonME0SegmentsValidation") << "+++ Info : finish to get geometry information from ES.\n";
0029 
0030   me0_simsegment_eta = ibooker.book1D("me0_simsegment_eta", "SimSegment Eta Distribution; #eta; entries", 8, 2.0, 2.8);
0031   me0_simsegment_pt = ibooker.book1D("me0_simsegment_pt", "SimSegment pT Distribution; p_{T}; entries", 20, 0.0, 100.0);
0032   me0_simsegment_phi =
0033       ibooker.book1D("me0_simsegment_phi", "SimSegments phi Distribution; #phi; entries", 18, -M_PI, +M_PI);
0034 
0035   me0_matchedsimsegment_eta =
0036       ibooker.book1D("me0_matchedsimsegment_eta", "Matched SimSegment Eta Distribution; #eta; entries", 8, 2.0, 2.8);
0037   me0_matchedsimsegment_pt =
0038       ibooker.book1D("me0_matchedsimsegment_pt", "Matched SimSegment pT Distribution; p_{T}; entries", 20, 0.0, 100.0);
0039   me0_matchedsimsegment_phi = ibooker.book1D(
0040       "me0_matchedsimsegment_phi", "Matched SimSegments phi Distribution; #phi; entries", 18, -M_PI, +M_PI);
0041 
0042   me0_segment_chi2 = ibooker.book1D("me0_seg_Chi2", "#chi^{2}; #chi^{2}; # Segments", 100, 0, 100);
0043   me0_segment_redchi2 = ibooker.book1D("me0_seg_ReducedChi2", "#chi^{2}/ndof; #chi^{2}/ndof; # Segments", 100, 0, 5);
0044   me0_segment_ndof = ibooker.book1D("me0_seg_ndof", "ndof; ndof; #Segments", 50, 0, 50);
0045   me0_segment_numRH =
0046       ibooker.book1D("me0_seg_NumberRH", "Number of fitted RecHits; # RecHits; entries", 11, -0.5, 10.5);
0047   me0_segment_numRHSig =
0048       ibooker.book1D("me0_seg_NumberRHSig", "Number of fitted Signal RecHits; # RecHits; entries", 11, -0.5, 10.5);
0049   me0_segment_numRHBkg =
0050       ibooker.book1D("me0_seg_NumberRHBkg", "Number of fitted BKG RecHits; # RecHits; entries", 11, -0.5, 10.5);
0051   // me0_segment_EtaRH   = ibooker.book1D("me0_specRH_globalEta","Fitted RecHits
0052   // Eta Distribution; #eta; entries",200,-4.0,4.0); me0_segment_PhiRH   =
0053   // ibooker.book1D("me0_specRH_globalPhi","Fitted RecHits Phi Distribution;
0054   // #eta; entries",18,-3.14,3.14);
0055   me0_segment_time = ibooker.book1D("me0_seg_time", "Segment Timing; ns; entries", 300, -150, 150);
0056   me0_segment_timeErr = ibooker.book1D("me0_seg_timErr", "Segment Timing Error; ns; entries", 50, 0, 0.5);
0057   me0_segment_size =
0058       ibooker.book1D("me0_seg_size", "Segment Multiplicity; Number of ME0 segments; entries", 200, 0, 200);
0059 
0060   for (unsigned int region_num = 0; region_num < nregion; region_num++) {
0061     me0_specRH_zr[region_num] = BookHistZR(ibooker, "me0_specRH_tot", "Segment RecHits", region_num);
0062     for (unsigned int layer_num = 0; layer_num < 6; layer_num++) {
0063       // me0_strip_dg_zr[region_num][layer_num] =
0064       // BookHistZR(ibooker,"me0_strip_dg","SimHit",region_num,layer_num);
0065       me0_specRH_xy[region_num][layer_num] =
0066           BookHistXY(ibooker, "me0_specRH", "Segment RecHits", region_num, layer_num);
0067       // me0_rh_xy_Muon[region_num][layer_num] =
0068       // BookHistXY(ibooker,"me0_rh","RecHit Muon",region_num,layer_num);
0069 
0070       std::string histo_name_DeltaX =
0071           std::string("me0_specRH_DeltaX_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0072       std::string histo_name_DeltaY =
0073           std::string("me0_specRH_DeltaY_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0074       std::string histo_label_DeltaX = "Segment RecHits Delta X : region" + regionLabel[region_num] + " layer " +
0075                                        layerLabel[layer_num] + " " + " ; x_{SimHit} - x_{Segment RecHits} ; entries";
0076       std::string histo_label_DeltaY = "Segment RecHits Delta Y : region" + regionLabel[region_num] + " layer " +
0077                                        layerLabel[layer_num] + " " + " ; y_{SimHit} - y_{Segment RecHit} ; entries";
0078 
0079       me0_specRH_DeltaX[region_num][layer_num] =
0080           ibooker.book1D(histo_name_DeltaX.c_str(), histo_label_DeltaX.c_str(), 100, -10, 10);
0081       me0_specRH_DeltaY[region_num][layer_num] =
0082           ibooker.book1D(histo_name_DeltaY.c_str(), histo_label_DeltaY.c_str(), 100, -10, 10);
0083 
0084       std::string histo_name_PullX =
0085           std::string("me0_specRH_PullX_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0086       std::string histo_name_PullY =
0087           std::string("me0_specRH_PullY_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0088       std::string histo_label_PullX = "Segment RecHits Pull X : region" + regionLabel[region_num] + " layer " +
0089                                       layerLabel[layer_num] + " " +
0090                                       " ; #frac{x_{SimHit} - x_{Segment "
0091                                       "RecHit}}{#sigma_{x,RecHit}} ; entries";
0092       std::string histo_label_PullY = "Segment RecHits Pull Y : region" + regionLabel[region_num] + " layer " +
0093                                       layerLabel[layer_num] + " " +
0094                                       " ; #frac{y_{SimHit} - y_{Segment "
0095                                       "RecHit}}{#sigma_{y,RecHit}} ; entries";
0096 
0097       me0_specRH_PullX[region_num][layer_num] =
0098           ibooker.book1D(histo_name_PullX.c_str(), histo_label_DeltaX.c_str(), 100, -10, 10);
0099       me0_specRH_PullY[region_num][layer_num] =
0100           ibooker.book1D(histo_name_PullY.c_str(), histo_label_DeltaY.c_str(), 100, -10, 10);
0101     }
0102   }
0103 }
0104 
0105 ME0SegmentsValidation::~ME0SegmentsValidation() {}
0106 
0107 void ME0SegmentsValidation::analyze(const edm::Event &e, const edm::EventSetup &iSetup) {
0108   const ME0Geometry *ME0Geometry_ = &iSetup.getData(geomToken_);
0109 
0110   edm::Handle<edm::PSimHitContainer> ME0Hits;
0111   e.getByToken(InputTagToken_, ME0Hits);
0112 
0113   edm::Handle<edm::SimTrackContainer> simTracks;
0114   e.getByToken(InputTagTokenST_, simTracks);
0115 
0116   edm::Handle<ME0SegmentCollection> ME0Segments;
0117   e.getByToken(InputTagToken_Segments, ME0Segments);
0118 
0119   edm::Handle<ME0DigiPreRecoCollection> ME0Digis;
0120   e.getByToken(InputTagToken_Digis, ME0Digis);
0121 
0122   if (!ME0Digis.isValid()) {
0123     edm::LogError("ME0SegmentsValidation") << "Cannot get ME0Digis by Token InputTagToken";
0124     return;
0125   }
0126 
0127   if (!ME0Segments.isValid()) {
0128     edm::LogError("ME0SegmentsValidation") << "Cannot get ME0RecHits/ME0Segments by Token InputTagToken";
0129     return;
0130   }
0131 
0132   if (!ME0Hits.isValid()) {
0133     edm::LogError("ME0HitsValidation") << "Cannot get ME0Hits by Token simInputTagToken";
0134     return;
0135   }
0136 
0137   MapTypeSim myMap;
0138   MapTypeSeg myMapSeg;
0139 
0140   int countST = 0;
0141 
0142   edm::SimTrackContainer::const_iterator simTrack;
0143   for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0144     edm::PSimHitContainer selectedME0Hits;
0145 
0146     if (!isSimTrackGood(simTrack))
0147       continue;
0148     int count = 0;
0149 
0150     for (edm::PSimHitContainer::const_iterator itHit = ME0Hits->begin(); itHit != ME0Hits->end(); ++itHit) {
0151       int particleType_sh = itHit->particleType();
0152       int evtId_sh = itHit->eventId().event();
0153       int bx_sh = itHit->eventId().bunchCrossing();
0154       int procType_sh = itHit->processType();
0155       if (!(abs(particleType_sh) == 13 && evtId_sh == 0 && bx_sh == 0 && procType_sh == 0))
0156         continue;
0157 
0158       if (isSimMatched(simTrack, itHit)) {
0159         ++count;
0160         selectedME0Hits.push_back(*itHit);
0161         ;
0162       }
0163 
0164     }  // End loop SHs
0165 
0166     if (selectedME0Hits.size() >= 3) {
0167       myMap.insert(MapTypeSim::value_type(simTrack, selectedME0Hits));
0168       me0_simsegment_eta->Fill(std::abs((*simTrack).momentum().eta()));
0169       me0_simsegment_pt->Fill((*simTrack).momentum().pt());
0170       me0_simsegment_phi->Fill((*simTrack).momentum().phi());
0171       ++countST;
0172     }
0173   }
0174 
0175   me0_segment_size->Fill(ME0Segments->size());
0176 
0177   for (auto me0s = ME0Segments->begin(); me0s != ME0Segments->end(); me0s++) {
0178     // The ME0 Ensamble DetId refers to layer = 1
0179     ME0DetId id = me0s->me0DetId();
0180     auto chamber = ME0Geometry_->chamber(id);
0181     auto segLP = me0s->localPosition();
0182     auto segLD = me0s->localDirection();
0183     auto me0rhs = me0s->specificRecHits();
0184 
0185     //   float localX = segLP.x();
0186     //   float localY = segLP.y();
0187     //   float dirTheta = segLD.theta();
0188     //   float dirPhi = segLD.phi();
0189     int numberRH = me0rhs.size();
0190     float chi2 = (float)me0s->chi2();
0191     float ndof = me0s->degreesOfFreedom();
0192     double time = me0s->time();
0193     double timeErr = me0s->timeErr();
0194 
0195     float reducedChi2 = chi2 / ndof;
0196 
0197     me0_segment_chi2->Fill(chi2);
0198     me0_segment_redchi2->Fill(reducedChi2);
0199     me0_segment_ndof->Fill(ndof);
0200     me0_segment_numRH->Fill(numberRH);
0201 
0202     me0_segment_time->Fill(time);
0203     me0_segment_timeErr->Fill(timeErr);
0204 
0205     int numberRHSig = 0;
0206     int numberRHBkg = 0;
0207     std::vector<ME0RecHit> selectedME0RecHits;
0208 
0209     for (auto rh = me0rhs.begin(); rh != me0rhs.end(); rh++) {
0210       auto me0id = rh->me0Id();
0211       auto rhr = ME0Geometry_->etaPartition(me0id);
0212       auto rhLP = rh->localPosition();
0213 
0214       auto result = isMatched(me0id, rhLP, ME0Digis);
0215       if (result.second == 1) {
0216         ++numberRHSig;
0217         selectedME0RecHits.push_back(*rh);
0218 
0219       } else
0220         ++numberRHBkg;
0221 
0222       auto erhLEP = rh->localPositionError();
0223       auto rhGP = rhr->toGlobal(rhLP);
0224       auto rhLPSegm = chamber->toLocal(rhGP);
0225       float xe = segLP.x() + segLD.x() * rhLPSegm.z() / segLD.z();
0226       float ye = segLP.y() + segLD.y() * rhLPSegm.z() / segLD.z();
0227       float ze = rhLPSegm.z();
0228       LocalPoint extrPoint(xe, ye, ze);                           // in segment rest frame
0229       auto extSegm = rhr->toLocal(chamber->toGlobal(extrPoint));  // in layer restframe
0230 
0231       int region = me0id.region();
0232       int layer = me0id.layer();
0233       //     int chamber = me0id.chamber();
0234 
0235       float x = rhLP.x();
0236       float xErr = erhLEP.xx();
0237       float y = rhLP.y();
0238       float yErr = erhLEP.yy();
0239 
0240       float globalR = rhGP.perp();
0241       float globalX = rhGP.x();
0242       float globalY = rhGP.y();
0243       float globalZ = rhGP.z();
0244 
0245       float xExt = extSegm.x();
0246       float yExt = extSegm.y();
0247 
0248       float pull_x = (x - xExt) / sqrt(xErr);
0249       float pull_y = (y - yExt) / sqrt(yErr);
0250 
0251       int region_num = 0;
0252       if (region == -1)
0253         region_num = 0;
0254       else if (region == 1)
0255         region_num = 1;
0256       int layer_num = layer - 1;
0257 
0258       me0_specRH_xy[region_num][layer_num]->Fill(globalX, globalY);
0259       me0_specRH_zr[region_num]->Fill(globalZ, globalR);
0260 
0261       me0_specRH_DeltaX[region_num][layer_num]->Fill(x - xExt);
0262       me0_specRH_DeltaY[region_num][layer_num]->Fill(y - yExt);
0263       me0_specRH_PullX[region_num][layer_num]->Fill(pull_x);
0264       me0_specRH_PullY[region_num][layer_num]->Fill(pull_y);
0265     }
0266 
0267     me0_segment_numRHSig->Fill(numberRHSig);
0268     me0_segment_numRHBkg->Fill(numberRHBkg);
0269     myMapSeg.insert(MapTypeSeg::value_type(me0s, selectedME0RecHits));
0270   }
0271 
0272   //------------------- SimToReco -------------------
0273 
0274   for (auto const &st : myMap) {  // loop over the signal simTracks
0275 
0276     int num_sh = st.second.size();
0277     bool isThereOneSegmentMatched = false;
0278 
0279     for (auto const &seg : myMapSeg) {  // loop over the reconstructed me0 segments
0280 
0281       int num_sh_matched = 0;
0282       if (seg.second.empty())
0283         continue;
0284 
0285       for (auto const &sh : st.second) {  // loop over the me0 simHits left by
0286                                           // the signal simTracks
0287 
0288         for (auto const &rh : seg.second) {  // loop over the tracking recHits
0289                                              // already matched to signal digis
0290 
0291           auto me0id = rh.me0Id();
0292           int region_rh = (int)me0id.region();
0293           int layer_rh = (int)me0id.layer();
0294           int chamber_rh = (int)me0id.chamber();
0295           int roll_rh = (int)me0id.roll();
0296 
0297           const ME0DetId id(sh.detUnitId());
0298           int region_sh = id.region();
0299           int layer_sh = id.layer();
0300           int chamber_sh = id.chamber();
0301           int roll_sh = id.roll();
0302 
0303           if (!(region_sh == region_rh && chamber_sh == chamber_rh && layer_sh == layer_rh && roll_sh == roll_rh))
0304             continue;
0305 
0306           LocalPoint lp_sh = sh.localPosition();
0307           LocalPoint lp_rh = rh.localPosition();
0308 
0309           GlobalPoint gp_sh = ME0Geometry_->idToDet(id)->surface().toGlobal(lp_sh);
0310           GlobalPoint gp = ME0Geometry_->idToDet((rh).me0Id())->surface().toGlobal(lp_rh);
0311           float dphi_glob = gp_sh.phi() - gp.phi();
0312           float deta_glob = gp_sh.eta() - gp.eta();
0313 
0314           if (fabs(dphi_glob) < 3 * sigma_x_ && fabs(deta_glob) < 3 * sigma_y_)
0315             ++num_sh_matched;
0316 
0317         }  // End loop over RHs
0318 
0319       }  // End loop over SHs
0320 
0321       float quality = 0;
0322       if (num_sh != 0)
0323         quality = num_sh_matched / (1.0 * num_sh);
0324       if (quality > 0)
0325         isThereOneSegmentMatched = true;
0326 
0327     }  // End loop over segments
0328 
0329     // Fill hsitograms
0330     if (isThereOneSegmentMatched) {
0331       me0_matchedsimsegment_eta->Fill(std::abs((*(st.first)).momentum().eta()));
0332       me0_matchedsimsegment_pt->Fill((*(st.first)).momentum().pt());
0333       me0_matchedsimsegment_phi->Fill((*(st.first)).momentum().phi());
0334     }
0335 
0336   }  // End loop over STs
0337 }
0338 
0339 std::pair<int, int> ME0SegmentsValidation::isMatched(ME0DetId me0id,
0340                                                      LocalPoint rhLP,
0341                                                      edm::Handle<ME0DigiPreRecoCollection> ME0Digis) {
0342   int region_rh = (int)me0id.region();
0343   int layer_rh = (int)me0id.layer();
0344   int roll_rh = (int)me0id.roll();
0345   int chamber_rh = (int)me0id.chamber();
0346 
0347   float l_x_rh = rhLP.x();
0348   float l_y_rh = rhLP.y();
0349 
0350   int particleType = 0;
0351   int isPrompt = -1;
0352 
0353   for (ME0DigiPreRecoCollection::DigiRangeIterator cItr = ME0Digis->begin(); cItr != ME0Digis->end(); cItr++) {
0354     ME0DetId id = (*cItr).first;
0355 
0356     int region_dg = (int)id.region();
0357     int layer_dg = (int)id.layer();
0358     int roll_dg = (int)id.roll();
0359     int chamber_dg = (int)id.chamber();
0360 
0361     if (region_rh != region_dg)
0362       continue;
0363     if (layer_rh != layer_dg)
0364       continue;
0365     if (chamber_rh != chamber_dg)
0366       continue;
0367     if (roll_rh != roll_dg)
0368       continue;
0369 
0370     ME0DigiPreRecoCollection::const_iterator digiItr;
0371     for (digiItr = (*cItr).second.first; digiItr != (*cItr).second.second; ++digiItr) {
0372       float l_x_dg = digiItr->x();
0373       float l_y_dg = digiItr->y();
0374 
0375       if (l_x_rh != l_x_dg)
0376         continue;
0377       if (l_y_rh != l_y_dg)
0378         continue;
0379 
0380       particleType = digiItr->pdgid();
0381       isPrompt = digiItr->prompt();
0382     }
0383   }
0384 
0385   std::pair<int, int> result;
0386   result = std::make_pair(particleType, isPrompt);
0387 
0388   return result;
0389 }
0390 
0391 bool ME0SegmentsValidation::isSimTrackGood(edm::SimTrackContainer::const_iterator t) {
0392   if ((*t).noVertex() && !isMuonGun_)
0393     return false;
0394   if ((*t).noGenpart() && !isMuonGun_)
0395     return false;
0396   if (std::abs((*t).type()) != 13)
0397     return false;  // only interested in direct muon simtracks
0398   if ((*t).momentum().pt() < pt_min_)
0399     return false;
0400   const float eta(std::abs((*t).momentum().eta()));
0401   if (eta < eta_min_ || eta > eta_max_)
0402     return false;  // no GEMs could be in such eta
0403   return true;
0404 }
0405 
0406 bool ME0SegmentsValidation::isSimMatched(edm::SimTrackContainer::const_iterator simTrack,
0407                                          edm::PSimHitContainer::const_iterator itHit) {
0408   bool result = false;
0409   int trackId = simTrack->trackId();
0410   int trackId_sim = itHit->trackId();
0411   if (trackId == trackId_sim)
0412     result = true;
0413   return result;
0414 }