File indexing completed on 2023-01-16 23:37:22
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include <cmath>
0017 #include <iomanip>
0018 #include <iostream>
0019 #include <memory>
0020 #include <vector>
0021 #include <string>
0022
0023
0024 #include "SimG4Core/Notification/interface/Observer.h"
0025 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0026 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0027 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0028 #include "SimG4Core/Watcher/interface/SimProducer.h"
0029 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0030
0031
0032 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0033 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
0034 #include "SimDataFormats/HcalTestBeam/interface/PHcalTB04Info.h"
0035 #include "SimG4CMS/Calo/interface/ECalSD.h"
0036 #include "SimG4CMS/Calo/interface/HCalSD.h"
0037 #include "SimG4CMS/Calo/interface/HcalQie.h"
0038 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0039 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0040 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0041 #include "DataFormats/Math/interface/Point3D.h"
0042
0043 #include "FWCore/Framework/interface/Event.h"
0044 #include "FWCore/Framework/interface/EventSetup.h"
0045 #include "FWCore/Framework/interface/MakerMacros.h"
0046 #include "FWCore/PluginManager/interface/ModuleDef.h"
0047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0048
0049 #include "SimG4CMS/HcalTestBeam/interface/HcalTBNumberingScheme.h"
0050 #include "SimG4CMS/HcalTestBeam/interface/HcalTB04Histo.h"
0051 #include "SimG4CMS/HcalTestBeam/interface/HcalTB04XtalNumberingScheme.h"
0052
0053 #include "G4SDManager.hh"
0054 #include "G4Step.hh"
0055 #include "G4Track.hh"
0056 #include "G4ThreeVector.hh"
0057 #include "G4VProcess.hh"
0058 #include "G4HCofThisEvent.hh"
0059
0060 #include <CLHEP/Random/RandGaussQ.h>
0061 #include <CLHEP/Random/Randomize.h>
0062 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0063 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0064
0065 #include <cstdint>
0066
0067
0068
0069 namespace CLHEP {
0070 class HepRandomEngine;
0071 }
0072
0073 class HcalTB04Analysis : public SimProducer,
0074 public Observer<const BeginOfRun*>,
0075 public Observer<const BeginOfEvent*>,
0076 public Observer<const EndOfEvent*>,
0077 public Observer<const G4Step*> {
0078 public:
0079 HcalTB04Analysis(const edm::ParameterSet& p);
0080 HcalTB04Analysis(const HcalTB04Analysis&) = delete;
0081 const HcalTB04Analysis& operator=(const HcalTB04Analysis&) = delete;
0082 ~HcalTB04Analysis() override;
0083
0084 void produce(edm::Event&, const edm::EventSetup&) override;
0085
0086 private:
0087 void init();
0088
0089
0090 void update(const BeginOfRun* run) override;
0091 void update(const BeginOfEvent* evt) override;
0092 void update(const G4Step* step) override;
0093 void update(const EndOfEvent* evt) override;
0094
0095
0096 void fillBuffer(const EndOfEvent* evt);
0097 void qieAnalysis(CLHEP::HepRandomEngine*);
0098 void xtalAnalysis(CLHEP::HepRandomEngine*);
0099 void finalAnalysis();
0100 void fillEvent(PHcalTB04Info&);
0101
0102 void clear();
0103 int unitID(uint32_t id);
0104 double scale(int det, int layer);
0105 double timeOfFlight(int det, int layer, double eta);
0106
0107 private:
0108
0109 const edm::ParameterSet m_Anal;
0110 const bool hcalOnly;
0111 const int mode, type;
0112 const double ecalNoise, beamOffset;
0113 const double scaleHB0, scaleHB16, scaleHO, scaleHE0;
0114 const std::vector<std::string> names;
0115
0116 HcalQie* myQie;
0117 HcalTB04Histo* histo;
0118
0119 int iceta, icphi;
0120 G4RotationMatrix* beamline_RM;
0121
0122
0123 int count;
0124 int nTower, nCrystal;
0125 std::vector<int> idHcal, idXtal;
0126 std::vector<uint32_t> idTower, idEcal;
0127
0128
0129 int nPrimary, particleType;
0130 double pInit, etaInit, phiInit;
0131 std::vector<CaloHit> ecalHitCache;
0132 std::vector<CaloHit> hcalHitCache, hcalHitLayer;
0133 std::vector<double> esimh, eqie, esime, enois;
0134 std::vector<double> eseta, eqeta, esphi, eqphi, eslay, eqlay;
0135 double etots, eecals, ehcals, etotq, eecalq, ehcalq;
0136
0137 bool pvFound;
0138 int evNum, pvType;
0139 G4ThreeVector pvPosition, pvMomentum, pvUVW;
0140 std::vector<int> secTrackID, secPartID;
0141 std::vector<G4ThreeVector> secMomentum;
0142 std::vector<double> secEkin;
0143 std::vector<int> shortLivedSecondaries;
0144 };
0145
0146
0147
0148
0149
0150 HcalTB04Analysis::HcalTB04Analysis(const edm::ParameterSet& p)
0151 : m_Anal(p.getParameter<edm::ParameterSet>("HcalTB04Analysis")),
0152 hcalOnly(m_Anal.getParameter<bool>("HcalOnly")),
0153 mode(m_Anal.getParameter<int>("Mode")),
0154 type(m_Anal.getParameter<int>("Type")),
0155 ecalNoise(m_Anal.getParameter<double>("EcalNoise")),
0156 beamOffset(-m_Anal.getParameter<double>("BeamPosition") * CLHEP::cm),
0157 scaleHB0(m_Anal.getParameter<double>("ScaleHB0")),
0158 scaleHB16(m_Anal.getParameter<double>("ScaleHB16")),
0159 scaleHO(m_Anal.getParameter<double>("ScaleHO")),
0160 scaleHE0(m_Anal.getParameter<double>("ScaleHE0")),
0161 names(m_Anal.getParameter<std::vector<std::string> >("Names")),
0162 myQie(nullptr),
0163 histo(nullptr) {
0164 double fMinEta = m_Anal.getParameter<double>("MinEta");
0165 double fMaxEta = m_Anal.getParameter<double>("MaxEta");
0166 double fMinPhi = m_Anal.getParameter<double>("MinPhi");
0167 double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
0168 double beamEta = (fMaxEta + fMinEta) / 2.;
0169 double beamPhi = (fMaxPhi + fMinPhi) / 2.;
0170 double beamThet = 2 * atan(exp(-beamEta));
0171 if (beamPhi < 0)
0172 beamPhi += twopi;
0173 iceta = static_cast<int>(beamEta / 0.087) + 1;
0174 icphi = static_cast<int>(std::fabs(beamPhi) / 0.087) + 5;
0175 if (icphi > 72)
0176 icphi -= 73;
0177
0178 produces<PHcalTB04Info>();
0179
0180 beamline_RM = new G4RotationMatrix;
0181 beamline_RM->rotateZ(-beamPhi);
0182 beamline_RM->rotateY(-beamThet);
0183
0184 edm::LogVerbatim("HcalTBSim")
0185 << "HcalTB04:: Initialised as observer of BeginOf Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent with Parameter "
0186 "values:\n \thcalOnly = "
0187 << hcalOnly << "\tecalNoise = " << ecalNoise << "\n\tMode = " << mode << " (0: HB2 Standard; 1:HB2 Segmented)"
0188 << "\tType = " << type << " (0: HB; 1 HE; 2 HB+HE)\n\tbeamOffset = " << beamOffset << "\ticeta = " << iceta
0189 << "\ticphi = " << icphi << "\n\tbeamline_RM = " << *beamline_RM;
0190
0191 init();
0192
0193 myQie = new HcalQie(p);
0194 histo = new HcalTB04Histo(m_Anal);
0195 }
0196
0197 HcalTB04Analysis::~HcalTB04Analysis() {
0198 #ifdef EDM_ML_DEBUG
0199 edm::LogVerbatim("HcalTBSim") << "\n --------> Total number of selected entries : " << count << "\nPointers:: QIE "
0200 << myQie << " Histo " << histo;
0201 #endif
0202 if (myQie) {
0203 delete myQie;
0204 myQie = nullptr;
0205 }
0206 if (histo) {
0207 delete histo;
0208 histo = nullptr;
0209 }
0210 }
0211
0212
0213
0214
0215
0216 void HcalTB04Analysis::produce(edm::Event& e, const edm::EventSetup&) {
0217 std::unique_ptr<PHcalTB04Info> product(new PHcalTB04Info);
0218 fillEvent(*product);
0219 e.put(std::move(product));
0220 }
0221
0222 void HcalTB04Analysis::init() {
0223 idTower = HcalTBNumberingScheme::getUnitIDs(type, mode);
0224 nTower = idTower.size();
0225 #ifdef EDM_ML_DEBUG
0226 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save information from " << nTower << " HCal towers";
0227 #endif
0228 idHcal.reserve(nTower);
0229 for (int i = 0; i < nTower; i++) {
0230 int id = unitID(idTower[i]);
0231 idHcal.push_back(id);
0232 #ifdef EDM_ML_DEBUG
0233 edm::LogVerbatim("HcalTBSim") << "\tTower[" << i << "] Original " << std::hex << idTower[i] << " Stored "
0234 << idHcal[i] << std::dec;
0235 #endif
0236 }
0237
0238 if (!hcalOnly) {
0239 int det = 10;
0240 uint32_t id1;
0241 nCrystal = 0;
0242 for (int lay = 1; lay < 8; lay++) {
0243 for (int icr = 1; icr < 8; icr++) {
0244 id1 = HcalTestNumbering::packHcalIndex(det, 0, 1, icr, lay, 1);
0245 int id = unitID(id1);
0246 idEcal.push_back(id1);
0247 idXtal.push_back(id);
0248 nCrystal++;
0249 }
0250 }
0251 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save information from " << nCrystal << " ECal Crystals";
0252 #ifdef EDM_ML_DEBUG
0253 for (int i = 0; i < nCrystal; i++) {
0254 edm::LogVerbatim("HcalTBSim") << "\tCrystal[" << i << "] Original " << std::hex << idEcal[i] << " Stored "
0255 << idXtal[i] << std::dec;
0256 }
0257 #endif
0258 }
0259
0260 eseta.reserve(5);
0261 eqeta.reserve(5);
0262 esphi.reserve(3);
0263 eqphi.reserve(3);
0264 eslay.reserve(20);
0265 eqlay.reserve(20);
0266 for (int i = 0; i < 5; i++) {
0267 eseta.push_back(0.);
0268 eqeta.push_back(0.);
0269 }
0270 for (int i = 0; i < 3; i++) {
0271 esphi.push_back(0.);
0272 eqphi.push_back(0.);
0273 }
0274 for (int i = 0; i < 20; i++) {
0275 eslay.push_back(0.);
0276 eqlay.push_back(0.);
0277 }
0278
0279
0280 count = 0;
0281 evNum = 0;
0282 clear();
0283 }
0284
0285 void HcalTB04Analysis::update(const BeginOfRun* run) {
0286 int irun = (*run)()->GetRunID();
0287 edm::LogVerbatim("HcalTBSim") << " =====> Begin of Run = " << irun;
0288
0289 G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
0290 if (sd != nullptr) {
0291 std::string sdname = names[0];
0292 G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
0293 if (aSD == nullptr) {
0294 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
0295 << " with name " << sdname << " in this "
0296 << "Setup";
0297 } else {
0298 HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
0299 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0300 << " in this Setup";
0301 HcalNumberingScheme* org = new HcalTestNumberingScheme(false);
0302 theCaloSD->setNumberingScheme(org);
0303 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a new numbering scheme";
0304 }
0305 if (!hcalOnly) {
0306 sdname = names[1];
0307 aSD = sd->FindSensitiveDetector(sdname);
0308 if (aSD == nullptr) {
0309 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
0310 << " with name " << sdname << " in this "
0311 << "Setup";
0312 } else {
0313 ECalSD* theCaloSD = dynamic_cast<ECalSD*>(aSD);
0314 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0315 << " in this Setup";
0316 EcalNumberingScheme* org = new HcalTB04XtalNumberingScheme();
0317 theCaloSD->setNumberingScheme(org);
0318 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a new numbering scheme";
0319 }
0320 }
0321 } else {
0322 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Could "
0323 << "not get SD Manager!";
0324 }
0325 }
0326
0327 void HcalTB04Analysis::update(const BeginOfEvent* evt) {
0328 clear();
0329 #ifdef EDM_ML_DEBUG
0330 evNum = (*evt)()->GetEventID();
0331 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis: =====> Begin of event = " << evNum;
0332 #endif
0333 }
0334
0335 void HcalTB04Analysis::update(const G4Step* aStep) {
0336 if (aStep != nullptr) {
0337
0338 G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
0339 G4ThreeVector thePostStepPoint;
0340
0341
0342 G4Track* aTrack = aStep->GetTrack();
0343 int trackID = aTrack->GetTrackID();
0344 int parentID = aTrack->GetParentID();
0345 const G4ThreeVector& position = aTrack->GetPosition();
0346 G4ThreeVector momentum = aTrack->GetMomentum();
0347 G4String partType = aTrack->GetDefinition()->GetParticleType();
0348 G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
0349 int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
0350 #ifdef EDM_ML_DEBUG
0351 bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
0352 #endif
0353 double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
0354 double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
0355
0356 if (!pvFound) {
0357 double stepDeltaEnergy = aStep->GetDeltaEnergy();
0358 double kinEnergy = aTrack->GetKineticEnergy();
0359
0360
0361 if (trackID == 1 && parentID == 0 && ((kinEnergy == 0.) || (std::fabs(stepDeltaEnergy / kinEnergy) > 0.1))) {
0362 pvType = -1;
0363 if (kinEnergy == 0.) {
0364 pvType = 0;
0365 } else {
0366 if (std::fabs(stepDeltaEnergy / kinEnergy) > 0.1)
0367 pvType = 1;
0368 }
0369 pvFound = true;
0370 pvPosition = position;
0371 pvMomentum = momentum;
0372
0373 pvUVW = (*beamline_RM) * (pvPosition);
0374
0375
0376 G4String thePostPVname = "NoName";
0377 G4StepPoint* thePostPoint = aStep->GetPostStepPoint();
0378 if (thePostPoint) {
0379 thePostStepPoint = thePostPoint->GetPosition();
0380 G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
0381 if (thePostPV)
0382 thePostPVname = thePostPV->GetName();
0383 }
0384 #ifdef EDM_ML_DEBUG
0385 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: V1 found at: " << thePostStepPoint
0386 << " G4VPhysicalVolume: " << thePostPVname;
0387 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::fill_v1Pos: Primary Track momentum: " << pvMomentum
0388 << " psoition " << pvPosition << " u/v/w " << pvUVW;
0389 #endif
0390 }
0391 } else {
0392
0393 if ((trackID != 1 && parentID == 1 && (aTrack->GetCurrentStepNumber() == 1) && (thePreStepPoint == pvPosition)) ||
0394 (trackID == 1 && thePreStepPoint == pvPosition)) {
0395 #ifdef EDM_ML_DEBUG
0396 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::A secondary... PDG:" << partPDGEncoding
0397 << " TrackID:" << trackID << " ParentID:" << parentID
0398 << " stable: " << isPDGStable << " Tau: " << pDGlifetime
0399 << " cTauGamma=" << c_light * pDGlifetime * gammaFactor * 1000.
0400 << "um GammaFactor: " << gammaFactor;
0401 #endif
0402 secTrackID.push_back(trackID);
0403 secPartID.push_back(partPDGEncoding);
0404 secMomentum.push_back(momentum);
0405 secEkin.push_back(aTrack->GetKineticEnergy());
0406
0407
0408 double ctaugamma_um = CLHEP::c_light * pDGlifetime * gammaFactor * 1000.;
0409 if ((ctaugamma_um > 0.) && (ctaugamma_um < 100.)) {
0410 shortLivedSecondaries.push_back(trackID);
0411 } else {
0412
0413 }
0414 }
0415
0416
0417 if (aTrack->GetCurrentStepNumber() == 1) {
0418 if (!shortLivedSecondaries.empty()) {
0419 int pid = parentID;
0420 std::vector<int>::iterator pos1 = shortLivedSecondaries.begin();
0421 std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
0422 std::vector<int>::iterator pos;
0423 for (pos = pos1; pos != pos2; pos++) {
0424 if (*pos == pid) {
0425
0426 #ifdef EDM_ML_DEBUG
0427 edm::LogVerbatim("HcalTBSim")
0428 << "HcalTB04Analysis:: A tertiary... PDG:" << partPDGEncoding << " TrackID:" << trackID
0429 << " ParentID:" << parentID << " stable: " << isPDGStable << " Tau: " << pDGlifetime
0430 << " cTauGamma=" << c_light * pDGlifetime * gammaFactor * 1000. << "um GammaFactor: " << gammaFactor;
0431 #endif
0432 }
0433 }
0434 }
0435 }
0436 }
0437 }
0438 }
0439
0440 void HcalTB04Analysis::update(const EndOfEvent* evt) {
0441 count++;
0442
0443
0444 #ifdef EDM_ML_DEBUG
0445 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Fill event " << (*evt)()->GetEventID();
0446 #endif
0447 fillBuffer(evt);
0448
0449
0450 #ifdef EDM_ML_DEBUG
0451 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Do QIE analysis with " << hcalHitCache.size() << " hits";
0452 #endif
0453 CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
0454 qieAnalysis(engine);
0455
0456
0457 if (!hcalOnly) {
0458 #ifdef EDM_ML_DEBUG
0459 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Do Xtal analysis with " << ecalHitCache.size() << " hits";
0460 #endif
0461 xtalAnalysis(engine);
0462 }
0463
0464
0465 #ifdef EDM_ML_DEBUG
0466 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Final analysis";
0467 #endif
0468 finalAnalysis();
0469
0470 int iEvt = (*evt)()->GetEventID();
0471 if (iEvt < 10)
0472 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0473 else if ((iEvt < 100) && (iEvt % 10 == 0))
0474 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0475 else if ((iEvt < 1000) && (iEvt % 100 == 0))
0476 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0477 else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0478 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0479 }
0480
0481 void HcalTB04Analysis::fillBuffer(const EndOfEvent* evt) {
0482 std::vector<CaloHit> hhits, hhitl;
0483 int idHC, j;
0484 CaloG4HitCollection* theHC;
0485 std::map<int, float, std::less<int> > primaries;
0486 #ifdef EDM_ML_DEBUG
0487 double etot1 = 0, etot2 = 0;
0488 #endif
0489
0490
0491 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0492 std::string sdName = names[0];
0493 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0494 theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0495 #ifdef EDM_ML_DEBUG
0496 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0497 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0498 #endif
0499 int thehc_entries = theHC->entries();
0500 if (idHC >= 0 && theHC != nullptr) {
0501 hhits.reserve(theHC->entries());
0502 hhitl.reserve(theHC->entries());
0503 for (j = 0; j < thehc_entries; j++) {
0504 CaloG4Hit* aHit = (*theHC)[j];
0505 double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0506 double time = aHit->getTimeSlice();
0507 math::XYZPoint pos = aHit->getEntry();
0508 unsigned int id = aHit->getUnitID();
0509 double theta = pos.theta();
0510 double eta = -std::log(std::tan(theta * 0.5));
0511 double phi = pos.phi();
0512 int det, z, group, ieta, iphi, layer;
0513 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0514 double jitter = time - timeOfFlight(det, layer, eta);
0515 if (jitter < 0)
0516 jitter = 0;
0517 if (e < 0 || e > 1.)
0518 e = 0;
0519 double escl = e * scale(det, layer);
0520 unsigned int idx = HcalTBNumberingScheme::getUnitID(id, mode);
0521 CaloHit hit(det, layer, escl, eta, phi, jitter, idx);
0522 hhits.push_back(hit);
0523 CaloHit hitl(det, layer, escl, eta, phi, jitter, id);
0524 hhitl.push_back(hitl);
0525 primaries[aHit->getTrackID()] += e;
0526 #ifdef EDM_ML_DEBUG
0527 etot1 += escl;
0528 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j << " ID 0x" << std::hex << id << " 0x"
0529 << idx << std::dec << " time " << std::setw(6) << time << " " << std::setw(6)
0530 << jitter << " theta " << std::setw(8) << theta << " eta " << std::setw(8) << eta
0531 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << e << " "
0532 << std::setw(8) << escl;
0533 #endif
0534 }
0535 }
0536
0537
0538 std::vector<CaloHit>::iterator itr;
0539 int nHit = hhits.size();
0540 std::vector<CaloHit*> hits(nHit);
0541 for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
0542 hits[j] = &hhits[j];
0543 }
0544 sort(hits.begin(), hits.end(), CaloHitIdMore());
0545 std::vector<CaloHit*>::iterator k1, k2;
0546 int nhit = 0;
0547 for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0548 int det = (**k1).det();
0549 int layer = (**k1).layer();
0550 double ehit = (**k1).e();
0551 double eta = (**k1).eta();
0552 double phi = (**k1).phi();
0553 double jitter = (**k1).t();
0554 uint32_t unitID = (**k1).id();
0555 int jump = 0;
0556 for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0557 ehit += (**k2).e();
0558 jump++;
0559 }
0560 nhit++;
0561 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0562 hcalHitCache.push_back(hit);
0563 k1 += jump;
0564 #ifdef EDM_ML_DEBUG
0565 etot2 += ehit;
0566 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit << " ID 0x" << std::hex << unitID
0567 << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0568 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0569 #endif
0570 }
0571 #ifdef EDM_ML_DEBUG
0572 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits from " << nHit
0573 << " input hits E(Hcal) " << etot1 << " " << etot2;
0574 #endif
0575
0576 for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
0577 hits[j] = &hhitl[j];
0578 }
0579 sort(hits.begin(), hits.end(), CaloHitIdMore());
0580 int nhitl = 0;
0581 #ifdef EDM_ML_DEBUG
0582 double etotl = 0;
0583 #endif
0584 for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0585 int det = (**k1).det();
0586 int layer = (**k1).layer();
0587 double ehit = (**k1).e();
0588 double eta = (**k1).eta();
0589 double phi = (**k1).phi();
0590 double jitter = (**k1).t();
0591 uint32_t unitID = (**k1).id();
0592 int jump = 0;
0593 for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0594 ehit += (**k2).e();
0595 jump++;
0596 }
0597 nhitl++;
0598 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0599 hcalHitLayer.push_back(hit);
0600 #ifdef EDM_ML_DEBUG
0601 etotl += ehit;
0602 #endif
0603 k1 += jump;
0604 #ifdef EDM_ML_DEBUG
0605 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl << " ID 0x" << std::hex << unitID
0606 << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0607 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0608 #endif
0609 }
0610 #ifdef EDM_ML_DEBUG
0611 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal hits from " << nHit
0612 << " input hits E(Hcal) " << etot1 << " " << etotl;
0613 #endif
0614
0615 std::vector<CaloHit> ehits;
0616 sdName = names[1];
0617 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0618 theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0619 #ifdef EDM_ML_DEBUG
0620 etot1 = etot2 = 0;
0621 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0622 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0623 #endif
0624 if (idHC >= 0 && theHC != nullptr) {
0625 thehc_entries = theHC->entries();
0626 ehits.reserve(theHC->entries());
0627 for (j = 0; j < thehc_entries; j++) {
0628 CaloG4Hit* aHit = (*theHC)[j];
0629 double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0630 double time = aHit->getTimeSlice();
0631 if (e < 0 || e > 100000.)
0632 e = 0;
0633 if (e > 0) {
0634 math::XYZPoint pos = aHit->getEntry();
0635 unsigned int id = aHit->getUnitID();
0636 double theta = pos.theta();
0637 double eta = -std::log(std::tan(theta * 0.5));
0638 double phi = pos.phi();
0639 int det, z, group, ieta, iphi, layer;
0640 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0641 CaloHit hit(det, 0, e, eta, phi, time, id);
0642 ehits.push_back(hit);
0643 primaries[aHit->getTrackID()] += e;
0644 #ifdef EDM_ML_DEBUG
0645 etot1 += e;
0646 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j << " ID 0x" << std::hex << id
0647 << std::dec << " time " << std::setw(6) << time << " theta " << std::setw(8)
0648 << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
0649 << " e " << std::setw(8) << e;
0650 #endif
0651 }
0652 }
0653 }
0654
0655
0656 nHit = ehits.size();
0657 std::vector<CaloHit*> hite(nHit);
0658 for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
0659 hite[j] = &ehits[j];
0660 }
0661 sort(hite.begin(), hite.end(), CaloHitIdMore());
0662 nhit = 0;
0663 for (k1 = hite.begin(); k1 != hite.end(); k1++) {
0664 int det = (**k1).det();
0665 int layer = (**k1).layer();
0666 double ehit = (**k1).e();
0667 double eta = (**k1).eta();
0668 double phi = (**k1).phi();
0669 double jitter = (**k1).t();
0670 uint32_t unitID = (**k1).id();
0671 int jump = 0;
0672 for (k2 = k1 + 1; k2 != hite.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0673 ehit += (**k2).e();
0674 jump++;
0675 }
0676 nhit++;
0677 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0678 ecalHitCache.push_back(hit);
0679 k1 += jump;
0680 #ifdef EDM_ML_DEBUG
0681 etot2 += ehit;
0682 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit << " ID 0x" << std::hex << unitID
0683 << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0684 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0685 #endif
0686 }
0687 #ifdef EDM_ML_DEBUG
0688 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits from " << nHit
0689 << " input hits E(Ecal) " << etot1 << " " << etot2;
0690 #endif
0691
0692 nPrimary = static_cast<int>(primaries.size());
0693 int trackID = 0;
0694 G4PrimaryParticle* thePrim = nullptr;
0695 int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0696 #ifdef EDM_ML_DEBUG
0697 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex << " verteices";
0698 #endif
0699 if (nvertex <= 0)
0700 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no vertex found for event " << evNum;
0701 for (int i = 0; i < nvertex; i++) {
0702 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0703 if (avertex == nullptr) {
0704 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer to vertex = 0 for event " << evNum;
0705 } else {
0706 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " " << avertex->GetPosition();
0707 int npart = avertex->GetNumberOfParticle();
0708 if (npart == 0)
0709 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: no primary!";
0710 if (thePrim == nullptr)
0711 thePrim = avertex->GetPrimary(trackID);
0712 }
0713 }
0714
0715 if (thePrim != nullptr) {
0716 double px = thePrim->GetPx();
0717 double py = thePrim->GetPy();
0718 double pz = thePrim->GetPz();
0719 double p = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0720 pInit = p / CLHEP::GeV;
0721 if (p == 0)
0722 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: primary has p=0 ";
0723 else {
0724 double costheta = pz / p;
0725 double theta = acos(std::min(std::max(costheta, -1.), 1.));
0726 etaInit = -std::log(std::tan(theta / 2));
0727 if (px != 0 || py != 0)
0728 phiInit = std::atan2(py, px);
0729 }
0730 particleType = thePrim->GetPDGcode();
0731 } else
0732 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could not find primary";
0733 }
0734
0735 void HcalTB04Analysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
0736 int hittot = hcalHitCache.size();
0737 if (hittot <= 0)
0738 hittot = 1;
0739 std::vector<CaloHit> hits(hittot);
0740 std::vector<int> todo(nTower, 0);
0741
0742 #ifdef EDM_ML_DEBUG
0743 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size " << hits.size() << " " << todo.size() << " "
0744 << idTower.size() << " " << esimh.size() << " " << eqie.size();
0745 #endif
0746
0747 for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
0748 CaloHit hit = hcalHitCache[k1];
0749 uint32_t id = hit.id();
0750 int nhit = 0;
0751 double esim = hit.e();
0752 hits[nhit] = hit;
0753 for (unsigned int k2 = k1 + 1; k2 < hcalHitCache.size(); k2++) {
0754 hit = hcalHitCache[k2];
0755 if (hit.id() == id) {
0756 nhit++;
0757 hits[nhit] = hit;
0758 esim += hit.e();
0759 }
0760 }
0761 k1 += nhit;
0762 nhit++;
0763 std::vector<int> cd = myQie->getCode(nhit, hits, engine);
0764 double eq = myQie->getEnergy(cd);
0765 #ifdef EDM_ML_DEBUG
0766 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id << std::dec << " registers " << esim
0767 << " energy from " << nhit << " hits starting with hit # " << k1
0768 << " energy with noise " << eq;
0769 #endif
0770 for (int k2 = 0; k2 < nTower; k2++) {
0771 if (id == idTower[k2]) {
0772 todo[k2] = 1;
0773 esimh[k2] = esim;
0774 eqie[k2] = eq;
0775 }
0776 }
0777 }
0778
0779
0780 for (int k2 = 0; k2 < nTower; k2++) {
0781 if (todo[k2] == 0) {
0782 std::vector<int> cd = myQie->getCode(0, hits, engine);
0783 double eq = myQie->getEnergy(cd);
0784 esimh[k2] = 0;
0785 eqie[k2] = eq;
0786 #ifdef EDM_ML_DEBUG
0787 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << idTower[k2] << std::dec
0788 << " registers " << esimh[k2] << " energy from hits and energy after QIE analysis "
0789 << eqie[k2];
0790 #endif
0791 }
0792 }
0793 }
0794
0795 void HcalTB04Analysis::xtalAnalysis(CLHEP::HepRandomEngine* engine) {
0796 CLHEP::RandGaussQ randGauss(*engine);
0797
0798
0799 std::vector<int> iok(nCrystal, 0);
0800 #ifdef EDM_ML_DEBUG
0801 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " << iok.size() << " " << idEcal.size() << " "
0802 << esime.size() << " " << enois.size();
0803 #endif
0804 for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
0805 uint32_t id = ecalHitCache[k1].id();
0806 int nhit = 0;
0807 double esim = ecalHitCache[k1].e();
0808 for (unsigned int k2 = k1 + 1; k2 < ecalHitCache.size(); k2++) {
0809 if (ecalHitCache[k2].id() == id) {
0810 nhit++;
0811 esim += ecalHitCache[k2].e();
0812 }
0813 }
0814 k1 += nhit;
0815 nhit++;
0816 double eq = esim + randGauss.fire(0., ecalNoise);
0817 #ifdef EDM_ML_DEBUG
0818 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id << std::dec << " registers " << esim
0819 << " energy from " << nhit << " hits starting with hit # " << k1
0820 << " energy with noise " << eq;
0821 #endif
0822 for (int k2 = 0; k2 < nCrystal; k2++) {
0823 if (id == idEcal[k2]) {
0824 iok[k2] = 1;
0825 esime[k2] = esim;
0826 enois[k2] = eq;
0827 }
0828 }
0829 }
0830
0831
0832 for (int k2 = 0; k2 < nCrystal; k2++) {
0833 if (iok[k2] == 0) {
0834 esime[k2] = 0;
0835 enois[k2] = randGauss.fire(0., ecalNoise);
0836 #ifdef EDM_ML_DEBUG
0837 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << idEcal[k2] << std::dec
0838 << " registers " << esime[k2] << " energy from hits and energy from noise "
0839 << enois[k2];
0840 #endif
0841 }
0842 }
0843 }
0844
0845 void HcalTB04Analysis::finalAnalysis() {
0846
0847 histo->fillPrimary(pInit, etaInit, phiInit);
0848
0849
0850 eecals = ehcals = eecalq = ehcalq = 0.;
0851 for (int i = 0; i < nTower; i++) {
0852 ehcals += esimh[i];
0853 ehcalq += eqie[i];
0854 }
0855 for (int i = 0; i < nCrystal; i++) {
0856 eecals += esime[i];
0857 eecalq += enois[i];
0858 }
0859 etots = eecals + ehcals;
0860 etotq = eecalq + ehcalq;
0861 #ifdef EDM_ML_DEBUG
0862 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level (Total) " << etots << " (ECal) "
0863 << eecals << " (HCal) " << ehcals
0864 << "\nHcalTB04Analysis:: Energy deposit at Qie Level (Total) " << etotq << " (ECal) "
0865 << eecalq << " (HCal) " << ehcalq;
0866 #endif
0867 histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0868
0869
0870 for (int i = 0; i < 5; i++) {
0871 eseta[i] = 0.;
0872 eqeta[i] = 0.;
0873 }
0874 for (int i = 0; i < 3; i++) {
0875 esphi[i] = 0.;
0876 eqphi[i] = 0.;
0877 }
0878 double e1 = 0, e2 = 0;
0879 unsigned int id;
0880 for (int i = 0; i < nTower; i++) {
0881 int det, z, group, ieta, iphi, layer;
0882 id = idTower[i];
0883 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0884 iphi -= (icphi - 1);
0885 if (icphi > 4) {
0886 if (ieta == 0)
0887 ieta = 2;
0888 else
0889 ieta = -1;
0890 } else {
0891 ieta = ieta - iceta + 2;
0892 }
0893 if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
0894 eseta[ieta] += esimh[i];
0895 esphi[iphi] += esimh[i];
0896 e1 += esimh[i];
0897 eqeta[ieta] += eqie[i];
0898 eqphi[iphi] += eqie[i];
0899 e2 += eqie[i];
0900 }
0901 }
0902 for (int i = 0; i < 3; i++) {
0903 if (e1 > 0)
0904 esphi[i] /= e1;
0905 if (e2 > 0)
0906 eqphi[i] /= e2;
0907 }
0908 for (int i = 0; i < 5; i++) {
0909 if (e1 > 0)
0910 eseta[i] /= e1;
0911 if (e2 > 0)
0912 eqeta[i] /= e2;
0913 }
0914 #ifdef EDM_ML_DEBUG
0915 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and Phi (Sim/Qie)";
0916 for (int i = 0; i < 5; i++)
0917 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = " << eseta[i] << " Qie = " << eqeta[i]
0918 << " Phi Sim = " << esphi[i] << " Qie = " << eqphi[i];
0919 #endif
0920 histo->fillTrnsProf(eseta, eqeta, esphi, eqphi);
0921
0922
0923 for (int i = 0; i < 20; i++) {
0924 eslay[i] = 0.;
0925 eqlay[i] = 0.;
0926 }
0927 e1 = 0;
0928 e2 = 0;
0929 for (int i = 0; i < nTower; i++) {
0930 int det, z, group, ieta, iphi, layer;
0931 id = idTower[i];
0932 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0933 iphi -= (icphi - 1);
0934 layer -= 1;
0935 if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
0936 eslay[layer] += esimh[i];
0937 e1 += esimh[i];
0938 eqlay[layer] += eqie[i];
0939 e2 += eqie[i];
0940 }
0941 }
0942 for (int i = 0; i < 20; i++) {
0943 if (e1 > 0)
0944 eslay[i] /= e1;
0945 if (e2 > 0)
0946 eqlay[i] /= e2;
0947 }
0948 #ifdef EDM_ML_DEBUG
0949 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
0950 for (int i = 0; i < 20; i++)
0951 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = " << eslay[i] << " Qie = " << eqlay[i];
0952 #endif
0953 histo->fillLongProf(eslay, eqlay);
0954 }
0955
0956 void HcalTB04Analysis::fillEvent(PHcalTB04Info& product) {
0957
0958 product.setIDs(idHcal, idXtal);
0959
0960
0961 product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
0962
0963
0964 product.setEdepHcal(esimh, eqie);
0965 product.setEdepHcal(esime, enois);
0966
0967
0968 product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0969
0970
0971 product.setTrnsProf(eseta, eqeta, esphi, eqphi);
0972
0973
0974 product.setLongProf(eslay, eqlay);
0975
0976
0977 int i, nhit = 0;
0978 std::vector<CaloHit>::iterator itr;
0979 for (i = 0, itr = ecalHitCache.begin(); itr != ecalHitCache.end(); i++, itr++) {
0980 uint32_t id = itr->id();
0981 int det, z, group, ieta, iphi, lay;
0982 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
0983 product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
0984 nhit++;
0985 #ifdef EDM_ML_DEBUG
0986 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
0987 << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
0988 << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
0989 << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
0990 #endif
0991 }
0992 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from Crystals";
0993 int hit = nhit;
0994 nhit = 0;
0995
0996 for (i = hit, itr = hcalHitCache.begin(); itr != hcalHitCache.end(); i++, itr++) {
0997 uint32_t id = itr->id();
0998 int det, z, group, ieta, iphi, lay;
0999 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1000 product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
1001 nhit++;
1002 #ifdef EDM_ML_DEBUG
1003 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
1004 << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
1005 << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
1006 << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
1007 #endif
1008 }
1009 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from HCal";
1010
1011
1012 product.setVtxPrim(evNum,
1013 pvType,
1014 pvPosition.x(),
1015 pvPosition.y(),
1016 pvPosition.z(),
1017 pvUVW.x(),
1018 pvUVW.y(),
1019 pvUVW.z(),
1020 pvMomentum.x(),
1021 pvMomentum.y(),
1022 pvMomentum.z());
1023 for (unsigned int i = 0; i < secTrackID.size(); i++) {
1024 product.setVtxSec(
1025 secTrackID[i], secPartID[i], secMomentum[i].x(), secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
1026 }
1027 }
1028
1029 void HcalTB04Analysis::clear() {
1030 pvFound = false;
1031 pvType = -2;
1032 pvPosition = G4ThreeVector();
1033 pvMomentum = G4ThreeVector();
1034 pvUVW = G4ThreeVector();
1035 secTrackID.clear();
1036 secPartID.clear();
1037 secMomentum.clear();
1038 secEkin.clear();
1039 shortLivedSecondaries.clear();
1040
1041 ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
1042 hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
1043 hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
1044 nPrimary = particleType = 0;
1045 pInit = etaInit = phiInit = 0;
1046
1047 esimh.clear();
1048 eqie.clear();
1049 esimh.reserve(nTower);
1050 eqie.reserve(nTower);
1051 for (int i = 0; i < nTower; i++) {
1052 esimh.push_back(0.);
1053 eqie.push_back(0.);
1054 }
1055 esime.clear();
1056 enois.clear();
1057 esime.reserve(nCrystal);
1058 enois.reserve(nCrystal);
1059 for (int i = 0; i < nCrystal; i++) {
1060 esime.push_back(0.);
1061 enois.push_back(0.);
1062 }
1063 }
1064
1065 int HcalTB04Analysis::unitID(uint32_t id) {
1066 int det, z, group, ieta, iphi, lay;
1067 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1068 group = (det & 15) << 20;
1069 group += ((lay - 1) & 31) << 15;
1070 group += (z & 1) << 14;
1071 group += (ieta & 127) << 7;
1072 group += (iphi & 127);
1073 return group;
1074 }
1075
1076 double HcalTB04Analysis::scale(int det, int layer) {
1077 double tmp = 1.;
1078 if (det == static_cast<int>(HcalBarrel)) {
1079 if (layer == 1)
1080 tmp = scaleHB0;
1081 else if (layer == 17)
1082 tmp = scaleHB16;
1083 else if (layer > 17)
1084 tmp = scaleHO;
1085 } else {
1086 if (layer <= 2)
1087 tmp = scaleHE0;
1088 }
1089 return tmp;
1090 }
1091
1092 double HcalTB04Analysis::timeOfFlight(int det, int layer, double eta) {
1093 double theta = 2.0 * std::atan(std::exp(-eta));
1094 double dist = beamOffset;
1095 if (det == static_cast<int>(HcalBarrel)) {
1096 const double rLay[19] = {1836.0,
1097 1902.0,
1098 1962.0,
1099 2022.0,
1100 2082.0,
1101 2142.0,
1102 2202.0,
1103 2262.0,
1104 2322.0,
1105 2382.0,
1106 2448.0,
1107 2514.0,
1108 2580.0,
1109 2646.0,
1110 2712.0,
1111 2776.0,
1112 2862.5,
1113 3847.0,
1114 4052.0};
1115 if (layer > 0 && layer <= 19)
1116 dist += rLay[layer - 1] * mm / sin(theta);
1117 } else {
1118 const double zLay[19] = {4034.0,
1119 4032.0,
1120 4123.0,
1121 4210.0,
1122 4297.0,
1123 4384.0,
1124 4471.0,
1125 4558.0,
1126 4645.0,
1127 4732.0,
1128 4819.0,
1129 4906.0,
1130 4993.0,
1131 5080.0,
1132 5167.0,
1133 5254.0,
1134 5341.0,
1135 5428.0,
1136 5515.0};
1137 if (layer > 0 && layer <= 19)
1138 dist += zLay[layer - 1] * mm / cos(theta);
1139 }
1140
1141 double tmp = dist / c_light / ns;
1142 #ifdef EDM_ML_DEBUG
1143 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
1144 << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
1145 #endif
1146 return tmp;
1147 }
1148
1149 DEFINE_SIMWATCHER(HcalTB04Analysis);