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/one/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::one::EDAnalyzer<edm::one::SharedResources> {
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   usesResource(TFileService::kSharedResource);
0112 }
0113
0114 RealCosmicDataAnalyzer::~RealCosmicDataAnalyzer() {
0115   if (theService)
0116     delete theService;
0117 }
0118
0119 void RealCosmicDataAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0120   theService->update(iSetup);
0121
0122   nEvent++;
0123   std::cout << "reading event " << nEvent << std::endl;
0124
0125   Handle<reco::TrackCollection> muons;
0126   iEvent.getByLabel(trackLabel_, muons);
0127   cout << "cosmic Muon: " << muons->size() << endl;
0128
0129   //  if (muons->empty()) return;
0130   if (!muons->empty())
0131     successR++;
0132
0133   // Get Segment collections from the Event
0134   edm::Handle<DTRecSegment4DCollection> dtSegments;
0135   iEvent.getByLabel("dt4DSegments", dtSegments);
0136
0137   edm::Handle<CSCSegmentCollection> cscSegments;
0138   iEvent.getByLabel("cscSegments", cscSegments);
0139
0140   int nSeg = dtSegments->size() + cscSegments->size();
0141   cout << "cscSegments: " << cscSegments->size() << endl;
0142
0143   htotalSeg->Fill(nSeg);
0144
0145   int n4D = 0;
0146
0147   for (DTRecSegment4DCollection::const_iterator idt = dtSegments->begin(); idt != dtSegments->end(); ++idt) {
0148     if (idt->dimension() < 4)
0149       continue;
0150     bool sameChamber = false;
0151     DetId thisId = idt->geographicalId();
0152
0153     for (DTRecSegment4DCollection::const_iterator idt2 = dtSegments->begin(); idt2 != idt; ++idt2) {
0154       if (idt2->geographicalId() == thisId)
0155         sameChamber = true;
0156     }
0157     if (!sameChamber)
0158       n4D++;
0159   }
0160
0161   for (CSCSegmentCollection::const_iterator icsc = cscSegments->begin(); icsc != cscSegments->end(); ++icsc) {
0162     if (icsc->dimension() < 4)
0163       continue;
0164     bool sameChamber = false;
0165     DetId thisId = icsc->geographicalId();
0166
0167     for (CSCSegmentCollection::const_iterator icsc2 = cscSegments->begin(); icsc2 != icsc; ++icsc2) {
0168       if (icsc2->geographicalId() == thisId)
0169         sameChamber = true;
0170     }
0171
0172     if (!sameChamber)
0173       n4D++;
0174   }
0175
0176   htotal4D->Fill(n4D);
0177   if (n4D < 2)
0178     nNoSignal++;
0179
0180   if (muons->empty())
0181     return;
0182
0183   for (reco::TrackCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0184     cout << "cosmic Muon Track: "
0185          << " mom: (" << muon->px() << "," << muon->py() << "," << muon->pz() << ")" << endl;
0186     if (fabs(muon->p()) < 1e-5)
0187       return;  //prevent those failed to extrapolation to vertex
0188
0189     math::XYZVector innerMo = muon->innerMomentum();
0190
0191     float ptreco = muon->pt();
0192
0193     h_pt_rec->Fill(ptreco);
0194
0195     hnhit->Fill(muon->recHitsSize());
0196
0197     GlobalVector im(innerMo.x(), innerMo.y(), innerMo.z());
0198
0199     h_theta->Fill(im.theta());
0200     h_phi->Fill(im.phi());
0201
0202     math::XYZPoint innerPo = muon->innerPosition();
0203     GlobalPoint ip(innerPo.x(), innerPo.y(), innerPo.z());
0204
0205     h_innerPosXY->Fill(ip.x(), ip.y());
0206     h_innerPosEP->Fill(ip.eta(), Geom::Phi<float>(ip.phi()));
0207
0208     math::XYZPoint outerPo = muon->outerPosition();
0209     GlobalPoint op(outerPo.x(), outerPo.y(), outerPo.z());
0210
0211     h_outerPosXY->Fill(op.x(), op.y());
0212     h_outerPosEP->Fill(op.eta(), Geom::Phi<float>(op.phi()));
0213     h_inOutDis->Fill((ip - op).perp());
0214   }
0215 }
0216
0217 void RealCosmicDataAnalyzer::beginJob() {
0218   edm::Service<TFileService> fs;
0219
0220   h_theta = fs->make<TH1F>("h_theta", "theta angle ", 50, -0.1, 0.1);
0221   h_phi = fs->make<TH1F>("h_phi", "phi angle ", 50, -2.0, 2.0);
0222
0223   hnhit = fs->make<TH1F>("hnhit", "Number of Hits in Track by Cos", 60, 0.0, 60.0);
0224
0225   h_innerPosXY = fs->make<TH2F>("h_innerPosXY", "inner x-y", 100, -700.0, 700.0, 100, -700.0, 700.0);
0226   h_innerPosEP = fs->make<TH2F>("h_innerPosEP", "inner #eta-#phi", 100, -2.4, 2.4, 100, -3.3, 3.3);
0227
0228   h_outerPosXY = fs->make<TH2F>("h_outerPosXY", "outer x-y", 100, -700.0, 700.0, 100, -700.0, 700.0);
0229   h_outerPosEP = fs->make<TH2F>("h_outerPosEP", "outer #eta-#phi", 100, -2.4, 2.4, 100, -3.3, 3.3);
0230
0231   h_pt_rec = fs->make<TH1F>("h_pt_rec", "distribution of P_{T}", 100, 0.0, 100.0);
0232   h_inOutDis = fs->make<TH1F>("h_inOutDis", "distribution of inner outer distance", 200, 0.0, 1500.0);
0233   htotal4D = fs->make<TH1F>("htotal4D", "# of Segments", 15, 0.0, 15.0);
0234   htotalSeg = fs->make<TH1F>("htotalSeg", "# of Segments", 15, 0.0, 15.0);
0235 }
0236
0237 void RealCosmicDataAnalyzer::endJob() {
0238   h_innerPosXY->SetXTitle("X");
0239   h_innerPosXY->SetYTitle("Y");
0240   h_innerPosXY->SetMarkerColor(2);
0241   h_innerPosXY->SetMarkerStyle(24);
0242
0243   h_outerPosXY->SetXTitle("X");
0244   h_outerPosXY->SetYTitle("Y");
0245   h_outerPosXY->SetMarkerColor(2);
0246   h_outerPosXY->SetMarkerStyle(24);
0247
0248   h_innerPosEP->SetXTitle("#eta");
0249   h_innerPosEP->SetYTitle("#phi");
0250   h_innerPosEP->SetMarkerColor(2);
0251   h_innerPosEP->SetMarkerStyle(24);
0252
0253   h_outerPosEP->SetXTitle("#eta");
0254   h_outerPosEP->SetYTitle("#phi");
0255   h_outerPosEP->SetMarkerColor(2);
0256   h_outerPosEP->SetMarkerStyle(24);
0257
0258   hnhit->SetXTitle("N_{RecHits}");
0259   hnhit->SetLineWidth(2);
0260   hnhit->SetLineColor(2);
0261   hnhit->SetLineStyle(1);
0262
0263   h_pt_rec->SetLineWidth(2);
0264   h_pt_rec->SetLineColor(2);
0265   h_pt_rec->SetLineStyle(2);
0266
0267   htotal4D->SetXTitle("No. of Segments");
0268   htotal4D->SetLineWidth(2);
0269   htotal4D->SetLineColor(2);
0270
0271   htotalSeg->SetLineWidth(2);
0272   htotalSeg->SetLineColor(4);
0273
0274   htotal4D->SetLineStyle(1);
0275   htotalSeg->SetLineStyle(2);
0276
0277   int theDrawOption = 1;
0278
0279   if (theDrawOption == 0) {
0280     TCanvas* c2 = new TCanvas("innerTSOSXY", "XY", 10, 10, 800, 600);
0281     c2->SetFillColor(0);
0282     c2->SetGrid(1);
0283     c2->SetRightMargin(0.03);
0284     c2->SetTopMargin(0.02);
0285     c2->cd();
0286     h_innerPosXY->Draw("SCAT");
0287     c2->Update();
0288     c2->Write();
0289
0290     TCanvas* c2a = new TCanvas("outerTSOSXY", "Outer XY", 10, 10, 800, 600);
0291     c2a->SetFillColor(0);
0292     c2a->SetGrid(1);
0293     c2a->SetRightMargin(0.03);
0294     c2a->SetTopMargin(0.02);
0295     c2a->cd();
0296     h_outerPosXY->Draw("SCAT");
0297     c2a->Update();
0298     c2a->Write();
0299
0300     TCanvas* c3 = new TCanvas("innerEtaPhi", "Inner #eta #phi", 10, 10, 800, 600);
0301     c3->SetFillColor(0);
0302     c3->SetGrid(1);
0303     c3->SetRightMargin(0.03);
0304     c3->SetTopMargin(0.02);
0305     c3->cd();
0306     h_innerPosEP->Draw("SCAT");
0307     c3->Update();
0308     c3->Write();
0309
0310     TCanvas* c3a = new TCanvas("outerEtaPhi", "Outer #eta #phi", 10, 10, 800, 600);
0311     c3a->SetFillColor(0);
0312     c3a->SetGrid(1);
0313     c3a->SetRightMargin(0.03);
0314     c3a->SetTopMargin(0.02);
0315     c3a->cd();
0316     h_outerPosEP->Draw("SCAT");
0317     c3a->Update();
0318     c3a->Write();
0319
0320     TCanvas* cRes2 = new TCanvas("TrackTheta", "Theta", 10, 10, 800, 600);
0321     cRes2->SetFillColor(0);
0322     cRes2->SetGrid(1);
0323     cRes2->SetRightMargin(0.03);
0324     cRes2->SetTopMargin(0.02);
0325     cRes2->cd();
0326     h_theta->DrawCopy("HE");
0327     cRes2->Update();
0328     cRes2->Write();
0329
0330     TCanvas* cRes3 = new TCanvas("TrackPhi", "phi", 10, 10, 800, 600);
0331     cRes3->SetFillColor(0);
0332     cRes3->SetGrid(1);
0333     cRes3->SetRightMargin(0.03);
0334     cRes3->SetTopMargin(0.02);
0335     cRes3->cd();
0336     h_phi->DrawCopy("HE");
0337     cRes3->Update();
0338     cRes3->Write();
0339
0340     TCanvas* c7 = new TCanvas("nHits", "Number of RecHits in Track", 10, 10, 800, 600);
0341     c7->SetFillColor(0);
0342     c7->SetGrid(1);
0343     c7->SetRightMargin(0.03);
0344     c7->SetTopMargin(0.02);
0345     c7->cd();
0346     hnhit->DrawCopy("HE");
0347     c7->Update();
0348     c7->Write();
0349
0350     TCanvas* cRes8 = new TCanvas("PtDis", "Distribution of P_{T} at SimHit", 10, 10, 800, 600);
0351     cRes8->SetFillColor(0);
0352     cRes8->SetGrid(1);
0353     cRes8->SetRightMargin(0.03);
0354     cRes8->SetTopMargin(0.02);
0355     cRes8->cd();
0356
0357     h_pt_rec->DrawCopy("HE");
0358     TLegend* legend8 = new TLegend(0.7, 0.6, 0.9, 0.8);
0359     legend8->SetTextAlign(32);
0360     legend8->SetTextColor(1);
0361     legend8->SetTextSize(0.04);
0362     legend8->AddEntry("h_pt_rec", "By Cos", "l");
0363     legend8->Draw();
0364     cRes8->Update();
0365     cRes8->Write();
0366
0367     TCanvas* csegs = new TCanvas("csegs", "Total & 4D Segments", 10, 10, 800, 600);
0368
0369     csegs->SetFillColor(0);
0370     csegs->SetGrid(1);
0371     csegs->SetRightMargin(0.03);
0372     csegs->SetTopMargin(0.02);
0373     csegs->cd();
0374
0375     htotal4D->DrawCopy("HE");
0376     htotalSeg->DrawCopy("HEsame");
0377
0378     TLegend* legendseg = new TLegend(0.6, 0.2, 0.9, 0.4);
0379     legendseg->SetTextAlign(32);
0380     legendseg->SetTextColor(1);
0381     legendseg->SetTextSize(0.04);
0382     legendseg->AddEntry("htotal4D", "4D Segments", "l");
0383     legendseg->AddEntry("htotalSeg", "total Segments", "l");
0384     legendseg->Draw();
0385
0386     csegs->Update();
0387     csegs->Write();
0388
0389   } else {
0390     TCanvas* cpdf = new TCanvas("cpdf", "", 0, 1, 500, 700);
0391     cpdf->SetTicks();
0392
0393     TPostScript* pdf = new TPostScript("cosmicValidation.ps", 111);
0394
0395     const int NUM_PAGES = 7;
0396     TPad* pad[NUM_PAGES];
0397     for (int i_page = 0; i_page < NUM_PAGES; i_page++)
0398       pad[i_page] = new TPad("", "", .05, .05, .95, .93);
0399
0400     ostringstream page_print;
0401     int page = 0;
0402
0403     TLatex ttl;
0404
0405     pdf->NewPage();
0406     cpdf->Clear();
0407     cpdf->cd(0);
0408     ttl.DrawLatex(.4, .95, "inner and outer state (x,y) ");
0409     page_print.str("");
0410     page_print << page + 1;
0411     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0412     pad[page]->Draw();
0413     pad[page]->Divide(1, 2);
0414     pad[page]->cd(1);
0415     h_innerPosXY->Draw("SCAT");
0416     pad[page]->cd(2);
0417     h_outerPosXY->Draw("SCAT");
0418
0419     page++;
0420     cpdf->Update();
0421
0422     pdf->NewPage();
0423     cpdf->Clear();
0424     cpdf->cd(0);
0425     ttl.DrawLatex(.4, .95, "inner and outer state (eta, phi) ");
0426     page_print.str("");
0427     page_print << page + 1;
0428     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0429     pad[page]->Draw();
0430     pad[page]->Divide(1, 2);
0431     pad[page]->cd(1);
0432     h_innerPosEP->Draw("SCAT");
0433     pad[page]->cd(2);
0434     h_outerPosEP->Draw("SCAT");
0435
0436     page++;
0437     cpdf->Update();
0438
0439     pdf->NewPage();
0440     cpdf->Clear();
0441     cpdf->cd(0);
0442
0443     pad[page]->cd(3);
0444     h_phi->Draw("HE");
0445
0446     page++;
0447     cpdf->Update();
0448
0449     pdf->NewPage();
0450     cpdf->Clear();
0451     cpdf->cd(0);
0452
0453     page++;
0454     cpdf->Update();
0455
0456     pdf->NewPage();
0457     cpdf->Clear();
0458     cpdf->cd(0);
0459
0460     ttl.DrawLatex(.4, .95, "number of hits ");
0461     page_print.str("");
0462     page_print << page + 1;
0463     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0464     pad[page]->Draw();
0465     pad[page]->cd();
0466     hnhit->Draw("HE");
0467
0468     page++;
0469     cpdf->Update();
0470
0471     pdf->NewPage();
0472     cpdf->Clear();
0473     cpdf->cd(0);
0474
0475     ttl.DrawLatex(.4, .95, "pt distribution ");
0476     page_print.str("");
0477     page_print << page + 1;
0478     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0479     pad[page]->Draw();
0480     pad[page]->Divide(1, 2);
0481     pad[page]->cd(1);
0482     h_pt_rec->Draw("HE");
0483
0484     page++;
0485     cpdf->Update();
0486
0487     pdf->NewPage();
0488     cpdf->Clear();
0489     cpdf->cd(0);
0490
0491     ttl.DrawLatex(.4, .95, "Number of hits ");
0492     page_print.str("");
0493     page_print << page + 1;
0494     ttl.DrawLatex(.9, .02, page_print.str().c_str());
0495     pad[page]->Draw();
0496     pad[page]->Divide(1, 2);
0497     pad[page]->cd(1);
0498     htotal4D->Draw("HE");
0499     pad[page]->cd(2);
0500     htotalSeg->Draw("HE");
0501
0502     pdf->Close();
0503   }
0504 }
0505
0506 edm::ESHandle<Propagator> RealCosmicDataAnalyzer::propagator() const {
0507   return theService->propagator("SteppingHelixPropagatorAny");
0508 }
0509
0510 //define this as a plug-in
0511
0512 DEFINE_FWK_MODULE(RealCosmicDataAnalyzer);