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