Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-10 23:56:23 

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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!