Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:44

0001 /** \class CosmicMuonValidator
0002  *
0003  *  compare reconstructed and simulated cosmic muons
0004  *
0005  *  the validator assumes single muon events
0006  *
0007  *  \author Chang Liu   -  Purdue University <Chang.Liu@cern.ch>
0008  */
0009 
0010 #include <memory>
0011 
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/EDAnalyzer.h"
0015 
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 
0021 #include "SimDataFormats/Track/interface/SimTrack.h"
0022 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0023 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0024 
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0029 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0030 #include "MagneticField/Engine/interface/MagneticField.h"
0031 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0032 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0033 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0034 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0035 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0036 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0037 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0038 #include "FWCore/ServiceRegistry/interface/Service.h"
0039 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0040 
0041 #include <iostream>
0042 
0043 #include <TH1D.h>
0044 #include <TSystem.h>
0045 #include <TH1.h>
0046 #include <TH2.h>
0047 #include <TLegend.h>
0048 #include <TStyle.h>
0049 #include <TCanvas.h>
0050 #include <TGraph.h>
0051 #include <TFrame.h>
0052 #include <TMath.h>
0053 #include <TF1.h>
0054 #include <TPostScript.h>
0055 #include <TPad.h>
0056 #include <TText.h>
0057 #include <TLatex.h>
0058 
0059 using namespace std;
0060 using namespace edm;
0061 
0062 class CosmicMuonValidator : public edm::EDAnalyzer {
0063 public:
0064   explicit CosmicMuonValidator(const edm::ParameterSet&);
0065   ~CosmicMuonValidator();
0066 
0067 private:
0068   virtual void beginJob();
0069 
0070   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0071 
0072   virtual void endJob();
0073 
0074   reco::Track bestTrack(const reco::TrackCollection&) const;
0075 
0076   PSimHitContainer matchedHit(const GlobalPoint&, const PSimHitContainer&) const;
0077 
0078   TrajectoryStateOnSurface updatedState(const TrajectoryStateOnSurface&, const PSimHit&) const;
0079 
0080   edm::ESHandle<Propagator> propagator() const;
0081 
0082   edm::InputTag trackLabel_;
0083   edm::InputTag simTrackLabel_;
0084 
0085   MuonServiceProxy* theService;
0086 
0087   int theDrawOption;
0088 
0089   int nEvent;
0090   int successR;
0091   int nNoSignal;
0092 
0093   TH2F* h_innerPosXY;
0094   TH2F* h_innerPosEP;
0095 
0096   TH2F* h_outerPosXY;
0097   TH2F* h_outerPosEP;
0098 
0099   TH1F* h_res;
0100   TH1F* h_theta;
0101   TH1F* h_phi;
0102   TH1F* h_pt_rec_sim;
0103   TH1F* h_phi_rec_sim;
0104   TH1F* h_theta_rec_sim;
0105   TH1F* h_Pres_inv_sim;
0106   TH1F* h_pt_sim;
0107   TH1F* h_pt_rec;
0108   TH1F* hnhit;
0109 
0110   TH1F* htotal4D;
0111   TH1F* htotalSeg;
0112 };
0113 
0114 CosmicMuonValidator::CosmicMuonValidator(const edm::ParameterSet& iConfig) {
0115   trackLabel_ = iConfig.getParameter<edm::InputTag>("TrackLabel");
0116   simTrackLabel_ = iConfig.getParameter<edm::InputTag>("SimTrackLabel");
0117   theDrawOption = iConfig.getUntrackedParameter<int>("DrawOption", 1);
0118 
0119   // service parameters
0120   edm::ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0121   // the services
0122   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0123 
0124   nEvent = 0;
0125   successR = 0;
0126   nNoSignal = 0;
0127 }
0128 
0129 CosmicMuonValidator::~CosmicMuonValidator() {
0130   if (theService)
0131     delete theService;
0132 }
0133 
0134 void CosmicMuonValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0135   theService->update(iSetup);
0136 
0137   nEvent++;
0138   std::cout << "reading event " << nEvent << std::endl;
0139 
0140   Handle<reco::TrackCollection> muons;
0141   iEvent.getByLabel(trackLabel_, muons);
0142   cout << "cosmic Muon: " << muons->size() << endl;
0143 
0144   //  if (muons->empty()) return;
0145   if (!muons->empty())
0146     successR++;
0147 
0148   float ptsim = 0;
0149   float simC = 0;
0150   float thetasim = 0;
0151   float phisim = 0;
0152 
0153   Handle<edm::SimTrackContainer> simTracks;
0154   iEvent.getByLabel(simTrackLabel_, simTracks);
0155   cout << "simTracks: " << simTracks->size() << endl;
0156 
0157   for (SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0158     if (abs((*simTrack).type()) == 13) {
0159       cout << "MC Muon: mom (" << simTrack->momentum().x() << "," << simTrack->momentum().y() << ","
0160            << simTrack->momentum().z() << ")" << endl;
0161       thetasim = simTrack->momentum().theta();
0162       phisim = simTrack->momentum().phi();
0163 
0164       simC = -simTrack->type() / 13.;
0165       ptsim = simTrack->momentum().pt();
0166     }
0167   }
0168 
0169   // Get Segment collections from the Event
0170   edm::Handle<DTRecSegment4DCollection> dtSegments;
0171   iEvent.getByLabel("dt4DSegments", dtSegments);
0172 
0173   edm::Handle<CSCSegmentCollection> cscSegments;
0174   iEvent.getByLabel("cscSegments", cscSegments);
0175 
0176   int nSeg = dtSegments->size() + cscSegments->size();
0177   cout << "cscSegments: " << cscSegments->size() << endl;
0178 
0179   htotalSeg->Fill(nSeg);
0180 
0181   int n4D = 0;
0182 
0183   for (DTRecSegment4DCollection::const_iterator idt = dtSegments->begin(); idt != dtSegments->end(); ++idt) {
0184     if (idt->dimension() < 4)
0185       continue;
0186     bool sameChamber = false;
0187     DetId thisId = idt->geographicalId();
0188 
0189     for (DTRecSegment4DCollection::const_iterator idt2 = dtSegments->begin(); idt2 != idt; ++idt2) {
0190       if (idt2->geographicalId() == thisId)
0191         sameChamber = true;
0192     }
0193     if (!sameChamber)
0194       n4D++;
0195   }
0196 
0197   for (CSCSegmentCollection::const_iterator icsc = cscSegments->begin(); icsc != cscSegments->end(); ++icsc) {
0198     if (icsc->dimension() < 4)
0199       continue;
0200     bool sameChamber = false;
0201     DetId thisId = icsc->geographicalId();
0202 
0203     for (CSCSegmentCollection::const_iterator icsc2 = cscSegments->begin(); icsc2 != icsc; ++icsc2) {
0204       if (icsc2->geographicalId() == thisId)
0205         sameChamber = true;
0206     }
0207 
0208     if (!sameChamber)
0209       n4D++;
0210   }
0211 
0212   htotal4D->Fill(n4D);
0213   if (n4D < 2)
0214     nNoSignal++;
0215 
0216   if (muons->empty())
0217     return;
0218   reco::Track muon = bestTrack(*muons);
0219 
0220   cout << "cosmic Muon Track: "
0221        << " mom: (" << muon.px() << "," << muon.py() << "," << muon.pz() << ")" << endl;
0222 
0223   if (fabs(muon.p()) < 1e-5)
0224     return;  //prevent those failed to extrapolation to vertex
0225 
0226   math::XYZVector innerMo = muon.innerMomentum();
0227 
0228   float ptreco = muon.pt();
0229   int qreco = muon.charge();
0230 
0231   h_pt_rec->Fill(ptreco);
0232   hnhit->Fill(muon.recHitsSize());
0233 
0234   cout << "resolution " << (qreco / ptreco - simC / ptsim) * ptsim / simC << endl;
0235 
0236   h_res->Fill((qreco / ptreco - simC / ptsim) * ptsim / simC);
0237 
0238   GlobalVector im(innerMo.x(), innerMo.y(), innerMo.z());
0239   float thetareco = im.theta();
0240   float phireco = im.phi();
0241   h_theta->Fill(Geom::Theta<float>(thetareco - thetasim));
0242   h_phi->Fill(Geom::Phi<float>(phireco - phisim));
0243 
0244   math::XYZPoint innerPo = muon.innerPosition();
0245   GlobalPoint ip(innerPo.x(), innerPo.y(), innerPo.z());
0246 
0247   h_innerPosXY->Fill(ip.x(), ip.y());
0248   h_innerPosEP->Fill(ip.eta(), Geom::Phi<float>(ip.phi()));
0249 
0250   math::XYZPoint outerPo = muon.outerPosition();
0251   GlobalPoint op(outerPo.x(), outerPo.y(), outerPo.z());
0252 
0253   h_outerPosXY->Fill(op.x(), op.y());
0254   h_outerPosEP->Fill(op.eta(), Geom::Phi<float>(op.phi()));
0255 
0256   //Now compare innermost state with associated sim hit
0257 
0258   Handle<PSimHitContainer> dtSimHits;
0259   Handle<PSimHitContainer> cscSimHits;
0260   Handle<PSimHitContainer> rpcSimHits;
0261 
0262   iEvent.getByLabel("g4SimHits", "MuonDTHits", dtSimHits);
0263   iEvent.getByLabel("g4SimHits", "MuonCSCHits", cscSimHits);
0264   iEvent.getByLabel("g4SimHits", "MuonRPCHits", rpcSimHits);
0265 
0266   cout << "DT simHits collections: " << dtSimHits->size() << endl;
0267   cout << "CSC simHits collections: " << cscSimHits->size() << endl;
0268   cout << "RPC simHits collections: " << rpcSimHits->size() << endl;
0269 
0270   PSimHitContainer allSimHits = *dtSimHits;
0271   allSimHits.insert(allSimHits.end(), (cscSimHits)->begin(), (cscSimHits)->end());
0272   allSimHits.insert(allSimHits.end(), (rpcSimHits)->begin(), (rpcSimHits)->end());
0273 
0274   cout << "allSimHits " << allSimHits.size() << endl;
0275 
0276   if (allSimHits.empty())
0277     return;
0278 
0279   PSimHitContainer msimh = matchedHit(ip, allSimHits);
0280 
0281   if (!msimh.empty()) {
0282     DetId idSim(msimh.front().detUnitId());
0283 
0284     GlobalVector simmom =
0285         theService->trackingGeometry()->idToDet(idSim)->surface().toGlobal(msimh.front().momentumAtEntry());
0286 
0287     TrajectoryStateOnSurface innerTSOS = trajectoryStateTransform::innerStateOnSurface(
0288         muon, *theService->trackingGeometry(), &*theService->magneticField());
0289 
0290     TrajectoryStateOnSurface stateAH = updatedState(innerTSOS, msimh.front());
0291     if (!stateAH.isValid())
0292       return;
0293     im = stateAH.globalMomentum();
0294 
0295     cout << "sim Momentum: " << simmom << endl;
0296     cout << "track Mom here: " << im << endl;
0297 
0298     h_pt_rec_sim->Fill(((double)simmom.perp()) - im.perp());
0299     h_phi_rec_sim->Fill(((Geom::Phi<float>(simmom.phi())) - Geom::Phi<float>(im.phi())) * 180 / acos(-1.));
0300 
0301     h_Pres_inv_sim->Fill((1 / im.perp() - 1 / ((double)simmom.perp())) / (1 / ((double)simmom.perp())));
0302 
0303     h_pt_sim->Fill(((double)simmom.perp()));
0304 
0305     h_theta_rec_sim->Fill((((double)simmom.theta()) - im.theta()) * 180 / acos(-1.));
0306   }
0307 }
0308 
0309 void CosmicMuonValidator::beginJob() {
0310   cout << "Prepare histograms "
0311        << "\n";
0312   edm::Service<TFileService> fs;
0313 
0314   h_res = fs->make<TH1F>("h_res", "resolution of P_{T}", 50, -5.0, 5.0);
0315 
0316   h_theta = fs->make<TH1F>("h_theta", "theta angle ", 50, -0.1, 0.1);
0317   h_phi = fs->make<TH1F>("h_phi", "phi angle ", 50, -2.0, 2.0);
0318 
0319   hnhit = fs->make<TH1F>("hnhit", "Number of Hits in Track by Cos", 60, 0.0, 60.0);
0320 
0321   h_innerPosXY = fs->make<TH2F>("h_innerPosXY", "inner x-y", 100, -700.0, 700.0, 100, -700.0, 700.0);
0322   h_innerPosEP = fs->make<TH2F>("h_innerPosEP", "inner #eta-#phi", 100, -2.4, 2.4, 100, -3.3, 3.3);
0323 
0324   h_outerPosXY = fs->make<TH2F>("h_outerPosXY", "outer x-y", 100, -700.0, 700.0, 100, -700.0, 700.0);
0325   h_outerPosEP = fs->make<TH2F>("h_outerPosEP", "outer #eta-#phi", 100, -2.4, 2.4, 100, -3.3, 3.3);
0326 
0327   h_pt_rec_sim = fs->make<TH1F>("h_pt_res_sim", "diff of P_{T} at SimHit", 50, -2.0, 2.0);
0328   h_phi_rec_sim = fs->make<TH1F>("h_phi_res_sim", "diff of #phi at SimHit", 50, -2.0, 2.0);
0329   h_theta_rec_sim = fs->make<TH1F>("h_theta_res_sim", "diff of #theta at SimHit", 50, -2.0, 2.0);
0330   h_Pres_inv_sim = fs->make<TH1F>("h_Pres_inv_sim", "resolution of P_{T} at SimHit", 70, -1.0, 1.0);
0331 
0332   h_pt_sim = fs->make<TH1F>("h_pt_sim", "distribution of P_{T} at SimHit", 100, 0.0, 100.0);
0333   h_pt_rec = fs->make<TH1F>("h_pt_rec", "distribution of P_{T} at SimHit", 100, 0.0, 100.0);
0334   htotal4D = fs->make<TH1F>("htotal4D", "# of Segments", 15, 0.0, 15.0);
0335   htotalSeg = fs->make<TH1F>("htotalSeg", "# of Segments", 15, 0.0, 15.0);
0336 }
0337 
0338 void CosmicMuonValidator::endJob() {
0339   float eff = (float)successR / ((float)nEvent - (float)nNoSignal) * 100;
0340 
0341   std::cout << "++++++++++++++++++++++++++++++++++++++++" << std::endl;
0342   std::cout << "Analyzed " << nEvent << " events, " << std::endl;
0343   std::cout << successR << " events are successfully reconstructed. " << std::endl;
0344   std::cout << nNoSignal << " events do not have good enough signals. " << std::endl;
0345   std::cout << "Reconstruction efficiency is approximately " << eff << "%. " << std::endl;
0346   std::cout << "++++++++++++++++++++++++++++++++++++++++" << std::endl;
0347 
0348   h_innerPosXY->SetXTitle("X");
0349   h_innerPosXY->SetYTitle("Y");
0350   h_innerPosXY->SetMarkerColor(2);
0351   h_innerPosXY->SetMarkerStyle(24);
0352 
0353   h_outerPosXY->SetXTitle("X");
0354   h_outerPosXY->SetYTitle("Y");
0355   h_outerPosXY->SetMarkerColor(2);
0356   h_outerPosXY->SetMarkerStyle(24);
0357 
0358   h_innerPosEP->SetXTitle("#eta");
0359   h_innerPosEP->SetYTitle("#phi");
0360   h_innerPosEP->SetMarkerColor(2);
0361   h_innerPosEP->SetMarkerStyle(24);
0362 
0363   h_outerPosEP->SetXTitle("#eta");
0364   h_outerPosEP->SetYTitle("#phi");
0365   h_outerPosEP->SetMarkerColor(2);
0366   h_outerPosEP->SetMarkerStyle(24);
0367 
0368   h_res->SetXTitle("(q^{reco}/P^{reco}_{T}-q^{simTrack}/P^{simTrack}_{T})/(q^{simTrack}/P^{simTrack}_{T})");
0369   h_res->SetLineWidth(2);
0370   h_res->SetLineColor(2);
0371   h_res->SetLineStyle(1);
0372 
0373   h_theta->SetXTitle("#theta^{reco}-#theta^{simTrack}");
0374   h_theta->SetLineWidth(2);
0375   h_theta->SetLineColor(2);
0376   h_theta->SetLineStyle(1);
0377 
0378   h_phi->SetXTitle("#phi^{reco}-#phi^{simTrack}");
0379   h_phi->SetLineWidth(2);
0380   h_phi->SetLineColor(2);
0381   h_phi->SetLineStyle(1);
0382 
0383   h_pt_rec_sim->SetXTitle("P^{simHit}_{T}-P^{reco}_{T}");
0384   h_pt_rec_sim->SetLineWidth(2);
0385   h_pt_rec_sim->SetLineColor(2);
0386   h_pt_rec_sim->SetLineStyle(1);
0387 
0388   h_phi_rec_sim->SetXTitle("#phi^{simHit}-#phi^{reco}");
0389   h_phi_rec_sim->SetLineWidth(2);
0390   h_phi_rec_sim->SetLineColor(2);
0391   h_phi_rec_sim->SetLineStyle(1);
0392 
0393   h_theta_rec_sim->SetXTitle("#theta^{simHit}-#theta^{reco}");
0394   h_theta_rec_sim->SetLineWidth(2);
0395   h_theta_rec_sim->SetLineColor(2);
0396   h_theta_rec_sim->SetLineStyle(1);
0397 
0398   h_Pres_inv_sim->SetXTitle("(q^{reco}/P^{reco}_{T}-q^{simHit}/P^{simHit}_{T})/(q^{simHit}/P^{simHit}_{T})");
0399   h_Pres_inv_sim->SetLineWidth(2);
0400   h_Pres_inv_sim->SetLineColor(2);
0401   h_Pres_inv_sim->SetLineStyle(1);
0402 
0403   hnhit->SetXTitle("N_{RecHits}");
0404   hnhit->SetLineWidth(2);
0405   hnhit->SetLineColor(2);
0406   hnhit->SetLineStyle(1);
0407 
0408   h_pt_sim->SetXTitle("P_{t}");
0409   h_pt_sim->SetLineWidth(2);
0410   h_pt_sim->SetLineColor(4);
0411   h_pt_rec->SetLineWidth(2);
0412   h_pt_rec->SetLineColor(2);
0413   h_pt_sim->SetLineStyle(1);
0414   h_pt_rec->SetLineStyle(2);
0415 
0416   htotal4D->SetXTitle("No. of Segments");
0417   htotal4D->SetLineWidth(2);
0418   htotal4D->SetLineColor(2);
0419 
0420   htotalSeg->SetLineWidth(2);
0421   htotalSeg->SetLineColor(4);
0422 
0423   htotal4D->SetLineStyle(1);
0424   htotalSeg->SetLineStyle(2);
0425 
0426   int theDrawOption = 1;
0427 
0428   if (theDrawOption == 0) {
0429     TCanvas* c2 = new TCanvas("innerTSOSXY", "XY", 10, 10, 800, 600);
0430     c2->SetFillColor(0);
0431     c2->SetGrid(1);
0432     c2->SetRightMargin(0.03);
0433     c2->SetTopMargin(0.02);
0434     c2->cd();
0435     h_innerPosXY->Draw("SCAT");
0436     c2->Update();
0437     c2->Write();
0438 
0439     TCanvas* c2a = new TCanvas("outerTSOSXY", "Outer XY", 10, 10, 800, 600);
0440     c2a->SetFillColor(0);
0441     c2a->SetGrid(1);
0442     c2a->SetRightMargin(0.03);
0443     c2a->SetTopMargin(0.02);
0444     c2a->cd();
0445     h_outerPosXY->Draw("SCAT");
0446     c2a->Update();
0447     c2a->Write();
0448 
0449     TCanvas* c3 = new TCanvas("innerEtaPhi", "Inner #eta #phi", 10, 10, 800, 600);
0450     c3->SetFillColor(0);
0451     c3->SetGrid(1);
0452     c3->SetRightMargin(0.03);
0453     c3->SetTopMargin(0.02);
0454     c3->cd();
0455     h_innerPosEP->Draw("SCAT");
0456     c3->Update();
0457     c3->Write();
0458 
0459     TCanvas* c3a = new TCanvas("outerEtaPhi", "Outer #eta #phi", 10, 10, 800, 600);
0460     c3a->SetFillColor(0);
0461     c3a->SetGrid(1);
0462     c3a->SetRightMargin(0.03);
0463     c3a->SetTopMargin(0.02);
0464     c3a->cd();
0465     h_outerPosEP->Draw("SCAT");
0466     c3a->Update();
0467     c3a->Write();
0468 
0469     TCanvas* cRes1 = new TCanvas("TrackPtRes", "Resolution of P_{T} wrt SimTrack", 10, 10, 800, 600);
0470     cRes1->SetFillColor(0);
0471     cRes1->SetGrid(1);
0472     cRes1->SetRightMargin(0.03);
0473     cRes1->SetTopMargin(0.02);
0474     cRes1->cd();
0475     h_res->DrawCopy("HE");
0476     cRes1->Update();
0477     cRes1->Write();
0478 
0479     TCanvas* cRes2 = new TCanvas("TrackTheta", "Resolution of Theta wrt SimTrack", 10, 10, 800, 600);
0480     cRes2->SetFillColor(0);
0481     cRes2->SetGrid(1);
0482     cRes2->SetRightMargin(0.03);
0483     cRes2->SetTopMargin(0.02);
0484     cRes2->cd();
0485     h_theta->DrawCopy("HE");
0486     cRes2->Update();
0487     cRes2->Write();
0488 
0489     TCanvas* cRes3 = new TCanvas("TrackPhi", "Resolution of phi wrt SimTrack", 10, 10, 800, 600);
0490     cRes3->SetFillColor(0);
0491     cRes3->SetGrid(1);
0492     cRes3->SetRightMargin(0.03);
0493     cRes3->SetTopMargin(0.02);
0494     cRes3->cd();
0495     h_phi->DrawCopy("HE");
0496     cRes3->Update();
0497     cRes3->Write();
0498 
0499     TCanvas* cRes4 = new TCanvas("inTsosPtDiff", "Resolution of P_{T} at SimHit", 10, 10, 800, 600);
0500     cRes4->SetFillColor(0);
0501     cRes4->SetGrid(1);
0502     cRes4->SetRightMargin(0.03);
0503     cRes4->SetTopMargin(0.02);
0504     cRes4->cd();
0505     h_pt_rec_sim->DrawCopy("HE");
0506     cRes4->Update();
0507     cRes4->Write();
0508 
0509     TCanvas* cRes5 = new TCanvas("inTsosPhi", "Resolution of Phi at SimHit", 10, 10, 800, 600);
0510     cRes5->SetFillColor(0);
0511     cRes5->SetGrid(1);
0512     cRes5->SetRightMargin(0.03);
0513     cRes5->SetTopMargin(0.02);
0514     cRes5->cd();
0515     h_phi_rec_sim->DrawCopy("HE");
0516     cRes5->Update();
0517     cRes5->Write();
0518 
0519     TCanvas* cRes6 = new TCanvas("inTsosTheta", "Resolution of #theta at SimHit", 10, 10, 800, 600);
0520     cRes6->SetFillColor(0);
0521     cRes6->SetGrid(1);
0522     cRes6->SetRightMargin(0.03);
0523     cRes6->SetTopMargin(0.02);
0524     cRes6->cd();
0525     h_theta_rec_sim->DrawCopy("HE");
0526     cRes6->Update();
0527     cRes6->Write();
0528 
0529     TCanvas* cRes7 = new TCanvas("inTsosPtRes", "Resolution of P_{T} at SimHit", 10, 10, 800, 600);
0530     cRes7->SetFillColor(0);
0531     cRes7->SetGrid(1);
0532     cRes7->SetRightMargin(0.03);
0533     cRes7->SetTopMargin(0.02);
0534     cRes7->cd();
0535     h_Pres_inv_sim->DrawCopy("HE");
0536     TF1* g2 = new TF1("g2", "gaus", -0.5, 0.5);
0537     g2->SetLineColor(4);
0538     h_Pres_inv_sim->Fit("g2", "R");
0539 
0540     cRes7->Update();
0541     cRes7->Write();
0542 
0543     TCanvas* c7 = new TCanvas("nHits", "Number of RecHits in Track", 10, 10, 800, 600);
0544     c7->SetFillColor(0);
0545     c7->SetGrid(1);
0546     c7->SetRightMargin(0.03);
0547     c7->SetTopMargin(0.02);
0548     c7->cd();
0549     hnhit->DrawCopy("HE");
0550     c7->Update();
0551     c7->Write();
0552 
0553     TCanvas* cRes8 = new TCanvas("PtDis", "Distribution of P_{T} at SimHit", 10, 10, 800, 600);
0554     cRes8->SetFillColor(0);
0555     cRes8->SetGrid(1);
0556     cRes8->SetRightMargin(0.03);
0557     cRes8->SetTopMargin(0.02);
0558     cRes8->cd();
0559 
0560     h_pt_sim->DrawCopy("HE");
0561     h_pt_rec->DrawCopy("HEsame");
0562     TLegend* legend8 = new TLegend(0.7, 0.6, 0.9, 0.8);
0563     legend8->SetTextAlign(32);
0564     legend8->SetTextColor(1);
0565     legend8->SetTextSize(0.04);
0566     legend8->AddEntry("h_pt_sim", "By Sim", "l");
0567     legend8->AddEntry("h_pt_rec", "By Cos", "l");
0568     legend8->Draw();
0569     cRes8->Update();
0570     cRes8->Write();
0571 
0572     TCanvas* csegs = new TCanvas("csegs", "Total & 4D Segments", 10, 10, 800, 600);
0573 
0574     csegs->SetFillColor(0);
0575     csegs->SetGrid(1);
0576     csegs->SetRightMargin(0.03);
0577     csegs->SetTopMargin(0.02);
0578     csegs->cd();
0579 
0580     htotal4D->DrawCopy("HE");
0581     htotalSeg->DrawCopy("HEsame");
0582 
0583     TLegend* legendseg = new TLegend(0.6, 0.2, 0.9, 0.4);
0584     legendseg->SetTextAlign(32);
0585     legendseg->SetTextColor(1);
0586     legendseg->SetTextSize(0.04);
0587     legendseg->AddEntry("htotal4D", "4D Segments", "l");
0588     legendseg->AddEntry("htotalSeg", "total Segments", "l");
0589     legendseg->Draw();
0590 
0591     csegs->Update();
0592     csegs->Write();
0593 
0594   } else {
0595     TCanvas* cpdf = new TCanvas("cpdf", "", 0, 1, 500, 700);
0596     cpdf->SetTicks();
0597 
0598     TPostScript* pdf = new TPostScript("cosmicValidation.ps", 111);
0599 
0600     const int NUM_PAGES = 7;
0601     TPad* pad[NUM_PAGES];
0602     for (int i_page = 0; i_page < NUM_PAGES; i_page++)
0603       pad[i_page] = new TPad("", "", .05, .05, .95, .93);
0604 
0605     ostringstream page_print;
0606     int page = 0;
0607 
0608     TLatex ttl;
0609 
0610     pdf->NewPage();
0611     cpdf->Clear();
0612     cpdf->cd(0);
0613     ttl.DrawLatex(.4, .95, "inner and outer state (x,y) ");
0614     page_print.str("");
0615     page_print << page + 1;
0616     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0617     pad[page]->Draw();
0618     pad[page]->Divide(1, 2);
0619     pad[page]->cd(1);
0620     h_innerPosXY->Draw("SCAT");
0621     pad[page]->cd(2);
0622     h_outerPosXY->Draw("SCAT");
0623 
0624     page++;
0625     cpdf->Update();
0626 
0627     pdf->NewPage();
0628     cpdf->Clear();
0629     cpdf->cd(0);
0630     ttl.DrawLatex(.4, .95, "inner and outer state (eta, phi) ");
0631     page_print.str("");
0632     page_print << page + 1;
0633     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0634     pad[page]->Draw();
0635     pad[page]->Divide(1, 2);
0636     pad[page]->cd(1);
0637     h_innerPosEP->Draw("SCAT");
0638     pad[page]->cd(2);
0639     h_outerPosEP->Draw("SCAT");
0640 
0641     page++;
0642     cpdf->Update();
0643 
0644     pdf->NewPage();
0645     cpdf->Clear();
0646     cpdf->cd(0);
0647 
0648     ttl.DrawLatex(.4, .95, "resolution wrt simTrack (pt, theta, phi) ");
0649     page_print.str("");
0650     page_print << page + 1;
0651     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0652 
0653     pad[page]->Draw();
0654     pad[page]->Divide(1, 3);
0655     pad[page]->cd(1);
0656     h_res->Draw("HE");
0657 
0658     pad[page]->cd(2);
0659     h_res->Draw("HE");
0660 
0661     pad[page]->cd(3);
0662     h_phi->Draw("HE");
0663 
0664     page++;
0665     cpdf->Update();
0666 
0667     pdf->NewPage();
0668     cpdf->Clear();
0669     cpdf->cd(0);
0670 
0671     ttl.DrawLatex(.4, .95, "resolution wrt simHit at innermost (pt, theta, phi) ");
0672     page_print.str("");
0673     page_print << page + 1;
0674     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0675     pad[page]->Draw();
0676     pad[page]->Divide(1, 3);
0677     pad[page]->cd(1);
0678     h_Pres_inv_sim->Draw("HE");
0679     pad[page]->cd(2);
0680     h_theta_rec_sim->Draw("HE");
0681     pad[page]->cd(3);
0682     h_phi_rec_sim->Draw("HE");
0683 
0684     page++;
0685     cpdf->Update();
0686 
0687     pdf->NewPage();
0688     cpdf->Clear();
0689     cpdf->cd(0);
0690 
0691     ttl.DrawLatex(.4, .95, "number of hits ");
0692     page_print.str("");
0693     page_print << page + 1;
0694     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0695     pad[page]->Draw();
0696     pad[page]->cd();
0697     hnhit->Draw("HE");
0698 
0699     page++;
0700     cpdf->Update();
0701 
0702     pdf->NewPage();
0703     cpdf->Clear();
0704     cpdf->cd(0);
0705 
0706     ttl.DrawLatex(.4, .95, "pt distribution ");
0707     page_print.str("");
0708     page_print << page + 1;
0709     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0710     pad[page]->Draw();
0711     pad[page]->Divide(1, 2);
0712     pad[page]->cd(1);
0713     h_pt_rec->Draw("HE");
0714     pad[page]->cd(2);
0715     h_pt_sim->Draw("HE");
0716 
0717     page++;
0718     cpdf->Update();
0719 
0720     pdf->NewPage();
0721     cpdf->Clear();
0722     cpdf->cd(0);
0723 
0724     ttl.DrawLatex(.4, .95, "Number of hits ");
0725     page_print.str("");
0726     page_print << page + 1;
0727     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0728     pad[page]->Draw();
0729     pad[page]->Divide(1, 2);
0730     pad[page]->cd(1);
0731     htotal4D->Draw("HE");
0732     pad[page]->cd(2);
0733     htotalSeg->Draw("HE");
0734 
0735     pdf->Close();
0736   }
0737 }
0738 
0739 reco::Track CosmicMuonValidator::bestTrack(const reco::TrackCollection& muons) const {
0740   reco::Track bestOne = muons.front();
0741 
0742   for (reco::TrackCollection::const_iterator muon = muons.begin() + 1; muon != muons.end(); ++muon) {
0743     if (((*muon).found() > bestOne.found()) ||
0744         (((*muon).found() == bestOne.found()) && ((*muon).chi2() < bestOne.chi2())))
0745       bestOne = (*muon);
0746   }
0747   return bestOne;
0748 }
0749 
0750 PSimHitContainer CosmicMuonValidator::matchedHit(const GlobalPoint& tp, const PSimHitContainer& simHs) const {
0751   float dcut = 3.0;
0752   PSimHitContainer result;
0753   PSimHit rs = simHs.front();
0754   bool hasMatched = false;
0755 
0756   if (simHs.empty())
0757     return result;
0758 
0759   for (PSimHitContainer::const_iterator ish = simHs.begin(); ish != simHs.end(); ish++) {
0760     if (abs((*ish).particleType()) != 13)
0761       continue;
0762 
0763     DetId idsim((*ish).detUnitId());
0764 
0765     GlobalPoint sp = theService->trackingGeometry()->idToDet(idsim)->surface().toGlobal(
0766         ish->entryPoint());  //entryPoint or localPosition??
0767     GlobalVector dist = sp - tp;
0768     float d = fabs(dist.y());
0769     if (d < dcut) {
0770       rs = (*ish);
0771       dcut = d;
0772       hasMatched = true;
0773     }
0774   }
0775   if (hasMatched) {
0776     result.push_back(rs);
0777     DetId idsim(rs.detUnitId());
0778     cout << "selected simhit: " << theService->trackingGeometry()->idToDet(idsim)->surface().toGlobal(rs.entryPoint())
0779          << endl;
0780     cout << "matched with   : " << tp << endl;
0781   }
0782   return result;
0783 }
0784 
0785 TrajectoryStateOnSurface CosmicMuonValidator::updatedState(const TrajectoryStateOnSurface& tsos,
0786                                                            const PSimHit& hit) const {
0787   DetId idsim(hit.detUnitId());
0788 
0789   TrajectoryStateOnSurface result =
0790       propagator()->propagate(tsos, theService->trackingGeometry()->idToDet(idsim)->surface());
0791 
0792   return result;
0793 }
0794 
0795 edm::ESHandle<Propagator> CosmicMuonValidator::propagator() const {
0796   return theService->propagator("SteppingHelixPropagatorAny");
0797 }