File indexing completed on 2024-05-10 02:21:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <cmath>
0015 #include <iostream>
0016 #include <iomanip>
0017 #include <memory>
0018 #include <map>
0019 #include <vector>
0020 #include <string>
0021
0022
0023
0024 #include "SimDataFormats/HcalTestBeam/interface/HcalTB02HistoClass.h"
0025 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0026 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0027 #include "SimG4Core/Notification/interface/Observer.h"
0028 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0029 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0030 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0031 #include "SimG4Core/Watcher/interface/SimProducer.h"
0032 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0033
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/Framework/interface/MakerMacros.h"
0037 #include "FWCore/PluginManager/interface/ModuleDef.h"
0038
0039 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02HcalNumberingScheme.h"
0040 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02XtalNumberingScheme.h"
0041 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02Histo.h"
0042
0043 #include "G4HCofThisEvent.hh"
0044 #include "G4SDManager.hh"
0045 #include "G4Step.hh"
0046 #include "G4Track.hh"
0047 #include "G4ThreeVector.hh"
0048 #include "G4VProcess.hh"
0049
0050 #include <CLHEP/Random/RandGaussQ.h>
0051 #include <CLHEP/Random/Randomize.h>
0052 #include <CLHEP/Units/SystemOfUnits.h>
0053 #include <CLHEP/Units/PhysicalConstants.h>
0054
0055 using CLHEP::GeV;
0056 using CLHEP::twopi;
0057
0058
0059
0060 namespace CLHEP {
0061 class HepRandomEngine;
0062 }
0063
0064 class HcalTB02Analysis : public SimProducer, public Observer<const BeginOfEvent*>, public Observer<const EndOfEvent*> {
0065 public:
0066 HcalTB02Analysis(const edm::ParameterSet& p);
0067 HcalTB02Analysis(const HcalTB02Analysis&) = delete;
0068 const HcalTB02Analysis& operator=(const HcalTB02Analysis&) = delete;
0069 ~HcalTB02Analysis() override;
0070
0071 void produce(edm::Event&, const edm::EventSetup&) override;
0072
0073 private:
0074
0075 void update(const BeginOfEvent* evt) override;
0076 void update(const EndOfEvent* evt) override;
0077
0078 void fillEvent(HcalTB02HistoClass&);
0079 void clear();
0080 void finish();
0081
0082 private:
0083
0084 std::unique_ptr<HcalTB02Histo> histo;
0085
0086
0087 const edm::ParameterSet m_Anal;
0088 const bool hcalOnly;
0089 const std::vector<std::string> names;
0090
0091
0092 std::map<int, float> energyInScints, energyInCrystals;
0093 std::map<int, float> primaries;
0094 int particleType;
0095 double eta, phi, pInit, incidentEnergy;
0096 float SEnergy, E7x7Matrix, E5x5Matrix;
0097 float SEnergyN, E7x7MatrixN, E5x5MatrixN;
0098 int maxTime;
0099 double xIncidentEnergy;
0100 float xSEnergy, xSEnergyN;
0101 float xE3x3Matrix, xE5x5Matrix;
0102 float xE3x3MatrixN, xE5x5MatrixN;
0103 };
0104
0105
0106
0107
0108
0109 HcalTB02Analysis::HcalTB02Analysis(const edm::ParameterSet& p)
0110 : m_Anal(p.getParameter<edm::ParameterSet>("HcalTB02Analysis")),
0111 hcalOnly(m_Anal.getUntrackedParameter<bool>("HcalClusterOnly", true)),
0112 names(m_Anal.getParameter<std::vector<std::string> >("Names")) {
0113 produces<HcalTB02HistoClass>();
0114
0115 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis:: Initialised as observer of "
0116 << "BeginOfJob/BeginOfEvent/EndOfEvent with "
0117 << "Parameter values:\n \thcalOnly = " << hcalOnly;
0118
0119 histo = std::make_unique<HcalTB02Histo>(m_Anal);
0120 }
0121
0122 HcalTB02Analysis::~HcalTB02Analysis() { finish(); }
0123
0124
0125
0126
0127
0128 void HcalTB02Analysis::produce(edm::Event& e, const edm::EventSetup&) {
0129 std::unique_ptr<HcalTB02HistoClass> product(new HcalTB02HistoClass);
0130 fillEvent(*product);
0131 e.put(std::move(product));
0132 }
0133
0134 void HcalTB02Analysis::update(const BeginOfEvent* evt) {
0135 #ifdef EDM_ML_DEBUG
0136 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis: =====> Begin of event = " << (*evt)()->GetEventID();
0137 #endif
0138 clear();
0139 }
0140
0141 void HcalTB02Analysis::update(const EndOfEvent* evt) {
0142 CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
0143 CLHEP::RandGaussQ randGauss(*engine);
0144
0145
0146 #ifdef EDM_ML_DEBUG
0147 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis::Fill event " << (*evt)()->GetEventID();
0148 #endif
0149
0150 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0151 int ihit = 0;
0152
0153
0154 std::string sd = names[0];
0155 int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
0156 CaloG4HitCollection* theHCHC = (CaloG4HitCollection*)allHC->GetHC(HCHCid);
0157 #ifdef EDM_ML_DEBUG
0158 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: Hit Collection for " << sd << " of ID " << HCHCid
0159 << " is obtained at " << theHCHC;
0160 #endif
0161 int nentries = 0;
0162 nentries = theHCHC->entries();
0163 if (nentries == 0)
0164 return;
0165
0166 int xentries = 0;
0167 int XTALSid = 0;
0168 CaloG4HitCollection* theXTHC = nullptr;
0169
0170 if (!hcalOnly) {
0171
0172 sd = names[1];
0173 XTALSid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
0174 theXTHC = (CaloG4HitCollection*)allHC->GetHC(XTALSid);
0175 #ifdef EDM_ML_DEBUG
0176 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: Hit Collection for " << sd << " of ID " << XTALSid
0177 << " is obtained at " << theXTHC;
0178 #endif
0179 xentries = theXTHC->entries();
0180 if (xentries == 0)
0181 return;
0182 }
0183
0184 #ifdef EDM_ML_DEBUG
0185 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: There are " << nentries << " HCal hits, and" << xentries
0186 << " xtal hits";
0187 #endif
0188 float ETot = 0., xETot = 0.;
0189 int scintID = 0, xtalID = 0;
0190
0191
0192 std::unique_ptr<HcalTB02HcalNumberingScheme> org(new HcalTB02HcalNumberingScheme());
0193
0194 if (HCHCid >= 0 && theHCHC != nullptr) {
0195 for (ihit = 0; ihit < nentries; ihit++) {
0196 CaloG4Hit* aHit = (*theHCHC)[ihit];
0197 scintID = aHit->getUnitID();
0198 int layer = org->getlayerID(scintID);
0199 float enEm = aHit->getEM();
0200 float enhad = aHit->getHadr();
0201
0202 if (layer == 0) {
0203 enEm = enEm / 2.;
0204 enhad = enhad / 2.;
0205 }
0206
0207 energyInScints[scintID] += enEm + enhad;
0208 primaries[aHit->getTrackID()] += enEm + enhad;
0209 float time = aHit->getTimeSliceID();
0210 if (time > maxTime)
0211 maxTime = static_cast<int>(time);
0212 histo->fillAllTime(time);
0213 }
0214
0215
0216
0217
0218
0219 float TowerEne[8][18], TowerEneCF[8][18], LayerEne[19], EnRing[100];
0220 for (int iphi = 0; iphi < 8; iphi++) {
0221 for (int jeta = 0; jeta < 18; jeta++) {
0222 TowerEne[iphi][jeta] = 0.;
0223 TowerEneCF[iphi][jeta] = 0.;
0224 }
0225 }
0226
0227 for (int ilayer = 0; ilayer < 19; ilayer++)
0228 LayerEne[ilayer] = 0.;
0229 for (int iring = 0; iring < 100; iring++)
0230 EnRing[iring] = 0.;
0231
0232 for (std::map<int, float>::iterator is = energyInScints.begin(); is != energyInScints.end(); is++) {
0233 ETot = (*is).second;
0234
0235 int layer = org->getlayerID((*is).first);
0236
0237 if ((layer != 17) && (layer != 18)) {
0238 float eta = org->getetaID((*is).first);
0239 float phi = org->getphiID((*is).first);
0240
0241 SEnergy += ETot;
0242 int iphi = static_cast<int>(phi);
0243 int ieta = static_cast<int>(eta);
0244 TowerEne[iphi][ieta] += ETot;
0245
0246 TowerEneCF[iphi][ieta] += ETot * (1. + 0.1 * randGauss.fire());
0247 double dR = 0.08727 * std::sqrt((eta - 8.) * (eta - 8.) + (phi - 3.) * (phi - 3.));
0248 EnRing[static_cast<int>(dR / 0.01)] += ETot;
0249 }
0250
0251 LayerEne[layer] += ETot;
0252 }
0253
0254 for (int ilayer = 0; ilayer < 19; ilayer++) {
0255 histo->fillProfile(ilayer, LayerEne[ilayer] / GeV);
0256 }
0257
0258 for (int iring = 0; iring < 100; iring++)
0259 histo->fillTransProf(0.01 * iring + 0.005, EnRing[iring] / GeV);
0260
0261 for (int iphi = 0; iphi < 8; iphi++) {
0262 for (int jeta = 0; jeta < 18; jeta++) {
0263 SEnergyN += TowerEneCF[iphi][jeta] + 3. * randGauss.fire();
0264
0265 double Rand = 3. * randGauss.fire();
0266
0267 if ((iphi >= 0) && (iphi < 7)) {
0268 if ((jeta >= 5) && (jeta < 12)) {
0269 E7x7Matrix += TowerEne[iphi][jeta];
0270 E7x7MatrixN += TowerEneCF[iphi][jeta] + Rand;
0271
0272 if ((iphi >= 1) && (iphi < 6)) {
0273 if ((jeta >= 6) && (jeta < 11)) {
0274 E5x5Matrix += TowerEne[iphi][jeta];
0275 E5x5MatrixN += TowerEneCF[iphi][jeta] + Rand;
0276 }
0277 }
0278 }
0279 }
0280 }
0281 }
0282
0283
0284
0285
0286 int trackID = 0;
0287 G4PrimaryParticle* thePrim = nullptr;
0288 G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0289 #ifdef EDM_ML_DEBUG
0290 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: Event has " << nvertex << " vertex";
0291 #endif
0292 if (nvertex == 0)
0293 edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
0294 << "ERROR: no vertex";
0295
0296 for (int i = 0; i < nvertex; i++) {
0297 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0298 if (avertex == nullptr) {
0299 edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
0300 << "ERROR: pointer to vertex = 0";
0301 } else {
0302 #ifdef EDM_ML_DEBUG
0303 int npart = avertex->GetNumberOfParticle();
0304 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis::Vertex number :" << i << " with " << npart << " particles";
0305 #endif
0306 if (thePrim == nullptr)
0307 thePrim = avertex->GetPrimary(trackID);
0308 }
0309 }
0310
0311 double px = 0., py = 0., pz = 0.;
0312
0313 if (thePrim != nullptr) {
0314 px = thePrim->GetPx();
0315 py = thePrim->GetPy();
0316 pz = thePrim->GetPz();
0317 pInit = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0318 if (pInit == 0) {
0319 edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
0320 << " ERROR: primary has p=0 ";
0321 } else {
0322 float costheta = pz / pInit;
0323 float theta = acos(std::min(std::max(costheta, float(-1.)), float(1.)));
0324 eta = -std::log(std::tan(theta / 2));
0325 if (px != 0)
0326 phi = std::atan(py / px);
0327 }
0328 particleType = thePrim->GetPDGcode();
0329 } else {
0330 #ifdef EDM_ML_DEBUG
0331 edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis:: End Of Event ERROR: could not find primary ";
0332 #endif
0333 }
0334
0335 CaloG4Hit* firstHit = (*theHCHC)[0];
0336 incidentEnergy = (firstHit->getIncidentEnergy() / GeV);
0337
0338 }
0339
0340 if (!hcalOnly) {
0341
0342
0343 if (XTALSid >= 0 && theXTHC != nullptr) {
0344 for (int xihit = 0; xihit < xentries; xihit++) {
0345 CaloG4Hit* xaHit = (*theXTHC)[xihit];
0346
0347 float xenEm = xaHit->getEM();
0348 float xenhad = xaHit->getHadr();
0349 xtalID = xaHit->getUnitID();
0350
0351 energyInCrystals[xtalID] += xenEm + xenhad;
0352 }
0353
0354 float xCrysEne[7][7];
0355 for (int irow = 0; irow < 7; irow++) {
0356 for (int jcol = 0; jcol < 7; jcol++) {
0357 xCrysEne[irow][jcol] = 0.;
0358 }
0359 }
0360
0361 for (std::map<int, float>::iterator is = energyInCrystals.begin(); is != energyInCrystals.end(); is++) {
0362 int xtalID = (*is).first;
0363 xETot = (*is).second;
0364
0365 int irow = static_cast<int>(xtalID / 100.);
0366 int jcol = static_cast<int>(xtalID - 100. * irow);
0367
0368 xSEnergy += xETot;
0369 xCrysEne[irow][jcol] = xETot;
0370
0371 float dR = std::sqrt(0.01619 * 0.01619 * (jcol - 3) * (jcol - 3) + 0.01606 * 0.01606 * (irow - 3) * (irow - 3));
0372 histo->fillTransProf(dR, xETot * 1.05);
0373
0374 if ((irow > 0) && (irow < 6)) {
0375 if ((jcol > 0) && (jcol < 6)) {
0376 xE5x5Matrix += xCrysEne[irow][jcol];
0377 xE5x5MatrixN += xCrysEne[irow][jcol] + 108.5 * randGauss.fire();
0378
0379 if ((irow > 1) && (irow < 5)) {
0380 if ((jcol > 1) && (jcol < 5)) {
0381 xE3x3Matrix += xCrysEne[irow][jcol];
0382 xE3x3MatrixN += xCrysEne[irow][jcol] + 108.5 * randGauss.fire();
0383 }
0384 }
0385 }
0386 }
0387 }
0388
0389 if (!hcalOnly) {
0390
0391 if (theXTHC != nullptr) {
0392 CaloG4Hit* xfirstHit = (*theXTHC)[0];
0393 xIncidentEnergy = xfirstHit->getIncidentEnergy() / GeV;
0394 }
0395 }
0396
0397 }
0398 }
0399
0400 int iEvt = (*evt)()->GetEventID();
0401 if (iEvt < 10)
0402 edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0403 else if ((iEvt < 100) && (iEvt % 10 == 0))
0404 edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0405 else if ((iEvt < 1000) && (iEvt % 100 == 0))
0406 edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0407 else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0408 edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0409 }
0410
0411 void HcalTB02Analysis::fillEvent(HcalTB02HistoClass& product) {
0412
0413 product.set_Nprim(float(primaries.size()));
0414 product.set_partType(particleType);
0415 product.set_Einit(pInit / GeV);
0416 product.set_eta(eta);
0417 product.set_phi(phi);
0418 product.set_Eentry(incidentEnergy);
0419
0420
0421 product.set_ETot(SEnergy / GeV);
0422 product.set_E7x7(E7x7Matrix / GeV);
0423 product.set_E5x5(E5x5Matrix / GeV);
0424 product.set_ETotN(SEnergyN / GeV);
0425 product.set_E7x7N(E7x7MatrixN / GeV);
0426 product.set_E5x5N(E5x5MatrixN / GeV);
0427 product.set_NUnit(float(energyInScints.size()));
0428 product.set_Ntimesli(float(maxTime));
0429
0430
0431 product.set_xEentry(xIncidentEnergy);
0432 product.set_xNUnit(float(energyInCrystals.size()));
0433 product.set_xETot(xSEnergy / GeV);
0434 product.set_xETotN(xSEnergyN / GeV);
0435 product.set_xE5x5(xE5x5Matrix / GeV);
0436 product.set_xE3x3(xE3x3Matrix / GeV);
0437 product.set_xE5x5N(xE5x5MatrixN / GeV);
0438 product.set_xE3x3N(xE3x3MatrixN / GeV);
0439 }
0440
0441 void HcalTB02Analysis::clear() {
0442 primaries.clear();
0443 particleType = 0;
0444 pInit = eta = phi = incidentEnergy = 0;
0445
0446 SEnergy = E7x7Matrix = E5x5Matrix = SEnergyN = 0;
0447 E7x7MatrixN = E5x5MatrixN = 0;
0448 energyInScints.clear();
0449 maxTime = 0;
0450
0451 xIncidentEnergy = 0;
0452 energyInCrystals.clear();
0453 xSEnergy = xSEnergyN = xE5x5Matrix = xE3x3Matrix = 0;
0454 xE5x5MatrixN = xE3x3MatrixN = 0;
0455 }
0456
0457 void HcalTB02Analysis::finish() {
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491 }
0492
0493 DEFINE_SIMWATCHER(HcalTB02Analysis);