CSCSegment/test/CSCSegmentVisualise.cc

0001 /** \file CSCSegmentVisualise
0002  *
0003  *  \author D. Fortin - UC Riverside
0004  */
0005
0006 #include <RecoLocalMuon/CSCSegment/test/CSCSegmentVisualise.h>
0007
0008 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0009 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0010 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
0011 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
0012 #include <DataFormats/CSCRecHit/interface/CSCRangeMapAccessor.h>
0013
0014 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0015 #include <Geometry/CSCGeometry/interface/CSCChamber.h>
0016 #include <Geometry/CSCGeometry/interface/CSCChamberSpecs.h>
0017 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
0018 #include <Geometry/CSCGeometry/interface/CSCLayerGeometry.h>
0019 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0020
0021 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0022
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028
0029 #include "TFile.h"
0030
0031 // Constructor
0032 CSCSegmentVisualise::CSCSegmentVisualise(const edm::ParameterSet& pset) {
0033   filename = pset.getUntrackedParameter<std::string>("RootFileName");
0034   minRechitChamber = pset.getUntrackedParameter<int>("minRechitPerChamber");
0035   maxRechitChamber = pset.getUntrackedParameter<int>("maxRechitPerChamber");
0036
0037   geomToken_ = esConsumes();
0038   simHitsToken_ = consumes(edm::InputTag("g4SimHits", "MuonCSCHits"));
0039   recHitsToken_ = consumes(edm::InputTag("csc2DRecHits"));
0040   segmentsToken_ = consumes(edm::InputTag("cscSegments"));
0041
0042   file = new TFile(filename.c_str(), "RECREATE");
0043
0044   if (file->IsOpen())
0045     std::cout << "file open!" << std::endl;
0046   else
0047     std::cout << "*** Error in opening file ***" << std::endl;
0048
0049   idxHisto = 0;
0050   ievt = 0;
0051 }
0052
0053 // Destructor
0054 CSCSegmentVisualise::~CSCSegmentVisualise() {
0055   file->cd();
0056
0057   for (int i = 0; i < idxHisto; i++) {
0058     hxvsz[i]->Write();
0059     hyvsz[i]->Write();
0060
0061     hxvszSeg[i]->Write();
0062     hyvszSeg[i]->Write();
0063
0064     hxvszSegP[i]->Write();
0065     hyvszSegP[i]->Write();
0066
0067     hxvszE[i]->Write();
0068     hyvszE[i]->Write();
0069   }
0070   file->Close();
0071 }
0072
0073 // The real analysis
0074 void CSCSegmentVisualise::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0075   ievt++;
0076
0077   if (idxHisto > 99)
0078     return;
0079
0080   const CSCGeometry& geom_ = eventSetup.getData(geomToken_);
0081   const edm::PSimHitContainer& simHits = event.get(simHitsToken_);
0082   const CSCRecHit2DCollection& recHits = event.get(recHitsToken_);
0083   const CSCSegmentCollection& segments = event.get(segmentsToken_);
0084
0085   std::vector<CSCDetId> chambers;
0086   std::vector<CSCDetId>::const_iterator chIt;
0087
0088   // First, create vector of chambers with rechits
0089   for (CSCRecHit2DCollection::const_iterator it2 = recHits.begin(); it2 != recHits.end(); it2++) {
0090     bool insert = true;
0091     for (chIt = chambers.begin(); chIt != chambers.end(); ++chIt)
0092       if (((*it2).cscDetId().chamber() == (*chIt).chamber()) && ((*it2).cscDetId().station() == (*chIt).station()) &&
0093           ((*it2).cscDetId().ring() == (*chIt).ring()) && ((*it2).cscDetId().endcap() == (*chIt).endcap()))
0094         insert = false;
0095
0096     if (insert)
0097       chambers.push_back((*it2).cscDetId().chamberId());
0098   }
0099
0100   for (chIt = chambers.begin(); chIt != chambers.end(); ++chIt) {
0101     std::vector<const CSCRecHit2D*> cscRecHits;
0102     std::vector<const CSCSegment*> cscSegments;
0103     std::vector<const CSCRecHit2D*> eRecHits;
0104     const CSCChamber* chamber = geom_.chamber(*chIt);
0105
0106     CSCRangeMapAccessor acc;
0107     CSCRecHit2DCollection::range range = recHits.get(acc.cscChamber(*chIt));
0108
0109     for (CSCRecHit2DCollection::const_iterator rechit = range.first; rechit != range.second; rechit++) {
0110       cscRecHits.push_back(&(*rechit));
0111     }
0112
0113     if (int(cscRecHits.size()) < minRechitChamber)
0114       continue;
0115     if (int(cscRecHits.size()) > maxRechitChamber)
0116       continue;
0117     if (chamber->id().ring() == 4)
0118       continue;  // not interested in ME1a for now...
0119
0120     std::cout << "chamber hits satisfy criteria " << std::endl;
0121
0122     float xmin = -10.;
0123     float xmax = 10.;
0124     float ymin = -10.;
0125     float ymax = 10.;
0126
0127     // Fill histograms for rechits:
0128     int nHits = cscRecHits.size();
0129     for (int i = 0; i < nHits; ++i) {
0130       const CSCRecHit2D* rhit = cscRecHits[i];
0131       LocalPoint lphit = (*rhit).localPosition();
0132       if (i == 0) {
0133         xmin = lphit.x();
0134         ymin = lphit.y();
0135         xmax = lphit.x();
0136         ymax = lphit.y();
0137       } else {
0138         if (lphit.x() < xmin)
0139           xmin = lphit.x();
0140         if (lphit.y() < ymin)
0141           ymin = lphit.y();
0142         if (lphit.x() > xmax)
0143           xmax = lphit.x();
0144         if (lphit.y() > ymax)
0145           ymax = lphit.y();
0146       }
0147     }
0148
0149     xmin = xmin - 5.;
0150     xmax = xmax + 5.;
0151     ymin = ymin - 5.;
0152     ymax = ymax + 5.;
0153
0154     char a[14];
0155     char evt[10];
0156
0157     sprintf(evt, "Event %d", ievt);
0158
0159     // Create histograms on the fly
0160
0161     // rechits
0162     sprintf(a, "h%d", idxHisto + 100);
0163     hxvsz[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, xmin, xmax);
0164     sprintf(a, "h%d", idxHisto + 200);
0165     hyvsz[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, ymin, ymax);
0166
0167     // Hits on segment
0168     sprintf(a, "h%d", idxHisto + 300);
0169     hxvszSeg[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, xmin, xmax);
0170     sprintf(a, "h%d", idxHisto + 400);
0171     hyvszSeg[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, ymin, ymax);
0172
0173     // Segment projection
0174     sprintf(a, "h%d", idxHisto + 500);
0175     hxvszSegP[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, xmin, xmax);
0176     sprintf(a, "h%d", idxHisto + 600);
0177     hyvszSegP[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, ymin, ymax);
0178
0179     // recHits from electrons/delta rays
0180     sprintf(a, "h%d", idxHisto + 700);
0181     hxvszE[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, xmin, xmax);
0182     sprintf(a, "h%d", idxHisto + 800);
0183     hyvszE[idxHisto] = new TH2F(a, evt, 100, -10., 10., 100, ymin, ymax);
0184
0185     std::cout << "done creating histograms" << std::endl;
0186
0187     // Fill histograms for rechits:
0188     for (int i = 0; i < nHits; ++i) {
0189       const CSCRecHit2D* rhit = cscRecHits[i];
0190       CSCDetId id = (CSCDetId)(*rhit).cscDetId();
0191       const CSCLayer* csclayer = geom_.layer(id);
0192       LocalPoint lphit = (*rhit).localPosition();
0193       GlobalPoint gphit = csclayer->toGlobal(lphit);
0194       LocalPoint lphitChamber = chamber->toLocal(gphit);
0195       float xhit = lphitChamber.x();
0196       float yhit = lphitChamber.y();
0197       float zhit = lphitChamber.z();
0198
0199       hxvsz[idxHisto]->Fill(zhit, xhit);
0200       hyvsz[idxHisto]->Fill(zhit, yhit);
0201     }
0202
0203     std::cout << "done filling rechit histos " << std::endl;
0204
0205     // Sort rechits which comes from electron hits per chamber type:
0206
0207     for (CSCRecHit2DCollection::const_iterator rec_it = range.first; rec_it != range.second; rec_it++) {
0208       bool isElec = false;
0209       float r_closest = 9999;
0210
0211       CSCDetId idrec = (CSCDetId)(*rec_it).cscDetId();
0212       LocalPoint rhitlocal = (*rec_it).localPosition();
0213
0214       for (edm::PSimHitContainer::const_iterator sim_it = simHits.begin(); sim_it != simHits.end(); sim_it++) {
0215         CSCDetId idsim = (CSCDetId)(*sim_it).detUnitId();
0216
0217         if (idrec.endcap() == idsim.endcap() && idrec.station() == idsim.station() && idrec.ring() == idsim.ring() &&
0218             idrec.chamber() == idsim.chamber() && idrec.layer() == idsim.layer()) {
0219           LocalPoint shitlocal = (*sim_it).localPosition();
0220
0221           float dx2 = (rhitlocal.x() - shitlocal.x()) * (rhitlocal.x() - shitlocal.x());
0222           float dy2 = (rhitlocal.y() - shitlocal.y()) * (rhitlocal.y() - shitlocal.y());
0223           float dr2 = dx2 + dy2;
0224           if (dr2 < r_closest) {
0225             r_closest = dr2;
0226             if (abs((*sim_it).particleType()) != 13) {
0227               isElec = true;
0228             } else {
0229               isElec = false;
0230             }
0231           }
0232         }
0233       }
0234       if (isElec)
0235         eRecHits.push_back(&(*rec_it));
0236     }
0237
0238     nHits = eRecHits.size();
0239
0240     // Fill histograms for rechits from electrons:
0241     for (int i = 0; i < nHits; ++i) {
0242       const CSCRecHit2D* rhit = eRecHits[i];
0243       CSCDetId id = (CSCDetId)(*rhit).cscDetId();
0244       const CSCLayer* csclayer = geom_.layer(id);
0245       LocalPoint lphit = (*rhit).localPosition();
0246       GlobalPoint gphit = csclayer->toGlobal(lphit);
0247       LocalPoint lphitChamber = chamber->toLocal(gphit);
0248       float xhit = lphitChamber.x();
0249       float yhit = lphitChamber.y();
0250       float zhit = lphitChamber.z();
0251
0252       hxvszE[idxHisto]->Fill(zhit, xhit);
0253       hyvszE[idxHisto]->Fill(zhit, yhit);
0254     }
0255
0256     // Then, sort segments per chamber type as well:
0257
0258     CSCSegmentCollection::range range3 = segments.get(acc.cscChamber(*chIt));
0259     for (CSCSegmentCollection::const_iterator segments = range3.first; segments != range3.second; segments++) {
0260       cscSegments.push_back(&(*segments));
0261     }
0262
0263     // Fill histograms for segments
0264     int nSegs = cscSegments.size();
0265     for (int j = 0; j < nSegs; ++j) {
0266       const CSCSegment* segs = cscSegments[j];
0267
0268       LocalVector vec = (*segs).localDirection();
0269       LocalPoint ori = (*segs).localPosition();
0270
0271       double dxdz = vec.x() / vec.z();
0272       double dydz = vec.y() / vec.z();
0273
0274       float zstart = -10.;
0275       float zstop = 10.;
0276
0277       float step_size = fabs(zstop - zstart) / 100.;
0278
0279       std::cout << "Projected Segment" << std::endl;
0280
0281       for (int i = 0; i < 100; ++i) {
0282         float z_proj = zstart + (i * step_size);
0283         float x_proj = dxdz * z_proj + ori.x();
0284         float y_proj = dydz * z_proj + ori.y();
0285
0286         std::cout << i << " " << x_proj << " " << y_proj << " " << z_proj << std::endl;
0287
0288         hxvszSegP[idxHisto]->Fill(z_proj, x_proj);
0289         hyvszSegP[idxHisto]->Fill(z_proj, y_proj);
0290       }
0291
0292       // Loop over hits used in segment
0293       const std::vector<CSCRecHit2D>& rhseg = (*segs).specificRecHits();
0294       std::vector<CSCRecHit2D>::const_iterator rh_i;
0295
0296       for (rh_i = rhseg.begin(); rh_i != rhseg.end(); ++rh_i) {
0297         CSCDetId id = (CSCDetId)(*rh_i).cscDetId();
0298         const CSCLayer* csclayer = geom_.layer(id);
0299         LocalPoint lphit = (*rh_i).localPosition();
0300         GlobalPoint gphit = csclayer->toGlobal(lphit);
0301         LocalPoint lphitChamber = chamber->toLocal(gphit);
0302         float xhit = lphitChamber.x();
0303         float yhit = lphitChamber.y();
0304         float zhit = lphitChamber.z();
0305
0306         hxvszSeg[idxHisto]->Fill(zhit, xhit);
0307         hyvszSeg[idxHisto]->Fill(zhit, yhit);
0308       }
0309     }
0310
0311     idxHisto++;
0312   }
0313 }
0314
0315 DEFINE_FWK_MODULE(CSCSegmentVisualise);