Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-08-29 02:16:44

0001 #include "DQM/GEM/interface/GEMDQMBase.h"
0002 #include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
0003 
0004 using namespace std;
0005 using namespace edm;
0006 
0007 GEMDQMBase::GEMDQMBase(const edm::ParameterSet& cfg) : geomToken_(esConsumes<edm::Transition::BeginRun>()) {
0008   std::string strRunType = cfg.getUntrackedParameter<std::string>("runType");
0009 
0010   nRunType_ = GEMDQM_RUNTYPE_ONLINE;
0011 
0012   if (strRunType == "online") {
0013     nRunType_ = GEMDQM_RUNTYPE_ONLINE;
0014   } else if (strRunType == "offline") {
0015     nRunType_ = GEMDQM_RUNTYPE_OFFLINE;
0016   } else if (strRunType == "relval") {
0017     nRunType_ = GEMDQM_RUNTYPE_RELVAL;
0018   } else if (strRunType == "allplots") {
0019     nRunType_ = GEMDQM_RUNTYPE_ALLPLOTS;
0020   } else {
0021     edm::LogError(log_category_) << "+++ Error : GEM geometry is unavailable on event loop. +++\n";
0022   }
0023 
0024   log_category_ = cfg.getUntrackedParameter<std::string>("logCategory");
0025 }
0026 
0027 int GEMDQMBase::initGeometry(edm::EventSetup const& iSetup) {
0028   GEMGeometry_ = nullptr;
0029   if (auto handle = iSetup.getHandle(geomToken_)) {
0030     GEMGeometry_ = handle.product();
0031   } else {
0032     edm::LogError(log_category_) << "+++ Error : GEM geometry is unavailable on event loop. +++\n";
0033     return -1;
0034   }
0035 
0036   return 0;
0037 }
0038 
0039 // Borrowed from DQM/GEM/src/GEMOfflineDQMBase.cc
0040 int GEMDQMBase::getNumEtaPartitions(const GEMStation* station) {
0041   const auto&& superchambers = station->superChambers();
0042   if (not checkRefs(superchambers)) {
0043     edm::LogError(log_category_) << "failed to get a valid vector of GEMSuperChamber ptrs" << std::endl;
0044     return 0;
0045   }
0046 
0047   const auto& chambers = superchambers.front()->chambers();
0048   if (not checkRefs(chambers)) {
0049     edm::LogError(log_category_) << "failed to get a valid vector of GEMChamber ptrs" << std::endl;
0050     return 0;
0051   }
0052 
0053   return chambers.front()->nEtaPartitions();
0054 }
0055 
0056 int GEMDQMBase::loadChambers() {
0057   if (GEMGeometry_ == nullptr)
0058     return -1;
0059   listChamberId_.clear();
0060   mapEtaPartition_.clear();
0061   for (const GEMRegion* region : GEMGeometry_->regions()) {
0062     for (const GEMStation* station : region->stations()) {
0063       for (auto sch : station->superChambers()) {
0064         for (auto pchamber : sch->chambers()) {
0065           GEMDetId gid = pchamber->id();
0066           listChamberId_.push_back(pchamber->id());
0067           for (auto iEta : pchamber->etaPartitions()) {
0068             mapEtaPartition_[gid].push_back(iEta);
0069           }
0070         }
0071       }
0072     }
0073   }
0074 
0075   // Borrwed from DQM/GEM/src/GEMOfflineMonitor.cc
0076   nMaxNumCh_ = 0;
0077   for (const GEMRegion* region : GEMGeometry_->regions()) {
0078     const int region_number = region->region();
0079 
0080     for (const GEMStation* station : region->stations()) {
0081       const auto&& superchambers = station->superChambers();
0082 
0083       const int station_number = station->station();
0084       const int num_superchambers = (station_number == 1 ? 36 : 18);
0085       const int num_mod = getNumModule(station->station());
0086       const int max_vfat = getMaxVFAT(station->station());  // the number of VFATs per GEMEtaPartition
0087       const int num_etas = getNumEtaPartitions(station);    // the number of eta partitions per GEMChamber
0088       const int num_vfat = num_etas * max_vfat;             // the number of VFATs per GEMChamber
0089       const int strip1st = (station_number == 2 ? 1 : 0);   // the index of the first strip
0090       const int num_digi = GEMeMap::maxChan_;               // the number of digis (channels) per VFAT
0091 
0092       nMaxNumCh_ = std::max(nMaxNumCh_, num_superchambers);
0093 
0094       Int_t nMinIdxChamber = 1048576;
0095       Int_t nMaxIdxChamber = -1048576;
0096       for (auto sch : superchambers) {
0097         auto nIdxChamber = sch->chambers().front()->id().chamber();
0098         if (nMinIdxChamber > nIdxChamber)
0099           nMinIdxChamber = nIdxChamber;
0100         if (nMaxIdxChamber < nIdxChamber)
0101           nMaxIdxChamber = nIdxChamber;
0102       }
0103 
0104       const auto& chambers = superchambers.front()->chambers();
0105 
0106       for (auto pchamber : chambers) {
0107         int layer_number = pchamber->id().layer();
0108         ME3IdsKey key3(region_number, station_number, layer_number);
0109         mapStationInfo_[key3] = MEStationInfo(region_number,
0110                                               station_number,
0111                                               layer_number,
0112                                               num_superchambers,
0113                                               num_mod,
0114                                               num_etas,
0115                                               num_vfat,
0116                                               strip1st,
0117                                               num_digi,
0118                                               nMinIdxChamber,
0119                                               nMaxIdxChamber);
0120         readGeometryRadiusInfoChamber(station, mapStationInfo_[key3]);
0121         readGeometryPhiInfoChamber(station, mapStationInfo_[key3]);
0122       }
0123     }
0124   }
0125 
0126   return 0;
0127 }
0128 
0129 int GEMDQMBase::SortingLayers(std::vector<ME4IdsKey>& listLayers) {
0130   std::sort(listLayers.begin(), listLayers.end(), [](ME4IdsKey key1, ME4IdsKey key2) {
0131     Int_t re1 = std::get<0>(key1), re2 = std::get<0>(key2);
0132     Int_t st1 = std::get<1>(key1), st2 = std::get<1>(key2);
0133     Int_t la1 = std::get<2>(key1), la2 = std::get<2>(key2);
0134     Int_t mo1 = std::get<3>(key1), mo2 = std::get<3>(key2);
0135     if (re1 < 0 && re2 > 0)
0136       return false;
0137     if (re1 > 0 && re2 < 0)
0138       return true;
0139     Bool_t bRes = (re1 < 0);  // == re2 < 0
0140     Int_t sum1 = 4096 * std::abs(re1) + 256 * st1 + 16 * la1 + mo1;
0141     Int_t sum2 = 4096 * std::abs(re2) + 256 * st2 + 16 * la2 + mo2;
0142     if (sum1 <= sum2)
0143       return bRes;
0144     return !bRes;
0145   });
0146 
0147   return 0;
0148 }
0149 
0150 dqm::impl::MonitorElement* GEMDQMBase::CreateSummaryHist(DQMStore::IBooker& ibooker, TString strName) {
0151   std::vector<ME4IdsKey> listLayers;
0152   for (auto const& [key3, stationInfo] : mapStationInfo_) {
0153     for (int module_number = 1; module_number <= stationInfo.nNumModules_; module_number++) {
0154       ME4IdsKey key4{keyToRegion(key3), keyToStation(key3), keyToLayer(key3), module_number};
0155       listLayers.push_back(key4);  // Note: Not only count layers but also modules
0156     }
0157   }
0158   SortingLayers(listLayers);
0159   for (Int_t i = 0; i < (Int_t)listLayers.size(); i++)
0160     mapStationToIdx_[listLayers[i]] = i + 1;
0161 
0162   auto h2Res =
0163       ibooker.book2D(strName, "", nMaxNumCh_, 0.5, nMaxNumCh_ + 0.5, listLayers.size(), 0.5, listLayers.size() + 0.5);
0164   h2Res->setXTitle("Chamber");
0165   h2Res->setYTitle("Layer");
0166 
0167   if (h2Res == nullptr)
0168     return nullptr;
0169 
0170   for (Int_t i = 1; i <= nMaxNumCh_; i++)
0171     h2Res->setBinLabel(i, Form("%i", i), 1);
0172   for (Int_t i = 1; i <= (Int_t)listLayers.size(); i++) {
0173     auto key = listLayers[i - 1];
0174     ME3IdsKey key3 = key4Tokey3(key);
0175 
0176     auto region = keyToRegion(key);
0177     auto strInfo = GEMUtils::getSuffixName(key3);  // NOTE: It starts with '_'
0178     if (mapStationInfo_[key3].nNumModules_ > 1) {
0179       strInfo += Form("-M%i", keyToModule(key));
0180     }
0181     auto label = Form("GE%+i1-%cL%i-M%i;%s",
0182                       region * keyToStation(key),
0183                       (region > 0 ? 'P' : 'M'),
0184                       keyToLayer(key),
0185                       keyToModule(key),
0186                       strInfo.Data());
0187     h2Res->setBinLabel(i, label, 2);
0188     Int_t nNumCh = mapStationInfo_[key3].nNumChambers_;
0189     h2Res->setBinContent(0, i, nNumCh);
0190   }
0191 
0192   h2Res->setBinContent(0, 0, 1.0);
0193 
0194   return h2Res;
0195 }
0196 
0197 int GEMDQMBase::GenerateMEPerChamber(DQMStore::IBooker& ibooker) {
0198   MEMap2Check_.clear();
0199   MEMap3Check_.clear();
0200   MEMap4Check_.clear();
0201   MEMap5Check_.clear();
0202   MEMap2WithEtaCheck_.clear();
0203   MEMap2AbsReWithEtaCheck_.clear();
0204 
0205   MEMap2WithChCheck_.clear();
0206 
0207   MEMap4WithChCheck_.clear();
0208   MEMap5WithChCheck_.clear();
0209 
0210   MEMap2WithEtaChCheck_.clear();
0211   for (auto gid : listChamberId_) {
0212     ME2IdsKey key2{gid.region(), gid.station()};
0213     ME3IdsKey key3{gid.region(), gid.station(), gid.layer()};
0214     /*******************/
0215     ME3IdsKey key2WithChamber{gid.region(), gid.station(), gid.chamber()};
0216     /******************/
0217     const auto num_mod = mapStationInfo_[key3].nNumModules_;
0218     for (int module_number = 1; module_number <= num_mod; module_number++) {
0219       ME4IdsKey key4{gid.region(), gid.station(), gid.layer(), module_number};
0220       ME4IdsKey key4WithChamber{gid.region(), gid.station(), gid.layer(), gid.chamber()};
0221       ME5IdsKey key5WithChamber{gid.region(), gid.station(), gid.layer(), module_number, gid.chamber()};
0222       if (!MEMap2Check_[key2]) {
0223         auto strSuffixName = GEMUtils::getSuffixName(key2);
0224         auto strSuffixTitle = GEMUtils::getSuffixTitle(key2);
0225         BookingHelper bh2(ibooker, strSuffixName, strSuffixTitle);
0226         ProcessWithMEMap2(bh2, key2);
0227         MEMap2Check_[key2] = true;
0228       }
0229       if (!MEMap2WithChCheck_[key2WithChamber]) {
0230         Int_t nCh = gid.chamber();
0231         //Int_t nLa = gid.layer();
0232         char cLS = (nCh % 2 == 0 ? 'L' : 'S');  // FIXME: Is it general enough?
0233         auto strSuffixName = GEMUtils::getSuffixName(key2) + Form("-%02i-%c", nCh, cLS);
0234         auto strSuffixTitle = GEMUtils::getSuffixTitle(key2) + Form("-%02i-%c", nCh, cLS);
0235         BookingHelper bh2Ch(ibooker, strSuffixName, strSuffixTitle);
0236         ProcessWithMEMap2WithChamber(bh2Ch, key2WithChamber);
0237         MEMap2WithChCheck_[key2WithChamber] = true;
0238       }
0239       if (!MEMap3Check_[key3]) {
0240         auto strSuffixName = GEMUtils::getSuffixName(key3);
0241         auto strSuffixTitle = GEMUtils::getSuffixTitle(key3);
0242         BookingHelper bh3(ibooker, strSuffixName, strSuffixTitle);
0243         ProcessWithMEMap3(bh3, key3);
0244         MEMap3Check_[key3] = true;
0245       }
0246       if (!MEMap4Check_[key4]) {
0247         Int_t nLa = gid.layer();
0248         TString strSuffixCh = Form("-L%i", nLa);
0249         if (mapStationInfo_[key3].nNumModules_ > 1)
0250           strSuffixCh = Form("-L%i-M%i", nLa, module_number);
0251         auto strSuffixName = GEMUtils::getSuffixName(key2) + strSuffixCh;
0252         auto strSuffixTitle = GEMUtils::getSuffixTitle(key2) + strSuffixCh;
0253         BookingHelper bh4(ibooker, strSuffixName, strSuffixTitle);
0254         ProcessWithMEMap4(bh4, key4);
0255         MEMap4Check_[key4] = true;
0256       }
0257       if (!MEMap4WithChCheck_[key4WithChamber]) {
0258         Int_t nCh = gid.chamber();
0259         Int_t nLa = gid.layer();
0260         char cLS = (nCh % 2 == 0 ? 'L' : 'S');  // FIXME: Is it general enough?
0261         TString strSuffixCh = Form("-%02iL%i-%c", nCh, nLa, cLS);
0262         auto strSuffixName = GEMUtils::getSuffixName(key2) + strSuffixCh;
0263         auto strSuffixTitle = GEMUtils::getSuffixTitle(key2) + strSuffixCh;
0264         BookingHelper bh4Ch(ibooker, strSuffixName, strSuffixTitle);
0265         ProcessWithMEMap4WithChamber(bh4Ch, key4WithChamber);
0266         MEMap4WithChCheck_[key4WithChamber] = true;
0267       }
0268       if (!MEMap5WithChCheck_[key5WithChamber]) {
0269         Int_t nCh = gid.chamber();
0270         Int_t nLa = gid.layer();
0271         char cLS = (nCh % 2 == 0 ? 'L' : 'S');  // FIXME: Is it general enough?
0272         TString strSuffixCh = Form("-%02iL%i-%c", nCh, nLa, cLS);
0273         if (mapStationInfo_[key3].nNumModules_ > 1)
0274           strSuffixCh = Form("-%02iL%i-M%i-%c", nCh, nLa, module_number, cLS);
0275         auto strSuffixName = GEMUtils::getSuffixName(key2) + strSuffixCh;
0276         auto strSuffixTitle = GEMUtils::getSuffixTitle(key2) + strSuffixCh;
0277         BookingHelper bh5Ch(ibooker, strSuffixName, strSuffixTitle);
0278         ProcessWithMEMap5WithChamber(bh5Ch, key5WithChamber);
0279         MEMap5WithChCheck_[key5WithChamber] = true;
0280       }
0281       for (auto iEta : mapEtaPartition_[gid]) {
0282         GEMDetId eId = iEta->id();
0283         ME5IdsKey key5{gid.region(), gid.station(), gid.layer(), module_number, eId.ieta()};
0284         /*******************/
0285         ME4IdsKey key2WithEtaCh{gid.region(), gid.station(), eId.ieta(), gid.chamber()};
0286         /******************/
0287         ME3IdsKey key2WithEta{gid.region(), gid.station(), eId.ieta()};
0288         ME3IdsKey key2AbsReWithEta{std::abs(gid.region()), gid.station(), eId.ieta()};
0289         if (!MEMap5Check_[key5]) {
0290           auto strSuffixName = GEMUtils::getSuffixName(key3) + Form("-E%02i", eId.ieta());
0291           auto strSuffixTitle = GEMUtils::getSuffixTitle(key3) + Form("-E%02i", eId.ieta());
0292           BookingHelper bh5(ibooker, strSuffixName, strSuffixTitle);
0293           ProcessWithMEMap5(bh5, key5);
0294           MEMap5Check_[key5] = true;
0295         }
0296         if (!MEMap2WithEtaChCheck_[key2WithEtaCh]) {
0297           Int_t nCh = gid.chamber();
0298           //Int_t nLa = gid.layer();
0299           char cLS = (nCh % 2 == 0 ? 'L' : 'S');
0300           auto strSuffixName = GEMUtils::getSuffixName(key2) + Form("-%02i-%c-E%02i", nCh, cLS, eId.ieta());
0301           auto strSuffixTitle = GEMUtils::getSuffixTitle(key2) + Form("-%02i-%c-E%02i", nCh, cLS, eId.ieta());
0302           BookingHelper bh4(ibooker, strSuffixName, strSuffixTitle);
0303           ProcessWithMEMap2WithEtaCh(bh4, key2WithEtaCh);
0304           MEMap2WithEtaChCheck_[key2WithEtaCh] = true;
0305         }
0306         if (!MEMap2WithEtaCheck_[key2WithEta]) {
0307           auto strSuffixName = GEMUtils::getSuffixName(key2) + Form("-E%02i", eId.ieta());
0308           auto strSuffixTitle = GEMUtils::getSuffixTitle(key2) + Form("-E%02i", eId.ieta());
0309           BookingHelper bh3(ibooker, strSuffixName, strSuffixTitle);
0310           ProcessWithMEMap2WithEta(bh3, key2WithEta);
0311           MEMap2WithEtaCheck_[key2WithEta] = true;
0312         }
0313         if (!MEMap2AbsReWithEtaCheck_[key2AbsReWithEta]) {
0314           auto strSuffixName = Form("_GE%d1-E%02i", gid.station(), eId.ieta());
0315           auto strSuffixTitle = Form(" GE%d1-E%02i", gid.station(), eId.ieta());
0316           BookingHelper bh3(ibooker, strSuffixName, strSuffixTitle);
0317           ProcessWithMEMap2AbsReWithEta(bh3, key2AbsReWithEta);
0318           MEMap2AbsReWithEtaCheck_[key2AbsReWithEta] = true;
0319         }
0320       }
0321     }
0322   }
0323   return 0;
0324 }
0325 
0326 int GEMDQMBase::readGeometryRadiusInfoChamber(const GEMStation* station, MEStationInfo& stationInfo) {
0327   auto listSuperChambers = station->superChambers();
0328 
0329   Bool_t bDoneEven = false, bDoneOdd = false;
0330 
0331   // Obtaining radius intervals of even/odd chambers
0332   for (auto superchamber : listSuperChambers) {
0333     Int_t chamberNo = superchamber->id().chamber();
0334     if (chamberNo % 2 == 0 && bDoneEven)
0335       continue;
0336     if (chamberNo % 2 != 0 && bDoneOdd)
0337       continue;
0338 
0339     auto& etaPartitions = superchamber->chambers().front()->etaPartitions();
0340 
0341     // A little of additional procedures to list up the radius intervals
0342     // It would be independent to local direction of chambers and the order of eta partitions
0343     //   1. Obtain the radius of the middle top/bottom points of the trapezoid
0344     //   2. Sort these two values and determine which one is the lower/upper one
0345     //   3. Keep them all and then sort them
0346     //   4. The intermediate radii are set as the mean of the corresponding values of upper/lowers.
0347     std::vector<Float_t> listRadiusLower, listRadiusUpper;
0348     for (auto iEta : etaPartitions) {
0349       const GEMStripTopology& stripTopology = dynamic_cast<const GEMStripTopology&>(iEta->specificTopology());
0350       Float_t fHeight = stripTopology.stripLength();
0351       LocalPoint lp1(0.0, -0.5 * fHeight), lp2(0.0, 0.5 * fHeight);
0352       auto& surface = iEta->surface();
0353       GlobalPoint gp1 = surface.toGlobal(lp1), gp2 = surface.toGlobal(lp2);
0354       Float_t fR1 = gp1.perp(), fR2 = gp2.perp();
0355       Float_t fRL = std::min(fR1, fR2), fRH = std::max(fR1, fR2);
0356       listRadiusLower.push_back(fRL);
0357       listRadiusUpper.push_back(fRH);
0358       // For a future usage
0359       //std::cout << "GEO_RADIUS: " << iEta->id().chamber() << " " << iEta->id().ieta() << " "
0360       //  << fRL << " " << fRH << std::endl;
0361     }
0362 
0363     std::sort(listRadiusLower.begin(), listRadiusLower.end());
0364     std::sort(listRadiusUpper.begin(), listRadiusUpper.end());
0365 
0366     std::vector<Float_t>& listR =
0367         (chamberNo % 2 == 0 ? stationInfo.listRadiusEvenChamber_ : stationInfo.listRadiusOddChamber_);
0368     listR.clear();
0369     listR.push_back(listRadiusLower.front());
0370     for (int i = 1; i < (int)listRadiusLower.size(); i++) {
0371       listR.push_back(0.5 * (listRadiusLower[i] + listRadiusUpper[i - 1]));
0372     }
0373     listR.push_back(listRadiusUpper.back());
0374 
0375     if (chamberNo % 2 == 0)
0376       bDoneEven = true;
0377     if (chamberNo % 2 != 0)
0378       bDoneOdd = true;
0379 
0380     if (bDoneEven && bDoneOdd)
0381       break;
0382   }
0383 
0384   return 0;
0385 }
0386 
0387 int GEMDQMBase::readGeometryPhiInfoChamber(const GEMStation* station, MEStationInfo& stationInfo) {
0388   auto listSuperChambers = station->superChambers();
0389   Int_t nNumStripEta = stationInfo.nNumDigi_ * (stationInfo.nMaxVFAT_ / stationInfo.nNumEtaPartitions_);
0390 
0391   std::vector<std::pair<Int_t, std::pair<std::pair<Float_t, Float_t>, Bool_t>>> listDivPhi;
0392 
0393   // Obtaining phi intervals of chambers
0394   for (auto superchamber : listSuperChambers) {
0395     auto iEta = superchamber->chambers().front()->etaPartitions().front();
0396 
0397     // What is the index of the first strip? Rather than to ask to someone, let's calculate it!
0398     Float_t fWidthStrip = std::abs(iEta->centreOfStrip((Int_t)1).x() - iEta->centreOfStrip((Int_t)0).x());
0399     LocalPoint lpRef(-fWidthStrip / 3.0, 0.0);
0400     Int_t nStripMid = (Int_t)iEta->strip(lpRef);
0401     Int_t nFirstStrip = 1 - ((nNumStripEta / 2) - nStripMid);
0402     Int_t nLastStrip = nFirstStrip + nNumStripEta - 1;
0403 
0404     auto& surface = iEta->surface();
0405     LocalPoint lpF = iEta->centreOfStrip((Float_t)(nFirstStrip - 0.5));  // To avoid the round error(?)
0406     LocalPoint lpL = iEta->centreOfStrip((Float_t)(nLastStrip + 0.5));   // To avoid the round error(?)
0407     GlobalPoint gpF = surface.toGlobal(lpF);
0408     GlobalPoint gpL = surface.toGlobal(lpL);
0409 
0410     Float_t fPhiF = gpF.phi();
0411     Float_t fPhiL = gpL.phi();
0412     if (fPhiF * fPhiL < 0 && std::abs(fPhiF) > 0.5 * 3.14159265359) {
0413       if (fPhiF < 0)
0414         fPhiF += 2 * 3.14159265359;
0415       if (fPhiL < 0)
0416         fPhiL += 2 * 3.14159265359;
0417     }
0418     Bool_t bFlipped = fPhiF > fPhiL;
0419     Float_t fPhiMin = std::min(fPhiF, fPhiL);
0420     Float_t fPhiMax = std::max(fPhiF, fPhiL);
0421 
0422     listDivPhi.emplace_back();
0423     listDivPhi.back().first = iEta->id().chamber();
0424     listDivPhi.back().second.first.first = fPhiMin;
0425     listDivPhi.back().second.first.second = fPhiMax;
0426     listDivPhi.back().second.second = bFlipped;
0427   }
0428 
0429   stationInfo.fMinPhi_ = 0.0;
0430   for (auto p : listDivPhi) {
0431     if (p.first == 1) {
0432       stationInfo.fMinPhi_ = p.second.first.first;
0433       break;
0434     }
0435   }
0436 
0437   // For a future usage
0438   //for ( auto p : listDivPhi ) {
0439   //  std::cout << "GEO_PHI: " << p.first << " "
0440   //    << p.second.first.first << " " << p.second.first.second << " " << p.second.second << std::endl;
0441   //}
0442 
0443   return 0;
0444 }