File indexing completed on 2023-10-25 10:01:31
0001
0002
0003
0004
0005
0006
0007 #include "RecoMuon/StandAloneMuonProducer/test/STAMuonAnalyzer.h"
0008
0009
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014
0015 #include "MagneticField/Engine/interface/MagneticField.h"
0016 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0017 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0018 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0019 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0020
0021 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0022 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0025
0026 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0027
0028 #include "TFile.h"
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031
0032 using namespace std;
0033 using namespace edm;
0034
0035
0036 STAMuonAnalyzer::STAMuonAnalyzer(const ParameterSet& pset) {
0037 theSTAMuonLabel = pset.getUntrackedParameter<string>("StandAloneTrackCollectionLabel");
0038 theSeedCollectionLabel = pset.getUntrackedParameter<string>("MuonSeedCollectionLabel");
0039
0040 theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0041
0042 theDataType = pset.getUntrackedParameter<string>("DataType");
0043
0044 if (theDataType != "RealData" && theDataType != "SimData")
0045 cout << "Error in Data Type!!" << endl;
0046
0047 numberOfSimTracks = 0;
0048 numberOfRecTracks = 0;
0049
0050 theFieldToken = esConsumes();
0051 theGeomToken = esConsumes();
0052 }
0053
0054
0055 STAMuonAnalyzer::~STAMuonAnalyzer() {}
0056
0057 void STAMuonAnalyzer::beginJob() {
0058
0059 theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0060 theFile->cd();
0061
0062 hPtRec = new TH1F("pTRec", "p_{T}^{rec}", 250, 0, 120);
0063 hPtSim = new TH1F("pTSim", "p_{T}^{gen} ", 250, 0, 120);
0064
0065 hPTDiff = new TH1F("pTDiff", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0066 hPTDiff2 = new TH1F("pTDiff2", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0067
0068 hPTDiffvsEta = new TH2F("PTDiffvsEta", "p_{T}^{rec} - p_{T}^{gen} VS #eta", 100, -2.5, 2.5, 250, -120, 120);
0069 hPTDiffvsPhi = new TH2F("PTDiffvsPhi", "p_{T}^{rec} - p_{T}^{gen} VS #phi", 100, -6, 6, 250, -120, 120);
0070
0071 hPres = new TH1F("pTRes", "pT Resolution", 100, -2, 2);
0072 h1_Pres = new TH1F("invPTRes", "1/pT Resolution", 100, -2, 2);
0073 }
0074
0075 void STAMuonAnalyzer::endJob() {
0076 if (theDataType == "SimData") {
0077 cout << endl << endl << "Number of Sim tracks: " << numberOfSimTracks << endl;
0078 }
0079
0080 cout << "Number of Reco tracks: " << numberOfRecTracks << endl << endl;
0081
0082
0083 theFile->cd();
0084 hPtRec->Write();
0085 hPtSim->Write();
0086 hPres->Write();
0087 h1_Pres->Write();
0088 hPTDiff->Write();
0089 hPTDiff2->Write();
0090 hPTDiffvsEta->Write();
0091 hPTDiffvsPhi->Write();
0092 theFile->Close();
0093 }
0094
0095 void STAMuonAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0096 cout << "Run: " << event.id().run() << " Event: " << event.id().event() << endl;
0097 MuonPatternRecoDumper debug;
0098
0099
0100 Handle<reco::TrackCollection> staTracks;
0101 event.getByLabel(theSTAMuonLabel, staTracks);
0102
0103 ESHandle<MagneticField> theMGField = eventSetup.getHandle(theFieldToken);
0104 ESHandle<GlobalTrackingGeometry> theTrackingGeometry = eventSetup.getHandle(theGeomToken);
0105
0106 double recPt = 0.;
0107 double simPt = 0.;
0108
0109
0110 if (theDataType == "SimData") {
0111 Handle<SimTrackContainer> simTracks;
0112 event.getByLabel("g4SimHits", simTracks);
0113
0114 numberOfRecTracks += staTracks->size();
0115
0116 SimTrackContainer::const_iterator simTrack;
0117
0118 cout << "Simulated tracks: " << endl;
0119 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0120 if (abs((*simTrack).type()) == 13) {
0121 cout << "Sim pT: " << (*simTrack).momentum().pt() << endl;
0122 simPt = (*simTrack).momentum().pt();
0123 cout << "Sim Eta: " << (*simTrack).momentum().eta() << endl;
0124 numberOfSimTracks++;
0125 }
0126 }
0127 cout << endl;
0128 }
0129
0130 reco::TrackCollection::const_iterator staTrack;
0131
0132 cout << "Reconstructed tracks: " << staTracks->size() << endl;
0133
0134 for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack) {
0135 reco::TransientTrack track(*staTrack, &*theMGField, theTrackingGeometry);
0136
0137 cout << debug.dumpFTS(track.impactPointTSCP().theState());
0138
0139 recPt = track.impactPointTSCP().momentum().perp();
0140 cout << " p: " << track.impactPointTSCP().momentum().mag() << " pT: " << recPt << endl;
0141 cout << " chi2: " << track.chi2() << endl;
0142
0143 hPtRec->Fill(recPt);
0144
0145 TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0146 cout << "Inner TSOS:" << endl;
0147 cout << debug.dumpTSOS(innerTSOS);
0148 cout << " p: " << innerTSOS.globalMomentum().mag() << " pT: " << innerTSOS.globalMomentum().perp() << endl;
0149
0150 trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
0151 trackingRecHit_iterator rhend = staTrack->recHitsEnd();
0152
0153 cout << "RecHits:" << endl;
0154 for (trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit) {
0155 const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
0156 double r = geomDet->surface().position().perp();
0157 double z = geomDet->toGlobal((*recHit)->localPosition()).z();
0158 cout << "r: " << r << " z: " << z << endl;
0159 }
0160
0161 if (recPt && theDataType == "SimData") {
0162 hPres->Fill((recPt - simPt) / simPt);
0163 hPtSim->Fill(simPt);
0164
0165 hPTDiff->Fill(recPt - simPt);
0166
0167
0168 hPTDiffvsEta->Fill(track.impactPointTSCP().position().eta(), recPt - simPt);
0169 hPTDiffvsPhi->Fill(track.impactPointTSCP().position().phi(), recPt - simPt);
0170
0171 if (((recPt - simPt) / simPt) <= -0.4)
0172 cout << "Out of Res: " << (recPt - simPt) / simPt << endl;
0173 h1_Pres->Fill((1 / recPt - 1 / simPt) / (1 / simPt));
0174 }
0175 }
0176 cout << "---" << endl;
0177 }
0178
0179 DEFINE_FWK_MODULE(STAMuonAnalyzer);