File indexing completed on 2022-02-23 01:39:02
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 double etot1 = 0, etot2 = 0;
0487
0488
0489 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0490 std::string sdName = names[0];
0491 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0492 theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0493 #ifdef EDM_ML_DEBUG
0494 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0495 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0496 #endif
0497 int thehc_entries = theHC->entries();
0498 if (idHC >= 0 && theHC != nullptr) {
0499 hhits.reserve(theHC->entries());
0500 hhitl.reserve(theHC->entries());
0501 for (j = 0; j < thehc_entries; j++) {
0502 CaloG4Hit* aHit = (*theHC)[j];
0503 double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0504 double time = aHit->getTimeSlice();
0505 math::XYZPoint pos = aHit->getEntry();
0506 unsigned int id = aHit->getUnitID();
0507 double theta = pos.theta();
0508 double eta = -std::log(std::tan(theta * 0.5));
0509 double phi = pos.phi();
0510 int det, z, group, ieta, iphi, layer;
0511 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0512 double jitter = time - timeOfFlight(det, layer, eta);
0513 if (jitter < 0)
0514 jitter = 0;
0515 if (e < 0 || e > 1.)
0516 e = 0;
0517 double escl = e * scale(det, layer);
0518 unsigned int idx = HcalTBNumberingScheme::getUnitID(id, mode);
0519 CaloHit hit(det, layer, escl, eta, phi, jitter, idx);
0520 hhits.push_back(hit);
0521 CaloHit hitl(det, layer, escl, eta, phi, jitter, id);
0522 hhitl.push_back(hitl);
0523 primaries[aHit->getTrackID()] += e;
0524 etot1 += escl;
0525 #ifdef EDM_ML_DEBUG
0526 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j << " ID 0x" << std::hex << id << " 0x"
0527 << idx << std::dec << " time " << std::setw(6) << time << " " << std::setw(6)
0528 << jitter << " theta " << std::setw(8) << theta << " eta " << std::setw(8) << eta
0529 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << e << " "
0530 << std::setw(8) << escl;
0531 #endif
0532 }
0533 }
0534
0535
0536 std::vector<CaloHit>::iterator itr;
0537 int nHit = hhits.size();
0538 std::vector<CaloHit*> hits(nHit);
0539 for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
0540 hits[j] = &hhits[j];
0541 }
0542 sort(hits.begin(), hits.end(), CaloHitIdMore());
0543 std::vector<CaloHit*>::iterator k1, k2;
0544 int nhit = 0;
0545 for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0546 int det = (**k1).det();
0547 int layer = (**k1).layer();
0548 double ehit = (**k1).e();
0549 double eta = (**k1).eta();
0550 double phi = (**k1).phi();
0551 double jitter = (**k1).t();
0552 uint32_t unitID = (**k1).id();
0553 int jump = 0;
0554 for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0555 ehit += (**k2).e();
0556 jump++;
0557 }
0558 nhit++;
0559 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0560 hcalHitCache.push_back(hit);
0561 etot2 += ehit;
0562 k1 += jump;
0563 #ifdef EDM_ML_DEBUG
0564 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit << " ID 0x" << std::hex << unitID
0565 << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0566 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0567 #endif
0568 }
0569 #ifdef EDM_ML_DEBUG
0570 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits from " << nHit
0571 << " input hits E(Hcal) " << etot1 << " " << etot2;
0572 #endif
0573
0574 for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
0575 hits[j] = &hhitl[j];
0576 }
0577 sort(hits.begin(), hits.end(), CaloHitIdMore());
0578 int nhitl = 0;
0579 double etotl = 0;
0580 for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0581 int det = (**k1).det();
0582 int layer = (**k1).layer();
0583 double ehit = (**k1).e();
0584 double eta = (**k1).eta();
0585 double phi = (**k1).phi();
0586 double jitter = (**k1).t();
0587 uint32_t unitID = (**k1).id();
0588 int jump = 0;
0589 for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0590 ehit += (**k2).e();
0591 jump++;
0592 }
0593 nhitl++;
0594 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0595 hcalHitLayer.push_back(hit);
0596 etotl += ehit;
0597 k1 += jump;
0598 #ifdef EDM_ML_DEBUG
0599 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl << " ID 0x" << std::hex << unitID
0600 << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0601 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0602 #endif
0603 }
0604 #ifdef EDM_ML_DEBUG
0605 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal hits from " << nHit
0606 << " input hits E(Hcal) " << etot1 << " " << etotl;
0607 #endif
0608
0609 std::vector<CaloHit> ehits;
0610 sdName = names[1];
0611 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0612 theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0613 etot1 = etot2 = 0;
0614 #ifdef EDM_ML_DEBUG
0615 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0616 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0617 #endif
0618 if (idHC >= 0 && theHC != nullptr) {
0619 thehc_entries = theHC->entries();
0620 ehits.reserve(theHC->entries());
0621 for (j = 0; j < thehc_entries; j++) {
0622 CaloG4Hit* aHit = (*theHC)[j];
0623 double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0624 double time = aHit->getTimeSlice();
0625 if (e < 0 || e > 100000.)
0626 e = 0;
0627 if (e > 0) {
0628 math::XYZPoint pos = aHit->getEntry();
0629 unsigned int id = aHit->getUnitID();
0630 double theta = pos.theta();
0631 double eta = -std::log(std::tan(theta * 0.5));
0632 double phi = pos.phi();
0633 int det, z, group, ieta, iphi, layer;
0634 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0635 CaloHit hit(det, 0, e, eta, phi, time, id);
0636 ehits.push_back(hit);
0637 primaries[aHit->getTrackID()] += e;
0638 etot1 += e;
0639 #ifdef EDM_ML_DEBUG
0640 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j << " ID 0x" << std::hex << id
0641 << std::dec << " time " << std::setw(6) << time << " theta " << std::setw(8)
0642 << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
0643 << " e " << std::setw(8) << e;
0644 #endif
0645 }
0646 }
0647 }
0648
0649
0650 nHit = ehits.size();
0651 std::vector<CaloHit*> hite(nHit);
0652 for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
0653 hite[j] = &ehits[j];
0654 }
0655 sort(hite.begin(), hite.end(), CaloHitIdMore());
0656 nhit = 0;
0657 for (k1 = hite.begin(); k1 != hite.end(); k1++) {
0658 int det = (**k1).det();
0659 int layer = (**k1).layer();
0660 double ehit = (**k1).e();
0661 double eta = (**k1).eta();
0662 double phi = (**k1).phi();
0663 double jitter = (**k1).t();
0664 uint32_t unitID = (**k1).id();
0665 int jump = 0;
0666 for (k2 = k1 + 1; k2 != hite.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0667 ehit += (**k2).e();
0668 jump++;
0669 }
0670 nhit++;
0671 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0672 ecalHitCache.push_back(hit);
0673 etot2 += ehit;
0674 k1 += jump;
0675 #ifdef EDM_ML_DEBUG
0676 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit << " ID 0x" << std::hex << unitID
0677 << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0678 << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0679 #endif
0680 }
0681 #ifdef EDM_ML_DEBUG
0682 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits from " << nHit
0683 << " input hits E(Ecal) " << etot1 << " " << etot2;
0684 #endif
0685
0686 nPrimary = static_cast<int>(primaries.size());
0687 int trackID = 0;
0688 G4PrimaryParticle* thePrim = nullptr;
0689 int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0690 #ifdef EDM_ML_DEBUG
0691 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex << " verteices";
0692 #endif
0693 if (nvertex <= 0)
0694 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no vertex found for event " << evNum;
0695 for (int i = 0; i < nvertex; i++) {
0696 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0697 if (avertex == nullptr) {
0698 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer to vertex = 0 for event " << evNum;
0699 } else {
0700 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " " << avertex->GetPosition();
0701 int npart = avertex->GetNumberOfParticle();
0702 if (npart == 0)
0703 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: no primary!";
0704 if (thePrim == nullptr)
0705 thePrim = avertex->GetPrimary(trackID);
0706 }
0707 }
0708
0709 if (thePrim != nullptr) {
0710 double px = thePrim->GetPx();
0711 double py = thePrim->GetPy();
0712 double pz = thePrim->GetPz();
0713 double p = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0714 pInit = p / CLHEP::GeV;
0715 if (p == 0)
0716 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: primary has p=0 ";
0717 else {
0718 double costheta = pz / p;
0719 double theta = acos(std::min(std::max(costheta, -1.), 1.));
0720 etaInit = -std::log(std::tan(theta / 2));
0721 if (px != 0 || py != 0)
0722 phiInit = std::atan2(py, px);
0723 }
0724 particleType = thePrim->GetPDGcode();
0725 } else
0726 edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could not find primary";
0727 }
0728
0729 void HcalTB04Analysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
0730 int hittot = hcalHitCache.size();
0731 if (hittot <= 0)
0732 hittot = 1;
0733 std::vector<CaloHit> hits(hittot);
0734 std::vector<int> todo(nTower, 0);
0735
0736 #ifdef EDM_ML_DEBUG
0737 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size " << hits.size() << " " << todo.size() << " "
0738 << idTower.size() << " " << esimh.size() << " " << eqie.size();
0739 #endif
0740
0741 for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
0742 CaloHit hit = hcalHitCache[k1];
0743 uint32_t id = hit.id();
0744 int nhit = 0;
0745 double esim = hit.e();
0746 hits[nhit] = hit;
0747 for (unsigned int k2 = k1 + 1; k2 < hcalHitCache.size(); k2++) {
0748 hit = hcalHitCache[k2];
0749 if (hit.id() == id) {
0750 nhit++;
0751 hits[nhit] = hit;
0752 esim += hit.e();
0753 }
0754 }
0755 k1 += nhit;
0756 nhit++;
0757 std::vector<int> cd = myQie->getCode(nhit, hits, engine);
0758 double eq = myQie->getEnergy(cd);
0759 #ifdef EDM_ML_DEBUG
0760 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id << std::dec << " registers " << esim
0761 << " energy from " << nhit << " hits starting with hit # " << k1
0762 << " energy with noise " << eq;
0763 #endif
0764 for (int k2 = 0; k2 < nTower; k2++) {
0765 if (id == idTower[k2]) {
0766 todo[k2] = 1;
0767 esimh[k2] = esim;
0768 eqie[k2] = eq;
0769 }
0770 }
0771 }
0772
0773
0774 for (int k2 = 0; k2 < nTower; k2++) {
0775 if (todo[k2] == 0) {
0776 std::vector<int> cd = myQie->getCode(0, hits, engine);
0777 double eq = myQie->getEnergy(cd);
0778 esimh[k2] = 0;
0779 eqie[k2] = eq;
0780 #ifdef EDM_ML_DEBUG
0781 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << idTower[k2] << std::dec
0782 << " registers " << esimh[k2] << " energy from hits and energy after QIE analysis "
0783 << eqie[k2];
0784 #endif
0785 }
0786 }
0787 }
0788
0789 void HcalTB04Analysis::xtalAnalysis(CLHEP::HepRandomEngine* engine) {
0790 CLHEP::RandGaussQ randGauss(*engine);
0791
0792
0793 std::vector<int> iok(nCrystal, 0);
0794 #ifdef EDM_ML_DEBUG
0795 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " << iok.size() << " " << idEcal.size() << " "
0796 << esime.size() << " " << enois.size();
0797 #endif
0798 for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
0799 uint32_t id = ecalHitCache[k1].id();
0800 int nhit = 0;
0801 double esim = ecalHitCache[k1].e();
0802 for (unsigned int k2 = k1 + 1; k2 < ecalHitCache.size(); k2++) {
0803 if (ecalHitCache[k2].id() == id) {
0804 nhit++;
0805 esim += ecalHitCache[k2].e();
0806 }
0807 }
0808 k1 += nhit;
0809 nhit++;
0810 double eq = esim + randGauss.fire(0., ecalNoise);
0811 #ifdef EDM_ML_DEBUG
0812 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id << std::dec << " registers " << esim
0813 << " energy from " << nhit << " hits starting with hit # " << k1
0814 << " energy with noise " << eq;
0815 #endif
0816 for (int k2 = 0; k2 < nCrystal; k2++) {
0817 if (id == idEcal[k2]) {
0818 iok[k2] = 1;
0819 esime[k2] = esim;
0820 enois[k2] = eq;
0821 }
0822 }
0823 }
0824
0825
0826 for (int k2 = 0; k2 < nCrystal; k2++) {
0827 if (iok[k2] == 0) {
0828 esime[k2] = 0;
0829 enois[k2] = randGauss.fire(0., ecalNoise);
0830 #ifdef EDM_ML_DEBUG
0831 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << idEcal[k2] << std::dec
0832 << " registers " << esime[k2] << " energy from hits and energy from noise "
0833 << enois[k2];
0834 #endif
0835 }
0836 }
0837 }
0838
0839 void HcalTB04Analysis::finalAnalysis() {
0840
0841 histo->fillPrimary(pInit, etaInit, phiInit);
0842
0843
0844 eecals = ehcals = eecalq = ehcalq = 0.;
0845 for (int i = 0; i < nTower; i++) {
0846 ehcals += esimh[i];
0847 ehcalq += eqie[i];
0848 }
0849 for (int i = 0; i < nCrystal; i++) {
0850 eecals += esime[i];
0851 eecalq += enois[i];
0852 }
0853 etots = eecals + ehcals;
0854 etotq = eecalq + ehcalq;
0855 #ifdef EDM_ML_DEBUG
0856 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level (Total) " << etots << " (ECal) "
0857 << eecals << " (HCal) " << ehcals
0858 << "\nHcalTB04Analysis:: Energy deposit at Qie Level (Total) " << etotq << " (ECal) "
0859 << eecalq << " (HCal) " << ehcalq;
0860 #endif
0861 histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0862
0863
0864 for (int i = 0; i < 5; i++) {
0865 eseta[i] = 0.;
0866 eqeta[i] = 0.;
0867 }
0868 for (int i = 0; i < 3; i++) {
0869 esphi[i] = 0.;
0870 eqphi[i] = 0.;
0871 }
0872 double e1 = 0, e2 = 0;
0873 unsigned int id;
0874 for (int i = 0; i < nTower; i++) {
0875 int det, z, group, ieta, iphi, layer;
0876 id = idTower[i];
0877 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0878 iphi -= (icphi - 1);
0879 if (icphi > 4) {
0880 if (ieta == 0)
0881 ieta = 2;
0882 else
0883 ieta = -1;
0884 } else {
0885 ieta = ieta - iceta + 2;
0886 }
0887 if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
0888 eseta[ieta] += esimh[i];
0889 esphi[iphi] += esimh[i];
0890 e1 += esimh[i];
0891 eqeta[ieta] += eqie[i];
0892 eqphi[iphi] += eqie[i];
0893 e2 += eqie[i];
0894 }
0895 }
0896 for (int i = 0; i < 3; i++) {
0897 if (e1 > 0)
0898 esphi[i] /= e1;
0899 if (e2 > 0)
0900 eqphi[i] /= e2;
0901 }
0902 for (int i = 0; i < 5; i++) {
0903 if (e1 > 0)
0904 eseta[i] /= e1;
0905 if (e2 > 0)
0906 eqeta[i] /= e2;
0907 }
0908 #ifdef EDM_ML_DEBUG
0909 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and Phi (Sim/Qie)";
0910 for (int i = 0; i < 5; i++)
0911 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = " << eseta[i] << " Qie = " << eqeta[i]
0912 << " Phi Sim = " << esphi[i] << " Qie = " << eqphi[i];
0913 #endif
0914 histo->fillTrnsProf(eseta, eqeta, esphi, eqphi);
0915
0916
0917 for (int i = 0; i < 20; i++) {
0918 eslay[i] = 0.;
0919 eqlay[i] = 0.;
0920 }
0921 e1 = 0;
0922 e2 = 0;
0923 for (int i = 0; i < nTower; i++) {
0924 int det, z, group, ieta, iphi, layer;
0925 id = idTower[i];
0926 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0927 iphi -= (icphi - 1);
0928 layer -= 1;
0929 if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
0930 eslay[layer] += esimh[i];
0931 e1 += esimh[i];
0932 eqlay[layer] += eqie[i];
0933 e2 += eqie[i];
0934 }
0935 }
0936 for (int i = 0; i < 20; i++) {
0937 if (e1 > 0)
0938 eslay[i] /= e1;
0939 if (e2 > 0)
0940 eqlay[i] /= e2;
0941 }
0942 #ifdef EDM_ML_DEBUG
0943 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
0944 for (int i = 0; i < 20; i++)
0945 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = " << eslay[i] << " Qie = " << eqlay[i];
0946 #endif
0947 histo->fillLongProf(eslay, eqlay);
0948 }
0949
0950 void HcalTB04Analysis::fillEvent(PHcalTB04Info& product) {
0951
0952 product.setIDs(idHcal, idXtal);
0953
0954
0955 product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
0956
0957
0958 product.setEdepHcal(esimh, eqie);
0959 product.setEdepHcal(esime, enois);
0960
0961
0962 product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0963
0964
0965 product.setTrnsProf(eseta, eqeta, esphi, eqphi);
0966
0967
0968 product.setLongProf(eslay, eqlay);
0969
0970
0971 int i, nhit = 0;
0972 std::vector<CaloHit>::iterator itr;
0973 for (i = 0, itr = ecalHitCache.begin(); itr != ecalHitCache.end(); i++, itr++) {
0974 uint32_t id = itr->id();
0975 int det, z, group, ieta, iphi, lay;
0976 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
0977 product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
0978 nhit++;
0979 #ifdef EDM_ML_DEBUG
0980 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
0981 << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
0982 << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
0983 << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
0984 #endif
0985 }
0986 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from Crystals";
0987 int hit = nhit;
0988 nhit = 0;
0989
0990 for (i = hit, itr = hcalHitCache.begin(); itr != hcalHitCache.end(); i++, itr++) {
0991 uint32_t id = itr->id();
0992 int det, z, group, ieta, iphi, lay;
0993 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
0994 product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
0995 nhit++;
0996 #ifdef EDM_ML_DEBUG
0997 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
0998 << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
0999 << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
1000 << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
1001 #endif
1002 }
1003 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from HCal";
1004
1005
1006 product.setVtxPrim(evNum,
1007 pvType,
1008 pvPosition.x(),
1009 pvPosition.y(),
1010 pvPosition.z(),
1011 pvUVW.x(),
1012 pvUVW.y(),
1013 pvUVW.z(),
1014 pvMomentum.x(),
1015 pvMomentum.y(),
1016 pvMomentum.z());
1017 for (unsigned int i = 0; i < secTrackID.size(); i++) {
1018 product.setVtxSec(
1019 secTrackID[i], secPartID[i], secMomentum[i].x(), secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
1020 }
1021 }
1022
1023 void HcalTB04Analysis::clear() {
1024 pvFound = false;
1025 pvType = -2;
1026 pvPosition = G4ThreeVector();
1027 pvMomentum = G4ThreeVector();
1028 pvUVW = G4ThreeVector();
1029 secTrackID.clear();
1030 secPartID.clear();
1031 secMomentum.clear();
1032 secEkin.clear();
1033 shortLivedSecondaries.clear();
1034
1035 ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
1036 hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
1037 hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
1038 nPrimary = particleType = 0;
1039 pInit = etaInit = phiInit = 0;
1040
1041 esimh.clear();
1042 eqie.clear();
1043 esimh.reserve(nTower);
1044 eqie.reserve(nTower);
1045 for (int i = 0; i < nTower; i++) {
1046 esimh.push_back(0.);
1047 eqie.push_back(0.);
1048 }
1049 esime.clear();
1050 enois.clear();
1051 esime.reserve(nCrystal);
1052 enois.reserve(nCrystal);
1053 for (int i = 0; i < nCrystal; i++) {
1054 esime.push_back(0.);
1055 enois.push_back(0.);
1056 }
1057 }
1058
1059 int HcalTB04Analysis::unitID(uint32_t id) {
1060 int det, z, group, ieta, iphi, lay;
1061 HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1062 group = (det & 15) << 20;
1063 group += ((lay - 1) & 31) << 15;
1064 group += (z & 1) << 14;
1065 group += (ieta & 127) << 7;
1066 group += (iphi & 127);
1067 return group;
1068 }
1069
1070 double HcalTB04Analysis::scale(int det, int layer) {
1071 double tmp = 1.;
1072 if (det == static_cast<int>(HcalBarrel)) {
1073 if (layer == 1)
1074 tmp = scaleHB0;
1075 else if (layer == 17)
1076 tmp = scaleHB16;
1077 else if (layer > 17)
1078 tmp = scaleHO;
1079 } else {
1080 if (layer <= 2)
1081 tmp = scaleHE0;
1082 }
1083 return tmp;
1084 }
1085
1086 double HcalTB04Analysis::timeOfFlight(int det, int layer, double eta) {
1087 double theta = 2.0 * std::atan(std::exp(-eta));
1088 double dist = beamOffset;
1089 if (det == static_cast<int>(HcalBarrel)) {
1090 const double rLay[19] = {1836.0,
1091 1902.0,
1092 1962.0,
1093 2022.0,
1094 2082.0,
1095 2142.0,
1096 2202.0,
1097 2262.0,
1098 2322.0,
1099 2382.0,
1100 2448.0,
1101 2514.0,
1102 2580.0,
1103 2646.0,
1104 2712.0,
1105 2776.0,
1106 2862.5,
1107 3847.0,
1108 4052.0};
1109 if (layer > 0 && layer <= 19)
1110 dist += rLay[layer - 1] * mm / sin(theta);
1111 } else {
1112 const double zLay[19] = {4034.0,
1113 4032.0,
1114 4123.0,
1115 4210.0,
1116 4297.0,
1117 4384.0,
1118 4471.0,
1119 4558.0,
1120 4645.0,
1121 4732.0,
1122 4819.0,
1123 4906.0,
1124 4993.0,
1125 5080.0,
1126 5167.0,
1127 5254.0,
1128 5341.0,
1129 5428.0,
1130 5515.0};
1131 if (layer > 0 && layer <= 19)
1132 dist += zLay[layer - 1] * mm / cos(theta);
1133 }
1134
1135 double tmp = dist / c_light / ns;
1136 #ifdef EDM_ML_DEBUG
1137 edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
1138 << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
1139 #endif
1140 return tmp;
1141 }
1142
1143 DEFINE_SIMWATCHER(HcalTB04Analysis);