CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
BtlSimHitsValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Validation/MtdValidation
4 // Class: BtlSimHitsValidation
5 //
14 #include <string>
15 
20 
23 
27 
31 
36 
39 
41 
42 struct MTDHit {
43  float energy;
44  float time;
45  float x;
46  float y;
47  float z;
48 };
49 
51 public:
52  explicit BtlSimHitsValidation(const edm::ParameterSet&);
53  ~BtlSimHitsValidation() override;
54 
55  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
56 
57 private:
58  void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
59 
60  void analyze(const edm::Event&, const edm::EventSetup&) override;
61 
62  // ------------ member data ------------
63 
65  const float hitMinEnergy_;
66 
68 
71 
72  // --- histograms declaration
73 
75 
78 
82 
86 
88 
94 
102 };
103 
104 // ------------ constructor and destructor --------------
106  : folder_(iConfig.getParameter<std::string>("folder")),
107  hitMinEnergy_(iConfig.getParameter<double>("hitMinimumEnergy")) {
108  btlSimHitsToken_ = consumes<CrossingFrame<PSimHit> >(iConfig.getParameter<edm::InputTag>("inputTag"));
109  mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
110  mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>();
111 }
112 
114 
115 // ------------ method called for each event ------------
117  using namespace edm;
118  using namespace geant_units::operators;
119 
120  auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
121  const MTDGeometry* geom = geometryHandle.product();
122 
123  auto topologyHandle = iSetup.getTransientHandle(mtdtopoToken_);
124  const MTDTopology* topology = topologyHandle.product();
125 
126  auto btlSimHitsHandle = makeValid(iEvent.getHandle(btlSimHitsToken_));
127  MixCollection<PSimHit> btlSimHits(btlSimHitsHandle.product());
128 
129  std::unordered_map<uint32_t, MTDHit> m_btlHits;
130  std::unordered_map<uint32_t, std::set<int> > m_btlTrkPerCell;
131 
132  // --- Loop over the BLT SIM hits
133  for (auto const& simHit : btlSimHits) {
134  // --- Use only hits compatible with the in-time bunch-crossing
135  if (simHit.tof() < 0 || simHit.tof() > 25.)
136  continue;
137 
138  DetId id = simHit.detUnitId();
139 
140  m_btlTrkPerCell[id.rawId()].insert(simHit.trackId());
141 
142  auto simHitIt = m_btlHits.emplace(id.rawId(), MTDHit()).first;
143 
144  // --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
145  (simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
146 
147  // --- Get the time of the first SIM hit in the cell
148  if ((simHitIt->second).time == 0 || simHit.tof() < (simHitIt->second).time) {
149  (simHitIt->second).time = simHit.tof();
150 
151  auto hit_pos = simHit.entryPoint();
152  (simHitIt->second).x = hit_pos.x();
153  (simHitIt->second).y = hit_pos.y();
154  (simHitIt->second).z = hit_pos.z();
155  }
156 
157  } // simHit loop
158 
159  // ==============================================================================
160  // Histogram filling
161  // ==============================================================================
162 
163  if (!m_btlHits.empty())
164  meNhits_->Fill(log10(m_btlHits.size()));
165 
166  for (auto const& hit : m_btlTrkPerCell)
167  meNtrkPerCell_->Fill((hit.second).size());
168 
169  for (auto const& hit : m_btlHits) {
170  meHitLogEnergy_->Fill(log10((hit.second).energy));
171 
172  if ((hit.second).energy < hitMinEnergy_)
173  continue;
174 
175  // --- Get the SIM hit global position
176  BTLDetId detId(hit.first);
177  DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
178  const MTDGeomDet* thedet = geom->idToDet(geoId);
179  if (thedet == nullptr)
180  throw cms::Exception("BtlSimHitsValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
181  << detId.rawId() << ") is invalid!" << std::dec << std::endl;
182  const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
183  const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
184 
185  Local3DPoint local_point(
186  convertMmToCm((hit.second).x), convertMmToCm((hit.second).y), convertMmToCm((hit.second).z));
187 
188  local_point = topo.pixelToModuleLocalPoint(local_point, detId.row(topo.nrows()), detId.column(topo.nrows()));
189  const auto& global_point = thedet->toGlobal(local_point);
190 
191  // --- Fill the histograms
192  meHitEnergy_->Fill((hit.second).energy);
193  meHitTime_->Fill((hit.second).time);
194 
195  meHitXlocal_->Fill((hit.second).x);
196  meHitYlocal_->Fill((hit.second).y);
197  meHitZlocal_->Fill((hit.second).z);
198 
199  meOccupancy_->Fill(global_point.z(), global_point.phi());
200 
201  meHitX_->Fill(global_point.x());
202  meHitY_->Fill(global_point.y());
203  meHitZ_->Fill(global_point.z());
204  meHitPhi_->Fill(global_point.phi());
205  meHitEta_->Fill(global_point.eta());
206 
207  meHitTvsE_->Fill((hit.second).energy, (hit.second).time);
208  meHitEvsPhi_->Fill(global_point.phi(), (hit.second).energy);
209  meHitEvsEta_->Fill(global_point.eta(), (hit.second).energy);
210  meHitEvsZ_->Fill(global_point.z(), (hit.second).energy);
211  meHitTvsPhi_->Fill(global_point.phi(), (hit.second).time);
212  meHitTvsEta_->Fill(global_point.eta(), (hit.second).time);
213  meHitTvsZ_->Fill(global_point.z(), (hit.second).time);
214 
215  } // hit loop
216 
217  // --- This is to count the number of processed events, needed in the harvesting step
218  meNevents_->Fill(0.5);
219 }
220 
221 // ------------ method for histogram booking ------------
223  edm::Run const& run,
224  edm::EventSetup const& iSetup) {
225  ibook.setCurrentFolder(folder_);
226 
227  // --- histograms booking
228 
229  meNevents_ = ibook.book1D("BtlNevents", "Number of events", 1, 0., 1.);
230 
231  meNhits_ = ibook.book1D("BtlNhits", "Number of BTL cells with SIM hits;log_{10}(N_{BTL cells})", 100, 0., 5.25);
232  meNtrkPerCell_ = ibook.book1D("BtlNtrkPerCell", "Number of tracks per BTL cell;N_{trk}", 10, 0., 10.);
233 
234  meHitEnergy_ = ibook.book1D("BtlHitEnergy", "BTL SIM hits energy;E_{SIM} [MeV]", 100, 0., 20.);
235  meHitLogEnergy_ = ibook.book1D("BtlHitLogEnergy", "BTL SIM hits energy;log_{10}(E_{SIM} [MeV])", 200, -6., 3.);
236  meHitTime_ = ibook.book1D("BtlHitTime", "BTL SIM hits ToA;ToA_{SIM} [ns]", 100, 0., 25.);
237 
238  meHitXlocal_ = ibook.book1D("BtlHitXlocal", "BTL SIM local X;X_{SIM}^{LOC} [mm]", 100, -30., 30.);
239  meHitYlocal_ = ibook.book1D("BtlHitYlocal", "BTL SIM local Y;Y_{SIM}^{LOC} [mm]", 100, -1.65, 1.65);
240  meHitZlocal_ = ibook.book1D("BtlHitZlocal", "BTL SIM local z;z_{SIM}^{LOC} [mm]", 100, -2., 2.);
241 
242  meOccupancy_ = ibook.book2D(
243  "BtlOccupancy", "BTL SIM hits occupancy;z_{SIM} [cm];#phi_{SIM} [rad]", 130, -260., 260., 200, -3.15, 3.15);
244 
245  meHitX_ = ibook.book1D("BtlHitX", "BTL SIM hits X;X_{SIM} [cm]", 100, -120., 120.);
246  meHitY_ = ibook.book1D("BtlHitY", "BTL SIM hits Y;Y_{SIM} [cm]", 100, -120., 120.);
247  meHitZ_ = ibook.book1D("BtlHitZ", "BTL SIM hits Z;Z_{SIM} [cm]", 100, -260., 260.);
248  meHitPhi_ = ibook.book1D("BtlHitPhi", "BTL SIM hits #phi;#phi_{SIM} [rad]", 200, -3.15, 3.15);
249  meHitEta_ = ibook.book1D("BtlHitEta", "BTL SIM hits #eta;#eta_{SIM}", 100, -1.55, 1.55);
250 
251  meHitTvsE_ =
252  ibook.bookProfile("BtlHitTvsE", "BTL SIM time vs energy;E_{SIM} [MeV];T_{SIM} [ns]", 50, 0., 20., 0., 100.);
253  meHitEvsPhi_ = ibook.bookProfile(
254  "BtlHitEvsPhi", "BTL SIM energy vs #phi;#phi_{SIM} [rad];E_{SIM} [MeV]", 50, -3.15, 3.15, 0., 100.);
255  meHitEvsEta_ =
256  ibook.bookProfile("BtlHitEvsEta", "BTL SIM energy vs #eta;#eta_{SIM};E_{SIM} [MeV]", 50, -1.55, 1.55, 0., 100.);
257  meHitEvsZ_ =
258  ibook.bookProfile("BtlHitEvsZ", "BTL SIM energy vs Z;Z_{SIM} [cm];E_{SIM} [MeV]", 50, -260., 260., 0., 100.);
259  meHitTvsPhi_ = ibook.bookProfile(
260  "BtlHitTvsPhi", "BTL SIM time vs #phi;#phi_{SIM} [rad];T_{SIM} [ns]", 50, -3.15, 3.15, 0., 100.);
261  meHitTvsEta_ =
262  ibook.bookProfile("BtlHitTvsEta", "BTL SIM time vs #eta;#eta_{SIM};T_{SIM} [ns]", 50, -1.55, 1.55, 0., 100.);
263  meHitTvsZ_ =
264  ibook.bookProfile("BtlHitTvsZ", "BTL SIM time vs Z;Z_{SIM} [cm];T_{SIM} [ns]", 50, -260., 260., 0., 100.);
265 }
266 
267 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
270 
271  desc.add<std::string>("folder", "MTD/BTL/SimHits");
272  desc.add<edm::InputTag>("inputTag", edm::InputTag("mix", "g4SimHitsFastTimerHitsBarrel"));
273  desc.add<double>("hitMinimumEnergy", 1.); // [MeV]
274 
275  descriptions.add("btlSimHits", desc);
276 }
277 
MonitorElement * meHitEta_
MonitorElement * meHitZlocal_
MonitorElement * meHitPhi_
void analyze(const edm::Event &, const edm::EventSetup &) override
std::string folder_
MonitorElement * meHitXlocal_
MonitorElement * meNtrkPerCell_
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:32
virtual const Topology & topology() const
Definition: GeomDet.cc:67
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
constexpr NumType convertUnitsTo(double desiredUnits, NumType val)
Definition: GeantUnits.h:100
void Fill(long long x)
MonitorElement * meHitEvsPhi_
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Definition: Event.h:563
int iEvent
Definition: GenABIO.cc:224
BTLDetId geographicalId(CrysLayout lay) const
Definition: BTLDetId.cc:171
const std::string folder_
MonitorElement * meHitEnergy_
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:322
MonitorElement * meOccupancy_
virtual const PixelTopology & specificTopology() const
MonitorElement * meHitEvsEta_
int nrows() const override
edm::ESGetToken< MTDGeometry, MTDDigiGeometryRecord > mtdgeoToken_
ParameterDescriptionBase * add(U const &iLabel, T const &value)
LocalPoint pixelToModuleLocalPoint(const LocalPoint &plp, int row, int col) const
MonitorElement * meHitEvsZ_
MonitorElement * meHitYlocal_
MonitorElement * meHitTvsPhi_
MonitorElement * meNevents_
MonitorElement * meHitTvsE_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: DetId.h:17
edm::ESGetToken< MTDTopology, MTDTopologyRcd > mtdtopoToken_
MonitorElement * meHitLogEnergy_
BtlSimHitsValidation(const edm::ParameterSet &)
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:177
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
void add(std::string const &label, ParameterSetDescription const &psetDescription)
MonitorElement * meHitTime_
edm::EDGetTokenT< CrossingFrame< PSimHit > > btlSimHitsToken_
ESTransientHandle< T > getTransientHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:168
Detector identifier class for the Barrel Timing Layer. The crystal count must start from 0...
Definition: BTLDetId.h:18
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
constexpr NumType convertMmToCm(NumType millimeters)
Definition: GeantUnits.h:63
int column(unsigned nrows=16) const
Definition: BTLDetId.h:110
BTLDetId::CrysLayout crysLayoutFromTopoMode(const int &topoMode)
MonitorElement * meHitTvsEta_
auto makeValid(const U &iOtherHandleType) noexcept(false)
Definition: ValidHandle.h:52
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
Definition: Run.h:45
int row(unsigned nrows=16) const
Definition: BTLDetId.h:105