File indexing completed on 2024-10-08 05:12:11
0001
0002
0003
0004
0005 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0006 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0007 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0008 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0009 #include "SimG4Core/Notification/interface/Observer.h"
0010 #include "SimG4Core/Watcher/interface/SimProducer.h"
0011
0012
0013 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0014 #include "DataFormats/Math/interface/GeantUnits.h"
0015 #include "DataFormats/Math/interface/Point3D.h"
0016 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
0017 #include "SimDataFormats/ValidationFormats/interface/PValidationFormats.h"
0018
0019 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0020 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0021 #include "SimG4CMS/Calo/interface/HCalSD.h"
0022 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0023
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028
0029 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0030 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
0031 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
0032
0033 #include "Validation/HcalHits/interface/SimG4HcalHitCluster.h"
0034 #include "Validation/HcalHits/interface/SimG4HcalHitJetFinder.h"
0035
0036 #include <CLHEP/Units/SystemOfUnits.h>
0037
0038 #include "G4HCofThisEvent.hh"
0039 #include "G4SDManager.hh"
0040 #include "G4Step.hh"
0041 #include "G4Track.hh"
0042 #include "G4VProcess.hh"
0043
0044 #include <iostream>
0045 #include <memory>
0046 #include <string>
0047 #include <vector>
0048
0049 using namespace geant_units::operators;
0050 using CLHEP::GeV;
0051 using CLHEP::MeV;
0052
0053 class SimG4HcalValidation : public SimProducer,
0054 public Observer<const BeginOfRun *>,
0055 public Observer<const BeginOfEvent *>,
0056 public Observer<const EndOfEvent *>,
0057 public Observer<const G4Step *> {
0058 public:
0059 SimG4HcalValidation(const edm::ParameterSet &p);
0060 SimG4HcalValidation(const SimG4HcalValidation &) = delete;
0061 const SimG4HcalValidation &operator=(const SimG4HcalValidation &) = delete;
0062 ~SimG4HcalValidation() override;
0063
0064 void registerConsumes(edm::ConsumesCollector) override;
0065 void produce(edm::Event &, const edm::EventSetup &) override;
0066 void beginRun(edm::EventSetup const &) override;
0067
0068 private:
0069 void init();
0070
0071
0072 void update(const BeginOfRun *run) override;
0073 void update(const BeginOfEvent *evt) override;
0074 void update(const G4Step *step) override;
0075 void update(const EndOfEvent *evt) override;
0076
0077
0078 void fill(const EndOfEvent *ev);
0079 void layerAnalysis(PHcalValidInfoLayer &);
0080 void nxNAnalysis(PHcalValidInfoNxN &);
0081 void jetAnalysis(PHcalValidInfoJets &);
0082 void fetchHits(PHcalValidInfoLayer &);
0083 void clear();
0084 void collectEnergyRdir(const double, const double);
0085 double getHcalScale(std::string, int) const;
0086
0087 private:
0088 edm::ESGetToken<HcalDDDSimConstants, HcalSimNumberingRecord> ddconsToken_;
0089
0090
0091 std::unique_ptr<SimG4HcalHitJetFinder> jetf;
0092
0093
0094 std::unique_ptr<HcalNumberingFromDDD> numberingFromDDD;
0095
0096
0097 std::unique_ptr<HcalTestNumberingScheme> org;
0098
0099
0100 std::vector<CaloHit> hitcache;
0101
0102
0103 std::vector<float> scaleHB;
0104 std::vector<float> scaleHE;
0105 std::vector<float> scaleHF;
0106
0107
0108 std::vector<std::string> names;
0109 double coneSize, ehitThreshold, hhitThreshold;
0110 float timeLowlim, timeUplim, eta0, phi0, jetThreshold;
0111 bool applySampling, hcalOnly;
0112 int infolevel;
0113 std::string labelLayer, labelNxN, labelJets;
0114
0115
0116 std::vector<double> dEta;
0117 std::vector<double> dPhi;
0118
0119
0120 unsigned int count;
0121 double edepEB, edepEE, edepHB, edepHE, edepHO;
0122 double edepd[5], edepl[20];
0123 double een, hen, hoen;
0124 double vhitec, vhithc, enEcal, enHcal;
0125 };
0126
0127 SimG4HcalValidation::SimG4HcalValidation(const edm::ParameterSet &p) {
0128 edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("SimG4HcalValidation");
0129 infolevel = m_Anal.getParameter<int>("InfoLevel");
0130 hcalOnly = m_Anal.getParameter<bool>("HcalClusterOnly");
0131 applySampling = m_Anal.getParameter<bool>("HcalSampling");
0132 coneSize = m_Anal.getParameter<double>("ConeSize");
0133 ehitThreshold = m_Anal.getParameter<double>("EcalHitThreshold");
0134 hhitThreshold = m_Anal.getParameter<double>("HcalHitThreshold");
0135 timeLowlim = m_Anal.getParameter<double>("TimeLowLimit");
0136 timeUplim = m_Anal.getParameter<double>("TimeUpLimit");
0137 jetThreshold = m_Anal.getParameter<double>("JetThreshold");
0138 eta0 = m_Anal.getParameter<double>("Eta0");
0139 phi0 = m_Anal.getParameter<double>("Phi0");
0140 names = m_Anal.getParameter<std::vector<std::string>>("Names");
0141 labelLayer = m_Anal.getUntrackedParameter<std::string>("LabelLayerInfo", "HcalInfoLayer");
0142 labelNxN = m_Anal.getUntrackedParameter<std::string>("LabelNxNInfo", "HcalInfoNxN");
0143 labelJets = m_Anal.getUntrackedParameter<std::string>("LabelJetsInfo", "HcalInfoJets");
0144
0145 produces<PHcalValidInfoLayer>(labelLayer);
0146 if (infolevel > 0)
0147 produces<PHcalValidInfoNxN>(labelNxN);
0148 if (infolevel > 1)
0149 produces<PHcalValidInfoJets>(labelJets);
0150
0151 edm::LogVerbatim("ValidHcal") << "HcalTestAnalysis:: Initialized as observer of begin/end events and "
0152 << "of G4step with Parameter values: \n\tInfoLevel = " << infolevel
0153 << "\n\thcalOnly = " << hcalOnly << "\n\tapplySampling = " << applySampling
0154 << "\n\tconeSize = " << coneSize << "\n\tehitThreshold = " << ehitThreshold
0155 << "\n\thhitThreshold = " << hhitThreshold << "\n\tttimeLowlim = " << timeLowlim
0156 << "\n\tttimeUplim = " << timeUplim << "\n\teta0 = " << eta0
0157 << "\n\tphi0 = " << phi0 << "\nLabels (Layer): " << labelLayer
0158 << " (NxN): " << labelNxN << " (Jets): " << labelJets;
0159
0160 init();
0161 }
0162
0163 SimG4HcalValidation::~SimG4HcalValidation() {
0164 edm::LogVerbatim("ValidHcal") << "\n --------> Total number of selected entries"
0165 << " : " << count;
0166 }
0167
0168 void SimG4HcalValidation::registerConsumes(edm::ConsumesCollector cc) {
0169 ddconsToken_ = cc.esConsumes<HcalDDDSimConstants, HcalSimNumberingRecord, edm::Transition::BeginRun>();
0170 edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::Initialize ESGetToken for HcalDDDSimConstants";
0171 }
0172
0173 void SimG4HcalValidation::produce(edm::Event &e, const edm::EventSetup &) {
0174 std::unique_ptr<PHcalValidInfoLayer> productLayer(new PHcalValidInfoLayer);
0175 layerAnalysis(*productLayer);
0176 e.put(std::move(productLayer), labelLayer);
0177
0178 if (infolevel > 0) {
0179 std::unique_ptr<PHcalValidInfoNxN> productNxN(new PHcalValidInfoNxN);
0180 nxNAnalysis(*productNxN);
0181 e.put(std::move(productNxN), labelNxN);
0182 }
0183
0184 if (infolevel > 1) {
0185 std::unique_ptr<PHcalValidInfoJets> productJets(new PHcalValidInfoJets);
0186 jetAnalysis(*productJets);
0187 e.put(std::move(productJets), labelJets);
0188 }
0189 }
0190
0191
0192 void SimG4HcalValidation::init() {
0193 float sHB[4] = {117., 117., 178., 217.};
0194 float sHE[3] = {178., 178., 178.};
0195 float sHF[3] = {2.84, 2.09, 0.};
0196
0197 float deta[4] = {0.0435, 0.1305, 0.2175, 0.3045};
0198 float dphi[4] = {0.0436, 0.1309, 0.2182, 0.3054};
0199
0200 int i = 0;
0201 for (i = 0; i < 4; i++) {
0202 scaleHB.push_back(sHB[i]);
0203 }
0204 for (i = 0; i < 3; i++) {
0205 scaleHE.push_back(sHE[i]);
0206 }
0207 for (i = 0; i < 3; i++) {
0208 scaleHF.push_back(sHF[i]);
0209 }
0210
0211
0212 for (i = 0; i < 4; i++) {
0213 dEta.push_back(deta[i]);
0214 dPhi.push_back(dphi[i]);
0215 }
0216
0217
0218 jetf = std::make_unique<SimG4HcalHitJetFinder>(coneSize);
0219
0220
0221 count = 0;
0222 }
0223
0224 void SimG4HcalValidation::beginRun(edm::EventSetup const &es) {
0225
0226 const HcalDDDSimConstants *hcons = &es.getData(ddconsToken_);
0227 edm::LogVerbatim("ValidHcal") << "HcalTestAnalysis:: Initialise HcalNumberingFromDDD";
0228 numberingFromDDD = std::make_unique<HcalNumberingFromDDD>(hcons);
0229
0230
0231 org = std::make_unique<HcalTestNumberingScheme>(false);
0232 }
0233
0234 void SimG4HcalValidation::update(const BeginOfRun *run) {
0235 int irun = (*run)()->GetRunID();
0236
0237 edm::LogVerbatim("ValidHcal") << " =====> Begin of Run = " << irun;
0238
0239 std::string sdname = names[0];
0240 G4SDManager *sd = G4SDManager::GetSDMpointerIfExist();
0241 if (sd != nullptr) {
0242 G4VSensitiveDetector *aSD = sd->FindSensitiveDetector(sdname);
0243 if (aSD == nullptr) {
0244 edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: No SD"
0245 << " with name " << sdname << " in this "
0246 << "Setup";
0247 } else {
0248 HCalSD *theCaloSD = dynamic_cast<HCalSD *>(aSD);
0249 edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0250 << " in this Setup";
0251 if (org.get()) {
0252 theCaloSD->setNumberingScheme(org.get());
0253 edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::beginOfRun: set a new numbering scheme";
0254 }
0255 }
0256 } else {
0257 edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: Could "
0258 << "not get SD Manager!";
0259 }
0260 }
0261
0262
0263 void SimG4HcalValidation::update(const BeginOfEvent *evt) {
0264 int i = 0;
0265 edepEB = edepEE = edepHB = edepHE = edepHO = 0.;
0266 for (i = 0; i < 5; i++)
0267 edepd[i] = 0.;
0268 for (i = 0; i < 20; i++)
0269 edepl[i] = 0.;
0270 vhitec = vhithc = enEcal = enHcal = 0;
0271
0272 clear();
0273
0274 LogDebug("ValidHcal") << "SimG4HcalValidation: =====> Begin of event = " << (*evt)()->GetEventID();
0275 }
0276
0277
0278 void SimG4HcalValidation::update(const G4Step *aStep) {
0279 if (aStep != nullptr) {
0280 G4VPhysicalVolume *curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
0281 G4String name = curPV->GetName();
0282 name.assign(name, 0, 3);
0283 double edeposit = aStep->GetTotalEnergyDeposit();
0284 int layer = -1, depth = -1;
0285 if (name == "EBR") {
0286 depth = 0;
0287 edepEB += edeposit;
0288 } else if (name == "EFR") {
0289 depth = 0;
0290 edepEE += edeposit;
0291 } else if (name == "HBS") {
0292 layer = (curPV->GetCopyNo() / 10) % 100;
0293 depth = (curPV->GetCopyNo()) % 10 + 1;
0294 if (depth > 0 && depth < 4 && layer >= 0 && layer < 17) {
0295 edepHB += edeposit;
0296 } else {
0297 edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
0298 depth = -1;
0299 layer = -1;
0300 }
0301 } else if (name == "HES") {
0302 layer = (curPV->GetCopyNo() / 10) % 100;
0303 depth = (curPV->GetCopyNo()) % 10 + 1;
0304 if (depth > 0 && depth < 3 && layer >= 0 && layer < 19) {
0305 edepHE += edeposit;
0306 } else {
0307 edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
0308 depth = -1;
0309 layer = -1;
0310 }
0311 } else if (name == "HTS") {
0312 layer = (curPV->GetCopyNo() / 10) % 100;
0313 depth = (curPV->GetCopyNo()) % 10 + 1;
0314 if (depth > 3 && depth < 5 && layer >= 17 && layer < 20) {
0315 edepHO += edeposit;
0316 } else {
0317 edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
0318 depth = -1;
0319 layer = -1;
0320 }
0321 }
0322 if (depth >= 0 && depth < 5)
0323 edepd[depth] += edeposit;
0324 if (layer >= 0 && layer < 20)
0325 edepl[layer] += edeposit;
0326
0327 if (layer >= 0 && layer < 20) {
0328 LogDebug("ValidHcal") << "SimG4HcalValidation:: G4Step: " << name << " Layer " << std::setw(3) << layer
0329 << " Depth " << std::setw(2) << depth << " Edep " << std::setw(6) << edeposit / MeV
0330 << " MeV";
0331 }
0332 }
0333 }
0334
0335
0336 void SimG4HcalValidation::update(const EndOfEvent *evt) {
0337 count++;
0338
0339
0340 fill(evt);
0341 LogDebug("ValidHcal") << "SimG4HcalValidation:: --- after Fill ";
0342 }
0343
0344
0345 void SimG4HcalValidation::fill(const EndOfEvent *evt) {
0346 LogDebug("ValidHcal") << "SimG4HcalValidation:Fill event " << (*evt)()->GetEventID();
0347
0348
0349 G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
0350
0351 int nhc = 0, j = 0;
0352
0353
0354 int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[0]);
0355 CaloG4HitCollection *theHCHC = (CaloG4HitCollection *)allHC->GetHC(HCHCid);
0356 LogDebug("ValidHcal") << "SimG4HcalValidation :: Hit Collection for " << names[0] << " of ID " << HCHCid
0357 << " is obtained at " << theHCHC;
0358 if (HCHCid >= 0 && theHCHC != nullptr) {
0359 int nhits = theHCHC->entries();
0360 for (j = 0; j < nhits; j++) {
0361 CaloG4Hit *aHit = (*theHCHC)[j];
0362
0363 double e = aHit->getEnergyDeposit() / GeV;
0364 double time = aHit->getTimeSlice();
0365
0366 math::XYZPoint pos = aHit->getPosition();
0367 double theta = pos.theta();
0368 double eta = -log(tan(theta * 0.5));
0369 double phi = pos.phi();
0370
0371 uint32_t unitID = aHit->getUnitID();
0372 int subdet, zside, layer, etaIndex, phiIndex, lay;
0373 org->unpackHcalIndex(unitID, subdet, zside, layer, etaIndex, phiIndex, lay);
0374
0375
0376 layer--;
0377 std::string det = "HB";
0378 if (subdet == static_cast<int>(HcalForward)) {
0379 det = "HF";
0380 uint16_t depth = aHit->getDepth();
0381 if (depth != 0) {
0382 layer += 2;
0383 lay += 2;
0384 }
0385 if (layer == 1)
0386 vhithc += e;
0387 else if (layer == 0)
0388 vhitec += e;
0389 else
0390 edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::HitPMT " << subdet << " " << (2 * zside - 1) * etaIndex
0391 << " " << phiIndex << " " << layer + 1 << " R " << pos.rho() << " Phi "
0392 << convertRadToDeg(phi) << " Edep " << e << " Time " << time;
0393 } else if (subdet == static_cast<int>(HcalEndcap)) {
0394 if (etaIndex <= 20) {
0395 det = "HES";
0396 } else {
0397 det = "HED";
0398 }
0399 }
0400 LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Hcal " << det << " layer " << std::setw(2) << layer
0401 << " lay " << std::setw(2) << lay << " time " << std::setw(6) << time << " theta "
0402 << std::setw(8) << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
0403 << " e " << std::setw(8) << e << " ID 0x" << std::hex << unitID << " ID dec " << std::dec
0404 << (int)unitID;
0405
0406
0407 if (applySampling)
0408 e *= getHcalScale(det, layer);
0409
0410
0411 if (time >= timeLowlim && time <= timeUplim && e > hhitThreshold) {
0412 enHcal += e;
0413 CaloHit ahit(subdet, lay, e, eta, phi, time, unitID);
0414 hitcache.push_back(ahit);
0415 ++nhc;
0416 }
0417 }
0418 }
0419 LogDebug("ValidHcal") << "SimG4HcalValidation:: HCAL hits : " << nhc;
0420
0421 if (!hcalOnly) {
0422 int ndets = names.size();
0423 if (ndets > 3)
0424 ndets = 3;
0425 for (int idty = 1; idty < ndets; idty++) {
0426 std::string det = "EB";
0427 if (idty == 2)
0428 det = "EE";
0429 else if (idty > 2)
0430 det = "ES";
0431
0432 int nec = 0;
0433 int ECHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[idty]);
0434 CaloG4HitCollection *theECHC = (CaloG4HitCollection *)allHC->GetHC(ECHCid);
0435 LogDebug("ValidHcal") << "SimG4HcalValidation:: Hit Collection for " << names[idty] << " of ID " << ECHCid
0436 << " is obtained at " << theECHC;
0437 int theechc_entries = theECHC->entries();
0438 if (ECHCid >= 0 && theECHC != nullptr) {
0439 for (j = 0; j < theechc_entries; j++) {
0440 CaloG4Hit *aHit = (*theECHC)[j];
0441
0442 double e = aHit->getEnergyDeposit() / GeV;
0443 double time = aHit->getTimeSlice();
0444
0445 math::XYZPoint pos = aHit->getPosition();
0446 double theta = pos.theta();
0447 double eta = -log(tan(theta / 2.));
0448 double phi = pos.phi();
0449 HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta, phi, 1, 1);
0450 uint32_t unitID = org->getUnitID(id);
0451 int subdet, zside, layer, ieta, iphi, lay;
0452 org->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0453 subdet = idty + 9;
0454 layer = 0;
0455 unitID = org->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
0456
0457
0458 if (time >= timeLowlim && time <= timeUplim && e > ehitThreshold) {
0459 enEcal += e;
0460 CaloHit ahit(subdet, lay, e, eta, phi, time, unitID);
0461 hitcache.push_back(ahit);
0462 ++nec;
0463 }
0464
0465 LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Ecal " << det << " layer " << std::setw(2) << layer
0466 << " lay " << std::setw(2) << lay << " time " << std::setw(6) << time << " theta "
0467 << std::setw(8) << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8)
0468 << phi << " e " << std::setw(8) << e << " ID 0x" << std::hex << unitID << " ID dec "
0469 << std::dec << (int)unitID;
0470 }
0471 }
0472
0473 LogDebug("ValidHcal") << "SimG4HcalValidation:: " << det << " hits : " << nec;
0474 }
0475 }
0476 }
0477
0478 void SimG4HcalValidation::layerAnalysis(PHcalValidInfoLayer &product) {
0479 int i = 0;
0480 LogDebug("ValidHcal") << "\n ===>>> SimG4HcalValidation: Energy deposit "
0481 << "in MeV\n at EB : " << std::setw(6) << edepEB / MeV << "\n at EE : " << std::setw(6)
0482 << edepEE / MeV << "\n at HB : " << std::setw(6) << edepHB / MeV
0483 << "\n at HE : " << std::setw(6) << edepHE / MeV << "\n at HO : " << std::setw(6)
0484 << edepHO / MeV << "\n ---- SimG4HcalValidation: Energy deposit in";
0485 for (i = 0; i < 5; i++)
0486 LogDebug("ValidHcal") << " Depth " << std::setw(2) << i << " E " << std::setw(8) << edepd[i] / MeV << " MeV";
0487 LogDebug("ValidHcal") << " ---- SimG4HcalValidation: Energy deposit in"
0488 << "layers";
0489 for (i = 0; i < 20; i++)
0490 LogDebug("ValidHcal") << " Layer " << std::setw(2) << i << " E " << std::setw(8) << edepl[i] / MeV << " MeV";
0491
0492 product.fillLayers(edepl, edepd, edepHO, edepHB + edepHE, edepEB + edepEE);
0493
0494
0495 product.fillHF(vhitec, vhithc, enEcal, enHcal);
0496 LogDebug("ValidHcal") << "SimG4HcalValidation::HF hits " << vhitec << " in EC and " << vhithc << " in HC\n"
0497 << " HB/HE " << enEcal << " in EC and " << enHcal << " in HC";
0498
0499
0500 fetchHits(product);
0501
0502 LogDebug("ValidHcal") << " layerAnalysis ===> after fetchHits";
0503 }
0504
0505
0506 void SimG4HcalValidation::nxNAnalysis(PHcalValidInfoNxN &product) {
0507 std::vector<CaloHit> *hits = &hitcache;
0508 std::vector<CaloHit>::iterator hit_itr;
0509
0510 LogDebug("ValidHcal") << "SimG4HcalValidation::NxNAnalysis : entrance ";
0511
0512 collectEnergyRdir(eta0, phi0);
0513
0514
0515 LogDebug("ValidHcal") << " NxNAnalysis : coolectEnergyRdir - Ecal " << een << " Hcal " << hen;
0516
0517 double etot = 0.;
0518 double ee = 0.;
0519 double he = 0.;
0520 double hoe = 0.;
0521
0522 int max = dEta.size();
0523
0524 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
0525 double e = hit_itr->e();
0526 double t = hit_itr->t();
0527 double eta = hit_itr->eta();
0528 double phi = hit_itr->phi();
0529 int type = hit_itr->det();
0530 int layer = hit_itr->layer();
0531
0532
0533
0534 if (fabs(eta0 - eta) <= dEta[max - 1] && fabs(phi0 - phi) <= dPhi[max - 1]) {
0535 etot += e;
0536 if (type == 10 || type == 11 || type == 12)
0537 ee += e;
0538 if (type == static_cast<int>(HcalBarrel) || type == static_cast<int>(HcalEndcap) ||
0539 type == static_cast<int>(HcalForward)) {
0540 he += e;
0541 if (type == static_cast<int>(HcalBarrel) && layer > 17)
0542 hoe += e;
0543
0544
0545 for (int i = 0; i < max; i++) {
0546 if ((eta0 - eta) <= dEta[i] && (eta0 - eta) >= -dEta[i] && (phi0 - phi) <= dPhi[i] &&
0547 (phi0 - phi) >= -dPhi[i]) {
0548 LogDebug("ValidHcal") << "SimG4HcalValidation:: nxNAnalysis eta0,"
0549 << " phi0 = " << eta0 << " " << phi0 << " type, layer = " << type << "," << layer
0550 << " eta, phi = " << eta << " " << phi;
0551
0552 product.fillTProfileNxN(e, i, t);
0553 break;
0554 }
0555 }
0556
0557 }
0558 }
0559 }
0560
0561 product.fillEcollectNxN(ee, he, hoe, etot);
0562 product.fillHvsE(een, hen, hoen, een + hen);
0563
0564 LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
0565 }
0566
0567
0568 void SimG4HcalValidation::jetAnalysis(PHcalValidInfoJets &product) {
0569 std::vector<CaloHit> *hhit = &hitcache;
0570
0571 jetf->setInput(hhit);
0572
0573
0574 std::vector<SimG4HcalHitCluster> *result = jetf->getClusters(hcalOnly);
0575
0576 std::vector<SimG4HcalHitCluster>::iterator clus_itr;
0577
0578 LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size() << " clusters ---------------";
0579 for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
0580 LogDebug("ValidHcal") << (*clus_itr);
0581
0582 std::vector<double> enevec, etavec, phivec;
0583
0584 if (!(*result).empty()) {
0585 sort((*result).begin(), (*result).end());
0586
0587 clus_itr = result->begin();
0588 double etac = clus_itr->eta();
0589 double phic = clus_itr->phi();
0590
0591 double ecal_collect = 0.;
0592 if (!hcalOnly) {
0593 [[clang::suppress]] ecal_collect = clus_itr->collectEcalEnergyR();
0594 } else {
0595 collectEnergyRdir(etac, phic);
0596 [[clang::suppress]] ecal_collect = een;
0597 }
0598 LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect = " << ecal_collect;
0599
0600
0601 double dist = jetf->rDist(eta0, phi0, etac, phic);
0602 LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation = " << dist;
0603 product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
0604
0605 LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
0606
0607 std::vector<CaloHit> *hits = clus_itr->getHits();
0608 std::vector<CaloHit>::iterator hit_itr;
0609
0610 double ee = 0., he = 0., hoe = 0., etot = 0.;
0611
0612
0613 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
0614 double e = hit_itr->e();
0615 double t = hit_itr->t();
0616 double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
0617
0618
0619 etot += e;
0620 if (hit_itr->det() == 10 || hit_itr->det() == 11 || hit_itr->det() == 12)
0621 ee += e;
0622 if (hit_itr->det() == static_cast<int>(HcalBarrel) || hit_itr->det() == static_cast<int>(HcalEndcap) ||
0623 hit_itr->det() == static_cast<int>(HcalForward)) {
0624 he += e;
0625 if (hit_itr->det() == static_cast<int>(HcalBarrel) && hit_itr->layer() > 17)
0626 hoe += e;
0627 }
0628
0629 if (hit_itr->det() == static_cast<int>(HcalBarrel) || hit_itr->det() == static_cast<int>(HcalEndcap) ||
0630 hit_itr->det() == static_cast<int>(HcalForward)) {
0631 product.fillTProfileJet(he, r, t);
0632 }
0633 }
0634
0635 product.fillEcollectJet(ee, he, hoe, etot);
0636
0637 LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
0638 << "ee/he/hoe/etot " << ee << "/" << he << "/" << hoe << "/" << etot;
0639
0640
0641 for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
0642 if (clus_itr->e() > jetThreshold) {
0643 enevec.push_back(clus_itr->e());
0644 etavec.push_back(clus_itr->eta());
0645 phivec.push_back(clus_itr->phi());
0646 }
0647 }
0648 product.fillJets(enevec, etavec, phivec);
0649
0650 LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
0651 << " JetAnalysis ===> (*result).size() " << (*result).size();
0652
0653
0654 if (etavec.size() > 1) {
0655 if (etavec[0] > -2.5 && etavec[0] < 2.5 && etavec[1] > -2.5 && etavec[1] < 2.5) {
0656 LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
0657 << " JetAnalysis ===> Di-jet vectors ";
0658 for (unsigned int i = 0; i < enevec.size(); i++)
0659 LogDebug("ValidHcal") << " e, eta, phi = " << enevec[i] << " " << etavec[i] << " " << phivec[i];
0660
0661 double et0 = enevec[0] / cosh(etavec[0]);
0662 double px0 = et0 * cos(phivec[0]);
0663 double py0 = et0 * sin(phivec[0]);
0664 double pz0 = et0 * sinh(etavec[0]);
0665 double et1 = enevec[1] / cosh(etavec[1]);
0666 double px1 = et1 * cos(phivec[1]);
0667 double py1 = et1 * sin(phivec[1]);
0668 double pz1 = et1 * sinh(etavec[1]);
0669
0670 double dijetmass2 = (enevec[0] + enevec[1]) * (enevec[0] + enevec[1]) - (px1 + px0) * (px1 + px0) -
0671 (py1 + py0) * (py1 + py0) - (pz1 + pz0) * (pz1 + pz0);
0672
0673 LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ " << dijetmass2;
0674
0675 double dijetmass;
0676 if (dijetmass2 >= 0.)
0677 dijetmass = sqrt(dijetmass2);
0678 else
0679 dijetmass = -sqrt(-1. * dijetmass2);
0680
0681 product.fillDiJets(dijetmass);
0682
0683 LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
0684 }
0685 }
0686 }
0687 }
0688
0689
0690 void SimG4HcalValidation::fetchHits(PHcalValidInfoLayer &product) {
0691 LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with " << hitcache.size() << " hits";
0692 int nHit = hitcache.size();
0693 int hit = 0;
0694 int i;
0695 std::vector<CaloHit>::iterator itr;
0696 std::vector<CaloHit *> lhits(nHit);
0697 for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
0698 uint32_t unitID = itr->id();
0699 int subdet, zside, group, ieta, iphi, lay;
0700 HcalTestNumbering::unpackHcalIndex(unitID, subdet, zside, group, ieta, iphi, lay);
0701 subdet = itr->det();
0702 lay = itr->layer();
0703 group = (subdet & 15) << 20;
0704 group += ((lay - 1) & 31) << 15;
0705 group += (zside & 1) << 14;
0706 group += (ieta & 127) << 7;
0707 group += (iphi & 127);
0708 itr->setId(group);
0709 lhits[i] = &hitcache[i];
0710 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i << " " << hitcache[i] << "\n"
0711 << "SimG4HcalValidation::fetchHits:Copied " << i << " " << *lhits[i];
0712 }
0713 sort(lhits.begin(), lhits.end(), CaloHitIdMore());
0714 std::vector<CaloHit *>::iterator k1, k2;
0715 for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
0716 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i << " " << **k1;
0717 int nHits = 0;
0718 for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
0719 double ehit = (**k1).e();
0720 double t = (**k1).t();
0721 uint32_t unitID = (**k1).id();
0722 int jump = 0;
0723 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i << " U/T/E"
0724 << " 0x" << std::hex << unitID << std::dec << " " << t << " " << ehit;
0725 for (k2 = k1 + 1; k2 != lhits.end() && (t - (**k2).t()) < 1 && (t - (**k2).t()) > -1 && unitID == (**k2).id();
0726 k2++) {
0727 ehit += (**k2).e();
0728 LogDebug("ValidHcal") << "\t + " << (**k2).e();
0729 jump++;
0730 }
0731 LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
0732
0733 double eta = (*k1)->eta();
0734 double phi = (*k1)->phi();
0735 int lay = ((unitID >> 15) & 31) + 1;
0736 int subdet = (unitID >> 20) & 15;
0737
0738
0739 product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
0740 nHits++;
0741
0742 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits << " " << i << " ID 0x" << std::hex
0743 << unitID << " det " << std::dec << subdet << " " << lay << " " << ((unitID >> 14) & 1)
0744 << " " << ((unitID >> 7) & 127) << " " << ((unitID) & 127) << " Time " << t << " E " << ehit;
0745
0746 i += jump;
0747 k1 += jump;
0748 }
0749
0750 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with " << nHit << " hits"
0751 << " and writes out " << nHits << '(' << hit << ") hits";
0752 }
0753
0754 void SimG4HcalValidation::clear() { hitcache.erase(hitcache.begin(), hitcache.end()); }
0755
0756
0757 void SimG4HcalValidation::collectEnergyRdir(const double eta0, const double phi0) {
0758 std::vector<CaloHit> *hits = &hitcache;
0759 std::vector<CaloHit>::iterator hit_itr;
0760
0761 double sume = 0., sumh = 0., sumho = 0.;
0762
0763 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
0764 double e = hit_itr->e();
0765 double eta = hit_itr->eta();
0766 double phi = hit_itr->phi();
0767
0768 int type = hit_itr->det();
0769
0770 double r = jetf->rDist(eta0, phi0, eta, phi);
0771 if (r < coneSize) {
0772 if (type == 10 || type == 11 || type == 12) {
0773 sume += e;
0774 } else {
0775 sumh += e;
0776 if (type == static_cast<int>(HcalBarrel) && hit_itr->layer() > 17)
0777 sumho += e;
0778 }
0779 }
0780 }
0781
0782 een = sume;
0783 hen = sumh;
0784 hoen = sumho;
0785 }
0786
0787
0788 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
0789 double tmp = 0.;
0790
0791 if (det == "HB") {
0792 tmp = scaleHB[layer];
0793 } else if (det == "HES" || det == "HED") {
0794 tmp = scaleHE[layer];
0795 } else if (det == "HF") {
0796 tmp = scaleHF[layer];
0797 }
0798
0799 return tmp;
0800 }
0801
0802 #include "FWCore/Framework/interface/MakerMacros.h"
0803 #include "FWCore/PluginManager/interface/ModuleDef.h"
0804 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0805
0806 DEFINE_SIMWATCHER(SimG4HcalValidation);