47 #include "TProfile2D.h"
73 static bool sortTime(
const std::pair<double, double>&
i,
const std::pair<double, double>&
j);
144 : ifEE_(iConfig.getParameter<
bool>(
"useEE")),
145 ifFH_(iConfig.getParameter<
bool>(
"useFH")),
146 ifBH_(iConfig.getParameter<
bool>(
"useBH")),
147 ifBeam_(iConfig.getParameter<
bool>(
"useBeam")),
148 doSimHits_(iConfig.getParameter<
bool>(
"doSimHits")),
149 doDigis_(iConfig.getParameter<
bool>(
"doDigis")),
150 doRecHits_(iConfig.getParameter<
bool>(
"doRecHits")),
151 doTree_(iConfig.getParameter<
bool>(
"doTree")),
152 doTreeCell_(iConfig.getParameter<
bool>(
"doTreeCell")),
153 doPassive_(iConfig.getParameter<
bool>(
"doPassive")),
154 doPassiveEE_(iConfig.getParameter<
bool>(
"doPassiveEE")),
155 doPassiveHE_(iConfig.getParameter<
bool>(
"doPassiveHE")),
156 doPassiveBH_(iConfig.getParameter<
bool>(
"doPassiveBH")),
157 addP_(iConfig.getParameter<
bool>(
"addP")),
158 doBeam_(iConfig.getParameter<
bool>(
"doBeam")),
159 detectorEE_(iConfig.getParameter<
std::
string>(
"detectorEE")),
160 detectorFH_(iConfig.getParameter<
std::
string>(
"detectorFH")),
161 detectorBH_(iConfig.getParameter<
std::
string>(
"detectorBH")),
162 detectorBeam_(iConfig.getParameter<
std::
string>(
"detectorBeam")),
163 zFrontEE_(iConfig.getParameter<double>(
"zFrontEE")),
164 zFrontFH_(iConfig.getParameter<double>(
"zFrontFH")),
165 zFrontBH_(iConfig.getParameter<double>(
"zFrontBH")),
166 sampleIndex_(iConfig.getParameter<
int>(
"sampleIndex")),
167 gev2mip200_(iConfig.getUntrackedParameter<double>(
"gev2mip200", 57.0
e-6)),
168 gev2mip300_(iConfig.getUntrackedParameter<double>(
"gev2mip300", 85.5
e-6)),
169 stoc_smear_time_200_(iConfig.getUntrackedParameter<double>(
"stoc_smear_time_200", 10.24)),
170 stoc_smear_time_300_(iConfig.getUntrackedParameter<double>(
"stoc_smear_time_300", 15.5)) {
171 usesResource(
"TFileService");
172 ahcalGeom_ = std::make_unique<AHCalGeometry>(iConfig);
198 tok_hepMC_ = consumes<edm::HepMCProduct>(tmp0);
201 edm::LogVerbatim(
"HGCSim") <<
"HGCalTBAnalyzer:: GeneratorSource = " << tmp0;
281 desc.add<
bool>(
"useEE",
true);
282 desc.add<
double>(
"zFrontEE", 0.0);
287 desc.add<
bool>(
"useFH",
false);
288 desc.add<
double>(
"zFrontFH", 0.0);
293 desc.add<
bool>(
"useBH",
false);
294 desc.add<
double>(
"zFrontBH", 0.0);
299 desc.add<
bool>(
"useBeam",
false);
301 std::vector<int> ids = {
302 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1011, 1012, 1013, 1014, 2001, 2002, 2003, 2004, 2005};
303 desc.add<std::vector<int>>(
"idBeams", ids);
311 desc.add<
bool>(
"doSimHits",
true);
312 desc.add<
bool>(
"doDigis",
true);
313 desc.add<
bool>(
"doRecHits",
true);
314 desc.add<
int>(
"sampleIndex", 0);
315 desc.add<
bool>(
"doTree",
true);
316 desc.add<
bool>(
"doTreeCell",
true);
317 desc.add<
bool>(
"doPassive",
false);
318 desc.add<
bool>(
"doPassiveEE",
false);
319 desc.add<
bool>(
"doPassiveHE",
false);
320 desc.add<
bool>(
"doPassiveBH",
false);
321 desc.add<
bool>(
"addP",
false);
322 desc.add<
bool>(
"doBeam",
false);
323 desc.addUntracked<
double>(
"gev2mip200", 57.0e-6);
324 desc.addUntracked<
double>(
"gev2mip300", 85.5e-6);
325 desc.addUntracked<
double>(
"stoc_smear_time_200", 10.24);
326 desc.addUntracked<
double>(
"stoc_smear_time_300", 15.5);
327 desc.addUntracked<
int>(
"maxDepth", 12);
328 desc.addUntracked<
double>(
"deltaX", 30.0);
329 desc.addUntracked<
double>(
"deltaY", 30.0);
330 desc.addUntracked<
double>(
"deltaZ", 81.0);
331 desc.addUntracked<
double>(
"zFirst", 17.6);
332 descriptions.
add(
"HGCalTBAnalyzer",
desc);
337 hBeam_ =
fs_->
make<TH1D>(
"BeamP",
"Beam Momentum", 1000, 0, 1000.0);
338 for (
int i = 0;
i < 3; ++
i) {
350 sprintf(
name,
"SimHitEn%s", det.c_str());
351 sprintf(
title,
"Sim Hit Energy for %s", det.c_str());
353 sprintf(
name,
"SimHitEnX%s", det.c_str());
354 sprintf(
title,
"Sim Hit Energy for %s", det.c_str());
356 sprintf(
name,
"SimHitTm%s", det.c_str());
357 sprintf(
title,
"Sim Hit Timing for %s", det.c_str());
359 sprintf(
name,
"SimHitLat%s", det.c_str());
360 sprintf(
title,
"Lateral Shower profile (Sim Hit) for %s", det.c_str());
362 sprintf(
name,
"SimHitLng%s", det.c_str());
363 sprintf(
title,
"Longitudinal Shower profile (Sim Hit) for %s", det.c_str());
365 sprintf(
name,
"SimHitLng1%s", det.c_str());
366 sprintf(
title,
"Longitudinal Shower profile (Layer) for %s", det.c_str());
368 sprintf(
name,
"SimHitLng2%s", det.c_str());
369 sprintf(
title,
"Longitudinal Shower profile (Layer) for %s", det.c_str());
374 sprintf(
name,
"DigiADC%s", det.c_str());
375 sprintf(
title,
"ADC at Digi level for %s", det.c_str());
377 sprintf(
name,
"DigiOcc%s", det.c_str());
378 sprintf(
title,
"Occupancy (Digi)for %s", det.c_str());
380 sprintf(
name,
"DigiLng%s", det.c_str());
381 sprintf(
title,
"Longitudinal Shower profile (Digi) for %s", det.c_str());
386 sprintf(
name,
"RecHitEn%s", det.c_str());
387 sprintf(
title,
"Rec Hit Energy for %s", det.c_str());
389 sprintf(
name,
"RecHitOcc%s", det.c_str());
390 sprintf(
title,
"Occupancy (Rec Hit)for %s", det.c_str());
392 sprintf(
name,
"RecHitLat%s", det.c_str());
393 sprintf(
title,
"Lateral Shower profile (Rec Hit) for %s", det.c_str());
395 sprintf(
name,
"RecHitLng%s", det.c_str());
396 sprintf(
title,
"Longitudinal Shower profile (Rec Hit) for %s", det.c_str());
398 sprintf(
name,
"RecHitLng1%s", det.c_str());
399 sprintf(
title,
"Longitudinal Shower profile vs Layer for %s", det.c_str());
500 sprintf(
title,
"Sim Hit Energy in SIM layer %d for %s",
l + 1,
detectorEE_.c_str());
504 sprintf(
title,
"Sim Hit Energy in layer %d for %s", (
l / 3 + 1),
detectorEE_.c_str());
530 sprintf(
title,
"Sim Hit Energy in layer %d for %s", (
l / 3 + 1),
detectorFH_.c_str());
555 for (
unsigned int l = 0;
l <
idBeams_.size(); ++
l) {
572 HepMC::FourVector pxyz(0, 0, 0, 0);
573 for (HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end();
575 edm::LogVerbatim(
"HGCSim") <<
"Particle [" <<
k <<
"] with p " << (*p)->momentum().rho() <<
" theta "
576 << (*p)->momentum().theta() <<
" phi " << (*p)->momentum().phi() <<
" pxyz ("
577 << (*p)->momentum().px() <<
", " << (*p)->momentum().py() <<
", "
578 << (*p)->momentum().pz() <<
")";
580 pxyz.setPx(pxyz.px() + (*p)->momentum().px());
581 pxyz.setPy(pxyz.py() + (*p)->momentum().py());
582 pxyz.setPz(pxyz.pz() + (*p)->momentum().pz());
583 pxyz.setE(pxyz.e() + (*p)->momentum().e());
584 }
else if (!
addP_ && (
k == 0)) {
585 pxyz = (*p)->momentum();
589 edm::LogVerbatim(
"HGCSim") <<
"Particle with p " << pxyz.rho() <<
" theta " << pxyz.theta() <<
" phi "
626 std::vector<PCaloHit> caloHits;
631 if (theCaloHitContainers.
isValid()) {
637 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
649 if (theCaloHitContainers.
isValid()) {
655 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
667 if (theCaloHitContainers.
isValid()) {
673 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
684 if (theCaloHitContainers.
isValid()) {
687 << theCaloHitContainers->size() <<
" hits";
690 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
753 if (theDigiContainers.
isValid()) {
758 for (
const auto& it : *theDigiContainers) {
761 uint16_t
adc = hgcSample.
data();
769 if (theDigiContainers.
isValid()) {
774 for (
const auto& it : *theDigiContainers) {
777 uint16_t
adc = hgcSample.
data();
789 if (theCaloHitContainers.
isValid()) {
792 << theCaloHitContainers->
size() <<
" hits";
803 if (theCaloHitContainers.
isValid()) {
806 << theCaloHitContainers->
size() <<
" hits";
820 std::map<uint32_t, double> map_hits, map_hitn;
821 std::map<uint32_t, double> map_hittime_firsthit, map_hittime_lasthit, map_hittime_15Mip;
822 std::map<int, double> map_hitDepth, map_hitWafer;
823 std::map<int, std::pair<uint32_t, double>> map_hitLayer, map_hitCell;
825 std::map<uint32_t, double>
nhits;
826 std::map<uint32_t, int>
ID,
Depth;
827 std::map<uint32_t, double> GeV2Mip;
828 std::map<uint32_t, double> StochTermTime;
829 std::vector<int> nSimLayers;
831 std::map<uint32_t, std::vector<std::pair<double, double>>> map_hitTimeEn;
834 for (
unsigned int i = 0;
i <
hits.size();
i++) {
837 uint32_t
id =
hits[
i].id();
852 }
else if (
type == 3) {
861 idx = sector * 1000 + cell;
865 <<
layer <<
":" << sector <<
":" << cell <<
" Position " <<
xy.first <<
":"
871 << sector <<
":" << subsector <<
":" << cell <<
":" <<
depth <<
" Energy " <<
energy
874 if (map_hits.count(
id) != 0) {
879 if (map_hitLayer.count(
layer) != 0) {
881 map_hitLayer[
layer] = std::make_pair(
id, ee);
885 if (map_hitWafer.count(sector) != 0)
886 map_hitWafer[sector] +=
energy;
888 map_hitWafer[sector] =
energy;
890 if (map_hitCell.count(
idx) != 0) {
891 double ee =
energy + map_hitCell[
idx].second;
892 map_hitCell[
idx] = std::make_pair(
id, ee);
894 map_hitCell[
idx] = std::make_pair(
id,
energy);
902 if (map_hitDepth.count(
depth) != 0) {
909 map_hitTimeEn[idn].push_back(std::make_pair(
time,
energy));
914 if (map_hitn.count(idn) != 0) {
926 edm::LogVerbatim(
"HGCSim") <<
"HGCalTAnalyzer:: " << map_hitWafer.size() <<
" wafers are hit in type " <<
type;
927 for (
auto itr = map_hitWafer.begin(); itr != map_hitWafer.end(); ++itr)
928 edm::LogVerbatim(
"HGCSim") <<
"Wafer: " << itr->first <<
" Deposited Energy " << itr->second;
930 for (
const auto& itr : map_hitTimeEn) {
931 uint32_t
id = itr.first;
940 <<
"\nID(sim) and id(reco) " << std::hex <<
ID[
id] <<
" " <<
id <<
std::dec;
953 double totEbeforeThreshold = 0.;
954 double timebeforeThreshold = 0.;
955 double timeAtThresohld = 0.;
956 for (
unsigned int ihit = 0; ihit < map_hitTimeEn[
id].size(); ihit++) {
957 double energy = (map_hitTimeEn[
id].at(ihit)).
second / GeV2Mip[
id];
961 edm::LogVerbatim(
"HGCSim") <<
"Tot E till now : time of that E : GeV2Mip[id] is " << totE <<
" " <<
time
962 <<
" " << GeV2Mip[
id];
964 totEbeforeThreshold = totE;
965 timebeforeThreshold =
time;
968 (
threshold - totEbeforeThreshold) * (
time - timebeforeThreshold) / (totE - totEbeforeThreshold) +
970 map_hittime_15Mip[
id] = timeAtThresohld;
972 edm::LogVerbatim(
"HGCSim") <<
"ihit : energyBefore : timeBefore : energyTot : timeTot : timeAt15MIP "
973 << ihit <<
" " << totEbeforeThreshold <<
" " << timebeforeThreshold <<
" "
974 << totE <<
" " <<
time <<
" " << map_hittime_15Mip[
id];
978 if (!map_hitTimeEn[
id].
empty()) {
979 map_hittime_firsthit[
id] = (map_hitTimeEn[
id].at(0)).
first;
980 map_hittime_lasthit[
id] = (map_hitTimeEn[
id].at(map_hitTimeEn[
id].
size() - 1)).
first;
981 if (map_hittime_15Mip[
id] < map_hittime_firsthit[
id])
982 map_hittime_15Mip[
id] = map_hittime_firsthit[
id];
992 map_hittime_15Mip[
id] = -99;
994 edm::LogVerbatim(
"HGCSim") <<
"id : first hit time : last hit time " <<
id <<
" " << map_hittime_firsthit[
id]
995 <<
" " << map_hittime_lasthit[
id] <<
"\nFinally for this cell, time is "
996 << map_hittime_15Mip[
id];
1001 for (
const auto& itr : map_hits) {
1007 for (
const auto& itr : map_hitLayer) {
1008 int layer = (
type == 2) ? itr.first : (itr.first - 1);
1027 }
else if (
type == 1) {
1032 }
else if (
type == 2) {
1040 for (
unsigned int k = 0;
k <
idBeams_.size(); ++
k) {
1049 for (
const auto& itr : map_hitDepth) {
1050 int layer = (
type == 2) ? itr.first : (itr.first - 1);
1051 double energy = itr.second;
1061 }
else if (
type == 1) {
1066 }
else if (
type == 2) {
1079 for (
const auto& itr : map_hitCell) {
1080 uint32_t
id = ((itr.second).
first);
1082 std::pair<float, float>
xy(0, 0);
1088 int subdet,
zside,
layer, sector, subsector, cell;
1092 xx = (zp < 0) ? -
xy.first :
xy.first;
1098 for (
const auto& itr : map_hitn) {
1099 uint32_t
id = itr.first;
1100 double energy = itr.second;
1102 double time_firsthit = map_hittime_firsthit[
id];
1103 double time15Mip = map_hittime_15Mip[
id];
1104 double time_lasthit = map_hittime_lasthit[
id];
1110 if (
debug && (
energy / GeV2Mip[
id] < 15) && (map_hittime_15Mip[
id] > 0))
1112 <<
" " << map_hittime_15Mip[
id];
1113 }
else if (
type == 1) {
1114 double time_firsthit = map_hittime_firsthit[
id];
1115 double time15Mip = map_hittime_15Mip[
id];
1116 double time_lasthit = map_hittime_lasthit[
id];
1122 }
else if (
type == 2) {
1126 int row = hid.
irow();
1136 }
else if (
type == 3) {
1159 std::vector<float> verX, verY, verZ;
1163 for (
const auto& simVtxItr : *SimVtx) {
1164 verX.push_back(simVtxItr.position().X());
1165 verY.push_back(simVtxItr.position().Y());
1166 verZ.push_back(simVtxItr.position().Z());
1171 HepMC::FourVector pxyz(0, 0, 0, 0);
1172 for (
const auto& simTrkItr : *SimTk) {
1173 if (
addP_ && !(simTrkItr.noGenpart())) {
1174 pxyz.setPx(pxyz.px() + simTrkItr.momentum().px());
1175 pxyz.setPy(pxyz.py() + simTrkItr.momentum().py());
1176 pxyz.setPz(pxyz.pz() + simTrkItr.momentum().pz());
1177 pxyz.setE(pxyz.e() + simTrkItr.momentum().e());
1179 edm::LogVerbatim(
"HGCSim") <<
"Track " << simTrkItr.trackId() <<
" Vertex " << simTrkItr.vertIndex() <<
" Type "
1180 << simTrkItr.type() <<
" Charge " << simTrkItr.charge() <<
" px "
1181 << simTrkItr.momentum().px() <<
" py " << simTrkItr.momentum().py() <<
" pz "
1182 << simTrkItr.momentum().pz() <<
" P " << simTrkItr.momentum().P() <<
" GenIndex "
1183 << simTrkItr.genpartIndex();
1185 <<
" position-> X: " << verX[simTrkItr.vertIndex()]
1186 <<
" Y: " << verY[simTrkItr.vertIndex()] <<
" Z: " << verZ[simTrkItr.vertIndex()];
1189 if (
doBeam_ && !(simTrkItr.noGenpart())) {
1192 xBeamMC_.push_back(verX[simTrkItr.vertIndex()]);
1193 yBeamMC_.push_back(verY[simTrkItr.vertIndex()]);
1194 zBeamMC_.push_back(verZ[simTrkItr.vertIndex()]);
1195 pxBeamMC_.push_back(simTrkItr.momentum().px());
1196 pyBeamMC_.push_back(simTrkItr.momentum().py());
1197 pzBeamMC_.push_back(simTrkItr.momentum().pz());
1198 pBeamMC_.push_back(simTrkItr.momentum().P());
1199 }
else if (!
addP_ && (vertIndex == -1)) {
1200 pxyz = simTrkItr.momentum();
1202 if (vertIndex == -1)
1203 vertIndex = simTrkItr.vertIndex();
1211 if (vertIndex != -1 && vertIndex < static_cast<int>(SimVtx->size())) {
1212 edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
1213 for (
int iv = 0;
iv < vertIndex;
iv++)
1215 edm::LogVerbatim(
"HGCSim") <<
"Vertex " << vertIndex <<
" position " << simVtxItr->position();
1232 std::map<int, double> map_hitLayer;
1233 std::map<int, std::pair<DetId, double>> map_hitCell;
1234 for (
const auto& it : *
hits) {
1235 DetId detId = it.id();
1237 double energy = it.energy();
1245 if (map_hitLayer.count(
layer) != 0) {
1250 if (map_hitCell.count(cell) != 0) {
1251 double ee =
energy + map_hitCell[cell].second;
1252 map_hitCell[cell] = std::pair<uint32_t, double>(detId, ee);
1254 map_hitCell[cell] = std::pair<uint32_t, double>(detId,
energy);
1262 for (
const auto& itr : map_hitLayer) {
1263 int layer = itr.first;
1264 double energy = itr.second;
1273 for (
const auto& itr : map_hitCell) {
1282 for (
const auto&
v : *hgcPH) {
1285 unsigned int id =
v.id();
1287 double time =
v.time();
1288 edm::LogVerbatim(
"HGCSim") <<
"HGCalTBAnalyzer::analyzePassiveHits:Energy:"
1296 }
else if (subdet == 2) {
1300 }
else if (subdet == 3) {
1304 }
else if (subdet == 4) {
1308 }
else if (subdet == 5) {
1317 return i.first <
j.first;