48 #include "TProfile2D.h"
74 static bool sortTime(
const std::pair<double, double>&
i,
const std::pair<double, double>&
j);
141 : ifEE_(iConfig.getParameter<
bool>(
"useEE")),
142 ifFH_(iConfig.getParameter<
bool>(
"useFH")),
143 ifBH_(iConfig.getParameter<
bool>(
"useBH")),
144 ifBeam_(iConfig.getParameter<
bool>(
"useBeam")),
145 doSimHits_(iConfig.getParameter<
bool>(
"doSimHits")),
146 doDigis_(iConfig.getParameter<
bool>(
"doDigis")),
147 doRecHits_(iConfig.getParameter<
bool>(
"doRecHits")),
148 doTree_(iConfig.getParameter<
bool>(
"doTree")),
149 doTreeCell_(iConfig.getParameter<
bool>(
"doTreeCell")),
150 doPassive_(iConfig.getParameter<
bool>(
"doPassive")),
151 doPassiveEE_(iConfig.getParameter<
bool>(
"doPassiveEE")),
152 doPassiveHE_(iConfig.getParameter<
bool>(
"doPassiveHE")),
153 doPassiveBH_(iConfig.getParameter<
bool>(
"doPassiveBH")),
154 addP_(iConfig.getParameter<
bool>(
"addP")),
155 doBeam_(iConfig.getParameter<
bool>(
"doBeam")),
156 detectorEE_(iConfig.getParameter<
std::
string>(
"detectorEE")),
157 detectorFH_(iConfig.getParameter<
std::
string>(
"detectorFH")),
158 detectorBH_(iConfig.getParameter<
std::
string>(
"detectorBH")),
159 detectorBeam_(iConfig.getParameter<
std::
string>(
"detectorBeam")),
160 zFrontEE_(iConfig.getParameter<double>(
"zFrontEE")),
161 zFrontFH_(iConfig.getParameter<double>(
"zFrontFH")),
162 zFrontBH_(iConfig.getParameter<double>(
"zFrontBH")),
163 sampleIndex_(iConfig.getParameter<
int>(
"sampleIndex")),
164 gev2mip200_(iConfig.getUntrackedParameter<double>(
"gev2mip200", 57.0
e-6)),
165 gev2mip300_(iConfig.getUntrackedParameter<double>(
"gev2mip300", 85.5
e-6)),
166 stoc_smear_time_200_(iConfig.getUntrackedParameter<double>(
"stoc_smear_time_200", 10.24)),
167 stoc_smear_time_300_(iConfig.getUntrackedParameter<double>(
"stoc_smear_time_300", 15.5)) {
168 usesResource(
"TFileService");
169 ahcalGeom_ = std::make_unique<AHCalGeometry>(iConfig);
195 tok_hepMC_ = consumes<edm::HepMCProduct>(tmp0);
198 edm::LogVerbatim(
"HGCSim") <<
"HGCalTBAnalyzer:: GeneratorSource = " << tmp0;
266 desc.add<
bool>(
"useEE",
true);
267 desc.add<
double>(
"zFrontEE", 0.0);
272 desc.add<
bool>(
"useFH",
false);
273 desc.add<
double>(
"zFrontFH", 0.0);
278 desc.add<
bool>(
"useBH",
false);
279 desc.add<
double>(
"zFrontBH", 0.0);
284 desc.add<
bool>(
"useBeam",
false);
286 std::vector<int> ids = {
287 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1011, 1012, 1013, 1014, 2001, 2002, 2003, 2004, 2005};
288 desc.add<std::vector<int>>(
"idBeams", ids);
296 desc.add<
bool>(
"doSimHits",
true);
297 desc.add<
bool>(
"doDigis",
true);
298 desc.add<
bool>(
"doRecHits",
true);
299 desc.add<
int>(
"sampleIndex", 0);
300 desc.add<
bool>(
"doTree",
true);
301 desc.add<
bool>(
"doTreeCell",
true);
302 desc.add<
bool>(
"doPassive",
false);
303 desc.add<
bool>(
"doPassiveEE",
false);
304 desc.add<
bool>(
"doPassiveHE",
false);
305 desc.add<
bool>(
"doPassiveBH",
false);
306 desc.add<
bool>(
"addP",
false);
307 desc.add<
bool>(
"doBeam",
false);
308 desc.addUntracked<
double>(
"gev2mip200", 57.0e-6);
309 desc.addUntracked<
double>(
"gev2mip300", 85.5e-6);
310 desc.addUntracked<
double>(
"stoc_smear_time_200", 10.24);
311 desc.addUntracked<
double>(
"stoc_smear_time_300", 15.5);
312 desc.addUntracked<
int>(
"maxDepth", 12);
313 desc.addUntracked<
double>(
"deltaX", 30.0);
314 desc.addUntracked<
double>(
"deltaY", 30.0);
315 desc.addUntracked<
double>(
"deltaZ", 81.0);
316 desc.addUntracked<
double>(
"zFirst", 17.6);
317 descriptions.
add(
"HGCalTBAnalyzer",
desc);
322 hBeam_ =
fs_->
make<TH1D>(
"BeamP",
"Beam Momentum", 1000, 0, 1000.0);
323 for (
int i = 0;
i < 3; ++
i) {
335 sprintf(
name,
"SimHitEn%s", det.c_str());
336 sprintf(
title,
"Sim Hit Energy for %s", det.c_str());
338 sprintf(
name,
"SimHitEnX%s", det.c_str());
339 sprintf(
title,
"Sim Hit Energy for %s", det.c_str());
341 sprintf(
name,
"SimHitTm%s", det.c_str());
342 sprintf(
title,
"Sim Hit Timing for %s", det.c_str());
344 sprintf(
name,
"SimHitLat%s", det.c_str());
345 sprintf(
title,
"Lateral Shower profile (Sim Hit) for %s", det.c_str());
347 sprintf(
name,
"SimHitLng%s", det.c_str());
348 sprintf(
title,
"Longitudinal Shower profile (Sim Hit) for %s", det.c_str());
350 sprintf(
name,
"SimHitLng1%s", det.c_str());
351 sprintf(
title,
"Longitudinal Shower profile (Layer) for %s", det.c_str());
353 sprintf(
name,
"SimHitLng2%s", det.c_str());
354 sprintf(
title,
"Longitudinal Shower profile (Layer) for %s", det.c_str());
359 sprintf(
name,
"DigiADC%s", det.c_str());
360 sprintf(
title,
"ADC at Digi level for %s", det.c_str());
362 sprintf(
name,
"DigiOcc%s", det.c_str());
363 sprintf(
title,
"Occupancy (Digi)for %s", det.c_str());
365 sprintf(
name,
"DigiLng%s", det.c_str());
366 sprintf(
title,
"Longitudinal Shower profile (Digi) for %s", det.c_str());
371 sprintf(
name,
"RecHitEn%s", det.c_str());
372 sprintf(
title,
"Rec Hit Energy for %s", det.c_str());
374 sprintf(
name,
"RecHitOcc%s", det.c_str());
375 sprintf(
title,
"Occupancy (Rec Hit)for %s", det.c_str());
377 sprintf(
name,
"RecHitLat%s", det.c_str());
378 sprintf(
title,
"Lateral Shower profile (Rec Hit) for %s", det.c_str());
380 sprintf(
name,
"RecHitLng%s", det.c_str());
381 sprintf(
title,
"Longitudinal Shower profile (Rec Hit) for %s", det.c_str());
383 sprintf(
name,
"RecHitLng1%s", det.c_str());
384 sprintf(
title,
"Longitudinal Shower profile vs Layer for %s", det.c_str());
489 sprintf(
title,
"Sim Hit Energy in SIM layer %d for %s",
l + 1,
detectorEE_.c_str());
493 sprintf(
title,
"Sim Hit Energy in layer %d for %s", (
l / 3 + 1),
detectorEE_.c_str());
523 sprintf(
title,
"Sim Hit Energy in layer %d for %s", (
l / 3 + 1),
detectorFH_.c_str());
548 for (
unsigned int l = 0;
l <
idBeams_.size(); ++
l) {
565 HepMC::FourVector pxyz(0, 0, 0, 0);
566 for (HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end();
568 edm::LogVerbatim(
"HGCSim") <<
"Particle [" <<
k <<
"] with p " << (*p)->momentum().rho() <<
" theta "
569 << (*p)->momentum().theta() <<
" phi " << (*p)->momentum().phi() <<
" pxyz ("
570 << (*p)->momentum().px() <<
", " << (*p)->momentum().py() <<
", "
571 << (*p)->momentum().pz() <<
")";
573 pxyz.setPx(pxyz.px() + (*p)->momentum().px());
574 pxyz.setPy(pxyz.py() + (*p)->momentum().py());
575 pxyz.setPz(pxyz.pz() + (*p)->momentum().pz());
576 pxyz.setE(pxyz.e() + (*p)->momentum().e());
577 }
else if (!
addP_ && (
k == 0)) {
578 pxyz = (*p)->momentum();
582 edm::LogVerbatim(
"HGCSim") <<
"Particle with p " << pxyz.rho() <<
" theta " << pxyz.theta() <<
" phi "
619 std::vector<PCaloHit> caloHits;
624 if (theCaloHitContainers.
isValid()) {
630 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
642 if (theCaloHitContainers.
isValid()) {
648 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
660 if (theCaloHitContainers.
isValid()) {
666 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
677 if (theCaloHitContainers.
isValid()) {
680 << theCaloHitContainers->size() <<
" hits";
683 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
746 if (theDigiContainers.
isValid()) {
751 for (
const auto& it : *theDigiContainers) {
754 uint16_t
adc = hgcSample.
data();
762 if (theDigiContainers.
isValid()) {
767 for (
const auto& it : *theDigiContainers) {
770 uint16_t
adc = hgcSample.
data();
782 if (theCaloHitContainers.
isValid()) {
785 << theCaloHitContainers->
size() <<
" hits";
796 if (theCaloHitContainers.
isValid()) {
799 << theCaloHitContainers->
size() <<
" hits";
813 std::map<uint32_t, double> map_hits, map_hitn;
814 std::map<uint32_t, double> map_hittime_firsthit, map_hittime_lasthit, map_hittime_15Mip;
815 std::map<int, double> map_hitDepth;
816 std::map<int, std::pair<uint32_t, double>> map_hitLayer, map_hitCell;
818 std::map<uint32_t, double>
nhits;
819 std::map<uint32_t, int>
ID,
Depth;
820 std::map<uint32_t, double> GeV2Mip;
821 std::map<uint32_t, double> StochTermTime;
822 std::vector<int> nSimLayers;
824 std::map<uint32_t, std::vector<std::pair<double, double>>> map_hitTimeEn;
827 for (
unsigned int i = 0;
i <
hits.size();
i++) {
830 uint32_t
id =
hits[
i].id();
832 int subdet,
zside, layer, sector, subsector(0), cell,
depth(0),
idx(0);
845 }
else if (
type == 3) {
849 idx = subdet * 1000 + layer;
854 idx = sector * 1000 + cell;
858 << sector <<
":" << cell <<
" Position " <<
xy.first <<
":" <<
xy.second <<
":"
863 edm::LogVerbatim(
"HGCSim") <<
"SimHit:Hit[" <<
i <<
"] Id " << subdet <<
":" <<
zside <<
":" << layer <<
":"
864 << sector <<
":" << subsector <<
":" << cell <<
":" <<
depth <<
" Energy " <<
energy
867 if (map_hits.count(
id) != 0) {
872 if (map_hitLayer.count(layer) != 0) {
873 double ee =
energy + map_hitLayer[layer].second;
874 map_hitLayer[layer] = std::make_pair(
id, ee);
876 map_hitLayer[layer] = std::make_pair(
id,
energy);
879 if (map_hitCell.count(
idx) != 0) {
880 double ee =
energy + map_hitCell[
idx].second;
881 map_hitCell[
idx] = std::make_pair(
id, ee);
883 map_hitCell[
idx] = std::make_pair(
id,
energy);
891 if (map_hitDepth.count(
depth) != 0) {
898 map_hitTimeEn[idn].push_back(std::make_pair(
time,
energy));
903 if (map_hitn.count(idn) != 0) {
916 for (
const auto&
itr : map_hitTimeEn) {
917 uint32_t
id =
itr.first;
925 <<
"\ntype : layer : wafer thickness " <<
type <<
" " << layer <<
" " <<
thickness
926 <<
"\nID(sim) and id(reco) " << std::hex <<
ID[
id] <<
" " <<
id <<
std::dec;
935 std::sort(map_hitTimeEn[
id].begin(), map_hitTimeEn[
id].
end(),
sortTime);
939 double totEbeforeThreshold = 0.;
940 double timebeforeThreshold = 0.;
941 double timeAtThresohld = 0.;
942 for (
unsigned int ihit = 0; ihit < map_hitTimeEn[
id].size(); ihit++) {
943 double energy = (map_hitTimeEn[
id].at(ihit)).
second / GeV2Mip[
id];
947 edm::LogVerbatim(
"HGCSim") <<
"Tot E till now : time of that E : GeV2Mip[id] is " << totE <<
" " <<
time
948 <<
" " << GeV2Mip[
id];
950 totEbeforeThreshold = totE;
951 timebeforeThreshold =
time;
954 (
threshold - totEbeforeThreshold) * (
time - timebeforeThreshold) / (totE - totEbeforeThreshold) +
956 map_hittime_15Mip[
id] = timeAtThresohld;
958 edm::LogVerbatim(
"HGCSim") <<
"ihit : energyBefore : timeBefore : energyTot : timeTot : timeAt15MIP "
959 << ihit <<
" " << totEbeforeThreshold <<
" " << timebeforeThreshold <<
" "
960 << totE <<
" " <<
time <<
" " << map_hittime_15Mip[
id];
964 if (!map_hitTimeEn[
id].
empty()) {
965 map_hittime_firsthit[
id] = (map_hitTimeEn[
id].at(0)).
first;
966 map_hittime_lasthit[
id] = (map_hitTimeEn[
id].at(map_hitTimeEn[
id].
size() - 1)).
first;
967 if (map_hittime_15Mip[
id] < map_hittime_firsthit[
id])
968 map_hittime_15Mip[
id] = map_hittime_firsthit[
id];
978 map_hittime_15Mip[
id] = -99;
980 edm::LogVerbatim(
"HGCSim") <<
"id : first hit time : last hit time " <<
id <<
" " << map_hittime_firsthit[
id]
981 <<
" " << map_hittime_lasthit[
id] <<
"\nFinally for this cell, time is "
982 << map_hittime_15Mip[
id];
987 for (
const auto&
itr : map_hits) {
993 for (
const auto&
itr : map_hitLayer) {
994 int layer = (
type == 2) ?
itr.first : (
itr.first - 1);
1002 edm::LogVerbatim(
"HGCSim") <<
"SimHit:Layer " << layer + 1 <<
" Z " << zp <<
":" << zp - zFront <<
" E " <<
energy;
1013 }
else if (
type == 1) {
1018 }
else if (
type == 2) {
1026 for (
unsigned int k = 0;
k <
idBeams_.size(); ++
k) {
1035 for (
const auto&
itr : map_hitDepth) {
1036 int layer = (
type == 2) ?
itr.first : (
itr.first - 1);
1047 }
else if (
type == 1) {
1052 }
else if (
type == 2) {
1065 for (
const auto&
itr : map_hitCell) {
1066 uint32_t
id = ((
itr.second).
first);
1068 std::pair<float, float>
xy(0, 0);
1074 int subdet,
zside, layer, sector, subsector, cell;
1078 xx = (zp < 0) ? -
xy.first :
xy.first;
1084 for (
const auto&
itr : map_hitn) {
1085 uint32_t
id =
itr.first;
1088 double time_firsthit = map_hittime_firsthit[
id];
1089 double time15Mip = map_hittime_15Mip[
id];
1090 double time_lasthit = map_hittime_lasthit[
id];
1096 if (
debug && (
energy / GeV2Mip[
id] < 15) && (map_hittime_15Mip[
id] > 0))
1098 <<
" " << map_hittime_15Mip[
id];
1099 }
else if (
type == 1) {
1100 double time_firsthit = map_hittime_firsthit[
id];
1101 double time15Mip = map_hittime_15Mip[
id];
1102 double time_lasthit = map_hittime_lasthit[
id];
1108 }
else if (
type == 2) {
1112 int row = hid.
irow();
1114 int layer = hid.
depth();
1122 }
else if (
type == 3) {
1145 std::vector<float> verX, verY, verZ;
1149 for (
const auto& simVtxItr : *SimVtx) {
1150 verX.push_back(simVtxItr.position().X());
1151 verY.push_back(simVtxItr.position().Y());
1152 verZ.push_back(simVtxItr.position().Z());
1157 HepMC::FourVector pxyz(0, 0, 0, 0);
1158 for (
const auto& simTrkItr : *SimTk) {
1159 if (
addP_ && !(simTrkItr.noGenpart())) {
1160 pxyz.setPx(pxyz.px() + simTrkItr.momentum().px());
1161 pxyz.setPy(pxyz.py() + simTrkItr.momentum().py());
1162 pxyz.setPz(pxyz.pz() + simTrkItr.momentum().pz());
1163 pxyz.setE(pxyz.e() + simTrkItr.momentum().e());
1165 edm::LogVerbatim(
"HGCSim") <<
"Track " << simTrkItr.trackId() <<
" Vertex " << simTrkItr.vertIndex() <<
" Type "
1166 << simTrkItr.type() <<
" Charge " << simTrkItr.charge() <<
" px "
1167 << simTrkItr.momentum().px() <<
" py " << simTrkItr.momentum().py() <<
" pz "
1168 << simTrkItr.momentum().pz() <<
" P " << simTrkItr.momentum().P() <<
" GenIndex "
1169 << simTrkItr.genpartIndex();
1171 <<
" position-> X: " << verX[simTrkItr.vertIndex()]
1172 <<
" Y: " << verY[simTrkItr.vertIndex()] <<
" Z: " << verZ[simTrkItr.vertIndex()];
1175 if (
doBeam_ && !(simTrkItr.noGenpart())) {
1178 xBeamMC_.push_back(verX[simTrkItr.vertIndex()]);
1179 yBeamMC_.push_back(verY[simTrkItr.vertIndex()]);
1180 zBeamMC_.push_back(verZ[simTrkItr.vertIndex()]);
1181 pxBeamMC_.push_back(simTrkItr.momentum().px());
1182 pyBeamMC_.push_back(simTrkItr.momentum().py());
1183 pzBeamMC_.push_back(simTrkItr.momentum().pz());
1184 pBeamMC_.push_back(simTrkItr.momentum().P());
1185 }
else if (!
addP_ && (vertIndex == -1)) {
1186 pxyz = simTrkItr.momentum();
1188 if (vertIndex == -1)
1189 vertIndex = simTrkItr.vertIndex();
1197 if (vertIndex != -1 && vertIndex < static_cast<int>(SimVtx->size())) {
1198 edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
1199 for (
int iv = 0; iv < vertIndex; iv++)
1201 edm::LogVerbatim(
"HGCSim") <<
"Vertex " << vertIndex <<
" position " << simVtxItr->position();
1218 std::map<int, double> map_hitLayer;
1219 std::map<int, std::pair<DetId, double>> map_hitCell;
1220 for (
const auto& it : *
hits) {
1221 DetId detId = it.id();
1223 double energy = it.energy();
1231 if (map_hitLayer.count(layer) != 0) {
1232 map_hitLayer[layer] +=
energy;
1234 map_hitLayer[layer] =
energy;
1236 if (map_hitCell.count(cell) != 0) {
1237 double ee =
energy + map_hitCell[cell].second;
1238 map_hitCell[cell] = std::pair<uint32_t, double>(detId, ee);
1240 map_hitCell[cell] = std::pair<uint32_t, double>(detId,
energy);
1243 edm::LogVerbatim(
"HGCSim") <<
"RecHit: " << layer <<
" " << global.
x() <<
" " << global.
y() <<
" " << global.
z()
1248 for (
const auto&
itr : map_hitLayer) {
1249 int layer =
itr.first;
1259 for (
const auto&
itr : map_hitCell) {
1268 for (
const auto&
v : *hgcPH) {
1271 unsigned int id =
v.id();
1273 double time =
v.time();
1274 edm::LogVerbatim(
"HGCSim") <<
"HGCalTBAnalyzer::analyzePassiveHits:Energy:"
1282 }
else if (subdet == 2) {
1286 }
else if (subdet == 3) {
1290 }
else if (subdet == 4) {
1294 }
else if (subdet == 5) {
1303 return i.first <
j.first;