44 #include "CLHEP/Geometry/Point3D.h"
45 #include "CLHEP/Geometry/Transform3D.h"
46 #include "CLHEP/Geometry/Vector3D.h"
47 #include "CLHEP/Units/GlobalSystemOfUnits.h"
48 #include "CLHEP/Units/GlobalPhysicalConstants.h"
55 for (
int i = 0;
i < 6; ++
i)
83 void fillHitsInfo(std::pair<hitsinfo, energysum> hit_,
unsigned int itimeslice,
double esum);
112 : nameDetector_(iConfig.getParameter<std::
string>(
"DetectorName")),
113 caloHitSource_(iConfig.getParameter<std::
string>(
"CaloHitSource")),
114 times_(iConfig.getParameter<std::
vector<double> >(
"TimeSlices")),
115 verbosity_(iConfig.getUntrackedParameter<int>(
"Verbosity", 0)),
116 fromDDD_(iConfig.getUntrackedParameter<bool>(
"fromDDD",
true)),
119 tok_cpv_(esConsumes<DDCompactView, IdealGeometryRecord, edm::Transition::BeginRun>()),
120 tok_cpvc_(esConsumes<cms::DDCompactView, IdealGeometryRecord, edm::Transition::BeginRun>()),
123 tok_hepMC_ = consumes<edm::HepMCProduct>(
edm::InputTag(
"generatorSmeared"));
124 tok_hits_ = consumes<edm::PCaloHitContainer>(
edm::InputTag(
"g4SimHits", caloHitSource_));
125 nTimes_ = (times_.size() > maxTime_) ? maxTime_ : times_.size();
130 std::vector<double> times = {25.0, 1000.0};
133 desc.
add<std::vector<double> >(
"TimeSlices", times);
137 descriptions.
add(
"hgcalSimHitValidationEE", desc);
150 for (HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end();
152 edm::LogVerbatim(
"HGCalValidation") <<
"Particle[" << k <<
"] with pt " << (*p)->momentum().perp() <<
" eta "
153 << (*p)->momentum().eta() <<
" phi " << (*p)->momentum().phi();
160 iEvent.
getByToken(tok_hits_, theCaloHitContainers);
161 if (theCaloHitContainers.
isValid()) {
163 edm::LogVerbatim(
"HGCalValidation") <<
" PcalohitItr = " << theCaloHitContainers->size();
164 std::vector<PCaloHit> caloHits;
165 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
173 std::map<int, int> OccupancyMap_plus, OccupancyMap_minus;
174 OccupancyMap_plus.clear();
175 OccupancyMap_minus.clear();
177 std::map<uint32_t, std::pair<hitsinfo, energysum> > map_hits;
182 unsigned int nused(0);
183 for (
unsigned int i = 0;
i < hits.size();
i++) {
184 double energy = hits[
i].energy();
185 double time = hits[
i].time();
186 uint32_t id_ = hits[
i].id();
188 int cell2(0),
type(0);
191 cell = detId.
cellU();
192 cell2 = detId.
cellV();
194 subsector = detId.
waferV();
196 layer = detId.
layer();
197 zside = detId.
zside();
204 layer = detId.
layer();
205 zside = detId.
zside();
213 <<
"Detector " <<
nameDetector_ <<
" zside = " << zside <<
" sector|wafer = " << sector <<
":" << subsector
214 <<
" type = " << type <<
" layer = " << layer <<
" cell = " << cell <<
":" << cell2 <<
" energy = " << energy
215 <<
" energyem = " << hits[
i].energyEM() <<
" energyhad = " << hits[
i].energyHad() <<
" time = " << time;
217 HepGeom::Point3D<float> gcoord;
218 std::pair<float, float>
xy;
229 float xp = (zp < 0) ? -xy.first : xy.first;
230 gcoord = HepGeom::Point3D<float>(xp, xy.second, zp);
231 double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
234 << std::hex << id_ <<
std::dec <<
" global coordinate " << gcoord <<
" time " << time <<
":" << tof;
239 if (map_hits.count(id_) != 0) {
240 hinfo = map_hits[id_].first;
241 esum = map_hits[id_].second;
243 hinfo.
x = gcoord.x();
244 hinfo.
y = gcoord.y();
245 hinfo.
z = gcoord.z();
252 hinfo.
phi = gcoord.getPhi();
253 hinfo.
eta = gcoord.getEta();
257 if (time > 0 && time <
times_[
k])
262 edm::LogVerbatim(
"HGCalValidation") <<
" ----------------------- gx = " << hinfo.
x <<
" gy = " << hinfo.
y
263 <<
" gz = " << hinfo.
z <<
" phi = " << hinfo.
phi <<
" eta = " << hinfo.
eta;
264 map_hits[id_] = std::pair<hitsinfo, energysum>(hinfo, esum);
268 <<
" detector elements being hit";
270 std::map<uint32_t, std::pair<hitsinfo, energysum> >::iterator itr;
271 for (itr = map_hits.begin(); itr != map_hits.end(); ++itr) {
277 for (
unsigned int itimeslice = 0; itimeslice <
nTimes_; itimeslice++) {
287 edm::LogVerbatim(
"HGCalValidation") <<
"With map:used:total " << hits.size() <<
"|" << nused <<
"|"
288 << map_hits.size() <<
" hits";
290 for (
auto const& itr : OccupancyMap_plus) {
291 int layer = itr.first;
292 int occupancy = itr.second;
295 for (
auto const& itr : OccupancyMap_minus) {
296 int layer = itr.first;
297 int occupancy = itr.second;
303 if (OccupancyMap.find(layer) != OccupancyMap.end()) {
304 ++OccupancyMap[
layer];
306 OccupancyMap[
layer] = 1;
311 unsigned int ilayer = hits.first.layer;
313 energy_[itimeslice].at(ilayer)->Fill(esum);
314 if (itimeslice == 0) {
315 EtaPhi_Plus_.at(ilayer)->Fill(hits.first.eta, hits.first.phi);
316 EtaPhi_Minus_.at(ilayer)->Fill(hits.first.eta, hits.first.phi);
321 <<
"Problematic Hit for " <<
nameDetector_ <<
" at sector " << hits.first.sector <<
":" << hits.first.sector2
322 <<
" layer " << hits.first.layer <<
" cell " << hits.first.cell <<
":" << hits.first.cell2 <<
" energy "
323 << hits.second.etotal;
344 int isd = (name.find(
nameDetector_) == std::string::npos) ? -1 : 1;
347 int nsiz =
static_cast<int>(copy.size());
348 int lay = (nsiz > 0) ? copy[nsiz - 1] : -1;
349 int sec = (nsiz > 1) ? copy[nsiz - 2] : -1;
350 int zp = (nsiz > 3) ? copy[nsiz - 4] : -1;
354 int subs = (trp.
alpha1() > 0 ? 1 : 0);
358 fv.
rotation().GetComponents(x, y, z);
359 const CLHEP::HepRep3x3
rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
360 const CLHEP::HepRotation hr(
rotation);
362 const HepGeom::Transform3D ht3d(hr, h3v);
363 transMap_.insert(std::make_pair(
id, ht3d));
390 int nsiz =
static_cast<int>(
copy.size());
391 int lay = (nsiz > 0) ?
copy[0] : -1;
392 int sec = (nsiz > 1) ?
copy[1] : -1;
393 int zp = (nsiz > 3) ?
copy[3] : -1;
397 int subs = (pars[6] > 0 ? 1 : 0);
398 symmDet_ = (pars[6] == 0 ?
true :
false);
401 fv.
rotation().GetComponents(x, y, z);
402 const CLHEP::HepRep3x3
rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
403 const CLHEP::HepRotation hr(
rotation);
405 const HepGeom::Transform3D ht3d(hr, h3v);
406 transMap_.insert(std::make_pair(
id, ht3d));
437 std::ostringstream histoname;
438 for (
unsigned int il = 0; il <
layers_; ++il) {
440 auto istr1 = std::to_string(ilayer);
441 while (istr1.size() < 2) {
442 istr1.insert(0,
"0");
445 histoname <<
"HitOccupancy_Plus_layer_" << istr1;
448 histoname <<
"HitOccupancy_Minus_layer_" << istr1;
452 histoname <<
"EtaPhi_Plus_"
453 <<
"layer_" << istr1;
456 histoname <<
"EtaPhi_Minus_"
457 <<
"layer_" << istr1;
461 for (
unsigned int itimeslice = 0; itimeslice <
nTimes_; itimeslice++) {
463 histoname <<
"energy_time_" << itimeslice <<
"_layer_" << istr1;
464 energy_[itimeslice].push_back(iB.
book1D(histoname.str().c_str(),
"energy_", 100, 0, 0.1));
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
MonitorElement * MeanHitOccupancy_Plus_
Log< level::Info, true > LogVerbatim
bool defineGeometry(const DDCompactView *ddViewH)
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >> DD3Vector
void fillOccupancyMap(std::map< int, int > &OccupancyMap, int layer)
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
std::vector< MonitorElement * > HitOccupancy_Plus_
void analyze(const edm::Event &, const edm::EventSetup &) override
const std::vector< double > times_
virtual void setCurrentFolder(std::string const &fullpath)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
nav_type copyNumbers() const
return the stack of copy numbers
const RotationMatrix rotation() const
#define DEFINE_FWK_MODULE(type)
edm::EDGetTokenT< edm::HepMCProduct > tok_hepMC_
const DDRotationMatrix & rotation() const
The absolute rotation of the current node.
int type() const
get/set the type
void analyzeHits(std::vector< PCaloHit > &hits)
int zside() const
get the z-side of the cell (1/-1)
std::vector< MonitorElement * > HitOccupancy_Minus_
const DDSolid & solid(void) const
Returns a reference object of the solid being the shape of this LogicalPart.
std::pair< float, float > locateCell(int cell, int lay, int type, bool reco) const
Compact representation of the geometrical detector hierarchy.
int cellU() const
get the cell #'s in u,v or in x,y
A DDSolid represents the shape of a part.
MonitorElement * MeanHitOccupancy_Minus_
constexpr std::array< uint8_t, layerIndexSize > layer
const Translation translation() const
bool getData(T &iHolder) const
const std::string nameDetector_
unsigned int layers(bool reco) const
bool next()
set current node to the next node in the filtered tree
int type() const
get the type
int layer() const
get the layer #
std::pair< float, float > locateCellTrap(int lay, int ieta, int iphi, bool reco) const
static const unsigned int maxTime_
bool tileTrapezoid() const
void dqmBeginRun(const edm::Run &, const edm::EventSetup &) override
const std::string fullname() const
Interface to a Trapezoid.
~HGCalSimHitValidation() override
int iphi() const
get the phi index
std::vector< MonitorElement * > EtaPhi_Plus_
std::string_view name() const
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool firstChild()
set the current node to the first child
HGCalSimHitValidation(const edm::ParameterSet &)
double waferZ(int layer, bool reco) const
const std::vector< int > copyNos() const
The list of the volume copy numbers.
Basic2DVector< T > xy() const
const std::string caloHitSource_
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
double alpha1(void) const
Angle with respect to the y axis from the centre of the side at y=-pDy1 to the centre at y=+pDy1 of t...
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hits_
void add(std::string const &label, ParameterSetDescription const &psetDescription)
std::vector< MonitorElement * > energy_[maxTime_]
int layer() const
get the layer #
int zside() const
get the z-side of the cell (1/-1)
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
std::map< uint32_t, HepGeom::Transform3D > transMap_
bool firstChild()
set the current node to the first child ...
const edm::ESGetToken< cms::DDCompactView, IdealGeometryRecord > tok_cpvc_
bool waferHexagon8() const
const HGCalDDDConstants * hgcons_
const std::vector< double > parameters() const
extract shape parameters
std::vector< MonitorElement * > EtaPhi_Minus_
const DDTranslation & translation() const
The absolute translation of the current node.
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
static uint32_t packSquareIndex(int z, int lay, int sec, int subsec, int cell)
const edm::ESGetToken< DDCompactView, IdealGeometryRecord > tok_cpv_
void fillHitsInfo(std::pair< hitsinfo, energysum > hit_, unsigned int itimeslice, double esum)
const edm::ESGetToken< HGCalDDDConstants, IdealGeometryRecord > tok_hgcal_
static void unpackHexagonIndex(const uint32_t &idx, int &subdet, int &z, int &lay, int &wafer, int &celltyp, int &cell)