File indexing completed on 2021-12-06 04:01:07
0001 #include "SimG4CMS/Calo/interface/HFGflash.h"
0002
0003 #include "G4VPhysicalVolume.hh"
0004 #include "G4Step.hh"
0005 #include "G4Track.hh"
0006 #include "G4Navigator.hh"
0007 #include "G4NavigationHistory.hh"
0008 #include "CLHEP/Units/PhysicalConstants.h"
0009 #include "CLHEP/Units/SystemOfUnits.h"
0010
0011 #include "Randomize.hh"
0012 #include "G4TransportationManager.hh"
0013 #include "G4VPhysicalVolume.hh"
0014 #include "G4LogicalVolume.hh"
0015 #include "G4VSensitiveDetector.hh"
0016 #include "G4EventManager.hh"
0017 #include "G4SteppingManager.hh"
0018 #include "G4FastTrack.hh"
0019 #include "G4ParticleTable.hh"
0020
0021 #include "CLHEP/GenericFunctions/IncompleteGamma.hh"
0022
0023 #include "SimG4Core/Application/interface/SteppingAction.h"
0024 #include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
0025 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
0026
0027 #include "FWCore/ServiceRegistry/interface/Service.h"
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029
0030 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0031 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0032
0033 #include <cmath>
0034
0035
0036
0037 HFGflash::HFGflash(edm::ParameterSet const& p) {
0038 edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFGflash");
0039 theBField = m_HF.getUntrackedParameter<double>("BField", 3.8);
0040 theWatcherOn = m_HF.getUntrackedParameter<bool>("WatcherOn", true);
0041 theFillHisto = m_HF.getUntrackedParameter<bool>("FillHisto", true);
0042 edm::LogVerbatim("HFShower") << "HFGFlash:: Set B-Field to " << theBField << ", WatcherOn to " << theWatcherOn
0043 << " and FillHisto to " << theFillHisto;
0044
0045 theHelix = std::make_unique<GflashTrajectory>();
0046 theGflashStep = std::make_unique<G4Step>();
0047 theGflashNavigator = std::make_unique<G4Navigator>();
0048 theGflashTouchableHandle = std::make_unique<G4TouchableHistory>();
0049
0050 #ifdef EDM_ML_DEBUG
0051 if (theFillHisto) {
0052 edm::Service<TFileService> tfile;
0053 if (tfile.isAvailable()) {
0054 TFileDirectory showerDir = tfile->mkdir("GflashEMShowerProfile");
0055
0056 em_incE = showerDir.make<TH1F>("em_incE", "Incoming energy (GeV)", 500, 0, 500.);
0057 em_ssp_rho =
0058 showerDir.make<TH1F>("em_ssp_rho", "Shower starting position;#rho (cm);Number of Events", 100, 100.0, 200.0);
0059 em_ssp_z =
0060 showerDir.make<TH1F>("em_ssp_z", "Shower starting position;z (cm);Number of Events", 2000, 0.0, 2000.0);
0061 em_long =
0062 showerDir.make<TH1F>("em_long", "Longitudinal Profile;Radiation Length;Number of Spots", 800, 800.0, 1600.0);
0063 em_lateral = showerDir.make<TH1F>("em_lateral", "Lateral Profile;Radiation Length;Moliere Radius", 100, 0.0, 5.0);
0064 em_2d = showerDir.make<TH2F>("em_2d",
0065 "Lateral Profile vs. Shower Depth;Radiation Length;Moliere Radius",
0066 800,
0067 800.0,
0068 1600.0,
0069 100,
0070 0.0,
0071 5.0);
0072 em_long_sd = showerDir.make<TH1F>("em_long_sd",
0073 "Longitudinal Profile in Sensitive Detector;Radiation Length;Number of Spots",
0074 800,
0075 800.0,
0076 1600.0);
0077 em_lateral_sd =
0078 showerDir.make<TH1F>("em_lateral_sd",
0079 "Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",
0080 100,
0081 0.0,
0082 5.0);
0083 em_2d_sd =
0084 showerDir.make<TH2F>("em_2d_sd",
0085 "Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",
0086 800,
0087 800.0,
0088 1600.0,
0089 100,
0090 0.0,
0091 5.0);
0092 em_ze = showerDir.make<TH2F>(
0093 "em_ze", "Profile vs. Energy;Radiation Length;Moliere Radius", 800, 800.0, 1600.0, 1000, 0.0, 1.0);
0094 em_ratio = showerDir.make<TH2F>(
0095 "em_ratio", "Profile vs. Energy;Radiation Length;Moliere Radius", 800, 800.0, 1600.0, 1000, 0.0, 100.0);
0096 em_ratio_selected = showerDir.make<TH2F>("em_ratio_selected",
0097 "Profile vs. Energy;Radiation Length;Moliere Radius",
0098 800,
0099 800.0,
0100 1600.0,
0101 1000,
0102 0.0,
0103 100.0);
0104 em_nSpots_sd =
0105 showerDir.make<TH1F>("em_nSpots_sd",
0106 "Number of Gflash Spots in Sensitive Detector;Number of Spots;Number of Events",
0107 100,
0108 0.0,
0109 100);
0110 em_ze_ratio = showerDir.make<TH1F>("em_ze_ratio", "Ratio of Energy and Z Position", 1000, 0.0, 0.001);
0111 } else {
0112 theFillHisto = false;
0113 edm::LogVerbatim("HFShower") << "HFGFlash::No file is available for saving"
0114 << " histos so the flag is set to false";
0115 }
0116 }
0117 #endif
0118 jCalorimeter = Gflash::kHF;
0119 }
0120
0121 HFGflash::~HFGflash() {}
0122
0123 std::vector<HFGflash::Hit> HFGflash::gfParameterization(const G4Step* aStep, bool onlyLong) {
0124 double tempZCalo = 26;
0125 double hfcriticalEnergy = 0.021;
0126
0127 std::vector<HFGflash::Hit> hit;
0128 HFGflash::Hit oneHit;
0129
0130 auto const preStepPoint = aStep->GetPreStepPoint();
0131 auto const track = aStep->GetTrack();
0132
0133
0134
0135
0136 const G4double energyCutoff = 1;
0137 const G4int maxNumberOfSpots = 10000000;
0138
0139 G4ThreeVector showerStartingPosition = track->GetPosition() / cm;
0140 G4ThreeVector showerMomentum = preStepPoint->GetMomentum() / GeV;
0141 jCalorimeter = Gflash::kHF;
0142
0143 const G4double invgev = 1.0 / CLHEP::GeV;
0144 G4double energy = preStepPoint->GetTotalEnergy() * invgev;
0145 G4double logEinc = std::log(energy);
0146
0147 G4double y = energy / hfcriticalEnergy;
0148 G4double logY = std::log(y);
0149
0150 G4double nSpots = 93.0 * std::log(tempZCalo) * energy;
0151 if (energy < 1.6)
0152 nSpots = 140.4 * std::log(tempZCalo) * std::pow(energy, 0.876);
0153
0154
0155 double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
0156 theHelix->initializeTrajectory(showerMomentum, showerStartingPosition, charge, theBField);
0157
0158 G4double pathLength0 = theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
0159 G4double pathLength = pathLength0;
0160
0161
0162
0163 G4double fluctuatedTmax = std::log(logY - 0.7157);
0164 G4double fluctuatedAlpha = std::log(0.7996 + (0.4581 + 1.8628 / tempZCalo) * logY);
0165
0166 G4double sigmaTmax = 1.0 / (-1.4 + 1.26 * logY);
0167 G4double sigmaAlpha = 1.0 / (-0.58 + 0.86 * logY);
0168 G4double rho = 0.705 - 0.023 * logY;
0169 G4double sqrtPL = std::sqrt((1.0 + rho) / 2.0);
0170 G4double sqrtLE = std::sqrt((1.0 - rho) / 2.0);
0171
0172 G4double norm1 = G4RandGauss::shoot();
0173 G4double norm2 = G4RandGauss::shoot();
0174 G4double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
0175 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);
0176
0177
0178 G4double tmax = std::exp(tempTmax);
0179 G4double alpha = std::exp(tempAlpha);
0180 G4double beta = (alpha - 1.0) / tmax;
0181
0182 if (!alpha)
0183 return hit;
0184 if (!beta)
0185 return hit;
0186 if (alpha < 0.00001)
0187 return hit;
0188 if (beta < 0.00001)
0189 return hit;
0190
0191
0192 G4double averageTmax = logY - 0.858;
0193 G4double averageAlpha = 0.21 + (0.492 + 2.38 / tempZCalo) * logY;
0194 G4double spotTmax = averageTmax * (0.698 + .00212 * tempZCalo);
0195 G4double spotAlpha = averageAlpha * (0.639 + .00334 * tempZCalo);
0196 G4double spotBeta = (spotAlpha - 1.0) / spotTmax;
0197
0198 if (!spotAlpha)
0199 return hit;
0200 if (!spotBeta)
0201 return hit;
0202 if (spotAlpha < 0.00001)
0203 return hit;
0204 if (spotBeta < 0.00001)
0205 return hit;
0206
0207 #ifdef EDM_ML_DEBUG
0208 edm::LogVerbatim("HFShower") << "HFGflash: Incoming energy = " << energy << " Position (rho,z) = ("
0209 << showerStartingPosition.rho() << ", " << showerStartingPosition.z() << ")";
0210
0211 if (theFillHisto) {
0212 em_incE->Fill(energy);
0213 em_ssp_rho->Fill(showerStartingPosition.rho());
0214 em_ssp_z->Fill(std::abs(showerStartingPosition.z()));
0215 }
0216 #endif
0217
0218 G4double z1 = 0.0251 + 0.00319 * logEinc;
0219 G4double z2 = 0.1162 - 0.000381 * tempZCalo;
0220
0221 G4double k1 = 0.659 - 0.00309 * tempZCalo;
0222 G4double k2 = 0.645;
0223 G4double k3 = -2.59;
0224 G4double k4 = 0.3585 + 0.0421 * logEinc;
0225
0226 G4double p1 = 2.623 - 0.00094 * tempZCalo;
0227 G4double p2 = 0.401 + 0.00187 * tempZCalo;
0228 G4double p3 = 1.313 - 0.0686 * logEinc;
0229
0230
0231
0232 const G4double e25Scale = 1.03551;
0233 z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790e+00 * std::log(energy) - 3.66237e+00);
0234 p1 *= 0.96;
0235
0236 G4double stepLengthLeft = 10000;
0237 G4int nSpots_sd = 0;
0238 G4double zInX0 = 0.0;
0239 G4double deltaZInX0 = 0.0;
0240 G4double deltaZ = 0.0;
0241 G4double stepLengthLeftInX0 = 0.0;
0242
0243 const G4double divisionStepInX0 = 0.1;
0244
0245 Genfun::IncompleteGamma gammaDist;
0246
0247 G4double energyInGamma = 0.0;
0248 G4double preEnergyInGamma = 0.0;
0249 G4double sigmaInGamma = 0.;
0250 G4double preSigmaInGamma = 0.0;
0251
0252
0253 G4double deltaEnergy = 0.0;
0254 G4int spotCounter = 0;
0255
0256
0257 G4double deltaStep = 0.0;
0258
0259
0260 G4double timeGlobal = preStepPoint->GetGlobalTime();
0261
0262 theGflashNavigator->SetWorldVolume(
0263 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0264
0265
0266
0267 #ifdef EDM_ML_DEBUG
0268 edm::LogVerbatim("HFShower") << "HFGflash: Energy = " << energy << " Step Length Left = " << stepLengthLeft;
0269 #endif
0270 while (energy > 0.0 && stepLengthLeft > 0.0) {
0271 stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
0272
0273 if (stepLengthLeftInX0 < divisionStepInX0) {
0274 deltaZInX0 = stepLengthLeftInX0;
0275 deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0276 stepLengthLeft = 0.0;
0277 } else {
0278 deltaZInX0 = divisionStepInX0;
0279 deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0280 stepLengthLeft -= deltaZ;
0281 }
0282
0283 zInX0 += deltaZInX0;
0284
0285 #ifdef EDM_ML_DEBUG
0286 edm::LogVerbatim("HFShower") << "HFGflash: zInX0 = " << zInX0 << " spotBeta*zInX0 = " << spotBeta * zInX0;
0287 #endif
0288 if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
0289 (!(zInX0 * beta != 0)) || (zInX0 * beta < 0.00001))
0290 return hit;
0291
0292 G4int nSpotsInStep = 0;
0293
0294 #ifdef EDM_ML_DEBUG
0295 edm::LogVerbatim("HFShower") << "HFGflash: Energy - Energy Cut off = " << energy - energyCutoff;
0296 #endif
0297
0298 if (energy > energyCutoff) {
0299 preEnergyInGamma = energyInGamma;
0300 gammaDist.a().setValue(alpha);
0301
0302 energyInGamma = gammaDist(beta * zInX0);
0303 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
0304 deltaEnergy = std::min(energy, energy * energyInDeltaZ);
0305
0306 preSigmaInGamma = sigmaInGamma;
0307 gammaDist.a().setValue(spotAlpha);
0308 sigmaInGamma = gammaDist(spotBeta * zInX0);
0309 nSpotsInStep = std::max(1, int(nSpots * (sigmaInGamma - preSigmaInGamma)));
0310 } else {
0311 deltaEnergy = energy;
0312 preSigmaInGamma = sigmaInGamma;
0313 nSpotsInStep = std::max(1, int(nSpots * (1.0 - preSigmaInGamma)));
0314 }
0315
0316 if (deltaEnergy > energy || (energy - deltaEnergy) < energyCutoff)
0317 deltaEnergy = energy;
0318
0319 energy -= deltaEnergy;
0320
0321 if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
0322 nSpotsInStep = maxNumberOfSpots - spotCounter;
0323 if (nSpotsInStep < 1)
0324 nSpotsInStep = 1;
0325 }
0326
0327
0328 deltaStep += 0.5 * deltaZ;
0329 pathLength += deltaStep;
0330 deltaStep = 0.5 * deltaZ;
0331
0332
0333 G4double tScale = tmax * alpha / (alpha - 1.0) * (1.0 - std::exp(-fluctuatedAlpha));
0334 G4double tau = std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
0335 G4double rCore = z1 + z2 * tau;
0336 G4double rTail = k1 * (std::exp(k3 * (tau - k2)) + std::exp(k4 * (tau - k2)));
0337 G4double p23 = (p2 - tau) / p3;
0338 G4double probabilityWeight = p1 * std::exp(p23 - std::exp(p23));
0339
0340
0341
0342
0343 G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);
0344
0345 #ifdef EDM_ML_DEBUG
0346 edm::LogVerbatim("HFShower") << "HFGflash: nSpotsInStep = " << nSpotsInStep;
0347 #endif
0348
0349 for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
0350 spotCounter++;
0351 G4double u1 = G4UniformRand();
0352 G4double u2 = G4UniformRand();
0353 G4double rInRM = 0.0;
0354
0355 if (u1 < probabilityWeight) {
0356 rInRM = rCore * std::sqrt(u2 / (1.0 - u2));
0357 } else {
0358 rInRM = rTail * std::sqrt(u2 / (1.0 - u2));
0359 }
0360
0361 G4double rShower = rInRM * Gflash::rMoliere[jCalorimeter];
0362
0363
0364 G4double azimuthalAngle = twopi * G4UniformRand();
0365
0366
0367
0368 G4double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
0369
0370
0371 GflashTrajectoryPoint trajectoryPoint;
0372 theHelix->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
0373
0374 G4ThreeVector SpotPosition0 = trajectoryPoint.getPosition() +
0375 rShower * std::cos(azimuthalAngle) * trajectoryPoint.getOrthogonalUnitVector() +
0376 rShower * std::sin(azimuthalAngle) * trajectoryPoint.getCrossUnitVector();
0377
0378
0379
0380 SpotPosition0 *= cm;
0381
0382
0383
0384
0385
0386
0387 timeGlobal += 0.0001 * nanosecond;
0388
0389
0390
0391 G4double zInX0_spot = std::abs(pathLength + incrementPath - pathLength0) / Gflash::radLength[jCalorimeter];
0392
0393 #ifdef EDM_ML_DEBUG
0394 edm::LogVerbatim("HFShower") << "HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot << " , "
0395 << emSpotEnergy / GeV << "emSpotEnergy/GeV =" << emSpotEnergy / GeV;
0396 #endif
0397
0398 if (zInX0_spot <= 0)
0399 continue;
0400 if (emSpotEnergy <= 0)
0401 continue;
0402 if (rShower / Gflash::rMoliere[jCalorimeter] <= 0)
0403 continue;
0404
0405 oneHit.depth = 1;
0406
0407 #ifdef EDM_ML_DEBUG
0408 if (theFillHisto) {
0409 em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy * invgev);
0410 em_lateral->Fill(rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
0411 em_2d->Fill(SpotPosition0.z() / cm, rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
0412 }
0413 #endif
0414
0415
0416
0417
0418 double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
0419
0420 if (emSpotEnergy / GeV < 0.0001)
0421 continue;
0422 if (energyratio > 80)
0423 continue;
0424
0425 double zshift = 0;
0426 if (SpotPosition0.z() > 0)
0427 zshift = 18;
0428 if (SpotPosition0.z() < 0)
0429 zshift = -18;
0430
0431 G4ThreeVector gfshift(0, 0, zshift * (pow(100, 0.1) / pow(energy, 0.1)));
0432
0433 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
0434
0435 double LengthWeight = std::fabs(std::pow(SpotPosition0.z() / 11370, 1));
0436 if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
0437 continue;
0438
0439 oneHit.position = SpotPosition;
0440 oneHit.time = timeGlobal;
0441 oneHit.edep = emSpotEnergy * invgev;
0442 hit.push_back(oneHit);
0443 nSpots_sd++;
0444
0445 }
0446
0447 }
0448 #ifdef EDM_ML_DEBUG
0449 if (theFillHisto) {
0450 em_nSpots_sd->Fill(nSpots_sd);
0451 }
0452 #endif
0453 return hit;
0454 }