File indexing completed on 2021-02-14 14:25:44
0001
0002
0003
0004
0005
0006 #include <memory>
0007
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/EDAnalyzer.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0020 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0021 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0022 #include "MagneticField/Engine/interface/MagneticField.h"
0023 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0024 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0027 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0028 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0029 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0032
0033 #include <iostream>
0034
0035 #include <TH1D.h>
0036 #include <TTree.h>
0037 #include <TROOT.h>
0038 #include <TSystem.h>
0039 #include <TFile.h>
0040 #include <TH1.h>
0041 #include <TH2.h>
0042 #include <TLegend.h>
0043 #include <TStyle.h>
0044 #include <TCanvas.h>
0045 #include <TGraph.h>
0046 #include <TFrame.h>
0047 #include <TMath.h>
0048 #include <TF1.h>
0049 #include <TPostScript.h>
0050 #include <TPad.h>
0051 #include <TText.h>
0052 #include <TLatex.h>
0053
0054 using namespace std;
0055 using namespace edm;
0056
0057 class RealCosmicDataAnalyzer : public edm::EDAnalyzer {
0058 public:
0059 explicit RealCosmicDataAnalyzer(const edm::ParameterSet&);
0060 ~RealCosmicDataAnalyzer();
0061
0062 private:
0063 virtual void beginJob();
0064
0065 virtual void analyze(const edm::Event&, const edm::EventSetup&);
0066
0067 virtual void endJob();
0068
0069 edm::ESHandle<Propagator> propagator() const;
0070
0071 edm::InputTag trackLabel_;
0072
0073 MuonServiceProxy* theService;
0074
0075 int theDrawOption;
0076
0077 int nEvent;
0078 int successR;
0079 int nNoSignal;
0080
0081 TH2F* h_innerPosXY;
0082 TH2F* h_innerPosEP;
0083
0084 TH2F* h_outerPosXY;
0085 TH2F* h_outerPosEP;
0086
0087 TH1F* h_inOutDis;
0088
0089 TH1F* h_theta;
0090 TH1F* h_phi;
0091 TH1F* h_pt_rec;
0092 TH1F* hnhit;
0093
0094 TH1F* htotal4D;
0095 TH1F* htotalSeg;
0096 };
0097
0098 RealCosmicDataAnalyzer::RealCosmicDataAnalyzer(const edm::ParameterSet& iConfig) {
0099 trackLabel_ = iConfig.getParameter<edm::InputTag>("TrackLabel");
0100 theDrawOption = iConfig.getUntrackedParameter<int>("DrawOption", 1);
0101
0102
0103 edm::ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0104
0105 theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0106
0107 nEvent = 0;
0108 successR = 0;
0109 nNoSignal = 0;
0110 }
0111
0112 RealCosmicDataAnalyzer::~RealCosmicDataAnalyzer() {
0113 if (theService)
0114 delete theService;
0115 }
0116
0117 void RealCosmicDataAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0118 theService->update(iSetup);
0119
0120 nEvent++;
0121 std::cout << "reading event " << nEvent << std::endl;
0122
0123 Handle<reco::TrackCollection> muons;
0124 iEvent.getByLabel(trackLabel_, muons);
0125 cout << "cosmic Muon: " << muons->size() << endl;
0126
0127
0128 if (!muons->empty())
0129 successR++;
0130
0131
0132 edm::Handle<DTRecSegment4DCollection> dtSegments;
0133 iEvent.getByLabel("dt4DSegments", dtSegments);
0134
0135 edm::Handle<CSCSegmentCollection> cscSegments;
0136 iEvent.getByLabel("cscSegments", cscSegments);
0137
0138 int nSeg = dtSegments->size() + cscSegments->size();
0139 cout << "cscSegments: " << cscSegments->size() << endl;
0140
0141 htotalSeg->Fill(nSeg);
0142
0143 int n4D = 0;
0144
0145 for (DTRecSegment4DCollection::const_iterator idt = dtSegments->begin(); idt != dtSegments->end(); ++idt) {
0146 if (idt->dimension() < 4)
0147 continue;
0148 bool sameChamber = false;
0149 DetId thisId = idt->geographicalId();
0150
0151 for (DTRecSegment4DCollection::const_iterator idt2 = dtSegments->begin(); idt2 != idt; ++idt2) {
0152 if (idt2->geographicalId() == thisId)
0153 sameChamber = true;
0154 }
0155 if (!sameChamber)
0156 n4D++;
0157 }
0158
0159 for (CSCSegmentCollection::const_iterator icsc = cscSegments->begin(); icsc != cscSegments->end(); ++icsc) {
0160 if (icsc->dimension() < 4)
0161 continue;
0162 bool sameChamber = false;
0163 DetId thisId = icsc->geographicalId();
0164
0165 for (CSCSegmentCollection::const_iterator icsc2 = cscSegments->begin(); icsc2 != icsc; ++icsc2) {
0166 if (icsc2->geographicalId() == thisId)
0167 sameChamber = true;
0168 }
0169
0170 if (!sameChamber)
0171 n4D++;
0172 }
0173
0174 htotal4D->Fill(n4D);
0175 if (n4D < 2)
0176 nNoSignal++;
0177
0178 if (muons->empty())
0179 return;
0180
0181 for (reco::TrackCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0182 cout << "cosmic Muon Track: "
0183 << " mom: (" << muon->px() << "," << muon->py() << "," << muon->pz() << ")" << endl;
0184 if (fabs(muon->p()) < 1e-5)
0185 return;
0186
0187 math::XYZVector innerMo = muon->innerMomentum();
0188
0189 float ptreco = muon->pt();
0190
0191 h_pt_rec->Fill(ptreco);
0192
0193 hnhit->Fill(muon->recHitsSize());
0194
0195 GlobalVector im(innerMo.x(), innerMo.y(), innerMo.z());
0196
0197 h_theta->Fill(im.theta());
0198 h_phi->Fill(im.phi());
0199
0200 math::XYZPoint innerPo = muon->innerPosition();
0201 GlobalPoint ip(innerPo.x(), innerPo.y(), innerPo.z());
0202
0203 h_innerPosXY->Fill(ip.x(), ip.y());
0204 h_innerPosEP->Fill(ip.eta(), Geom::Phi<float>(ip.phi()));
0205
0206 math::XYZPoint outerPo = muon->outerPosition();
0207 GlobalPoint op(outerPo.x(), outerPo.y(), outerPo.z());
0208
0209 h_outerPosXY->Fill(op.x(), op.y());
0210 h_outerPosEP->Fill(op.eta(), Geom::Phi<float>(op.phi()));
0211 h_inOutDis->Fill((ip - op).perp());
0212 }
0213 }
0214
0215 void RealCosmicDataAnalyzer::beginJob() {
0216 edm::Service<TFileService> fs;
0217
0218 h_theta = fs->make<TH1F>("h_theta", "theta angle ", 50, -0.1, 0.1);
0219 h_phi = fs->make<TH1F>("h_phi", "phi angle ", 50, -2.0, 2.0);
0220
0221 hnhit = fs->make<TH1F>("hnhit", "Number of Hits in Track by Cos", 60, 0.0, 60.0);
0222
0223 h_innerPosXY = fs->make<TH2F>("h_innerPosXY", "inner x-y", 100, -700.0, 700.0, 100, -700.0, 700.0);
0224 h_innerPosEP = fs->make<TH2F>("h_innerPosEP", "inner #eta-#phi", 100, -2.4, 2.4, 100, -3.3, 3.3);
0225
0226 h_outerPosXY = fs->make<TH2F>("h_outerPosXY", "outer x-y", 100, -700.0, 700.0, 100, -700.0, 700.0);
0227 h_outerPosEP = fs->make<TH2F>("h_outerPosEP", "outer #eta-#phi", 100, -2.4, 2.4, 100, -3.3, 3.3);
0228
0229 h_pt_rec = fs->make<TH1F>("h_pt_rec", "distribution of P_{T}", 100, 0.0, 100.0);
0230 h_inOutDis = fs->make<TH1F>("h_inOutDis", "distribution of inner outer distance", 200, 0.0, 1500.0);
0231 htotal4D = fs->make<TH1F>("htotal4D", "# of Segments", 15, 0.0, 15.0);
0232 htotalSeg = fs->make<TH1F>("htotalSeg", "# of Segments", 15, 0.0, 15.0);
0233 }
0234
0235 void RealCosmicDataAnalyzer::endJob() {
0236 h_innerPosXY->SetXTitle("X");
0237 h_innerPosXY->SetYTitle("Y");
0238 h_innerPosXY->SetMarkerColor(2);
0239 h_innerPosXY->SetMarkerStyle(24);
0240
0241 h_outerPosXY->SetXTitle("X");
0242 h_outerPosXY->SetYTitle("Y");
0243 h_outerPosXY->SetMarkerColor(2);
0244 h_outerPosXY->SetMarkerStyle(24);
0245
0246 h_innerPosEP->SetXTitle("#eta");
0247 h_innerPosEP->SetYTitle("#phi");
0248 h_innerPosEP->SetMarkerColor(2);
0249 h_innerPosEP->SetMarkerStyle(24);
0250
0251 h_outerPosEP->SetXTitle("#eta");
0252 h_outerPosEP->SetYTitle("#phi");
0253 h_outerPosEP->SetMarkerColor(2);
0254 h_outerPosEP->SetMarkerStyle(24);
0255
0256 hnhit->SetXTitle("N_{RecHits}");
0257 hnhit->SetLineWidth(2);
0258 hnhit->SetLineColor(2);
0259 hnhit->SetLineStyle(1);
0260
0261 h_pt_rec->SetLineWidth(2);
0262 h_pt_rec->SetLineColor(2);
0263 h_pt_rec->SetLineStyle(2);
0264
0265 htotal4D->SetXTitle("No. of Segments");
0266 htotal4D->SetLineWidth(2);
0267 htotal4D->SetLineColor(2);
0268
0269 htotalSeg->SetLineWidth(2);
0270 htotalSeg->SetLineColor(4);
0271
0272 htotal4D->SetLineStyle(1);
0273 htotalSeg->SetLineStyle(2);
0274
0275 int theDrawOption = 1;
0276
0277 if (theDrawOption == 0) {
0278 TCanvas* c2 = new TCanvas("innerTSOSXY", "XY", 10, 10, 800, 600);
0279 c2->SetFillColor(0);
0280 c2->SetGrid(1);
0281 c2->SetRightMargin(0.03);
0282 c2->SetTopMargin(0.02);
0283 c2->cd();
0284 h_innerPosXY->Draw("SCAT");
0285 c2->Update();
0286 c2->Write();
0287
0288 TCanvas* c2a = new TCanvas("outerTSOSXY", "Outer XY", 10, 10, 800, 600);
0289 c2a->SetFillColor(0);
0290 c2a->SetGrid(1);
0291 c2a->SetRightMargin(0.03);
0292 c2a->SetTopMargin(0.02);
0293 c2a->cd();
0294 h_outerPosXY->Draw("SCAT");
0295 c2a->Update();
0296 c2a->Write();
0297
0298 TCanvas* c3 = new TCanvas("innerEtaPhi", "Inner #eta #phi", 10, 10, 800, 600);
0299 c3->SetFillColor(0);
0300 c3->SetGrid(1);
0301 c3->SetRightMargin(0.03);
0302 c3->SetTopMargin(0.02);
0303 c3->cd();
0304 h_innerPosEP->Draw("SCAT");
0305 c3->Update();
0306 c3->Write();
0307
0308 TCanvas* c3a = new TCanvas("outerEtaPhi", "Outer #eta #phi", 10, 10, 800, 600);
0309 c3a->SetFillColor(0);
0310 c3a->SetGrid(1);
0311 c3a->SetRightMargin(0.03);
0312 c3a->SetTopMargin(0.02);
0313 c3a->cd();
0314 h_outerPosEP->Draw("SCAT");
0315 c3a->Update();
0316 c3a->Write();
0317
0318 TCanvas* cRes2 = new TCanvas("TrackTheta", "Theta", 10, 10, 800, 600);
0319 cRes2->SetFillColor(0);
0320 cRes2->SetGrid(1);
0321 cRes2->SetRightMargin(0.03);
0322 cRes2->SetTopMargin(0.02);
0323 cRes2->cd();
0324 h_theta->DrawCopy("HE");
0325 cRes2->Update();
0326 cRes2->Write();
0327
0328 TCanvas* cRes3 = new TCanvas("TrackPhi", "phi", 10, 10, 800, 600);
0329 cRes3->SetFillColor(0);
0330 cRes3->SetGrid(1);
0331 cRes3->SetRightMargin(0.03);
0332 cRes3->SetTopMargin(0.02);
0333 cRes3->cd();
0334 h_phi->DrawCopy("HE");
0335 cRes3->Update();
0336 cRes3->Write();
0337
0338 TCanvas* c7 = new TCanvas("nHits", "Number of RecHits in Track", 10, 10, 800, 600);
0339 c7->SetFillColor(0);
0340 c7->SetGrid(1);
0341 c7->SetRightMargin(0.03);
0342 c7->SetTopMargin(0.02);
0343 c7->cd();
0344 hnhit->DrawCopy("HE");
0345 c7->Update();
0346 c7->Write();
0347
0348 TCanvas* cRes8 = new TCanvas("PtDis", "Distribution of P_{T} at SimHit", 10, 10, 800, 600);
0349 cRes8->SetFillColor(0);
0350 cRes8->SetGrid(1);
0351 cRes8->SetRightMargin(0.03);
0352 cRes8->SetTopMargin(0.02);
0353 cRes8->cd();
0354
0355 h_pt_rec->DrawCopy("HE");
0356 TLegend* legend8 = new TLegend(0.7, 0.6, 0.9, 0.8);
0357 legend8->SetTextAlign(32);
0358 legend8->SetTextColor(1);
0359 legend8->SetTextSize(0.04);
0360 legend8->AddEntry("h_pt_rec", "By Cos", "l");
0361 legend8->Draw();
0362 cRes8->Update();
0363 cRes8->Write();
0364
0365 TCanvas* csegs = new TCanvas("csegs", "Total & 4D Segments", 10, 10, 800, 600);
0366
0367 csegs->SetFillColor(0);
0368 csegs->SetGrid(1);
0369 csegs->SetRightMargin(0.03);
0370 csegs->SetTopMargin(0.02);
0371 csegs->cd();
0372
0373 htotal4D->DrawCopy("HE");
0374 htotalSeg->DrawCopy("HEsame");
0375
0376 TLegend* legendseg = new TLegend(0.6, 0.2, 0.9, 0.4);
0377 legendseg->SetTextAlign(32);
0378 legendseg->SetTextColor(1);
0379 legendseg->SetTextSize(0.04);
0380 legendseg->AddEntry("htotal4D", "4D Segments", "l");
0381 legendseg->AddEntry("htotalSeg", "total Segments", "l");
0382 legendseg->Draw();
0383
0384 csegs->Update();
0385 csegs->Write();
0386
0387 } else {
0388 TCanvas* cpdf = new TCanvas("cpdf", "", 0, 1, 500, 700);
0389 cpdf->SetTicks();
0390
0391 TPostScript* pdf = new TPostScript("cosmicValidation.ps", 111);
0392
0393 const int NUM_PAGES = 7;
0394 TPad* pad[NUM_PAGES];
0395 for (int i_page = 0; i_page < NUM_PAGES; i_page++)
0396 pad[i_page] = new TPad("", "", .05, .05, .95, .93);
0397
0398 ostringstream page_print;
0399 int page = 0;
0400
0401 TLatex ttl;
0402
0403 pdf->NewPage();
0404 cpdf->Clear();
0405 cpdf->cd(0);
0406 ttl.DrawLatex(.4, .95, "inner and outer state (x,y) ");
0407 page_print.str("");
0408 page_print << page + 1;
0409 ttl.DrawLatex(.9, .02, page_print.str().c_str());
0410 pad[page]->Draw();
0411 pad[page]->Divide(1, 2);
0412 pad[page]->cd(1);
0413 h_innerPosXY->Draw("SCAT");
0414 pad[page]->cd(2);
0415 h_outerPosXY->Draw("SCAT");
0416
0417 page++;
0418 cpdf->Update();
0419
0420 pdf->NewPage();
0421 cpdf->Clear();
0422 cpdf->cd(0);
0423 ttl.DrawLatex(.4, .95, "inner and outer state (eta, phi) ");
0424 page_print.str("");
0425 page_print << page + 1;
0426 ttl.DrawLatex(.9, .02, page_print.str().c_str());
0427 pad[page]->Draw();
0428 pad[page]->Divide(1, 2);
0429 pad[page]->cd(1);
0430 h_innerPosEP->Draw("SCAT");
0431 pad[page]->cd(2);
0432 h_outerPosEP->Draw("SCAT");
0433
0434 page++;
0435 cpdf->Update();
0436
0437 pdf->NewPage();
0438 cpdf->Clear();
0439 cpdf->cd(0);
0440
0441 pad[page]->cd(3);
0442 h_phi->Draw("HE");
0443
0444 page++;
0445 cpdf->Update();
0446
0447 pdf->NewPage();
0448 cpdf->Clear();
0449 cpdf->cd(0);
0450
0451 page++;
0452 cpdf->Update();
0453
0454 pdf->NewPage();
0455 cpdf->Clear();
0456 cpdf->cd(0);
0457
0458 ttl.DrawLatex(.4, .95, "number of hits ");
0459 page_print.str("");
0460 page_print << page + 1;
0461 ttl.DrawLatex(.9, .02, page_print.str().c_str());
0462 pad[page]->Draw();
0463 pad[page]->cd();
0464 hnhit->Draw("HE");
0465
0466 page++;
0467 cpdf->Update();
0468
0469 pdf->NewPage();
0470 cpdf->Clear();
0471 cpdf->cd(0);
0472
0473 ttl.DrawLatex(.4, .95, "pt distribution ");
0474 page_print.str("");
0475 page_print << page + 1;
0476 ttl.DrawLatex(.9, .02, page_print.str().c_str());
0477 pad[page]->Draw();
0478 pad[page]->Divide(1, 2);
0479 pad[page]->cd(1);
0480 h_pt_rec->Draw("HE");
0481
0482 page++;
0483 cpdf->Update();
0484
0485 pdf->NewPage();
0486 cpdf->Clear();
0487 cpdf->cd(0);
0488
0489 ttl.DrawLatex(.4, .95, "Number of hits ");
0490 page_print.str("");
0491 page_print << page + 1;
0492 ttl.DrawLatex(.9, .02, page_print.str().c_str());
0493 pad[page]->Draw();
0494 pad[page]->Divide(1, 2);
0495 pad[page]->cd(1);
0496 htotal4D->Draw("HE");
0497 pad[page]->cd(2);
0498 htotalSeg->Draw("HE");
0499
0500 pdf->Close();
0501 }
0502 }
0503
0504 edm::ESHandle<Propagator> RealCosmicDataAnalyzer::propagator() const {
0505 return theService->propagator("SteppingHelixPropagatorAny");
0506 }
0507
0508
0509
0510 DEFINE_FWK_MODULE(RealCosmicDataAnalyzer);