28 #include "G4HCofThisEvent.hh"
29 #include "G4SDManager.hh"
32 #include "G4VProcess.hh"
37 : jetf(nullptr), numberingFromDDD(nullptr), org(nullptr) {
57 produces<PHcalValidInfoNxN>(
labelNxN);
61 edm::LogVerbatim(
"ValidHcal") <<
"HcalTestAnalysis:: Initialised as observer of begin/end events and "
62 <<
"of G4step with Parameter values: \n\tInfoLevel = " <<
infolevel
74 edm::LogVerbatim(
"ValidHcal") <<
"\n --------> Total number of selected entries"
75 <<
" : " <<
count <<
"\nPointers:: JettFinder " <<
jetf <<
", Numbering Scheme " <<
org
109 float sHB[4] = {117., 117., 178., 217.};
110 float sHE[3] = {178., 178., 178.};
111 float sHF[3] = {2.84, 2.09, 0.};
113 float deta[4] = {0.0435, 0.1305, 0.2175, 0.3045};
114 float dphi[4] = {0.0436, 0.1309, 0.2182, 0.3054};
117 for (
i = 0;
i < 4;
i++) {
120 for (
i = 0;
i < 3;
i++) {
123 for (
i = 0;
i < 3;
i++) {
128 for (
i = 0;
i < 4;
i++) {
129 dEta.push_back(deta[
i]);
130 dPhi.push_back(dphi[
i]);
143 (*job)()->get<HcalSimNumberingRecord>().get(hdc);
146 <<
"HcalNumberingFromDDD";
154 int irun = (*run)()->GetRunID();
159 G4SDManager *
sd = G4SDManager::GetSDMpointerIfExist();
161 G4VSensitiveDetector *aSD =
sd->FindSensitiveDetector(sdname);
162 if (aSD ==
nullptr) {
163 edm::LogWarning(
"ValidHcal") <<
"SimG4HcalValidation::beginOfRun: No SD"
164 <<
" with name " << sdname <<
" in this "
167 HCalSD *theCaloSD = dynamic_cast<HCalSD *>(aSD);
168 edm::LogVerbatim(
"ValidHcal") <<
"SimG4HcalValidation::beginOfRun: Finds SD with name " << theCaloSD->GetName()
172 edm::LogVerbatim(
"ValidHcal") <<
"SimG4HcalValidation::beginOfRun: set a new numbering scheme";
176 edm::LogWarning(
"ValidHcal") <<
"SimG4HcalValidation::beginOfRun: Could "
177 <<
"not get SD Manager!";
185 for (
i = 0;
i < 5;
i++)
187 for (
i = 0;
i < 20;
i++)
193 int iev = (*evt)()->GetEventID();
194 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation: =====> Begin of event = " << iev;
199 if (aStep !=
nullptr) {
200 G4VPhysicalVolume *curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
201 G4String
name = curPV->GetName();
203 double edeposit = aStep->GetTotalEnergyDeposit();
204 int layer = -1,
depth = -1;
208 }
else if (
name ==
"EFR") {
211 }
else if (
name ==
"HBS") {
212 layer = (curPV->GetCopyNo() / 10) % 100;
213 depth = (curPV->GetCopyNo()) % 10 + 1;
214 if (
depth > 0 && depth < 4 && layer >= 0 && layer < 17) {
217 edm::LogWarning(
"ValidHcal") <<
"SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
221 }
else if (
name ==
"HES") {
222 layer = (curPV->GetCopyNo() / 10) % 100;
223 depth = (curPV->GetCopyNo()) % 10 + 1;
224 if (
depth > 0 && depth < 3 && layer >= 0 && layer < 19) {
227 edm::LogWarning(
"ValidHcal") <<
"SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
231 }
else if (
name ==
"HTS") {
232 layer = (curPV->GetCopyNo() / 10) % 100;
233 depth = (curPV->GetCopyNo()) % 10 + 1;
234 if (
depth > 3 && depth < 5 && layer >= 17 && layer < 20) {
237 edm::LogWarning(
"ValidHcal") <<
"SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
244 if (layer >= 0 && layer < 20)
245 edepl[layer] += edeposit;
247 if (layer >= 0 && layer < 20) {
248 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:: G4Step: " <<
name <<
" Layer " << std::setw(3) << layer
249 <<
" Depth " << std::setw(2) <<
depth <<
" Edep " << std::setw(6) << edeposit /
MeV
261 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:: --- after Fill ";
266 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:Fill event " << (*evt)()->GetEventID();
269 G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
274 int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(
names[0]);
276 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation :: Hit Collection for " <<
names[0] <<
" of ID " << HCHCid
277 <<
" is obtained at " << theHCHC;
278 if (HCHCid >= 0 && theHCHC !=
nullptr) {
279 int nhits = theHCHC->entries();
292 int subdet,
zside, layer, etaIndex, phiIndex, lay;
310 edm::LogVerbatim(
"ValidHcal") <<
"SimG4HcalValidation::HitPMT " << subdet <<
" " << (2 *
zside - 1) * etaIndex
311 <<
" " << phiIndex <<
" " << layer + 1 <<
" R " <<
pos.rho() <<
" Phi "
313 }
else if (subdet == static_cast<int>(
HcalEndcap)) {
314 if (etaIndex <= 20) {
320 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::debugFill Hcal " << det <<
" layer " << std::setw(2) << layer
321 <<
" lay " << std::setw(2) << lay <<
" time " << std::setw(6) <<
time <<
" theta "
322 << std::setw(8) <<
theta <<
" eta " << std::setw(8) <<
eta <<
" phi " << std::setw(8) <<
phi
323 <<
" e " << std::setw(8) <<
e <<
" ID 0x" << std::hex << unitID <<
" ID dec " <<
std::dec
339 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:: HCAL hits : " << nhc;
342 int ndets =
names.size();
345 for (
int idty = 1; idty < ndets; idty++) {
353 int ECHCid = G4SDManager::GetSDMpointer()->GetCollectionID(
names[idty]);
355 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:: Hit Collection for " <<
names[idty] <<
" of ID " << ECHCid
356 <<
" is obtained at " << theECHC;
357 int theechc_entries = theECHC->entries();
358 if (ECHCid >= 0 && theECHC !=
nullptr) {
359 for (
j = 0;
j < theechc_entries;
j++) {
385 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::debugFill Ecal " << det <<
" layer " << std::setw(2) << layer
386 <<
" lay " << std::setw(2) << lay <<
" time " << std::setw(6) <<
time <<
" theta "
387 << std::setw(8) <<
theta <<
" eta " << std::setw(8) <<
eta <<
" phi " << std::setw(8)
388 <<
phi <<
" e " << std::setw(8) <<
e <<
" ID 0x" << std::hex << unitID <<
" ID dec "
393 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:: " << det <<
" hits : " << nec;
400 LogDebug(
"ValidHcal") <<
"\n ===>>> SimG4HcalValidation: Energy deposit "
401 <<
"in MeV\n at EB : " << std::setw(6) <<
edepEB /
MeV <<
"\n at EE : " << std::setw(6)
403 <<
"\n at HE : " << std::setw(6) <<
edepHE /
MeV <<
"\n at HO : " << std::setw(6)
404 <<
edepHO /
MeV <<
"\n ---- SimG4HcalValidation: Energy deposit in";
405 for (
i = 0;
i < 5;
i++)
406 LogDebug(
"ValidHcal") <<
" Depth " << std::setw(2) <<
i <<
" E " << std::setw(8) <<
edepd[
i] /
MeV <<
" MeV";
407 LogDebug(
"ValidHcal") <<
" ---- SimG4HcalValidation: Energy deposit in"
409 for (
i = 0;
i < 20;
i++)
410 LogDebug(
"ValidHcal") <<
" Layer " << std::setw(2) <<
i <<
" E " << std::setw(8) <<
edepl[
i] /
MeV <<
" MeV";
416 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::HF hits " <<
vhitec <<
" in EC and " <<
vhithc <<
" in HC\n"
417 <<
" HB/HE " <<
enEcal <<
" in EC and " <<
enHcal <<
" in HC";
422 LogDebug(
"ValidHcal") <<
" layerAnalysis ===> after fetchHits";
428 std::vector<CaloHit>::iterator hit_itr;
430 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::NxNAnalysis : entrance ";
435 LogDebug(
"ValidHcal") <<
" NxNAnalysis : coolectEnergyRdir - Ecal " <<
een <<
" Hcal " <<
hen;
444 for (hit_itr =
hits->begin(); hit_itr <
hits->end(); hit_itr++) {
445 double e = hit_itr->e();
446 double t = hit_itr->t();
447 double eta = hit_itr->eta();
448 double phi = hit_itr->phi();
449 int type = hit_itr->det();
450 int layer = hit_itr->layer();
465 for (
int i = 0;
i <
max;
i++) {
468 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation:: nxNAnalysis eta0,"
469 <<
" phi0 = " <<
eta0 <<
" " <<
phi0 <<
" type, layer = " <<
type <<
"," << layer
470 <<
" eta, phi = " <<
eta <<
" " <<
phi;
484 LogDebug(
"ValidHcal") <<
" nxNAnalysis ===> after fillHvsE";
489 std::vector<CaloHit> *hhit = &
hitcache;
496 std::vector<SimG4HcalHitCluster>::iterator clus_itr;
498 LogDebug(
"ValidHcal") <<
"\n ---------- Final list of " << (*result).size() <<
" clusters ---------------";
499 for (clus_itr =
result->begin(); clus_itr <
result->end(); clus_itr++)
500 LogDebug(
"ValidHcal") << (*clus_itr);
502 std::vector<double> enevec, etavec, phivec;
504 if (!(*result).empty()) {
505 sort((*result).begin(), (*result).end());
507 clus_itr =
result->begin();
508 double etac = clus_itr->eta();
509 double phic = clus_itr->phi();
511 double ecal_collect = 0.;
513 ecal_collect = clus_itr->collectEcalEnergyR();
518 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> ecal_collect = " << ecal_collect;
522 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> eta phi deviation = " << dist;
525 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> after fillEtaPhiProfileJet";
527 std::vector<CaloHit> *
hits = clus_itr->getHits();
528 std::vector<CaloHit>::iterator hit_itr;
530 double ee = 0.,
he = 0.,
hoe = 0., etot = 0.;
533 for (hit_itr =
hits->begin(); hit_itr <
hits->end(); hit_itr++) {
534 double e = hit_itr->e();
535 double t = hit_itr->t();
536 double r =
jetf->
rDist(&(*clus_itr), &(*hit_itr));
540 if (hit_itr->det() == 10 || hit_itr->det() == 11 || hit_itr->det() == 12)
542 if (hit_itr->det() == static_cast<int>(
HcalBarrel) || hit_itr->det() == static_cast<int>(
HcalEndcap) ||
545 if (hit_itr->det() == static_cast<int>(
HcalBarrel) && hit_itr->layer() > 17)
549 if (hit_itr->det() == static_cast<int>(
HcalBarrel) || hit_itr->det() == static_cast<int>(
HcalEndcap) ||
557 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> after fillEcollectJet: "
558 <<
"ee/he/hoe/etot " << ee <<
"/" <<
he <<
"/" <<
hoe <<
"/" << etot;
561 for (clus_itr =
result->begin(); clus_itr <
result->end(); clus_itr++) {
563 enevec.push_back(clus_itr->e());
564 etavec.push_back(clus_itr->eta());
565 phivec.push_back(clus_itr->phi());
568 product.
fillJets(enevec, etavec, phivec);
570 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> after fillJets\n"
571 <<
" JetAnalysis ===> (*result).size() " << (*result).size();
574 if (etavec.size() > 1) {
575 if (etavec[0] > -2.5 && etavec[0] < 2.5 && etavec[1] > -2.5 && etavec[1] < 2.5) {
576 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> Di-jet mass enter\n"
577 <<
" JetAnalysis ===> Di-jet vectors ";
578 for (
unsigned int i = 0;
i < enevec.size();
i++)
579 LogDebug(
"ValidHcal") <<
" e, eta, phi = " << enevec[
i] <<
" " << etavec[
i] <<
" " << phivec[
i];
581 double et0 = enevec[0] / cosh(etavec[0]);
582 double px0 = et0 *
cos(phivec[0]);
583 double py0 = et0 *
sin(phivec[0]);
584 double pz0 = et0 * sinh(etavec[0]);
585 double et1 = enevec[1] / cosh(etavec[1]);
586 double px1 = et1 *
cos(phivec[1]);
587 double py1 = et1 *
sin(phivec[1]);
588 double pz1 = et1 * sinh(etavec[1]);
590 double dijetmass2 = (enevec[0] + enevec[1]) * (enevec[0] + enevec[1]) - (px1 + px0) * (px1 + px0) -
591 (py1 + py0) * (py1 + py0) - (pz1 + pz0) * (pz1 + pz0);
593 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> Di-jet massSQ " << dijetmass2;
596 if (dijetmass2 >= 0.)
597 dijetmass =
sqrt(dijetmass2);
599 dijetmass = -
sqrt(-1. * dijetmass2);
603 LogDebug(
"ValidHcal") <<
" JetAnalysis ===> after fillDiJets";
611 LogDebug(
"ValidHcal") <<
"Enter SimG4HcalValidation::fetchHits with " <<
hitcache.size() <<
" hits";
615 std::vector<CaloHit>::iterator
itr;
616 std::vector<CaloHit *> lhits(nHit);
618 uint32_t unitID =
itr->id();
623 group = (subdet & 15) << 20;
624 group += ((lay - 1) & 31) << 15;
630 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::fetchHits:Original " <<
i <<
" " <<
hitcache[
i] <<
"\n"
631 <<
"SimG4HcalValidation::fetchHits:Copied " <<
i <<
" " << *lhits[
i];
634 std::vector<CaloHit *>::iterator k1, k2;
635 for (
i = 0, k1 = lhits.begin(); k1 != lhits.end();
i++, k1++)
636 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::fetchHits:Sorted " <<
i <<
" " << **k1;
638 for (
i = 0, k1 = lhits.begin(); k1 != lhits.end();
i++, k1++) {
639 double ehit = (**k1).e();
640 double t = (**k1).t();
641 uint32_t unitID = (**k1).id();
643 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::fetchHits:Start " <<
i <<
" U/T/E"
644 <<
" 0x" << std::hex << unitID <<
std::dec <<
" " <<
t <<
" " << ehit;
645 for (k2 = k1 + 1; k2 != lhits.end() && (
t - (**k2).t()) < 1 && (
t - (**k2).t()) > -1 && unitID == (**k2).id();
648 LogDebug(
"ValidHcal") <<
"\t + " << (**k2).e();
651 LogDebug(
"ValidHcal") <<
"\t = " << ehit <<
" in " << jump;
653 double eta = (*k1)->eta();
654 double phi = (*k1)->phi();
655 int lay = ((unitID >> 15) & 31) + 1;
656 int subdet = (unitID >> 20) & 15;
657 int zside = (unitID >> 14) & 1;
658 int ieta = (unitID >> 7) & 127;
659 int iphi = (unitID)&127;
665 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::fetchHits:Hit " << nHits <<
" " <<
i <<
" ID 0x" << std::hex
666 << unitID <<
" det " <<
std::dec << subdet <<
" " << lay <<
" " <<
zside <<
" " <<
ieta
667 <<
" " <<
iphi <<
" Time " <<
t <<
" E " << ehit;
673 LogDebug(
"ValidHcal") <<
"SimG4HcalValidation::fetchHits called with " << nHit <<
" hits"
674 <<
" and writes out " << nHits <<
'(' <<
hit <<
") hits";
682 std::vector<CaloHit>::iterator hit_itr;
684 double sume = 0., sumh = 0., sumho = 0.;
686 for (hit_itr =
hits->begin(); hit_itr <
hits->end(); hit_itr++) {
687 double e = hit_itr->e();
688 double eta = hit_itr->eta();
689 double phi = hit_itr->phi();
691 int type = hit_itr->det();
699 if (
type == static_cast<int>(
HcalBarrel) && hit_itr->layer() > 17)
716 }
else if (det ==
"HES" || det ==
"HED") {
718 }
else if (det ==
"HF") {