Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-15 04:21:52 

0001 /*
0002  * DataROOTDumper2.cc
0003  *
0004  *  Created on: Dec 11, 2019
0005  *      Author: kbunkow
0006  */
0007 
0008 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Tools/DataROOTDumper2.h"
0009 
0010 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0015 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0016 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0017 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0018 
0019 #include "TFile.h"
0020 #include "TTree.h"
0021 
0022 DataROOTDumper2::DataROOTDumper2(const edm::ParameterSet& edmCfg,
0023                                  const OMTFConfiguration* omtfConfig,
0024                                  CandidateSimMuonMatcher* candidateSimMuonMatcher)
0025     : EmulationObserverBase(edmCfg, omtfConfig), candidateSimMuonMatcher(candidateSimMuonMatcher) {
0026   edm::LogVerbatim("l1tOmtfEventPrint") << " omtfConfig->nTestRefHits() " << omtfConfig->nTestRefHits()
0027                                         << " event.omtfGpResultsPdfSum.num_elements() " << endl;
0028   initializeTTree();
0029 
0030   if (edmCfg.exists("dumpKilledOmtfCands"))
0031     if (edmCfg.getParameter<bool>("dumpKilledOmtfCands"))
0032       dumpKilledOmtfCands = true;
0033 
0034   if (edmCfg.exists("candidateSimMuonMatcherType")) {
0035     if (edmCfg.getParameter<std::string>("candidateSimMuonMatcherType") == "propagation")
0036       usePropagation = true;
0037     else if (edmCfg.getParameter<std::string>("candidateSimMuonMatcherType") == "matchSimple")
0038       usePropagation = false;
0039 
0040     edm::LogImportant("l1tOmtfEventPrint")
0041         << " CandidateSimMuonMatcher: candidateSimMuonMatcherType "
0042         << edmCfg.getParameter<std::string>("candidateSimMuonMatcherType") << std::endl;
0043   }
0044 
0045   edm::LogVerbatim("l1tOmtfEventPrint") << " DataROOTDumper2 created. dumpKilledOmtfCands " << dumpKilledOmtfCands
0046                                         << std::endl;
0047 }
0048 
0049 DataROOTDumper2::~DataROOTDumper2() {}
0050 
0051 void DataROOTDumper2::initializeTTree() {
0052   edm::Service<TFileService> fs;
0053 
0054   rootTree = fs->make<TTree>("OMTFHitsTree", "");
0055 
0056   rootTree->Branch("eventNum", &omtfEvent.eventNum);
0057   rootTree->Branch("muonEvent", &omtfEvent.muonEvent);
0058 
0059   rootTree->Branch("muonPt", &omtfEvent.muonPt);
0060   rootTree->Branch("muonEta", &omtfEvent.muonEta);
0061   rootTree->Branch("muonPhi", &omtfEvent.muonPhi);
0062   rootTree->Branch("muonPropEta", &omtfEvent.muonPropEta);
0063   rootTree->Branch("muonPropPhi", &omtfEvent.muonPropPhi);
0064   rootTree->Branch("muonCharge", &omtfEvent.muonCharge);
0065 
0066   rootTree->Branch("muonDxy", &omtfEvent.muonDxy);
0067   rootTree->Branch("muonRho", &omtfEvent.muonRho);
0068 
0069   rootTree->Branch("omtfPt", &omtfEvent.omtfPt);
0070   rootTree->Branch("omtfUPt", &omtfEvent.omtfUPt);
0071   rootTree->Branch("omtfEta", &omtfEvent.omtfEta);
0072   rootTree->Branch("omtfPhi", &omtfEvent.omtfPhi);
0073   rootTree->Branch("omtfCharge", &omtfEvent.omtfCharge);
0074 
0075   rootTree->Branch("omtfHwEta", &omtfEvent.omtfHwEta);
0076 
0077   rootTree->Branch("omtfProcessor", &omtfEvent.omtfProcessor);
0078   rootTree->Branch("omtfScore", &omtfEvent.omtfScore);
0079   rootTree->Branch("omtfQuality", &omtfEvent.omtfQuality);
0080   rootTree->Branch("omtfRefLayer", &omtfEvent.omtfRefLayer);
0081   rootTree->Branch("omtfRefHitNum", &omtfEvent.omtfRefHitNum);
0082 
0083   rootTree->Branch("omtfFiredLayers", &omtfEvent.omtfFiredLayers);  //<<<<<<<<<<<<<<<<<<<<<<!!!!TODOO
0084 
0085   rootTree->Branch("killed", &omtfEvent.killed);
0086 
0087   rootTree->Branch("hits", &omtfEvent.hits);
0088 
0089   rootTree->Branch("deltaEta", &omtfEvent.deltaEta);
0090   rootTree->Branch("deltaPhi", &omtfEvent.deltaPhi);
0091 
0092   ptGenPos = fs->make<TH1I>("ptGenPos", "ptGenPos, eta at vertex 0.8 - 1.24", 400, 0, 200);  //TODO
0093   ptGenNeg = fs->make<TH1I>("ptGenNeg", "ptGenNeg, eta at vertex 0.8 - 1.24", 400, 0, 200);
0094 }
0095 
0096 void DataROOTDumper2::observeProcesorEmulation(unsigned int iProcessor,
0097                                                l1t::tftype mtfType,
0098                                                const std::shared_ptr<OMTFinput>&,
0099                                                const AlgoMuons& algoCandidates,
0100                                                const AlgoMuons& gbCandidates,
0101                                                const std::vector<l1t::RegionalMuonCand>& candMuons) {}
0102 
0103 void DataROOTDumper2::observeEventEnd(const edm::Event& iEvent,
0104                                       std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
0105   /*
0106   int muonCharge = 0;
0107   if (simMuon) {
0108     if (std::abs(simMuon->momentum().eta()) < 0.8 || std::abs(simMuon->momentum().eta()) > 1.24)
0109       return;
0110 
0111     muonCharge = (std::abs(simMuon->type()) == 13) ? simMuon->type() / -13 : 0;
0112     if (muonCharge > 0)
0113       ptGenPos->Fill(simMuon->momentum().pt());
0114     else
0115       ptGenNeg->Fill(simMuon->momentum().pt());
0116   }
0117 
0118   if (simMuon == nullptr || !omtfCand->isValid())  //no sim muon or empty candidate
0119     return;
0120 
0121   omtfEvent.muonPt = simMuon->momentum().pt();
0122   omtfEvent.muonEta = simMuon->momentum().eta();
0123 
0124   //TODO add cut on ete if needed
0125     if(std::abs(event.muonEta) < 0.8 || std::abs(event.muonEta) > 1.24)
0126     return;
0127 
0128   omtfEvent.muonPhi = simMuon->momentum().phi();
0129   omtfEvent.muonCharge = muonCharge;  //TODO
0130    */
0131 
0132   std::vector<MatchingResult> matchingResults = candidateSimMuonMatcher->getMatchingResults();
0133   LogTrace("l1tOmtfEventPrint") << "\nDataROOTDumper2::observeEventEnd matchingResults.size() "
0134                                 << matchingResults.size() << std::endl;
0135 
0136   //candidateSimMuonMatcher should use the  trackingParticles, because the simTracks are not stored for the pile-up events
0137 
0138   //for some events there are more than one matchingResults,
0139   //Usually at least one them has  genPt 0, which means no simMuon was matched, so candidate is ghost (or fake)
0140   //so better is to to drop such event, as it is not sue if the correct simMuon was matched to the candidate.
0141   //So we assume here that when the propagation is not used it is a single mu sample and this filter has sense
0142   //the propagation is used for multi-muon sample, so then this fitler cannot be used
0143   //TODO add a flag to enable this filter? Disable it if not needed
0144   if (!usePropagation && matchingResults.size() > 1) {  //omtfConfig->cleanStubs() &&
0145     edm::LogVerbatim("l1tOmtfEventPrint")
0146         << "\nDataROOTDumper2::observeEventEnd matchingResults.size() " << matchingResults.size() << std::endl;
0147 
0148     for (auto& matchingResult : matchingResults) {
0149       edm::LogVerbatim("l1tOmtfEventPrint") << "matchingResult: genPt " << matchingResult.genPt;
0150       if (matchingResult.procMuon)
0151         edm::LogVerbatim("l1tOmtfEventPrint") << " procMuon.PtConstr " << matchingResult.procMuon->getPtConstr();
0152       else
0153         edm::LogVerbatim("l1tOmtfEventPrint") << " no procMuon" << std::endl;
0154     }
0155     edm::LogVerbatim("l1tOmtfEventPrint") << "dropping the event!!!\n" << std::endl;
0156     return;
0157   }
0158 
0159   for (auto& matchingResult : matchingResults) {
0160     omtfEvent.eventNum = iEvent.id().event();
0161 
0162     if (matchingResult.trackingParticle) {
0163       auto trackingParticle = matchingResult.trackingParticle;
0164 
0165       omtfEvent.muonEvent = trackingParticle->eventId().event();
0166 
0167       omtfEvent.muonPt = trackingParticle->pt();
0168       omtfEvent.muonEta = trackingParticle->momentum().eta();
0169       omtfEvent.muonPhi = trackingParticle->momentum().phi();
0170       omtfEvent.muonPropEta = matchingResult.propagatedEta;
0171       omtfEvent.muonPropPhi = matchingResult.propagatedPhi;
0172       omtfEvent.muonCharge = (std::abs(trackingParticle->pdgId()) == 13) ? trackingParticle->pdgId() / -13 : 0;
0173 
0174       if (trackingParticle->parentVertex().isNonnull()) {
0175         omtfEvent.muonDxy = trackingParticle->dxy();
0176         omtfEvent.muonRho = trackingParticle->parentVertex()->position().Rho();
0177       }
0178 
0179       omtfEvent.deltaEta = matchingResult.deltaEta;
0180       omtfEvent.deltaPhi = matchingResult.deltaPhi;
0181 
0182       LogTrace("l1tOmtfEventPrint") << "DataROOTDumper2::observeEventEnd trackingParticle: eventId "
0183                                     << trackingParticle->eventId().event() << " pdgId " << std::setw(3)
0184                                     << trackingParticle->pdgId() << " trackId "
0185                                     << trackingParticle->g4Tracks().at(0).trackId() << " pt " << std::setw(9)
0186                                     << trackingParticle->pt()  //<<" Beta "<<simMuon->momentum().Beta()
0187                                     << " eta " << std::setw(9) << trackingParticle->momentum().eta() << " phi "
0188                                     << std::setw(9) << trackingParticle->momentum().phi() << std::endl;
0189 
0190       if (std::abs(omtfEvent.muonEta) > 0.8 && std::abs(omtfEvent.muonEta) < 1.24) {
0191         if (omtfEvent.muonCharge > 0)
0192           ptGenPos->Fill(omtfEvent.muonPt);
0193         else
0194           ptGenNeg->Fill(omtfEvent.muonPt);
0195       }
0196     } else if (matchingResult.simTrack) {
0197       auto simTrack = matchingResult.simTrack;
0198 
0199       omtfEvent.muonEvent = simTrack->eventId().event();
0200 
0201       omtfEvent.muonPt = simTrack->momentum().pt();
0202       omtfEvent.muonEta = simTrack->momentum().eta();
0203       omtfEvent.muonPhi = simTrack->momentum().phi();
0204       omtfEvent.muonPropEta = matchingResult.propagatedEta;
0205       omtfEvent.muonPropPhi = matchingResult.propagatedPhi;
0206       omtfEvent.muonCharge = simTrack->charge();
0207 
0208       if (!simTrack->noVertex() && matchingResult.simVertex) {
0209         const math::XYZTLorentzVectorD& vtxPos = matchingResult.simVertex->position();
0210         omtfEvent.muonDxy = (-vtxPos.X() * simTrack->momentum().py() + vtxPos.Y() * simTrack->momentum().px()) /
0211                             simTrack->momentum().pt();
0212         omtfEvent.muonRho = vtxPos.Rho();
0213       }
0214 
0215       omtfEvent.deltaEta = matchingResult.deltaEta;
0216       omtfEvent.deltaPhi = matchingResult.deltaPhi;
0217 
0218       LogTrace("l1tOmtfEventPrint") << "DataROOTDumper2::observeEventEnd simTrack: eventId "
0219                                     << simTrack->eventId().event() << " pdgId " << std::setw(3)
0220                                     << simTrack->type()  //<< " trackId " << simTrack->g4Tracks().at(0).trackId()
0221                                     << " pt " << std::setw(9)
0222                                     << simTrack->momentum().pt()  //<<" Beta "<<simMuon->momentum().Beta()
0223                                     << " eta " << std::setw(9) << simTrack->momentum().eta() << " phi " << std::setw(9)
0224                                     << simTrack->momentum().phi() << std::endl;
0225 
0226       if (std::abs(omtfEvent.muonEta) > 0.8 && std::abs(omtfEvent.muonEta) < 1.24) {
0227         if (omtfEvent.muonCharge > 0)
0228           ptGenPos->Fill(omtfEvent.muonPt);
0229         else
0230           ptGenNeg->Fill(omtfEvent.muonPt);
0231       }
0232     } else {
0233       omtfEvent.muonEvent = -1;
0234 
0235       omtfEvent.muonPt = 0;
0236 
0237       omtfEvent.muonEta = 0;
0238       omtfEvent.muonPhi = 0;
0239 
0240       omtfEvent.muonPropEta = 0;
0241       omtfEvent.muonPropPhi = 0;
0242 
0243       omtfEvent.muonCharge = 0;  //TODO
0244 
0245       omtfEvent.muonDxy = 0;
0246       omtfEvent.muonRho = 0;
0247     }
0248 
0249     auto addOmtfCand = [&](AlgoMuonPtr& procMuon) {
0250       omtfEvent.omtfPt = omtfConfig->hwPtToGev(procMuon->getPtConstr());
0251       omtfEvent.omtfUPt = omtfConfig->hwUPtToGev(procMuon->getPtUnconstr());
0252       omtfEvent.omtfEta = omtfConfig->hwEtaToEta(procMuon->getEtaHw());
0253       omtfEvent.omtfPhi = procMuon->getPhi();
0254       omtfEvent.omtfCharge = procMuon->getChargeConstr();
0255       omtfEvent.omtfScore = procMuon->getPdfSum();
0256 
0257       omtfEvent.omtfHwEta = procMuon->getEtaHw();
0258 
0259       omtfEvent.omtfFiredLayers = procMuon->getFiredLayerBits();
0260       omtfEvent.omtfRefLayer = procMuon->getRefLayer();
0261       omtfEvent.omtfRefHitNum = procMuon->getRefHitNumber();
0262 
0263       omtfEvent.hits.clear();
0264 
0265       //TODO choose, which gpResult should be dumped
0266       //auto& gpResult = procMuon->getGpResultConstr();
0267       auto& gpResult = (procMuon->getGpResultUnconstr().getPdfSumUnconstr() > procMuon->getGpResultConstr().getPdfSum())
0268                            ? procMuon->getGpResultUnconstr()
0269                            : procMuon->getGpResultConstr();
0270 
0271       /*
0272         edm::LogVerbatim("l1tOmtfEventPrint")<<"DataROOTDumper2:;observeEventEnd muonPt "<<event.muonPt<<" muonCharge "<<event.muonCharge
0273             <<" omtfPt "<<event.omtfPt<<" RefLayer "<<event.omtfRefLayer<<" omtfPtCont "<<event.omtfPtCont
0274             <<std::endl;  */
0275 
0276       for (unsigned int iLogicLayer = 0; iLogicLayer < gpResult.getStubResults().size(); ++iLogicLayer) {
0277         auto& stubResult = gpResult.getStubResults()[iLogicLayer];
0278 
0279         //TODO it is to have the hit if it is below the quality cut
0280         /*if (omtfConfig->isBendingLayer(iLogicLayer) && !stubResult.getMuonStub()) {
0281           auto& stubResult = gpResult.getStubResults()[iLogicLayer-1];
0282         }*/
0283 
0284         if (stubResult.getMuonStub()) {  //&& stubResult.getValid() //TODO!!!!!!!!!!!!!!!!1
0285           OmtfEvent::Hit hit;
0286           hit.layer = iLogicLayer;
0287           hit.quality = stubResult.getMuonStub()->qualityHw;
0288           hit.eta = stubResult.getMuonStub()->etaHw;  //in which scale?
0289           hit.valid = stubResult.getValid();
0290 
0291           int hitPhi = stubResult.getMuonStub()->phiHw;
0292           unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[procMuon->getRefLayer()];
0293           int phiRefHit = gpResult.getStubResults()[refLayerLogicNum].getMuonStub()->phiHw;
0294 
0295           if (omtfConfig->isBendingLayer(iLogicLayer)) {
0296             hitPhi = stubResult.getMuonStub()->phiBHw;
0297             phiRefHit = 0;  //phi ref hit for the bending layer set to 0, since it should not be included in the phiDist
0298           }
0299 
0300           //phiDist = hitPhi - phiRefHit;
0301           hit.phiDist = hitPhi - phiRefHit;
0302 
0303           /* LogTrace("l1tOmtfEventPrint")<<" muonPt "<<event.muonPt<<" omtfPt "<<event.omtfPt<<" RefLayer "<<event.omtfRefLayer
0304                 <<" layer "<<int(hit.layer)<<" PdfBin "<<stubResult.getPdfBin()<<" hit.phiDist "<<hit.phiDist<<" valid "<<stubResult.getValid()<<" " //<<" phiDist "<<phiDist
0305                 <<" getDistPhiBitShift "<<procMuon->getGoldenPatern()->getDistPhiBitShift(iLogicLayer, procMuon->getRefLayer())
0306                 <<" meanDistPhiValue   "<<procMuon->getGoldenPatern()->meanDistPhiValue(iLogicLayer, procMuon->getRefLayer())//<<(phiDist != hit.phiDist? "!!!!!!!<<<<<" : "")
0307                 <<endl;*/
0308 
0309           /*if (hit.phiDist > 504 || hit.phiDist < -512) {
0310             edm::LogVerbatim("l1tOmtfEventPrint")
0311                 << " muonPt " << omtfEvent.muonPt << " omtfPt " << omtfEvent.omtfPt << " RefLayer "
0312                 << (int)omtfEvent.omtfRefLayer << " layer " << int(hit.layer) << " hit.phiDist " << hit.phiDist
0313                 << " valid " << stubResult.getValid() << " !!!!!!!!!!!!!!!!!!!!!!!!" << endl;
0314           }*/
0315 
0316           DetId detId(stubResult.getMuonStub()->detId);
0317           if (detId.subdetId() == MuonSubdetId::CSC) {
0318             CSCDetId cscId(detId);
0319             hit.z = cscId.chamber() % 2;
0320           }
0321 
0322           omtfEvent.hits.push_back(hit.rawData);
0323         }
0324       }
0325 
0326       LogTrace("l1tOmtfEventPrint") << "DataROOTDumper2::observeEventEnd adding omtfCand : " << std::endl;
0327       auto finalCandidate = matchingResult.muonCand;
0328       LogTrace("l1tOmtfEventPrint") << " hwPt " << finalCandidate->hwPt() << " hwSign " << finalCandidate->hwSign()
0329                                     << " hwQual " << finalCandidate->hwQual() << " hwEta " << std::setw(4)
0330                                     << finalCandidate->hwEta() << std::setw(4) << " hwPhi " << finalCandidate->hwPhi()
0331                                     << "    eta " << std::setw(9) << (finalCandidate->hwEta() * 0.010875)
0332                                     << " isKilled " << procMuon->isKilled() << " tRefLayer " << procMuon->getRefLayer()
0333                                     << " RefHitNumber " << procMuon->getRefHitNumber() << std::endl;
0334     };
0335 
0336     if (matchingResult.muonCand && matchingResult.procMuon->getPtConstr() > 0 &&
0337         matchingResult.muonCand->hwQual() >= 1) {
0338       //TODO set the quality, quality 0 has the candidates with eta > 1.3(?) EtaHw >= 121
0339       //&& matchingResult.genPt < 20
0340 
0341       omtfEvent.omtfQuality = matchingResult.muonCand->hwQual();  //procMuon->getQ();
0342       omtfEvent.killed = false;
0343       omtfEvent.omtfProcessor = matchingResult.muonCand->processor();
0344 
0345       if (matchingResult.muonCand->trackFinderType() == l1t::omtf_neg) {
0346         omtfEvent.omtfProcessor *= -1;
0347       }
0348 
0349       addOmtfCand(matchingResult.procMuon);
0350       rootTree->Fill();
0351 
0352       if (dumpKilledOmtfCands) {
0353         for (auto& killedCand : matchingResult.procMuon->getKilledMuons()) {
0354           omtfEvent.omtfQuality = 0;
0355           omtfEvent.killed = true;
0356           if (killedCand->isKilled() == false) {
0357             edm::LogVerbatim("l1tOmtfEventPrint") << " killedCand->isKilled() == false !!!!!!!!";
0358           }
0359           addOmtfCand(killedCand);
0360           rootTree->Fill();
0361         }
0362       }
0363     } else if (omtfEvent.muonPt > 0) {  //checking if there was a simMuon
0364       LogTrace("l1tOmtfEventPrint") << "DataROOTDumper2::observeEventEnd no matching omtfCand" << std::endl;
0365 
0366       omtfEvent.omtfPt = 0;
0367       omtfEvent.omtfUPt = 0;
0368       omtfEvent.omtfEta = 0;
0369       omtfEvent.omtfPhi = 0;
0370       omtfEvent.omtfCharge = 0;
0371       omtfEvent.omtfScore = 0;
0372 
0373       omtfEvent.omtfHwEta = 0;
0374 
0375       omtfEvent.omtfFiredLayers = 0;
0376       omtfEvent.omtfRefLayer = 0;
0377       omtfEvent.omtfRefHitNum = 0;
0378       omtfEvent.omtfProcessor = 10;
0379 
0380       omtfEvent.omtfQuality = 0;
0381       omtfEvent.killed = false;
0382 
0383       omtfEvent.hits.clear();
0384 
0385       rootTree->Fill();
0386     }
0387   }
0388   evntCnt++;
0389 }
0390 
0391 void DataROOTDumper2::endJob() { edm::LogVerbatim("l1tOmtfEventPrint") << " evntCnt " << evntCnt << endl; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!