L1TMuonOverlap/plugins/OMTFPatternMaker.cc

0001 #include <iostream>
0002
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Utilities/interface/Transition.h"
0006
0007 #include "L1Trigger/L1TMuonOverlap/plugins/OMTFPatternMaker.h"
0008 #include "L1Trigger/L1TMuonOverlap/interface/OMTFProcessor.h"
0009 #include "L1Trigger/L1TMuonOverlap/interface/OMTFinputMaker.h"
0010 #include "L1Trigger/L1TMuonOverlap/interface/OMTFinput.h"
0011 #include "L1Trigger/L1TMuonOverlap/interface/OMTFConfiguration.h"
0012 #include "L1Trigger/L1TMuonOverlap/interface/OMTFConfigMaker.h"
0013 #include "L1Trigger/L1TMuonOverlap/interface/XMLConfigWriter.h"
0014
0015 #include "SimDataFormats/Track/interface/SimTrack.h"
0016
0017 #include "Math/VectorUtil.h"
0018
0019 #include "L1Trigger/RPCTrigger/interface/RPCConst.h"
0020
0021 OMTFPatternMaker::OMTFPatternMaker(const edm::ParameterSet& cfg)
0022     : theConfig(cfg),
0023       g4SimTrackSrc(cfg.getParameter<edm::InputTag>("g4SimTrackSrc")),
0024       esTokenParams_(esConsumes<edm::Transition::BeginRun>()) {
0025   inputTokenDTPh = consumes<L1MuDTChambPhContainer>(theConfig.getParameter<edm::InputTag>("srcDTPh"));
0026   inputTokenDTTh = consumes<L1MuDTChambThContainer>(theConfig.getParameter<edm::InputTag>("srcDTTh"));
0027   inputTokenCSC = consumes<CSCCorrelatedLCTDigiCollection>(theConfig.getParameter<edm::InputTag>("srcCSC"));
0028   inputTokenRPC = consumes<RPCDigiCollection>(theConfig.getParameter<edm::InputTag>("srcRPC"));
0029   inputTokenSimHit = consumes<edm::SimTrackContainer>(theConfig.getParameter<edm::InputTag>("g4SimTrackSrc"));
0030
0031   edm::ConsumesCollector consumesColl(consumesCollector());
0032   myInputMaker = new OMTFinputMaker(consumesColl);
0033
0034   makeGoldenPatterns = theConfig.getParameter<bool>("makeGoldenPatterns");
0035   makeConnectionsMaps = theConfig.getParameter<bool>("makeConnectionsMaps");
0036   mergeXMLFiles = theConfig.getParameter<bool>("mergeXMLFiles");
0037
0038   myOMTFConfig = nullptr;
0039 }
0040 /////////////////////////////////////////////////////
0041 /////////////////////////////////////////////////////
0042 OMTFPatternMaker::~OMTFPatternMaker() {
0043   delete myOMTFConfig;
0044   delete myOMTFConfigMaker;
0045   delete myOMTF;
0046 }
0047 /////////////////////////////////////////////////////
0048 /////////////////////////////////////////////////////
0049 void OMTFPatternMaker::beginRun(edm::Run const& run, edm::EventSetup const& iSetup) {
0050   const L1TMuonOverlapParams* omtfParams = &iSetup.getData(esTokenParams_);
0051
0052   if (!omtfParams) {
0053     edm::LogError("L1TMuonOverlapTrackProducer") << "Could not retrieve parameters from Event Setup" << std::endl;
0054   }
0055
0056   ///Initialise XML writer with default pdf.
0057   myWriter = new XMLConfigWriter(myOMTFConfig);
0058
0059   ///For making the patterns use extended pdf width in phi, as pdf are later shifted by the mean value
0060   ///For low pt muons non shifted pdfs would go out of the default pdf range.
0061   L1TMuonOverlapParams omtfParamsMutable = *omtfParams;
0062   std::vector<int> generalParams = *omtfParamsMutable.generalParams();
0063   nPdfAddrBits = omtfParams->nPdfAddrBits();
0064
0065   if (!mergeXMLFiles)
0066     generalParams[L1TMuonOverlapParams::GENERAL_ADDRBITS] = 2 * nPdfAddrBits;
0067   omtfParamsMutable.setGeneralParams(generalParams);
0068
0069   myOMTFConfig->configure(&omtfParamsMutable);
0070   myOMTF->configure(myOMTFConfig, omtfParams);
0071   myOMTFConfigMaker = new OMTFConfigMaker(myOMTFConfig);
0072
0073   ///Clear existing GoldenPatterns
0074   if (!mergeXMLFiles) {
0075     const std::map<Key, GoldenPattern*>& theGPs = myOMTF->getPatterns();
0076     for (auto itGP : theGPs)
0077       itGP.second->reset();
0078   }
0079 }
0080 /////////////////////////////////////////////////////
0081 /////////////////////////////////////////////////////
0082 void OMTFPatternMaker::endRun(edm::Run const&, edm::EventSetup const&) {}
0083 /////////////////////////////////////////////////////
0084 /////////////////////////////////////////////////////
0085 void OMTFPatternMaker::beginJob() {
0086   myOMTFConfig = new OMTFConfiguration();
0087   myOMTF = new OMTFProcessor();
0088 }
0089 /////////////////////////////////////////////////////
0090 /////////////////////////////////////////////////////
0091 void OMTFPatternMaker::endJob() {
0092   if (makeGoldenPatterns && !makeConnectionsMaps) {
0093     myWriter->initialiseXMLDocument("OMTF");
0094     const std::map<Key, GoldenPattern*>& myGPmap = myOMTF->getPatterns();
0095     for (auto itGP : myGPmap) {
0096       if (!itGP.second->hasCounts())
0097         continue;
0098       itGP.second->normalise(nPdfAddrBits);
0099     }
0100
0101     GoldenPattern dummyGP(Key(0, 0, 0), myOMTFConfig);
0102     dummyGP.reset();
0103
0104     ///Put back default value of the pdf width.
0105     L1TMuonOverlapParams omtfParamsMutable = *myOMTFConfig->getRawParams();
0106     std::vector<int> generalParams = *omtfParamsMutable.generalParams();
0107     generalParams[L1TMuonOverlapParams::GENERAL_ADDRBITS] = nPdfAddrBits;
0108     omtfParamsMutable.setGeneralParams(generalParams);
0109     myOMTFConfig->configure(&omtfParamsMutable);
0110
0111     for (auto itGP : myGPmap) {
0112       ////
0113       unsigned int iPt = theConfig.getParameter<int>("ptCode") + 1;
0114       if (iPt > 31)
0115         iPt = 200 * 2 + 1;
0116       else
0117         iPt = RPCConst::ptFromIpt(iPt) * 2.0 +
0118               1;  //MicroGMT has 0.5 GeV step size, with lower bin edge  (uGMT_pt_code - 1)*step_size
0119       ////
0120       if (itGP.first.thePtCode == iPt && itGP.first.theCharge == theConfig.getParameter<int>("charge")) {
0121         std::cout << *itGP.second << std::endl;
0122         myWriter->writeGPData(*itGP.second, dummyGP, dummyGP, dummyGP);
0123       }
0124     }
0125     std::string fName = "GPs.xml";
0126     myWriter->finaliseXMLDocument(fName);
0127   }
0128
0129   if (makeConnectionsMaps && !makeGoldenPatterns) {
0130     myWriter->initialiseXMLDocument("OMTF");
0131     std::string fName = "Connections.xml";
0132     unsigned int iProcessor = 0;
0133     ///Order important: printPhiMap updates global vector in OMTFConfiguration
0134     myOMTFConfigMaker->printPhiMap(std::cout);
0135     myOMTFConfigMaker->printConnections(std::cout, iProcessor, 0);
0136     myOMTFConfigMaker->printConnections(std::cout, iProcessor, 1);
0137     myOMTFConfigMaker->printConnections(std::cout, iProcessor, 2);
0138     myOMTFConfigMaker->printConnections(std::cout, iProcessor, 3);
0139     myOMTFConfigMaker->printConnections(std::cout, iProcessor, 4);
0140     myOMTFConfigMaker->printConnections(std::cout, iProcessor, 5);
0141     myWriter->writeConnectionsData(myOMTFConfig->getMeasurements4D());
0142     myWriter->finaliseXMLDocument(fName);
0143   }
0144
0145   if (mergeXMLFiles) {
0146     GoldenPattern* dummy = new GoldenPattern(Key(0, 0, 0), myOMTFConfig);
0147     dummy->reset();
0148
0149     std::string fName = "OMTF";
0150     myWriter->initialiseXMLDocument(fName);
0151     const std::map<Key, GoldenPattern*>& myGPmap = myOMTF->getPatterns();
0152     for (auto itGP : myGPmap) {
0153       myWriter->writeGPData(*itGP.second, *dummy, *dummy, *dummy);
0154     }
0155     fName = "GPs.xml";
0156     myWriter->finaliseXMLDocument(fName);
0157     ///Write GPs merged by 4 above iPt=71, and by 2 below//
0158     //////////////////////////////////////////////////////
0159     ///4x merging
0160     fName = "OMTF";
0161     myWriter->initialiseXMLDocument(fName);
0162     myOMTF->averagePatterns(-1);
0163     myOMTF->averagePatterns(1);
0164     writeMergedGPs();
0165     fName = "GPs_4x.xml";
0166     myWriter->finaliseXMLDocument(fName);
0167   }
0168 }
0169 /////////////////////////////////////////////////////
0170 /////////////////////////////////////////////////////
0171 void OMTFPatternMaker::writeMergedGPs() {
0172   const std::map<Key, GoldenPattern*>& myGPmap = myOMTF->getPatterns();
0173
0174   GoldenPattern* dummy = new GoldenPattern(Key(0, 0, 0), myOMTFConfig);
0175   dummy->reset();
0176
0177   unsigned int iPtMin = 9;
0178   Key aKey = Key(0, iPtMin, 1);
0179   while (myGPmap.find(aKey) != myGPmap.end()) {
0180     GoldenPattern* aGP1 = myGPmap.find(aKey)->second;
0181     GoldenPattern* aGP2 = dummy;
0182     GoldenPattern* aGP3 = dummy;
0183     GoldenPattern* aGP4 = dummy;
0184
0185     ++aKey.thePtCode;
0186     while (myGPmap.find(aKey) == myGPmap.end() && aKey.thePtCode <= 401)
0187       ++aKey.thePtCode;
0188     if (aKey.thePtCode <= 401 && myGPmap.find(aKey) != myGPmap.end())
0189       aGP2 = myGPmap.find(aKey)->second;
0190
0191     if (aKey.thePtCode > 71) {
0192       ++aKey.thePtCode;
0193       while (myGPmap.find(aKey) == myGPmap.end() && aKey.thePtCode <= 401)
0194         ++aKey.thePtCode;
0195       if (aKey.thePtCode <= 401 && myGPmap.find(aKey) != myGPmap.end())
0196         aGP3 = myGPmap.find(aKey)->second;
0197
0198       ++aKey.thePtCode;
0199       while (myGPmap.find(aKey) == myGPmap.end() && aKey.thePtCode <= 401)
0200         ++aKey.thePtCode;
0201       if (aKey.thePtCode <= 401 && myGPmap.find(aKey) != myGPmap.end())
0202         aGP4 = myGPmap.find(aKey)->second;
0203     }
0204     ++aKey.thePtCode;
0205     while (myGPmap.find(aKey) == myGPmap.end() && aKey.thePtCode <= 401)
0206       ++aKey.thePtCode;
0207     myWriter->writeGPData(*aGP1, *aGP2, *aGP3, *aGP4);
0208
0209     ///Write the opposite charge.
0210     Key aTmpKey = aGP1->key();
0211     aTmpKey.theCharge = -1;
0212     if (myGPmap.find(aTmpKey) != myGPmap.end())
0213       aGP1 = myGPmap.find(aTmpKey)->second;
0214     else
0215       aGP1 = dummy;
0216
0217     aTmpKey = aGP2->key();
0218     aTmpKey.theCharge = -1;
0219     if (myGPmap.find(aTmpKey) != myGPmap.end())
0220       aGP2 = myGPmap.find(aTmpKey)->second;
0221     else
0222       aGP2 = dummy;
0223
0224     aTmpKey = aGP3->key();
0225     aTmpKey.theCharge = -1;
0226     if (myGPmap.find(aTmpKey) != myGPmap.end())
0227       aGP3 = myGPmap.find(aTmpKey)->second;
0228     else
0229       aGP3 = dummy;
0230
0231     aTmpKey = aGP4->key();
0232     aTmpKey.theCharge = -1;
0233     if (myGPmap.find(aTmpKey) != myGPmap.end())
0234       aGP4 = myGPmap.find(aTmpKey)->second;
0235     else
0236       aGP4 = dummy;
0237
0238     myWriter->writeGPData(*aGP1, *aGP2, *aGP3, *aGP4);
0239   }
0240 }
0241 /////////////////////////////////////////////////////
0242 /////////////////////////////////////////////////////
0243 void OMTFPatternMaker::analyze(const edm::Event& iEvent, const edm::EventSetup& evSetup) {
0244   if (mergeXMLFiles)
0245     return;
0246
0247   ///Get the simulated muon parameters
0248   const SimTrack* aSimMuon = findSimMuon(iEvent);
0249   if (!aSimMuon) {
0250     edm::LogError("OMTFPatternMaker") << "No SimMuon found in the event!";
0251     return;
0252   }
0253
0254   myInputMaker->initialize(evSetup, myOMTFConfig);
0255
0256   edm::Handle<L1MuDTChambPhContainer> dtPhDigis;
0257   edm::Handle<L1MuDTChambThContainer> dtThDigis;
0258   edm::Handle<CSCCorrelatedLCTDigiCollection> cscDigis;
0259   edm::Handle<RPCDigiCollection> rpcDigis;
0260
0261   ///Filter digis by dropping digis from selected (by cfg.py) subsystems
0262   if (!theConfig.getParameter<bool>("dropDTPrimitives")) {
0263     iEvent.getByToken(inputTokenDTPh, dtPhDigis);
0264     iEvent.getByToken(inputTokenDTTh, dtThDigis);
0265   }
0266   if (!theConfig.getParameter<bool>("dropRPCPrimitives"))
0267     iEvent.getByToken(inputTokenRPC, rpcDigis);
0268   if (!theConfig.getParameter<bool>("dropCSCPrimitives"))
0269     iEvent.getByToken(inputTokenCSC, cscDigis);
0270
0271   //l1t::tftype mtfType = l1t::tftype::bmtf;
0272   l1t::tftype mtfType = l1t::tftype::omtf_pos;
0273   //l1t::tftype mtfType = l1t::tftype::emtf_pos;
0274
0275   ///Loop over all processors, each covering 60 deg in phi
0276   for (unsigned int iProcessor = 0; iProcessor < 6; ++iProcessor) {
0277     ///Input data with phi ranges shifted for each processor, so it fits 11 bits range
0278     OMTFinput myInput = myInputMaker->buildInputForProcessor(
0279         dtPhDigis.product(), dtThDigis.product(), cscDigis.product(), rpcDigis.product(), iProcessor, mtfType);
0280
0281     ///Connections maps are made by hand. makeConnetionsMap method
0282     ///provides tables for checking their consistency.
0283     if (makeConnectionsMaps)
0284       myOMTFConfigMaker->makeConnetionsMap(iProcessor, myInput);
0285
0286     if (makeGoldenPatterns)
0287       myOMTF->fillCounts(iProcessor, myInput, aSimMuon);
0288   }
0289 }
0290 /////////////////////////////////////////////////////
0291 /////////////////////////////////////////////////////
0292 const SimTrack* OMTFPatternMaker::findSimMuon(const edm::Event& ev, const SimTrack* previous) {
0293   const SimTrack* result = nullptr;
0294   edm::Handle<edm::SimTrackContainer> simTks;
0295   ev.getByToken(inputTokenSimHit, simTks);
0296
0297   for (std::vector<SimTrack>::const_iterator it = simTks->begin(); it < simTks->end(); it++) {
0298     const SimTrack& aTrack = *it;
0299     if (!(aTrack.type() == 13 || aTrack.type() == -13))
0300       continue;
0301     if (previous && ROOT::Math::VectorUtil::DeltaR(aTrack.momentum(), previous->momentum()) < 0.07)
0302       continue;
0303     if (!result || aTrack.momentum().pt() > result->momentum().pt())
0304       result = &aTrack;
0305   }
0306   return result;
0307 }
0308 /////////////////////////////////////////////////////
0309 /////////////////////////////////////////////////////
0310 #include "FWCore/Framework/interface/MakerMacros.h"
0311 DEFINE_FWK_MODULE(OMTFPatternMaker);