File indexing completed on 2023-03-17 11:24:12
0001 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0002 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0003 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0004 #include "SimG4Core/Notification/interface/Observer.h"
0005 #include "SimG4Core/Watcher/interface/SimProducer.h"
0006
0007
0008 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0009 #include "DataFormats/Math/interface/Point3D.h"
0010 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
0011 #include "SimDataFormats/CaloTest/interface/HcalTestHistoClass.h"
0012
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015
0016 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
0017 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0018 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
0019
0020 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0021 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0022 #include "SimG4CMS/Calo/interface/HCalSD.h"
0023 #include "SimG4CMS/Calo/interface/HcalQie.h"
0024 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0025
0026 #include "G4SDManager.hh"
0027 #include "G4Step.hh"
0028 #include "G4Track.hh"
0029 #include "G4VProcess.hh"
0030 #include "G4HCofThisEvent.hh"
0031
0032 #include <CLHEP/Random/Randomize.h>
0033 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0034 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0035
0036 #include <cmath>
0037 #include <iomanip>
0038 #include <iostream>
0039 #include <string>
0040 #include <vector>
0041
0042 class HcalTestAnalysis : public SimProducer,
0043 public Observer<const BeginOfRun*>,
0044 public Observer<const BeginOfEvent*>,
0045 public Observer<const EndOfEvent*>,
0046 public Observer<const G4Step*> {
0047 public:
0048 HcalTestAnalysis(const edm::ParameterSet& p);
0049 ~HcalTestAnalysis() override;
0050
0051 void registerConsumes(edm::ConsumesCollector) override;
0052 void produce(edm::Event&, const edm::EventSetup&) override;
0053 void beginRun(edm::EventSetup const&) override;
0054
0055 private:
0056
0057 void update(const BeginOfRun* run) override;
0058 void update(const BeginOfEvent* evt) override;
0059 void update(const EndOfEvent* evt) override;
0060 void update(const G4Step* step) override;
0061
0062
0063 std::vector<int> layerGrouping(int);
0064 std::vector<int> towersToAdd(int centre, int nadd);
0065 void fill(const EndOfEvent* ev);
0066 void qieAnalysis(CLHEP::HepRandomEngine*);
0067 void layerAnalysis();
0068 double timeOfFlight(int det, int layer, double eta);
0069
0070 private:
0071
0072 std::unique_ptr<HcalQie> myqie_;
0073 int addTower_;
0074
0075
0076 std::unique_ptr<HcalTestHistoClass> tuples_;
0077
0078
0079 const edm::ParameterSet m_Anal;
0080 double eta0_, phi0_;
0081 const int laygroup_, centralTower_;
0082 const std::vector<std::string> names_;
0083
0084 const std::string fileName_;
0085
0086
0087 edm::ESGetToken<HcalDDDSimConstants, HcalSimNumberingRecord> ddconsToken_;
0088 std::unique_ptr<HcalNumberingFromDDD> numberingFromDDD_;
0089 const HcalDDDSimConstants* hcons_;
0090 std::unique_ptr<HcalTestNumberingScheme> org_;
0091
0092
0093 std::vector<CaloHit> caloHitCache_;
0094 std::vector<int> group_, tower_;
0095 int nGroup_, nTower_;
0096
0097
0098 unsigned int count_;
0099 double edepEB_, edepEE_, edepHB_, edepHE_;
0100 double edepHO_, edepl_[20];
0101 double mudist_[20];
0102 };
0103
0104 HcalTestAnalysis::HcalTestAnalysis(const edm::ParameterSet& p)
0105 : addTower_(3),
0106 m_Anal(p.getParameter<edm::ParameterSet>("HcalTestAnalysis")),
0107 eta0_(m_Anal.getParameter<double>("Eta0")),
0108 phi0_(m_Anal.getParameter<double>("Phi0")),
0109 laygroup_(m_Anal.getParameter<int>("LayerGrouping")),
0110 centralTower_(m_Anal.getParameter<int>("CentralTower")),
0111 names_(m_Anal.getParameter<std::vector<std::string> >("Names")),
0112 fileName_(m_Anal.getParameter<std::string>("FileName")),
0113 hcons_(nullptr) {
0114 org_.reset(nullptr);
0115
0116 tuples_.reset(nullptr);
0117 numberingFromDDD_.reset(nullptr);
0118 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Initialised as observer of begin/end events"
0119 << " and of G4step";
0120
0121 count_ = 0;
0122 group_ = layerGrouping(laygroup_);
0123 nGroup_ = 0;
0124 for (unsigned int i = 0; i < group_.size(); i++)
0125 if (group_[i] > nGroup_)
0126 nGroup_ = group_[i];
0127 tower_ = towersToAdd(centralTower_, addTower_);
0128 nTower_ = tower_.size() / 2;
0129
0130 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: initialised for " << nGroup_ << " Longitudinal groups and "
0131 << nTower_ << " towers";
0132
0133
0134 myqie_ = std::make_unique<HcalQie>(p);
0135
0136 produces<HcalTestHistoClass>();
0137 }
0138
0139 HcalTestAnalysis::~HcalTestAnalysis() {
0140 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: --------> Total number of selected entries : " << count_;
0141 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Pointers:: Numbering Scheme " << org_.get();
0142 }
0143
0144 void HcalTestAnalysis::registerConsumes(edm::ConsumesCollector cc) {
0145 ddconsToken_ = cc.esConsumes<HcalDDDSimConstants, HcalSimNumberingRecord, edm::Transition::BeginRun>();
0146 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Initialize ESGetToken for HcalDDDSimConstants";
0147 }
0148
0149 void HcalTestAnalysis::produce(edm::Event& e, const edm::EventSetup&) {
0150 e.put(std::move(tuples_));
0151 tuples_.reset(nullptr);
0152 }
0153
0154 std::vector<int> HcalTestAnalysis::layerGrouping(int group) {
0155 std::vector<int> temp(19);
0156 if (group <= 1) {
0157 int grp[19] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2};
0158 for (int i = 0; i < 19; i++)
0159 temp[i] = grp[i];
0160 } else if (group == 2) {
0161 int grp[19] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
0162 for (int i = 0; i < 19; i++)
0163 temp[i] = grp[i];
0164 } else if (group == 3) {
0165 int grp[19] = {1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7};
0166 for (int i = 0; i < 19; i++)
0167 temp[i] = grp[i];
0168 } else if (group == 4) {
0169 int grp[19] = {1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7};
0170 for (int i = 0; i < 19; i++)
0171 temp[i] = grp[i];
0172 } else {
0173 int grp[19] = {1, 1, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7};
0174 for (int i = 0; i < 19; i++)
0175 temp[i] = grp[i];
0176 }
0177
0178 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Layer Grouping ";
0179 for (int i = 0; i < 19; i++)
0180 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Group[" << i << "] = " << temp[i];
0181 return temp;
0182 }
0183
0184 std::vector<int> HcalTestAnalysis::towersToAdd(int centre, int nadd) {
0185 int etac = (centre / 100) % 100;
0186 int phic = (centre % 100);
0187 int etamin, etamax, phimin, phimax;
0188 if (etac > 0) {
0189 etamin = etac - nadd;
0190 etamax = etac + nadd;
0191 } else {
0192 etamin = etac;
0193 etamax = etac;
0194 }
0195 if (phic > 0) {
0196 phimin = phic - nadd;
0197 phimax = phic + nadd;
0198 } else {
0199 phimin = phic;
0200 phimax = phic;
0201 }
0202
0203 int nbuf, kount = 0;
0204 nbuf = (etamax - etamin + 1) * (phimax - phimin + 1);
0205 std::vector<int> temp(2 * nbuf);
0206 for (int eta = etamin; eta <= etamax; eta++) {
0207 for (int phi = phimin; phi <= phimax; phi++) {
0208 temp[kount] = (eta * 100 + phi);
0209 temp[kount + nbuf] = std::max(abs(eta - etac), abs(phi - phic));
0210 kount++;
0211 }
0212 }
0213
0214 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Towers to be considered for Central " << centre << " and " << nadd
0215 << " on either side";
0216 for (int i = 0; i < nbuf; i++)
0217 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Tower[" << std::setw(3) << i << "] " << temp[i] << " "
0218 << temp[nbuf + i];
0219 return temp;
0220 }
0221
0222
0223
0224 void HcalTestAnalysis::beginRun(edm::EventSetup const& es) {
0225
0226 hcons_ = &es.getData(ddconsToken_);
0227 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Initialise HcalNumberingFromDDD for " << names_[0];
0228 numberingFromDDD_ = std::make_unique<HcalNumberingFromDDD>(hcons_);
0229
0230
0231 org_ = std::make_unique<HcalTestNumberingScheme>(false);
0232 }
0233
0234
0235 void HcalTestAnalysis::update(const BeginOfRun* run) {
0236 int irun = (*run)()->GetRunID();
0237 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Begin of Run = " << irun;
0238
0239 bool loop = true, eta = true, phi = true;
0240 int etac = (centralTower_ / 100) % 100;
0241 if (etac == 0) {
0242 etac = 1;
0243 eta = false;
0244 }
0245 int phic = (centralTower_ % 100);
0246 if (phic == 0) {
0247 phic = 1;
0248 phi = false;
0249 }
0250 int idet = static_cast<int>(HcalBarrel);
0251 while (loop) {
0252 HcalCellType::HcalCell tmp = hcons_->cell(idet, 1, 1, etac, phic);
0253 if (tmp.ok) {
0254 if (eta)
0255 eta0_ = tmp.eta;
0256 if (phi)
0257 phi0_ = tmp.phi;
0258 loop = false;
0259 } else if (idet == static_cast<int>(HcalBarrel)) {
0260 idet = static_cast<int>(HcalEndcap);
0261 } else if (idet == static_cast<int>(HcalEndcap)) {
0262 idet = static_cast<int>(HcalForward);
0263 } else {
0264 loop = false;
0265 }
0266 }
0267
0268 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Central Tower " << centralTower_
0269 << " corresponds to eta0 = " << eta0_ << " phi0 = " << phi0_;
0270
0271 std::string sdname = names_[0];
0272 G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
0273 if (sd != nullptr) {
0274 G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
0275 if (aSD == nullptr) {
0276 edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: No SD with "
0277 << "name " << sdname << " in this Setup";
0278 } else {
0279 HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
0280 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0281 << " in this Setup";
0282 if (org_.get()) {
0283 theCaloSD->setNumberingScheme(org_.get());
0284 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::beginOfRun: set a new numbering scheme";
0285 }
0286 }
0287 } else {
0288 edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: Could not get"
0289 << " SD Manager!";
0290 }
0291 }
0292
0293
0294 void HcalTestAnalysis::update(const BeginOfEvent* evt) {
0295
0296 tuples_ = std::make_unique<HcalTestHistoClass>();
0297
0298 tuples_->setCounters();
0299
0300 int i = 0;
0301 edepEB_ = edepEE_ = edepHB_ = edepHE_ = edepHO_ = 0.;
0302 for (i = 0; i < 20; i++)
0303 edepl_[i] = 0.;
0304 for (i = 0; i < 20; i++)
0305 mudist_[i] = -1.;
0306
0307 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Begin of event = " << (*evt)()->GetEventID();
0308 }
0309
0310
0311 void HcalTestAnalysis::update(const G4Step* aStep) {
0312 if (aStep != nullptr) {
0313 G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
0314 G4String name = curPV->GetName();
0315 name.assign(name, 0, 3);
0316 double edeposit = aStep->GetTotalEnergyDeposit();
0317 int layer = -1;
0318 if (name == "EBR") {
0319 edepEB_ += edeposit;
0320 } else if (name == "EFR") {
0321 edepEE_ += edeposit;
0322 } else if (name == "HBS") {
0323 layer = (curPV->GetCopyNo() / 10) % 100;
0324 if (layer >= 0 && layer < 17) {
0325 edepHB_ += edeposit;
0326 } else {
0327 edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HB " << curPV->GetName() << curPV->GetCopyNo();
0328 layer = -1;
0329 }
0330 } else if (name == "HES") {
0331 layer = (curPV->GetCopyNo() / 10) % 100;
0332 if (layer >= 0 && layer < 19) {
0333 edepHE_ += edeposit;
0334 } else {
0335 edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HE " << curPV->GetName() << curPV->GetCopyNo();
0336 layer = -1;
0337 }
0338 } else if (name == "HTS") {
0339 layer = (curPV->GetCopyNo() / 10) % 100;
0340 if (layer >= 17 && layer < 20) {
0341 edepHO_ += edeposit;
0342 } else {
0343 edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HO " << curPV->GetName() << curPV->GetCopyNo();
0344 layer = -1;
0345 }
0346 }
0347 if (layer >= 0 && layer < 20) {
0348 edepl_[layer] += edeposit;
0349
0350
0351 G4String part = aStep->GetTrack()->GetDefinition()->GetParticleName();
0352 if ((part == "mu-" || part == "mu+") && mudist_[layer] < 0) {
0353 math::XYZPoint pos(aStep->GetPreStepPoint()->GetPosition().x(),
0354 aStep->GetPreStepPoint()->GetPosition().y(),
0355 aStep->GetPreStepPoint()->GetPosition().z());
0356 double theta = pos.theta();
0357 double eta = -log(tan(theta * 0.5));
0358 double phi = pos.phi();
0359 double dist = sqrt((eta - eta0_) * (eta - eta0_) + (phi - phi0_) * (phi - phi0_));
0360 mudist_[layer] = dist * std::sqrt(pos.perp2());
0361 }
0362 }
0363
0364 if (layer >= 0 && layer < 20) {
0365 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: G4Step: " << name << " Layer " << std::setw(3) << layer
0366 << " Edep " << std::setw(6) << edeposit / MeV << " MeV";
0367 }
0368 } else {
0369 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: G4Step: Null Step";
0370 }
0371 }
0372
0373
0374 void HcalTestAnalysis::update(const EndOfEvent* evt) {
0375 ++count_;
0376
0377 fill(evt);
0378 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after Fill";
0379
0380
0381 CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
0382 qieAnalysis(engine);
0383 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after QieAnalysis";
0384
0385
0386 layerAnalysis();
0387 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after LayerAnalysis";
0388 }
0389
0390
0391 void HcalTestAnalysis::fill(const EndOfEvent* evt) {
0392 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Fill event " << (*evt)()->GetEventID();
0393
0394
0395 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0396
0397 int nhc = 0, neb = 0, nef = 0, j = 0;
0398 caloHitCache_.erase(caloHitCache_.begin(), caloHitCache_.end());
0399
0400
0401 int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[0]);
0402 CaloG4HitCollection* theHCHC = (CaloG4HitCollection*)allHC->GetHC(HCHCid);
0403 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[0] << " of ID " << HCHCid
0404 << " is obtained at " << theHCHC;
0405 int hchc_entries = theHCHC->entries();
0406 if (HCHCid >= 0 && theHCHC != nullptr) {
0407 for (j = 0; j < hchc_entries; j++) {
0408 CaloG4Hit* aHit = (*theHCHC)[j];
0409
0410 double e = aHit->getEnergyDeposit() / GeV;
0411 double time = aHit->getTimeSlice();
0412
0413 math::XYZPoint pos = aHit->getPosition();
0414 double theta = pos.theta();
0415 double eta = -log(tan(theta * 0.5));
0416 double phi = pos.phi();
0417
0418 uint32_t unitID = aHit->getUnitID();
0419 int subdet, zside, layer, etaIndex, phiIndex, lay;
0420 org_->unpackHcalIndex(unitID, subdet, zside, layer, etaIndex, phiIndex, lay);
0421 double jitter = time - timeOfFlight(subdet, lay, eta);
0422 if (jitter < 0)
0423 jitter = 0;
0424 CaloHit hit(subdet, lay, e, eta, phi, jitter, unitID);
0425 caloHitCache_.push_back(hit);
0426 nhc++;
0427
0428 std::string det = "HB";
0429 if (subdet == static_cast<int>(HcalForward)) {
0430 det = "HF";
0431 } else if (subdet == static_cast<int>(HcalEndcap)) {
0432 if (etaIndex <= 20) {
0433 det = "HES";
0434 } else {
0435 det = "HED";
0436 }
0437 }
0438
0439 edm::LogVerbatim("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2) << layer << " time "
0440 << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
0441 << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
0442 << e;
0443 }
0444 }
0445 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::HCAL hits : " << nhc;
0446
0447
0448 int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[1]);
0449 CaloG4HitCollection* theEBHC = (CaloG4HitCollection*)allHC->GetHC(EBHCid);
0450 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[1] << " of ID " << EBHCid
0451 << " is obtained at " << theEBHC;
0452 int ebhc_entries = theEBHC->entries();
0453 if (EBHCid >= 0 && theEBHC != nullptr) {
0454 for (j = 0; j < ebhc_entries; j++) {
0455 CaloG4Hit* aHit = (*theEBHC)[j];
0456
0457 double e = aHit->getEnergyDeposit() / GeV;
0458 double time = aHit->getTimeSlice();
0459 std::string det = "EB";
0460
0461 math::XYZPoint pos = aHit->getPosition();
0462 double theta = pos.theta();
0463 double eta = -log(tan(theta / 2.));
0464 double phi = pos.phi();
0465
0466 HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta, phi, 1, 1);
0467 uint32_t unitID = org_->getUnitID(id);
0468 int subdet, zside, layer, ieta, iphi, lay;
0469 org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0470 subdet = 10;
0471 layer = 0;
0472 unitID = org_->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
0473 CaloHit hit(subdet, lay, e, eta, phi, time, unitID);
0474 caloHitCache_.push_back(hit);
0475 neb++;
0476 edm::LogVerbatim("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2) << layer << " time "
0477 << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
0478 << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
0479 << e;
0480 }
0481 }
0482 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::EB hits : " << neb;
0483
0484
0485 int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[2]);
0486 CaloG4HitCollection* theEEHC = (CaloG4HitCollection*)allHC->GetHC(EEHCid);
0487 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[2] << " of ID " << EEHCid
0488 << " is obtained at " << theEEHC;
0489 int eehc_entries = theEEHC->entries();
0490 if (EEHCid >= 0 && theEEHC != nullptr) {
0491 for (j = 0; j < eehc_entries; j++) {
0492 CaloG4Hit* aHit = (*theEEHC)[j];
0493
0494 double e = aHit->getEnergyDeposit() / GeV;
0495 double time = aHit->getTimeSlice();
0496 std::string det = "EE";
0497
0498 math::XYZPoint pos = aHit->getPosition();
0499 double theta = pos.theta();
0500 double eta = -log(tan(theta / 2.));
0501 double phi = pos.phi();
0502
0503 HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta, phi, 1, 1);
0504 uint32_t unitID = org_->getUnitID(id);
0505 int subdet, zside, layer, ieta, iphi, lay;
0506 org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0507 subdet = 11;
0508 layer = 0;
0509 unitID = org_->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
0510 CaloHit hit(subdet, lay, e, eta, phi, time, unitID);
0511 caloHitCache_.push_back(hit);
0512 nef++;
0513 edm::LogVerbatim("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2) << layer << " time "
0514 << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
0515 << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
0516 << e;
0517 }
0518 }
0519 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::EE hits : " << nef;
0520 }
0521
0522
0523 void HcalTestAnalysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
0524
0525 int hittot = caloHitCache_.size();
0526 tuples_->fillHits(caloHitCache_);
0527
0528
0529 HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta0_, phi0_, 1, 1);
0530 uint32_t unitID = org_->getUnitID(id);
0531 int subdet, zside, layer, ieta, iphi, lay;
0532 org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0533 int laymax = 0;
0534 std::string det = "Unknown";
0535 if (subdet == static_cast<int>(HcalBarrel)) {
0536 laymax = 4;
0537 det = "HB";
0538 } else if (subdet == static_cast<int>(HcalEndcap)) {
0539 laymax = 2;
0540 det = "HES";
0541 }
0542 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: " << det << " Eta " << ieta << " Phi " << iphi << " Laymax "
0543 << laymax << " Hits " << hittot;
0544
0545 if (laymax > 0 && hittot > 0) {
0546 std::vector<CaloHit> hits(hittot);
0547 std::vector<double> eqielay(80, 0.0), esimlay(80, 0.0), esimtot(4, 0.0);
0548 std::vector<double> eqietow(200, 0.0), esimtow(200, 0.0), eqietot(4, 0.0);
0549 int etac = (centralTower_ / 100) % 100;
0550 int phic = (centralTower_ % 100);
0551
0552 for (int layr = 0; layr < nGroup_; layr++) {
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567 for (int it = 0; it < nTower_; it++) {
0568 int nhit = 0;
0569 double esim = 0;
0570 for (int k1 = 0; k1 < hittot; k1++) {
0571 CaloHit hit = caloHitCache_[k1];
0572 int subdetc = hit.det();
0573 int layer = hit.layer();
0574 int group = 0;
0575 if (layer > 0 && layer < 20)
0576 group = group_[layer];
0577 if (subdetc == subdet && group == layr + 1) {
0578 int zsidec, ietac, iphic, idx;
0579 unitID = hit.id();
0580 org_->unpackHcalIndex(unitID, subdetc, zsidec, layer, ietac, iphic, lay);
0581 if (etac > 0 && phic > 0) {
0582 idx = ietac * 100 + iphic;
0583 } else if (etac > 0) {
0584 idx = ietac * 100;
0585 } else if (phic > 0) {
0586 idx = iphic;
0587 } else {
0588 idx = 0;
0589 }
0590 if (zsidec == zside && idx == tower_[it]) {
0591 hits[nhit] = hit;
0592 edm::LogVerbatim("HcalSim") << "HcalTest: Hit " << nhit << " " << hit;
0593 nhit++;
0594 esim += hit.e();
0595 }
0596 }
0597 }
0598
0599 std::vector<int> cd = myqie_->getCode(nhit, hits, engine);
0600 double eqie = myqie_->getEnergy(cd);
0601
0602 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: Energy in layer " << layr << " Sim " << esim
0603 << " After QIE " << eqie;
0604 for (int i = 0; i < 4; i++) {
0605 if (tower_[nTower_ + it] <= i) {
0606 esimtot[i] += esim;
0607 eqietot[i] += eqie;
0608 esimlay[20 * i + layr] += esim;
0609 eqielay[20 * i + layr] += eqie;
0610 esimtow[50 * i + it] += esim;
0611 eqietow[50 * i + it] += eqie;
0612 }
0613 }
0614 }
0615 }
0616 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: Total energy " << esimtot[3] << " (SimHit) " << eqietot[3]
0617 << " (After QIE)";
0618
0619 std::vector<double> latphi(10);
0620 int nt = 2 * addTower_ + 1;
0621 for (int it = 0; it < nt; it++)
0622 latphi[it] = it - addTower_;
0623 for (int i = 0; i < 4; i++) {
0624 double scals = 1, scalq = 1;
0625 std::vector<double> latfs(10, 0.), latfq(10, 0.), longs(20), longq(20);
0626 if (esimtot[i] > 0)
0627 scals = 1. / esimtot[i];
0628 if (eqietot[i] > 0)
0629 scalq = 1. / eqietot[i];
0630 for (int it = 0; it < nTower_; it++) {
0631 int phib = it % nt;
0632 latfs[phib] += scals * esimtow[50 * i + it];
0633 latfq[phib] += scalq * eqietow[50 * i + it];
0634 }
0635 for (int layr = 0; layr <= nGroup_; layr++) {
0636 longs[layr] = scals * esimlay[20 * i + layr];
0637 longq[layr] = scalq * eqielay[20 * i + layr];
0638 }
0639 tuples_->fillQie(i, esimtot[i], eqietot[i], nGroup_, longs, longq, nt, latphi, latfs, latfq);
0640 }
0641 }
0642 }
0643
0644
0645 void HcalTestAnalysis::layerAnalysis() {
0646 int i = 0;
0647 edm::LogVerbatim("HcalSim") << "\n ===>>> HcalTestAnalysis: Energy deposit "
0648 << "\n at EB : " << std::setw(6) << edepEB_ / MeV << "\n at EE : " << std::setw(6)
0649 << edepEE_ / MeV << "\n at HB : " << std::setw(6) << edepHB_ / MeV
0650 << "\n at HE : " << std::setw(6) << edepHE_ / MeV << "\n at HO : " << std::setw(6)
0651 << edepHO_ / MeV << "\n ---- HcalTestAnalysis: Energy deposit in Layers";
0652 for (i = 0; i < 20; i++)
0653 edm::LogVerbatim("HcalSim") << " Layer " << std::setw(2) << i << " E " << std::setw(8) << edepl_[i] / MeV << " MeV";
0654
0655 tuples_->fillLayers(edepl_, edepHO_, edepHB_ + edepHE_, mudist_);
0656 }
0657
0658
0659 double HcalTestAnalysis::timeOfFlight(int det, int layer, double eta) {
0660 double theta = 2.0 * atan(exp(-eta));
0661 double dist = 0.;
0662 if (det == static_cast<int>(HcalBarrel)) {
0663 const double rLay[19] = {1836.0,
0664 1902.0,
0665 1962.0,
0666 2022.0,
0667 2082.0,
0668 2142.0,
0669 2202.0,
0670 2262.0,
0671 2322.0,
0672 2382.0,
0673 2448.0,
0674 2514.0,
0675 2580.0,
0676 2646.0,
0677 2712.0,
0678 2776.0,
0679 2862.5,
0680 3847.0,
0681 4052.0};
0682 if (layer > 0 && layer < 20)
0683 dist += rLay[layer - 1] * mm / sin(theta);
0684 } else {
0685 const double zLay[19] = {4034.0,
0686 4032.0,
0687 4123.0,
0688 4210.0,
0689 4297.0,
0690 4384.0,
0691 4471.0,
0692 4558.0,
0693 4645.0,
0694 4732.0,
0695 4819.0,
0696 4906.0,
0697 4993.0,
0698 5080.0,
0699 5167.0,
0700 5254.0,
0701 5341.0,
0702 5428.0,
0703 5515.0};
0704 if (layer > 0 && layer < 20)
0705 dist += zLay[layer - 1] * mm / cos(theta);
0706 }
0707 double tmp = dist / c_light / ns;
0708 edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
0709 << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
0710 return tmp;
0711 }
0712
0713 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0714 #include "FWCore/PluginManager/interface/ModuleDef.h"
0715
0716 DEFINE_SIMWATCHER(HcalTestAnalysis);