CosmicMuonProducer/test/RealCosmicDataAnalyzer.cc

0001 /** \class RealCosmicDataAnalyzer
0002  *
0003  *  \author Chang Liu   -  Purdue University <Chang.Liu@cern.ch>
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   // service parameters
0103   edm::ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0104   // the services
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   //  if (muons->empty()) return;
0128   if (!muons->empty())
0129     successR++;
0130
0131   // Get Segment collections from the Event
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;  //prevent those failed to extrapolation to vertex
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 //define this as a plug-in
0509
0510 DEFINE_FWK_MODULE(RealCosmicDataAnalyzer);