File indexing completed on 2023-03-17 11:20:26
0001 #include <memory>
0002 #include <algorithm>
0003 #include "FWCore/Utilities/interface/InputTag.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "CondFormats/RecoMuonObjects/interface/DYTThrObject.h"
0010 #include "CondFormats/DataRecord/interface/DYTThrObjectRcd.h"
0011 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0014 #include "DataFormats/MuonReco/interface/DYTInfo.h"
0015 #include "DataFormats/DetId/interface/DetId.h"
0016 #include "DataFormats/Common/interface/ValueMap.h"
0017 #include "DataFormats/Common/interface/Handle.h"
0018 #include "DataFormats/MuonReco/interface/Muon.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0023 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0024 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0025 #include "FWCore/Framework/interface/ConsumesCollector.h"
0026
0027 #define MAX_THR 1e7
0028
0029 class DYTTuner : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::WatchLuminosityBlocks> {
0030 public:
0031 explicit DYTTuner(const edm::ParameterSet&);
0032 ~DYTTuner() override;
0033
0034 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0035
0036 private:
0037 void beginJob() override;
0038 void analyze(const edm::Event&, const edm::EventSetup&) override;
0039 void endJob() override;
0040 virtual double doIntegral(std::vector<double>&, DetId&);
0041 virtual void writePlots();
0042 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0043 void endRun(edm::Run const&, edm::EventSetup const&) override;
0044 void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0045 void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0046
0047 typedef edm::ValueMap<reco::DYTInfo> DYTestimators;
0048 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0049 double MinEnVal, MaxEnVal, IntegralCut, MaxEstVal;
0050 unsigned int MinNumValues, MaxValPlots, NBinsPlots;
0051 bool saveROOTfile;
0052 std::map<DetId, std::vector<double> > mapId;
0053 std::map<DetId, TH1F*> EstPlots;
0054 edm::EDGetTokenT<DYTestimators> dytInfoToken;
0055 edm::EDGetTokenT<reco::MuonCollection> muonsToken;
0056 };
0057
0058 DYTTuner::DYTTuner(const edm::ParameterSet& iConfig) {
0059 MinEnVal = iConfig.getParameter<double>("MinEnergyVal");
0060 MaxEnVal = iConfig.getParameter<double>("MaxEnergyVal");
0061 IntegralCut = iConfig.getParameter<double>("IntegralCut");
0062 MaxEstVal = iConfig.getParameter<double>("MaxEstVal");
0063 MinNumValues = iConfig.getParameter<unsigned int>("MinNumValues");
0064 saveROOTfile = iConfig.getParameter<bool>("writePlots");
0065 MaxValPlots = iConfig.getParameter<unsigned int>("MaxValPlots");
0066 NBinsPlots = iConfig.getParameter<unsigned int>("NBinsPlots");
0067
0068 edm::ConsumesCollector iC = consumesCollector();
0069 dytInfoToken = iC.consumes<DYTestimators>(edm::InputTag("tevMuons", "dytInfo"));
0070 muonsToken = iC.consumes<reco::MuonCollection>(edm::InputTag("muons"));
0071
0072 if (MaxEstVal == -1)
0073 MaxEstVal = MAX_THR;
0074
0075 if (!poolDbService->isNewTagRequest("DYTThrObjectRcd"))
0076 throw cms::Exception("NotAvailable") << "The output file already contains a valid \"DYTThrObjectRcd\" "
0077 "record.\nPlease provide a different file name or tag.";
0078 }
0079
0080 DYTTuner::~DYTTuner() {}
0081
0082 void DYTTuner::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0083 using namespace std;
0084 using namespace edm;
0085 using namespace reco;
0086
0087 Handle<DYTestimators> dytInfoH;
0088 iEvent.getByToken(dytInfoToken, dytInfoH);
0089 const DYTestimators& dytInfoC = *dytInfoH;
0090 Handle<MuonCollection> muons;
0091 iEvent.getByToken(muonsToken, muons);
0092 for (size_t i = 0; i != muons->size(); ++i) {
0093
0094 if (muons->at(i).pt() < MinEnVal || muons->at(i).pt() > MaxEnVal)
0095 continue;
0096
0097 const TrackRef& tkRef = muons->at(i).globalTrack();
0098 if (dytInfoC.contains(tkRef.id())) {
0099 DYTInfo dytInfo = dytInfoC[muons->at(i).globalTrack()];
0100 vector<double> estimators = dytInfo.DYTEstimators();
0101 vector<DetId> ids = dytInfo.IdChambers();
0102
0103 for (int j = 0; j < 4; j++) {
0104 if (ids[j].null())
0105 continue;
0106 DetId chamberId = ids[j];
0107 double estValue = estimators[j];
0108 if (estValue >= 0 && estValue <= MaxEstVal)
0109 mapId[chamberId].push_back(estValue);
0110 }
0111 } else {
0112 continue;
0113 }
0114 }
0115 }
0116
0117 void DYTTuner::beginJob() {}
0118
0119 void DYTTuner::endJob() {
0120 if (saveROOTfile)
0121 writePlots();
0122 DYTThrObject thresholds;
0123
0124
0125 std::map<DetId, std::vector<double> >::iterator it;
0126 std::map<int, std::vector<double> > estBarrel, estEndcap;
0127 for (it = mapId.begin(); it != mapId.end(); it++) {
0128 DetId id = (*it).first;
0129 std::vector<double> estValCh = (*it).second;
0130 if ((*it).first.subdetId() == MuonSubdetId::DT)
0131 for (unsigned int b = 0; b < estValCh.size(); b++) {
0132 int station = DTChamberId(id).station();
0133 estBarrel[station].push_back(estValCh[b]);
0134 }
0135 if ((*it).first.subdetId() == MuonSubdetId::CSC)
0136 for (unsigned int e = 0; e < estValCh.size(); e++) {
0137 int station = CSCDetId(id).station();
0138 estEndcap[station].push_back(estValCh[e]);
0139 }
0140 }
0141 DetId empty = DetId();
0142 double barrelCut[4], endcapCut[4];
0143 for (unsigned st = 1; st <= 4; st++) {
0144 barrelCut[st - 1] = doIntegral(estBarrel[st], empty);
0145 endcapCut[st - 1] = doIntegral(estEndcap[st], empty);
0146 }
0147
0148
0149 for (it = mapId.begin(); it != mapId.end(); it++) {
0150 DetId id = (*it).first;
0151 std::vector<double> estValCh = (*it).second;
0152 DYTThrObject::DytThrStruct obj;
0153 obj.id = id;
0154 if (estValCh.size() < MinNumValues) {
0155 if (id.subdetId() == MuonSubdetId::DT) {
0156 int station = DTChamberId(id).station();
0157 obj.thr = barrelCut[station - 1];
0158 }
0159 if (id.subdetId() == MuonSubdetId::CSC) {
0160 int station = CSCDetId(id).station();
0161 obj.thr = endcapCut[station - 1];
0162 }
0163 thresholds.thrsVec.push_back(obj);
0164 continue;
0165 }
0166 obj.thr = doIntegral(estValCh, id);
0167 thresholds.thrsVec.push_back(obj);
0168 }
0169
0170
0171 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0172 if (poolDbService.isAvailable()) {
0173 poolDbService->writeOneIOV(thresholds, poolDbService->beginOfTime(), "DYTThrObjectRcd");
0174 } else
0175 throw cms::Exception("NotAvailable") << "PoolDBOutputService is not available.";
0176 }
0177
0178 double DYTTuner::doIntegral(std::vector<double>& estValues, DetId& id) {
0179 double cutOnE = -1;
0180 int nPosVal = 0;
0181 sort(estValues.begin(), estValues.end());
0182 for (unsigned int j = 0; j < estValues.size(); j++)
0183 if (estValues[j] > 0 && estValues[j] < MaxEstVal)
0184 nPosVal++;
0185 double limit = nPosVal * IntegralCut;
0186 int nVal = 0;
0187 for (unsigned int j = 0; j < estValues.size(); j++) {
0188 if (estValues[j] < 0)
0189 continue;
0190 nVal++;
0191 if (nVal >= limit) {
0192 cutOnE = estValues[j - 1];
0193 break;
0194 }
0195 }
0196 std::cout << "Det Id: " << id.rawId() << " - Threshold:: " << cutOnE << std::endl;
0197 return cutOnE;
0198 }
0199
0200 void DYTTuner::writePlots() {
0201 edm::Service<TFileService> fs;
0202 std::map<DetId, std::vector<double> >::iterator it;
0203 for (it = mapId.begin(); it != mapId.end(); it++) {
0204 DetId id = (*it).first;
0205 std::vector<double> estValCh = (*it).second;
0206 char plotName[200];
0207 sprintf(plotName, "%i", id.rawId());
0208 TH1F* tmpPlot = new TH1F(plotName, plotName, NBinsPlots, 0., MaxValPlots);
0209 for (unsigned int i = 0; i < estValCh.size(); i++)
0210 tmpPlot->Fill(estValCh[i]);
0211 EstPlots[id] = fs->make<TH1F>(*tmpPlot);
0212 delete tmpPlot;
0213 }
0214 }
0215
0216 void DYTTuner::beginRun(edm::Run const&, edm::EventSetup const&) {}
0217
0218 void DYTTuner::endRun(edm::Run const&, edm::EventSetup const&) {}
0219
0220 void DYTTuner::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
0221
0222 void DYTTuner::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
0223
0224 void DYTTuner::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0225
0226
0227 edm::ParameterSetDescription desc;
0228 desc.setUnknown();
0229 descriptions.addDefault(desc);
0230 }
0231
0232 DEFINE_FWK_MODULE(DYTTuner);