MuonME0Validation/src/ME0RecHitsValidation.cc

0001 #include "Validation/MuonME0Validation/interface/ME0RecHitsValidation.h"
0002 #include <TMath.h>
0003
0004 ME0RecHitsValidation::ME0RecHitsValidation(const edm::ParameterSet &cfg) : ME0BaseValidation(cfg) {
0005   InputTagToken_ = consumes<edm::PSimHitContainer>(cfg.getParameter<edm::InputTag>("simInputLabel"));
0006   InputTagToken_RecHit = consumes<ME0RecHitCollection>(cfg.getParameter<edm::InputTag>("recHitInputLabel"));
0007 }
0008
0009 void ME0RecHitsValidation::bookHistograms(DQMStore::IBooker &ibooker,
0010                                           edm::Run const &Run,
0011                                           edm::EventSetup const &iSetup) {
0012   LogDebug("MuonME0RecHitsValidation") << "Info : Loading Geometry information\n";
0013   ibooker.setCurrentFolder("MuonME0RecHitsV/ME0RecHitsTask");
0014
0015   unsigned int nregion = 2;
0016
0017   edm::LogInfo("MuonME0RecHitsValidation") << "+++ Info : # of region : " << nregion << std::endl;
0018
0019   LogDebug("MuonME0RecHitsValidation") << "+++ Info : finish to get geometry information from ES.\n";
0020
0021   for (unsigned int region_num = 0; region_num < nregion; region_num++) {
0022     me0_rh_zr[region_num] = BookHistZR(ibooker, "me0_rh_tot", "Digi", region_num);
0023     for (unsigned int layer_num = 0; layer_num < 6; layer_num++) {
0024       me0_rh_xy[region_num][layer_num] = BookHistXY(ibooker, "me0_rh", "RecHit", region_num, layer_num);
0025
0026       std::string histo_name_DeltaX =
0027           std::string("me0_rh_DeltaX_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0028       std::string histo_name_DeltaY =
0029           std::string("me0_rh_DeltaY_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0030       std::string histo_label_DeltaX = "RecHit Delta X : region" + regionLabel[region_num] + " layer " +
0031                                        layerLabel[layer_num] + " " + " ; x_{SimHit} - x_{RecHit} ; entries";
0032       std::string histo_label_DeltaY = "RecHit Delta Y : region" + regionLabel[region_num] + " layer " +
0033                                        layerLabel[layer_num] + " " + " ; y_{SimHit} - y_{RecHit} ; entries";
0034
0035       me0_rh_DeltaX[region_num][layer_num] =
0036           ibooker.book1D(histo_name_DeltaX.c_str(), histo_label_DeltaX.c_str(), 100, -10, 10);
0037       me0_rh_DeltaY[region_num][layer_num] =
0038           ibooker.book1D(histo_name_DeltaY.c_str(), histo_label_DeltaY.c_str(), 100, -10, 10);
0039
0040       std::string histo_name_PullX =
0041           std::string("me0_rh_PullX_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0042       std::string histo_name_PullY =
0043           std::string("me0_rh_PullY_r") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0044       std::string histo_label_PullX = "RecHit Pull X : region" + regionLabel[region_num] + " layer " +
0045                                       layerLabel[layer_num] + " " +
0046                                       " ; #frac{x_{SimHit} - x_{RecHit}}{#sigma_{x,RecHit}} ; entries";
0047       std::string histo_label_PullY = "RecHit Pull Y : region" + regionLabel[region_num] + " layer " +
0048                                       layerLabel[layer_num] + " " +
0049                                       " ; #frac{y_{SimHit} - y_{RecHit}}{#sigma_{y,RecHit}} ; entries";
0050
0051       me0_rh_PullX[region_num][layer_num] =
0052           ibooker.book1D(histo_name_PullX.c_str(), histo_label_DeltaX.c_str(), 100, -10, 10);
0053       me0_rh_PullY[region_num][layer_num] =
0054           ibooker.book1D(histo_name_PullY.c_str(), histo_label_DeltaY.c_str(), 100, -10, 10);
0055     }
0056   }
0057 }
0058
0059 ME0RecHitsValidation::~ME0RecHitsValidation() {}
0060
0061 void ME0RecHitsValidation::analyze(const edm::Event &e, const edm::EventSetup &iSetup) {
0062   const ME0Geometry *ME0Geometry_ = &iSetup.getData(geomToken_);
0063
0064   edm::Handle<edm::PSimHitContainer> ME0Hits;
0065   e.getByToken(InputTagToken_, ME0Hits);
0066
0067   edm::Handle<ME0RecHitCollection> ME0RecHits;
0068   e.getByToken(InputTagToken_RecHit, ME0RecHits);
0069
0070   if (!ME0Hits.isValid() || !ME0RecHits.isValid()) {
0071     edm::LogError("ME0RecHitsValidation") << "Cannot get ME0Hits/ME0RecHits by Token simInputTagToken";
0072     return;
0073   }
0074
0075   for (auto hits = ME0Hits->begin(); hits != ME0Hits->end(); hits++) {
0076     const ME0DetId id(hits->detUnitId());
0077     Int_t sh_region = id.region();
0078     Int_t sh_layer = id.layer();
0079     Int_t sh_station = 0;
0080     Int_t sh_chamber = id.chamber();
0081     Int_t sh_roll = id.roll();
0082
0083     // Int_t even_odd = id.chamber()%2;
0084     if (ME0Geometry_->idToDet(hits->detUnitId()) == nullptr) {
0085       edm::LogInfo("ME0RecHitsValidation") << "simHit did not matched with GEMGeometry." << std::endl;
0086       continue;
0087     }
0088
0089     const LocalPoint hitLP(hits->localPosition());
0090
0091     for (ME0RecHitCollection::const_iterator recHit = ME0RecHits->begin(); recHit != ME0RecHits->end(); ++recHit) {
0092       Float_t x = recHit->localPosition().x();
0093       Float_t xErr = recHit->localPositionError().xx();
0094       Float_t y = recHit->localPosition().y();
0095       Float_t yErr = recHit->localPositionError().yy();
0096       // Float_t detId = (Short_t) (*recHit).me0Id();
0097       Float_t bx = 0;
0098       Float_t clusterSize = 0;
0099       Float_t firstClusterStrip = 0;
0100
0101       ME0DetId id((*recHit).me0Id());
0102
0103       Short_t region = (Short_t)id.region();
0104       Short_t station = 0;
0105       Short_t layer = (Short_t)id.layer();
0106       Short_t chamber = (Short_t)id.chamber();
0107       Short_t roll = (Short_t)id.roll();
0108
0109       LocalPoint rhLP = recHit->localPosition();
0110       GlobalPoint rhGP = ME0Geometry_->idToDet((*recHit).me0Id())->surface().toGlobal(rhLP);
0111
0112       Float_t globalR = rhGP.perp();
0113       Float_t globalX = rhGP.x();
0114       Float_t globalY = rhGP.y();
0115       Float_t globalZ = rhGP.z();
0116
0117       Float_t x_sim = hitLP.x();
0118       Float_t y_sim = hitLP.y();
0119       Float_t pullX = (x_sim - x) / sqrt(xErr);
0120       Float_t pullY = (y_sim - y) / sqrt(yErr);
0121
0122       if (bx != 0)
0123         continue;
0124
0125       std::vector<int> stripsFired;
0126       for (int i = firstClusterStrip; i < (firstClusterStrip + clusterSize); i++) {
0127         stripsFired.push_back(i);
0128       }
0129
0130       const bool cond1(sh_region == region and sh_layer == layer and sh_station == station);
0131       const bool cond2(sh_chamber == chamber and sh_roll == roll);
0132       const bool cond3(
0133           (std::find(stripsFired.begin(), stripsFired.end(), (firstClusterStrip + 1)) != stripsFired.end()) or
0134           clusterSize == 0);
0135
0136       int region_num = 0;
0137       if (region == -1)
0138         region_num = 0;
0139       else if (region == 1)
0140         region_num = 1;
0141       int layer_num = layer - 1;
0142
0143       if (cond1 and cond2 and cond3) {
0144         me0_rh_xy[region_num][layer_num]->Fill(globalX, globalY);
0145         me0_rh_zr[region_num]->Fill(globalZ, globalR);
0146
0147         me0_rh_DeltaX[region_num][layer_num]->Fill(x_sim - x);
0148         me0_rh_DeltaY[region_num][layer_num]->Fill(y_sim - y);
0149         me0_rh_PullX[region_num][layer_num]->Fill(pullX);
0150         me0_rh_PullY[region_num][layer_num]->Fill(pullY);
0151       }
0152     }
0153   }
0154 }