StandAloneMuonProducer/test/STAMuonAnalyzer.cc

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   theFieldToken = esConsumes();
0051   theGeomToken = esConsumes();
0052 }
0053
0054 /// Destructor
0055 STAMuonAnalyzer::~STAMuonAnalyzer() {}
0056
0057 void STAMuonAnalyzer::beginJob() {
0058   // Create the root file
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   // Write the histos to file
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   // Get the RecTrack collection from the event
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   // 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);