File indexing completed on 2024-04-06 12:11:16
0001
0002
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006
0007
0008
0009 #include "FastSimulation/Event/interface/FSimEvent.h"
0010 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0011 #include "FastSimulation/TrackerSetup/interface/TrackerLayer.h"
0012 #include "FastSimulation/MaterialEffects/interface/MaterialEffects.h"
0013
0014 #include "FastSimulation/MaterialEffects/interface/PairProductionSimulator.h"
0015 #include "FastSimulation/MaterialEffects/interface/MultipleScatteringSimulator.h"
0016 #include "FastSimulation/MaterialEffects/interface/BremsstrahlungSimulator.h"
0017 #include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
0018 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionSimulator.h"
0019 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionFTFSimulator.h"
0020 #include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
0021
0022 #include <list>
0023 #include <map>
0024 #include <string>
0025
0026 MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)
0027 : PairProduction(nullptr),
0028 Bremsstrahlung(nullptr),
0029 MuonBremsstrahlung(nullptr),
0030 MultipleScattering(nullptr),
0031 EnergyLoss(nullptr),
0032 NuclearInteraction(nullptr),
0033 pTmin(999.),
0034 use_hardcoded(true) {
0035
0036
0037 use_hardcoded = matEff.getParameter<bool>("use_hardcoded_geometry");
0038
0039 bool doPairProduction = matEff.getParameter<bool>("PairProduction");
0040 bool doBremsstrahlung = matEff.getParameter<bool>("Bremsstrahlung");
0041 bool doEnergyLoss = matEff.getParameter<bool>("EnergyLoss");
0042 bool doMultipleScattering = matEff.getParameter<bool>("MultipleScattering");
0043 bool doNuclearInteraction = matEff.getParameter<bool>("NuclearInteraction");
0044 bool doG4NuclInteraction = matEff.getParameter<bool>("G4NuclearInteraction");
0045 bool doMuonBremsstrahlung = matEff.getParameter<bool>("MuonBremsstrahlung");
0046
0047 double A = matEff.getParameter<double>("A");
0048 double Z = matEff.getParameter<double>("Z");
0049 double density = matEff.getParameter<double>("Density");
0050 double radLen = matEff.getParameter<double>("RadiationLength");
0051
0052
0053
0054 if (doPairProduction) {
0055 double photonEnergy = matEff.getParameter<double>("photonEnergy");
0056 PairProduction = new PairProductionSimulator(photonEnergy);
0057 }
0058
0059 if (doBremsstrahlung) {
0060 double bremEnergy = matEff.getParameter<double>("bremEnergy");
0061 double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
0062 Bremsstrahlung = new BremsstrahlungSimulator(bremEnergy, bremEnergyFraction);
0063 }
0064
0065 if (doMuonBremsstrahlung) {
0066 double bremEnergy = matEff.getParameter<double>("bremEnergy");
0067 double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
0068 MuonBremsstrahlung = new MuonBremsstrahlungSimulator(A, Z, density, radLen, bremEnergy, bremEnergyFraction);
0069 }
0070
0071
0072
0073 if (doEnergyLoss) {
0074 pTmin = matEff.getParameter<double>("pTmin");
0075 EnergyLoss = new EnergyLossSimulator(A, Z, density, radLen);
0076 }
0077
0078 if (doMultipleScattering) {
0079 MultipleScattering = new MultipleScatteringSimulator(A, Z, density, radLen);
0080 }
0081
0082 if (doNuclearInteraction) {
0083
0084 std::vector<double> hadronEnergies = matEff.getUntrackedParameter<std::vector<double> >("hadronEnergies");
0085
0086
0087 std::vector<int> hadronTypes = matEff.getUntrackedParameter<std::vector<int> >("hadronTypes");
0088
0089
0090 std::vector<std::string> hadronNames = matEff.getUntrackedParameter<std::vector<std::string> >("hadronNames");
0091
0092
0093 std::vector<double> hadronMasses = matEff.getUntrackedParameter<std::vector<double> >("hadronMasses");
0094
0095
0096 std::vector<double> hadronPMin = matEff.getUntrackedParameter<std::vector<double> >("hadronMinP");
0097
0098
0099 std::vector<double> lengthRatio = matEff.getParameter<std::vector<double> >("lengthRatio");
0100
0101
0102
0103
0104
0105
0106 theTECFudgeFactor = matEff.getParameter<double>("fudgeFactor");
0107
0108
0109 std::vector<double> theRatios = matEff.getUntrackedParameter<std::vector<double> >("ratios");
0110
0111
0112
0113
0114
0115
0116 std::vector<std::vector<double> > ratios;
0117 ratios.resize(hadronTypes.size());
0118 for (unsigned i = 0; i < hadronTypes.size(); ++i) {
0119 for (unsigned j = 0; j < hadronEnergies.size(); ++j) {
0120 ratios[i].push_back(theRatios[i * hadronEnergies.size() + j]);
0121 }
0122 }
0123
0124
0125 double pionEnergy = matEff.getParameter<double>("pionEnergy");
0126
0127
0128
0129 unsigned distAlgo = matEff.getParameter<unsigned>("distAlgo");
0130 double distCut = matEff.getParameter<double>("distCut");
0131
0132
0133
0134 std::string inputFile = matEff.getUntrackedParameter<std::string>("inputFile");
0135
0136
0137 std::map<int, int> idMap;
0138
0139 std::vector<int> idProtons = matEff.getUntrackedParameter<std::vector<int> >("protons");
0140 for (unsigned i = 0; i < idProtons.size(); ++i)
0141 idMap[idProtons[i]] = 2212;
0142
0143 std::vector<int> idAntiProtons = matEff.getUntrackedParameter<std::vector<int> >("antiprotons");
0144 for (unsigned i = 0; i < idAntiProtons.size(); ++i)
0145 idMap[idAntiProtons[i]] = -2212;
0146
0147 std::vector<int> idNeutrons = matEff.getUntrackedParameter<std::vector<int> >("neutrons");
0148 for (unsigned i = 0; i < idNeutrons.size(); ++i)
0149 idMap[idNeutrons[i]] = 2112;
0150
0151 std::vector<int> idAntiNeutrons = matEff.getUntrackedParameter<std::vector<int> >("antineutrons");
0152 for (unsigned i = 0; i < idAntiNeutrons.size(); ++i)
0153 idMap[idAntiNeutrons[i]] = -2112;
0154
0155 std::vector<int> idK0Ls = matEff.getUntrackedParameter<std::vector<int> >("K0Ls");
0156 for (unsigned i = 0; i < idK0Ls.size(); ++i)
0157 idMap[idK0Ls[i]] = 130;
0158
0159 std::vector<int> idKplusses = matEff.getUntrackedParameter<std::vector<int> >("Kplusses");
0160 for (unsigned i = 0; i < idKplusses.size(); ++i)
0161 idMap[idKplusses[i]] = 321;
0162
0163 std::vector<int> idKminusses = matEff.getUntrackedParameter<std::vector<int> >("Kminusses");
0164 for (unsigned i = 0; i < idKminusses.size(); ++i)
0165 idMap[idKminusses[i]] = -321;
0166
0167 std::vector<int> idPiplusses = matEff.getUntrackedParameter<std::vector<int> >("Piplusses");
0168 for (unsigned i = 0; i < idPiplusses.size(); ++i)
0169 idMap[idPiplusses[i]] = 211;
0170
0171 std::vector<int> idPiminusses = matEff.getUntrackedParameter<std::vector<int> >("Piminusses");
0172 for (unsigned i = 0; i < idPiminusses.size(); ++i)
0173 idMap[idPiminusses[i]] = -211;
0174
0175
0176 if (doG4NuclInteraction) {
0177 double elimit = matEff.getParameter<double>("EkinBertiniGeV") * CLHEP::GeV;
0178 double eth = matEff.getParameter<double>("EkinLimitGeV") * CLHEP::GeV;
0179 NuclearInteraction = new NuclearInteractionFTFSimulator(distAlgo, distCut, elimit, eth);
0180 } else {
0181 NuclearInteraction = new NuclearInteractionSimulator(hadronEnergies,
0182 hadronTypes,
0183 hadronNames,
0184 hadronMasses,
0185 hadronPMin,
0186 pionEnergy,
0187 lengthRatio,
0188 ratios,
0189 idMap,
0190 inputFile,
0191 distAlgo,
0192 distCut);
0193 }
0194 }
0195 }
0196
0197 MaterialEffects::~MaterialEffects() {
0198 if (PairProduction)
0199 delete PairProduction;
0200 if (Bremsstrahlung)
0201 delete Bremsstrahlung;
0202 if (EnergyLoss)
0203 delete EnergyLoss;
0204 if (MultipleScattering)
0205 delete MultipleScattering;
0206 if (NuclearInteraction)
0207 delete NuclearInteraction;
0208
0209 if (MuonBremsstrahlung)
0210 delete MuonBremsstrahlung;
0211 }
0212
0213 void MaterialEffects::interact(FSimEvent& mySimEvent,
0214 const TrackerLayer& layer,
0215 ParticlePropagator& myTrack,
0216 unsigned itrack,
0217 RandomEngineAndDistribution const* random) {
0218 MaterialEffectsSimulator::RHEP_const_iter DaughterIter;
0219 double radlen;
0220 theEnergyLoss = 0;
0221 theNormalVector = normalVector(layer, myTrack);
0222 radlen = radLengths(layer, myTrack);
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232 if (PairProduction && myTrack.particle().pid() == 22) {
0233
0234 PairProduction->updateState(myTrack, radlen, random);
0235
0236 if (PairProduction->nDaughters()) {
0237
0238 int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::PAIR_VERTEX);
0239
0240
0241 if (ivertex >= 0) {
0242
0243 for (DaughterIter = PairProduction->beginDaughters(); DaughterIter != PairProduction->endDaughters();
0244 ++DaughterIter) {
0245 mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0246 }
0247
0248 return;
0249 } else {
0250 edm::LogWarning("MaterialEffects")
0251 << " WARNING: A non valid vertex was found in photon conv. -> " << ivertex << std::endl;
0252 }
0253 }
0254 }
0255
0256 if (myTrack.particle().pid() == 22)
0257 return;
0258
0259
0260
0261
0262
0263 if (NuclearInteraction && abs(myTrack.particle().pid()) > 100 && abs(myTrack.particle().pid()) < 1000000) {
0264
0265 double factor = 1.0;
0266 if (use_hardcoded) {
0267 if (layer.layerNumber() >= 19 && layer.layerNumber() <= 27)
0268 factor = theTECFudgeFactor;
0269 }
0270 NuclearInteraction->updateState(myTrack, radlen * factor, random);
0271
0272
0273
0274 if (NuclearInteraction->nDaughters()) {
0275
0276 int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::NUCL_VERTEX);
0277
0278
0279
0280 if (ivertex >= 0) {
0281
0282 int idaugh = 0;
0283 for (DaughterIter = NuclearInteraction->beginDaughters(); DaughterIter != NuclearInteraction->endDaughters();
0284 ++DaughterIter) {
0285
0286 int daughId = mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0287
0288
0289 if (NuclearInteraction->closestDaughterId() == idaugh) {
0290 if (mySimEvent.track(itrack).vertex().position().Pt() < 4.0)
0291 mySimEvent.track(itrack).setClosestDaughterId(daughId);
0292 }
0293 ++idaugh;
0294 }
0295
0296 return;
0297 } else {
0298 edm::LogWarning("MaterialEffects")
0299 << " WARNING: A non valid vertex was found in nucl. int. -> " << ivertex << std::endl;
0300 }
0301 }
0302 }
0303
0304 if (myTrack.particle().charge() == 0)
0305 return;
0306
0307 if (!MuonBremsstrahlung && !Bremsstrahlung && !EnergyLoss && !MultipleScattering)
0308 return;
0309
0310
0311
0312
0313
0314 if (Bremsstrahlung && abs(myTrack.particle().pid()) == 11) {
0315 Bremsstrahlung->updateState(myTrack, radlen, random);
0316
0317 if (Bremsstrahlung->nDaughters()) {
0318
0319
0320 int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::BREM_VERTEX);
0321
0322
0323 if (ivertex >= 0) {
0324 for (DaughterIter = Bremsstrahlung->beginDaughters(); DaughterIter != Bremsstrahlung->endDaughters();
0325 ++DaughterIter) {
0326 mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0327 }
0328 } else {
0329 edm::LogWarning("MaterialEffects")
0330 << " WARNING: A non valid vertex was found in brem -> " << ivertex << std::endl;
0331 }
0332 }
0333 }
0334
0335
0336
0337
0338
0339 if (MuonBremsstrahlung && abs(myTrack.particle().pid()) == 13) {
0340 MuonBremsstrahlung->updateState(myTrack, radlen, random);
0341
0342 if (MuonBremsstrahlung->nDaughters()) {
0343
0344
0345 int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::BREM_VERTEX);
0346
0347
0348 if (ivertex >= 0) {
0349 for (DaughterIter = MuonBremsstrahlung->beginDaughters(); DaughterIter != MuonBremsstrahlung->endDaughters();
0350 ++DaughterIter) {
0351 mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0352 }
0353 } else {
0354 edm::LogWarning("MaterialEffects")
0355 << " WARNING: A non valid vertex was found in muon brem -> " << ivertex << std::endl;
0356 }
0357 }
0358 }
0359
0360
0361
0362
0363
0364 if (EnergyLoss) {
0365 theEnergyLoss = myTrack.particle().E();
0366 EnergyLoss->updateState(myTrack, radlen, random);
0367 theEnergyLoss -= myTrack.particle().E();
0368 }
0369
0370
0371
0372
0373
0374 if (MultipleScattering && myTrack.particle().Pt() > pTmin) {
0375
0376 MultipleScattering->setNormalVector(theNormalVector);
0377 MultipleScattering->updateState(myTrack, radlen, random);
0378 }
0379 }
0380
0381 double MaterialEffects::radLengths(const TrackerLayer& layer, ParticlePropagator& myTrack) {
0382
0383 theThickness = layer.surface().mediumProperties().radLen();
0384
0385 GlobalVector P(myTrack.particle().Px(), myTrack.particle().Py(), myTrack.particle().Pz());
0386
0387
0388
0389 double radlen = theThickness / fabs(P.dot(theNormalVector)) * P.mag();
0390
0391
0392
0393 double rad = myTrack.particle().R();
0394 double zed = fabs(myTrack.particle().Z());
0395
0396 double factor = 1;
0397
0398
0399 if (layer.fudgeNumber())
0400
0401
0402 for (unsigned int iLayer = 0; iLayer < layer.fudgeNumber(); ++iLayer) {
0403
0404 if ((layer.forward() && layer.fudgeMin(iLayer) < rad && rad < layer.fudgeMax(iLayer)) ||
0405 (!layer.forward() && layer.fudgeMin(iLayer) < zed && zed < layer.fudgeMax(iLayer))) {
0406 factor = layer.fudgeFactor(iLayer);
0407 break;
0408 }
0409 }
0410
0411 theThickness *= factor;
0412
0413 return radlen * factor;
0414 }
0415
0416 GlobalVector MaterialEffects::normalVector(const TrackerLayer& layer, ParticlePropagator& myTrack) const {
0417 return layer.forward() ? layer.disk()->normalVector()
0418 : GlobalVector(myTrack.particle().X(), myTrack.particle().Y(), 0.) / myTrack.particle().R();
0419 }
0420
0421 void MaterialEffects::save() {
0422
0423 if (NuclearInteraction)
0424 NuclearInteraction->save();
0425 }