File indexing completed on 2023-05-08 23:19:24
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010
0011
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019
0020 #include "MagneticField/Engine/interface/MagneticField.h"
0021 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0022 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0023 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0024 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0025
0026 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0027 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0028 #include "DataFormats/TrackReco/interface/Track.h"
0029 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0030
0031 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0032
0033 #include "TFile.h"
0034 #include "TH1F.h"
0035 #include "TH2F.h"
0036
0037 class GLBMuonAnalyzer : public edm::one::EDAnalyzer<> {
0038 public:
0039
0040 GLBMuonAnalyzer(const edm::ParameterSet &pset);
0041
0042
0043 ~GLBMuonAnalyzer() override;
0044
0045
0046
0047 void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override;
0048
0049 void beginJob() override;
0050 void endJob() override;
0051
0052 protected:
0053 private:
0054 edm::InputTag theGLBMuonLabel;
0055
0056 std::string theRootFileName;
0057 TFile *theFile;
0058
0059
0060 TH1F *hPtRec;
0061 TH1F *hPtSim;
0062 TH1F *hPres;
0063 TH1F *h1_Pres;
0064 TH1F *hPTDiff;
0065 TH1F *hPTDiff2;
0066 TH2F *hPTDiffvsEta;
0067 TH2F *hPTDiffvsPhi;
0068
0069
0070 int numberOfSimTracks;
0071 int numberOfRecTracks;
0072
0073 std::string theDataType;
0074
0075 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMGFieldToken;
0076 edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> theTrackingGeometryToken;
0077 };
0078
0079 using namespace std;
0080 using namespace edm;
0081
0082
0083 GLBMuonAnalyzer::GLBMuonAnalyzer(const ParameterSet &pset) {
0084 theMGFieldToken = esConsumes();
0085 theTrackingGeometryToken = esConsumes();
0086 theGLBMuonLabel = pset.getUntrackedParameter<edm::InputTag>("GlobalTrackCollectionLabel");
0087
0088 theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0089
0090 theDataType = pset.getUntrackedParameter<string>("DataType");
0091
0092 if (theDataType != "RealData" && theDataType != "SimData")
0093 LogTrace("Analyzer") << "Error in Data Type!!" << endl;
0094
0095 numberOfSimTracks = 0;
0096 numberOfRecTracks = 0;
0097 }
0098
0099
0100 GLBMuonAnalyzer::~GLBMuonAnalyzer() {}
0101
0102 void GLBMuonAnalyzer::beginJob() {
0103
0104 theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0105 theFile->cd();
0106
0107 hPtRec = new TH1F("pTRec", "p_{T}^{rec}", 250, 0, 120);
0108 hPtSim = new TH1F("pTSim", "p_{T}^{gen} ", 250, 0, 120);
0109
0110 hPTDiff = new TH1F("pTDiff", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0111 hPTDiff2 = new TH1F("pTDiff2", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0112
0113 hPTDiffvsEta = new TH2F("PTDiffvsEta", "p_{T}^{rec} - p_{T}^{gen} VS #eta", 100, -2.5, 2.5, 250, -120, 120);
0114 hPTDiffvsPhi = new TH2F("PTDiffvsPhi", "p_{T}^{rec} - p_{T}^{gen} VS #phi", 100, -6, 6, 250, -120, 120);
0115
0116 hPres = new TH1F("pTRes", "pT Resolution", 100, -2, 2);
0117 h1_Pres = new TH1F("invPTRes", "1/pT Resolution", 100, -2, 2);
0118 }
0119
0120 void GLBMuonAnalyzer::endJob() {
0121 if (theDataType == "SimData") {
0122 LogTrace("Analyzer") << endl << endl << "Number of Sim tracks: " << numberOfSimTracks << endl;
0123 }
0124
0125 LogTrace("Analyzer") << "Number of Reco tracks: " << numberOfRecTracks << endl << endl;
0126
0127
0128 theFile->cd();
0129 hPtRec->Write();
0130 hPtSim->Write();
0131 hPres->Write();
0132 h1_Pres->Write();
0133 hPTDiff->Write();
0134 hPTDiff2->Write();
0135 hPTDiffvsEta->Write();
0136 hPTDiffvsPhi->Write();
0137 theFile->Close();
0138 }
0139
0140 void GLBMuonAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
0141
0142 MuonPatternRecoDumper debug;
0143
0144
0145 Handle<reco::TrackCollection> staTracks;
0146 event.getByLabel(theGLBMuonLabel, staTracks);
0147
0148 ESHandle<MagneticField> theMGField = eventSetup.getHandle(theMGFieldToken);
0149
0150 ESHandle<GlobalTrackingGeometry> theTrackingGeometry = eventSetup.getHandle(theTrackingGeometryToken);
0151
0152 double recPt = 0.;
0153 double simPt = 0.;
0154
0155
0156 if (theDataType == "SimData") {
0157 Handle<SimTrackContainer> simTracks;
0158 event.getByLabel("g4SimHits", simTracks);
0159
0160 numberOfRecTracks += staTracks->size();
0161
0162 SimTrackContainer::const_iterator simTrack;
0163
0164
0165 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0166 if (abs((*simTrack).type()) == 13) {
0167
0168 simPt = (*simTrack).momentum().pt();
0169
0170 numberOfSimTracks++;
0171 }
0172 }
0173 LogTrace("Analyzer") << endl;
0174 }
0175
0176 reco::TrackCollection::const_iterator staTrack;
0177
0178
0179
0180 for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack) {
0181 reco::TransientTrack track(*staTrack, &*theMGField, theTrackingGeometry);
0182
0183
0184
0185 recPt = track.impactPointTSCP().momentum().perp();
0186
0187
0188
0189 hPtRec->Fill(recPt);
0190
0191 TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0192
0193
0194
0195
0196 trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
0197 trackingRecHit_iterator rhend = staTrack->recHitsEnd();
0198
0199 int tkHit = 0;
0200
0201
0202 for (trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit) {
0203
0204
0205
0206
0207 if ((*recHit)->geographicalId().det() == DetId::Tracker)
0208 ++tkHit;
0209 }
0210 if (tkHit == 0)
0211 LogTrace("GlobalMuonAnalyzer") << "+++++++++ This track does not contain TH hits +++++++++ ";
0212
0213 if (recPt && theDataType == "SimData") {
0214 hPres->Fill((recPt - simPt) / simPt);
0215 hPtSim->Fill(simPt);
0216
0217 hPTDiff->Fill(recPt - simPt);
0218
0219
0220 hPTDiffvsEta->Fill(track.impactPointTSCP().position().eta(), recPt - simPt);
0221 hPTDiffvsPhi->Fill(track.impactPointTSCP().position().phi(), recPt - simPt);
0222
0223 if (((recPt - simPt) / simPt) <= -0.4)
0224 LogTrace("Analyzer") << "Out of Res: " << (recPt - simPt) / simPt << endl;
0225 h1_Pres->Fill((1 / recPt - 1 / simPt) / (1 / simPt));
0226 }
0227 }
0228 LogTrace("Analyzer") << "---" << endl;
0229 }
0230
0231 DEFINE_FWK_MODULE(GLBMuonAnalyzer);