GlobalMuonProducer/test/GLBMuonAnalyzer.cc

0001 /** \class GLBMuonAnalyzer
0002  *  Analyzer of the Global muon tracks
0003  *
0004  *  \author R. Bellan  - INFN Torino       <riccardo.bellan@cern.ch>
0005  *  \author A. Everett - Purdue University <adam.everett@cern.ch>
0006  */
0007
0008 // Base Class Headers
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010
0011 // Collaborating Class Header
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   /// Constructor
0040   GLBMuonAnalyzer(const edm::ParameterSet &pset);
0041
0042   /// Destructor
0043   ~GLBMuonAnalyzer() override;
0044
0045   // Operations
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   // Histograms
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   // Counters
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 /// Constructor
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 /// Destructor
0100 GLBMuonAnalyzer::~GLBMuonAnalyzer() {}
0101
0102 void GLBMuonAnalyzer::beginJob() {
0103   // Create the root file
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   // Write the histos to file
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   //LogTrace("Analyzer") << "Run: " << event.id().run() << " Event: " << event.id().event() << endl;
0142   MuonPatternRecoDumper debug;
0143
0144   // Get the RecTrack collection from the event
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   // Get the SimTrack collection from the event
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     //LogTrace("Analyzer")<<"Simulated tracks: "<<endl;
0165     for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0166       if (abs((*simTrack).type()) == 13) {
0167         //LogTrace("Analyzer")<<"Sim pT: "<<(*simTrack).momentum().perp()<<endl;
0168         simPt = (*simTrack).momentum().pt();
0169         //LogTrace("Analyzer")<<"Sim Eta: "<<(*simTrack).momentum().eta()<<endl;
0170         numberOfSimTracks++;
0171       }
0172     }
0173     LogTrace("Analyzer") << endl;
0174   }
0175
0176   reco::TrackCollection::const_iterator staTrack;
0177
0178   //LogTrace("Analyzer")<<"Reconstructed tracks: " << staTracks->size() << endl;
0179
0180   for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack) {
0181     reco::TransientTrack track(*staTrack, &*theMGField, theTrackingGeometry);
0182
0183     //LogTrace("Analyzer") << debug.dumpFTS(track.impactPointTSCP().theState());
0184
0185     recPt = track.impactPointTSCP().momentum().perp();
0186     //LogTrace("Analyzer")<<" p: "<<track.impactPointTSCP().momentum().mag()<< " pT: "<<recPt<<endl;
0187     //LogTrace("Analyzer")<<" chi2: "<<track.chi2()<<endl;
0188
0189     hPtRec->Fill(recPt);
0190
0191     TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0192     //LogTrace("Analyzer") << "Inner TSOS:"<<endl;
0193     //LogTrace("Analyzer") << debug.dumpTSOS(innerTSOS);
0194     //LogTrace("Analyzer")<<" p: "<<innerTSOS.globalMomentum().mag()<< " pT: "<<innerTSOS.globalMomentum().perp()<<endl;
0195
0196     trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
0197     trackingRecHit_iterator rhend = staTrack->recHitsEnd();
0198
0199     int tkHit = 0;
0200
0201     //LogTrace("Analyzer")<<"RecHits:"<<endl;
0202     for (trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit) {
0203       //      const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
0204       //      double r = geomDet->surface().position().perp();
0205       //      double z = geomDet->toGlobal((*recHit)->localPosition()).z();
0206       //LogTrace("Analyzer")<<"r: "<< r <<" z: "<<z <<endl;
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       //      hPTDiff2->Fill(track.innermostMeasurementState().globalMomentum().perp()-simPt);
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);