Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:12

0001 /** \class STAMuonAnalyzer
0002  *  Analyzer of the StandAlone muon tracks
0003  *
0004  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0005  */
0006 
0007 #include "RecoMuon/StandAloneMuonProducer/test/STAMuonAnalyzer.h"
0008 
0009 // Collaborating Class Header
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 /// Constructor
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 
0051 /// Destructor
0052 STAMuonAnalyzer::~STAMuonAnalyzer() {}
0053 
0054 void STAMuonAnalyzer::beginJob() {
0055   // Create the root file
0056   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0057   theFile->cd();
0058 
0059   hPtRec = new TH1F("pTRec", "p_{T}^{rec}", 250, 0, 120);
0060   hPtSim = new TH1F("pTSim", "p_{T}^{gen} ", 250, 0, 120);
0061 
0062   hPTDiff = new TH1F("pTDiff", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0063   hPTDiff2 = new TH1F("pTDiff2", "p_{T}^{rec} - p_{T}^{gen} ", 250, -120, 120);
0064 
0065   hPTDiffvsEta = new TH2F("PTDiffvsEta", "p_{T}^{rec} - p_{T}^{gen} VS #eta", 100, -2.5, 2.5, 250, -120, 120);
0066   hPTDiffvsPhi = new TH2F("PTDiffvsPhi", "p_{T}^{rec} - p_{T}^{gen} VS #phi", 100, -6, 6, 250, -120, 120);
0067 
0068   hPres = new TH1F("pTRes", "pT Resolution", 100, -2, 2);
0069   h1_Pres = new TH1F("invPTRes", "1/pT Resolution", 100, -2, 2);
0070 }
0071 
0072 void STAMuonAnalyzer::endJob() {
0073   if (theDataType == "SimData") {
0074     cout << endl << endl << "Number of Sim tracks: " << numberOfSimTracks << endl;
0075   }
0076 
0077   cout << "Number of Reco tracks: " << numberOfRecTracks << endl << endl;
0078 
0079   // Write the histos to file
0080   theFile->cd();
0081   hPtRec->Write();
0082   hPtSim->Write();
0083   hPres->Write();
0084   h1_Pres->Write();
0085   hPTDiff->Write();
0086   hPTDiff2->Write();
0087   hPTDiffvsEta->Write();
0088   hPTDiffvsPhi->Write();
0089   theFile->Close();
0090 }
0091 
0092 void STAMuonAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0093   cout << "Run: " << event.id().run() << " Event: " << event.id().event() << endl;
0094   MuonPatternRecoDumper debug;
0095 
0096   // Get the RecTrack collection from the event
0097   Handle<reco::TrackCollection> staTracks;
0098   event.getByLabel(theSTAMuonLabel, staTracks);
0099 
0100   ESHandle<MagneticField> theMGField;
0101   eventSetup.get<IdealMagneticFieldRecord>().get(theMGField);
0102 
0103   ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
0104   eventSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
0105 
0106   double recPt = 0.;
0107   double simPt = 0.;
0108 
0109   // Get the SimTrack collection from the event
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       //      hPTDiff2->Fill(track.innermostMeasurementState().globalMomentum().perp()-simPt);
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);