GEMDigitizer/test/ME0DigiReader.cc

0001 /** \class ME0DigiReader
0002  *
0003  *  Dumps ME0 digis
0004  *
0005  *  \authors: Roumyana Hadjiiska
0006  */
0007
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013
0014 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0015 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0016 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0017 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0018 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0019 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0020 #include "DataFormats/GEMDigi/interface/ME0DigiCollection.h"
0021 #include "DataFormats/Common/interface/DetSetVector.h"
0022 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0023 #include <map>
0024 #include <set>
0025
0026 #include "DataFormats/Common/interface/DetSet.h"
0027 #include <iostream>
0028
0029 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
0030 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0031
0032 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0033 #include "TH1F.h"
0034 #include "TH2F.h"
0035 #include "TGraph.h"
0036 #include "TGraphErrors.h"
0037
0038 using namespace std;
0039
0040 class ME0DigiReader : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0041 public:
0042   explicit ME0DigiReader(const edm::ParameterSet &pset);
0043
0044   ~ME0DigiReader() override {}
0045
0046   void beginJob() override;
0047   void analyze(const edm::Event &, const edm::EventSetup &) override;
0048   void endJob() override;
0049
0050 private:
0051   const edm::EDGetTokenT<edm::PSimHitContainer> simhitToken_;
0052   const edm::EDGetTokenT<ME0DigiCollection> me0DigiToken_;
0053   const edm::ESGetToken<ME0Geometry, MuonGeometryRecord> geomToken_;
0054   const bool debug_;
0055
0056   TH1F *hProces;
0057   TH2F *hNstripEtaParts;
0058   TH1F *hBx;
0059   TH2F *hRadiusEtaPartVsNdigi;
0060   TH2F *hRadiusEtaPartVsNdigiOvTrArea;
0061   TH1F *hRadiusEtaPart;
0062   TH1F *hdeltaXEntryPointVsCentreStrip;
0063   TH1F *hResidualsSimPhi;
0064   TH1F *hResidualsDigiPhi;
0065   TH1F *hResidualsSimVsDigiPhi;
0066   TGraphErrors *grRatePerRoll;
0067
0068   int numbEvents;
0069   int ndigiEtaPart[8] = {0};
0070   double ndigiVsArea[8] = {0};
0071   double rollRadiusEtaPart[8] = {0};
0072 };
0073
0074 ME0DigiReader::ME0DigiReader(const edm::ParameterSet &pset)
0075     : simhitToken_(consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("simhitToken"))),
0076       me0DigiToken_(consumes<ME0DigiCollection>(pset.getParameter<edm::InputTag>("me0DigiToken"))),
0077       geomToken_(esConsumes<ME0Geometry, MuonGeometryRecord>()),
0078       debug_(pset.getParameter<bool>("debugFlag")) {
0079   usesResource("TFileService");
0080   edm::Service<TFileService> fs;
0081
0082   hProces = fs->make<TH1F>("hProces", "Process type for all the simHits", 20, 0, 20);
0083   hNstripEtaParts = fs->make<TH2F>("NstripEtaParts", "Nstrips in all EtaPartitions", 20, 0, 10, 770, 1, 770);
0084   hBx = fs->make<TH1F>("hBx", "bx from digi - for all #eta partiotions", 9, -5.5, 3.5);
0085   hRadiusEtaPartVsNdigi =
0086       fs->make<TH2F>("hRadiusEtaPartVsNdigi", "Radius Eta Partition vs Ndigi", 2500, 0., 250., 200, 0., 20.);
0087   hRadiusEtaPartVsNdigiOvTrArea = fs->make<TH2F>(
0088       "hRadiusEtaPartVsNdigiOvTrArea", "Ndigi/TrArea vs Radius Eta Partition", 2500, 0., 250., 1000, 0., 0.1);
0089   hRadiusEtaPart = fs->make<TH1F>("hRadiusEtaPart", "Radius Eta Partition", 200, 0., 200.);
0090   hdeltaXEntryPointVsCentreStrip = fs->make<TH1F>("deltaX", "delta X Residuals", 200, -10., 10.);
0091   hResidualsSimPhi = fs->make<TH1F>("ResidualsSimPhi", "Global SimMuon Phi", 200, -10., 10.);
0092   hResidualsDigiPhi = fs->make<TH1F>("ResidualsDigiPhi", "Global DigiMuon Phi", 200, -10., 10.);
0093   hResidualsSimVsDigiPhi = fs->make<TH1F>("ResidualsSimVsDigiPhi", "Residuals (SimMuon-Digi) Phi", 50000, -0.5, 0.5);
0094
0095   grRatePerRoll = fs->make<TGraphErrors>(8);
0096   grRatePerRoll->SetName("grRatePerRoll");
0097   grRatePerRoll->SetTitle("ME0 Rate vs Roll Radius - BKG model");
0098
0099   numbEvents = 0;
0100 }
0101
0102 void ME0DigiReader::beginJob() {}
0103
0104 void ME0DigiReader::analyze(const edm::Event &event, const edm::EventSetup &eventSetup) {
0105   const auto &pDD = eventSetup.getHandle(geomToken_);
0106
0107   const auto &simHits = event.getHandle(simhitToken_);
0108
0109   const auto &digis = event.getHandle(me0DigiToken_);
0110
0111   ME0DigiCollection::DigiRangeIterator detUnitIt;
0112
0113   int countRoll[8] = {0};
0114   for (detUnitIt = digis->begin(); detUnitIt != digis->end(); ++detUnitIt) {
0115     const ME0DetId &id = (*detUnitIt).first;
0116     const ME0EtaPartition *roll = pDD->etaPartition(id);
0117
0118     int ndigi = 0;
0119     double trArea(0.0);
0120     double trStripArea(0.0);
0121     Local3DPoint locMuonEntry(0., 0., 0.);
0122     GlobalPoint globMuonEntry(0., 0., 0.);
0123     double simMuPhi = -99.;
0124     double deltaPhi = -99.;
0125     Local3DPoint locDigi(0., 0., 0.);
0126     GlobalPoint pointDigiHit;
0127
0128     const TrapezoidalStripTopology *top_(dynamic_cast<const TrapezoidalStripTopology *>(&(roll->topology())));
0129
0130     const float rollRadius = top_->radius();
0131     const float striplength(top_->stripLength());
0132     const int nstrips = roll->nstrips();
0133     trStripArea = (roll->pitch()) * striplength;
0134     trArea = trStripArea * nstrips;
0135
0136     if (id.roll() > 0 && id.roll() < 9)
0137       countRoll[id.roll() - 1]++;
0138
0139     // Loop over the digis of this DetUnit
0140     const ME0DigiCollection::Range &range = (*detUnitIt).second;
0141     for (ME0DigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0142       hNstripEtaParts->Fill(id.roll(), digiIt->strip());
0143       if (id.roll() > 0 && id.roll() < 9)
0144         ndigiEtaPart[id.roll() - 1]++;
0145
0146       //bx
0147       hBx->Fill(digiIt->bx());
0148       ndigi++;
0149
0150       if (digiIt->strip() < 1 || digiIt->strip() > roll->nstrips()) {
0151         LogDebug("ME0DigiReader") << " XXXXXXXXXXXXX Problemt with " << id
0152                                   << "  a digi has strip# = " << digiIt->strip() << endl;
0153       }
0154
0155       for (const auto &simHit : *simHits) {
0156         hProces->Fill(simHit.processType());
0157
0158         ME0DetId me0Id(simHit.detUnitId());
0159
0160         if (me0Id == id)  //&& abs(simHit.particleType()) == 13)
0161         {
0162           locMuonEntry = simHit.entryPoint();
0163           globMuonEntry = roll->toGlobal(locMuonEntry);
0164           simMuPhi = globMuonEntry.phi();
0165           locDigi = roll->centreOfStrip(digiIt->strip());
0166           pointDigiHit = roll->toGlobal(locDigi);
0167           double digiPhi = pointDigiHit.phi();
0168           deltaPhi = simMuPhi - digiPhi;
0169           hdeltaXEntryPointVsCentreStrip->Fill((simHit.entryPoint().x() - roll->centreOfStrip(digiIt->strip()).x()));
0170           hResidualsSimPhi->Fill(simMuPhi);
0171           hResidualsDigiPhi->Fill(digiPhi);
0172           hResidualsSimVsDigiPhi->Fill(deltaPhi);
0173         }
0174       }
0175     }  // for digis in roll
0176
0177     hRadiusEtaPartVsNdigi->Fill(rollRadius, ndigi);
0178     hRadiusEtaPartVsNdigiOvTrArea->Fill(rollRadius, ndigi / trArea);
0179     hRadiusEtaPart->Fill(rollRadius);
0180
0181     if (id.roll() > 0 && id.roll() < 9) {
0182       ndigiVsArea[id.roll() - 1] = ndigiVsArea[id.roll() - 1] + ndigiEtaPart[id.roll() - 1] * 1. / trArea;
0183       rollRadiusEtaPart[id.roll() - 1] = rollRadius;
0184     }
0185   }  // for eta partitions (rolls)
0186
0187   for (int i = 0; i <= 8; ++i) {
0188     LogDebug("ME0DigiReader") << "roll " << i + 1 << " numbers = " << countRoll[i] << "\tndigi = " << ndigiEtaPart[i]
0189                               << std::endl;
0190   }
0191
0192   numbEvents++;
0193 }
0194
0195 void ME0DigiReader::endJob() {
0196   LogDebug("ME0DigiReader") << "number of events = " << numbEvents << std::endl;
0197   LogDebug("ME0DigiReader") << "--------------" << std::endl;
0198   //  hRadiusEtaPartVsNdigiOvTrArea->GetXProjections();
0199
0200   std::vector<double> myRadii, myRates;
0201
0202   for (int i = 0; i <= 8; ++i) {
0203     LogDebug("ME0DigiReader") << "ndigiVsArea" << i + 1 << " = " << ndigiVsArea[i];
0204     ndigiVsArea[0] = ndigiVsArea[i] / (numbEvents * 9 * 25 * 1.0e-9 * 2. * 18. * 6. * 1.5);
0205     LogDebug("ME0DigiReader") << "\tRate [Hz/cm2] = " << ndigiVsArea[i] << std::endl;
0206     myRadii.push_back(rollRadiusEtaPart[i]);
0207     myRates.push_back(ndigiVsArea[i]);
0208   }
0209
0210   LogDebug("ME0DigiReader") << "rollRadius[cm]\tRate[Hz/cm2]" << std::endl;
0211   for (int i = 0; i <= 8; ++i) {
0212     LogDebug("ME0DigiReader") << rollRadiusEtaPart[i] << "\t" << ndigiVsArea[i] << std::endl;
0213   }
0214
0215   for (unsigned int i = 0; i < myRadii.size(); i++) {
0216     LogDebug("ME0DigiReader") << "radius = " << myRadii[i] << "\tRate = " << myRates[i] << std::endl;
0217     grRatePerRoll->SetPoint(i, myRadii[i], myRates[i]);
0218   }
0219 }
0220
0221 #include "FWCore/Framework/interface/MakerMacros.h"
0222 DEFINE_FWK_MODULE(ME0DigiReader);