Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-03 02:57:00

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);