GlobalTrackingTools/test/DYTTuner.cc

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     // Energy range cut
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   // Full barrel/endcap computation
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   // Chamber by chamber computation
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   // Writing to DB
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   //The following says we do not know what parameters are allowed so do no validation
0226   // Please change this to state exactly what you do use, even if it is no parameters
0227   edm::ParameterSetDescription desc;
0228   desc.setUnknown();
0229   descriptions.addDefault(desc);
0230 }
0231
0232 DEFINE_FWK_MODULE(DYTTuner);