Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-20 02:15:05

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/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::EDAnalyzer {
0038 public:
0039   /// Constructor
0040   GLBMuonAnalyzer(const edm::ParameterSet &pset);
0041 
0042   /// Destructor
0043   virtual ~GLBMuonAnalyzer();
0044 
0045   // Operations
0046 
0047   void analyze(const edm::Event &event, const edm::EventSetup &eventSetup);
0048 
0049   virtual void beginJob();
0050   virtual void endJob();
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 
0076 using namespace std;
0077 using namespace edm;
0078 
0079 /// Constructor
0080 GLBMuonAnalyzer::GLBMuonAnalyzer(const ParameterSet &pset) {
0081   theGLBMuonLabel = pset.getUntrackedParameter<edm::InputTag>("GlobalTrackCollectionLabel");
0082 
0083   theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0084 
0085   theDataType = pset.getUntrackedParameter<string>("DataType");
0086 
0087   if (theDataType != "RealData" && theDataType != "SimData")
0088     LogTrace("Analyzer") << "Error in Data Type!!" << endl;
0089 
0090   numberOfSimTracks = 0;
0091   numberOfRecTracks = 0;
0092 }
0093 
0094 /// Destructor
0095 GLBMuonAnalyzer::~GLBMuonAnalyzer() {}
0096 
0097 void GLBMuonAnalyzer::beginJob() {
0098   // Create the root file
0099   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0100   theFile->cd();
0101 
0102   hPtRec = new TH1F("pTRec", "p_{T}^{rec}", 250, 0, 120);
0103   hPtSim = new TH1F("pTSim", "p_{T}^{gen} ", 250, 0, 120);
0104 
0105   hPTDiff = new TH1F("pTDiff", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0106   hPTDiff2 = new TH1F("pTDiff2", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0107 
0108   hPTDiffvsEta = new TH2F("PTDiffvsEta", "p_{T}^{rec} - p_{T}^{gen} VS #eta", 100, -2.5, 2.5, 250, -120, 120);
0109   hPTDiffvsPhi = new TH2F("PTDiffvsPhi", "p_{T}^{rec} - p_{T}^{gen} VS #phi", 100, -6, 6, 250, -120, 120);
0110 
0111   hPres = new TH1F("pTRes", "pT Resolution", 100, -2, 2);
0112   h1_Pres = new TH1F("invPTRes", "1/pT Resolution", 100, -2, 2);
0113 }
0114 
0115 void GLBMuonAnalyzer::endJob() {
0116   if (theDataType == "SimData") {
0117     LogTrace("Analyzer") << endl << endl << "Number of Sim tracks: " << numberOfSimTracks << endl;
0118   }
0119 
0120   LogTrace("Analyzer") << "Number of Reco tracks: " << numberOfRecTracks << endl << endl;
0121 
0122   // Write the histos to file
0123   theFile->cd();
0124   hPtRec->Write();
0125   hPtSim->Write();
0126   hPres->Write();
0127   h1_Pres->Write();
0128   hPTDiff->Write();
0129   hPTDiff2->Write();
0130   hPTDiffvsEta->Write();
0131   hPTDiffvsPhi->Write();
0132   theFile->Close();
0133 }
0134 
0135 void GLBMuonAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
0136   //LogTrace("Analyzer") << "Run: " << event.id().run() << " Event: " << event.id().event() << endl;
0137   MuonPatternRecoDumper debug;
0138 
0139   // Get the RecTrack collection from the event
0140   Handle<reco::TrackCollection> staTracks;
0141   event.getByLabel(theGLBMuonLabel, staTracks);
0142 
0143   ESHandle<MagneticField> theMGField;
0144   eventSetup.get<IdealMagneticFieldRecord>().get(theMGField);
0145 
0146   ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
0147   eventSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
0148 
0149   double recPt = 0.;
0150   double simPt = 0.;
0151 
0152   // Get the SimTrack collection from the event
0153   if (theDataType == "SimData") {
0154     Handle<SimTrackContainer> simTracks;
0155     event.getByLabel("g4SimHits", simTracks);
0156 
0157     numberOfRecTracks += staTracks->size();
0158 
0159     SimTrackContainer::const_iterator simTrack;
0160 
0161     //LogTrace("Analyzer")<<"Simulated tracks: "<<endl;
0162     for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0163       if (abs((*simTrack).type()) == 13) {
0164         //LogTrace("Analyzer")<<"Sim pT: "<<(*simTrack).momentum().perp()<<endl;
0165         simPt = (*simTrack).momentum().pt();
0166         //LogTrace("Analyzer")<<"Sim Eta: "<<(*simTrack).momentum().eta()<<endl;
0167         numberOfSimTracks++;
0168       }
0169     }
0170     LogTrace("Analyzer") << endl;
0171   }
0172 
0173   reco::TrackCollection::const_iterator staTrack;
0174 
0175   //LogTrace("Analyzer")<<"Reconstructed tracks: " << staTracks->size() << endl;
0176 
0177   for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack) {
0178     reco::TransientTrack track(*staTrack, &*theMGField, theTrackingGeometry);
0179 
0180     //LogTrace("Analyzer") << debug.dumpFTS(track.impactPointTSCP().theState());
0181 
0182     recPt = track.impactPointTSCP().momentum().perp();
0183     //LogTrace("Analyzer")<<" p: "<<track.impactPointTSCP().momentum().mag()<< " pT: "<<recPt<<endl;
0184     //LogTrace("Analyzer")<<" chi2: "<<track.chi2()<<endl;
0185 
0186     hPtRec->Fill(recPt);
0187 
0188     TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0189     //LogTrace("Analyzer") << "Inner TSOS:"<<endl;
0190     //LogTrace("Analyzer") << debug.dumpTSOS(innerTSOS);
0191     //LogTrace("Analyzer")<<" p: "<<innerTSOS.globalMomentum().mag()<< " pT: "<<innerTSOS.globalMomentum().perp()<<endl;
0192 
0193     trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
0194     trackingRecHit_iterator rhend = staTrack->recHitsEnd();
0195 
0196     int muHit = 0;
0197     int tkHit = 0;
0198 
0199     //LogTrace("Analyzer")<<"RecHits:"<<endl;
0200     for (trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit) {
0201       //      const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
0202       //      double r = geomDet->surface().position().perp();
0203       //      double z = geomDet->toGlobal((*recHit)->localPosition()).z();
0204       //LogTrace("Analyzer")<<"r: "<< r <<" z: "<<z <<endl;
0205       if ((*recHit)->geographicalId().det() == DetId::Muon)
0206         ++muHit;
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);