File indexing completed on 2023-10-25 10:04:28
0001 #include "SimG4Core/CustomPhysics/interface/G4ProcessHelper.h"
0002 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
0003 #include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
0004 #include "SimG4Core/CustomPhysics/interface/CustomParticleFactory.h"
0005
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ParameterSet/interface/FileInPath.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include "G4ParticleTable.hh"
0011 #include "Randomize.hh"
0012
0013 #include <iostream>
0014 #include <fstream>
0015 #include <string>
0016
0017 using namespace CLHEP;
0018
0019 G4ProcessHelper::G4ProcessHelper(const edm::ParameterSet& p, CustomParticleFactory* ptr) {
0020 fParticleFactory = ptr;
0021
0022 particleTable = G4ParticleTable::GetParticleTable();
0023
0024 theProton = particleTable->FindParticle("proton");
0025 theNeutron = particleTable->FindParticle("neutron");
0026
0027 G4String line;
0028
0029 edm::FileInPath fp = p.getParameter<edm::FileInPath>("processesDef");
0030 std::string processDefFilePath = fp.fullPath();
0031 std::ifstream process_stream(processDefFilePath.c_str());
0032
0033 resonant = p.getParameter<bool>("resonant");
0034 ek_0 = p.getParameter<double>("resonanceEnergy") * GeV;
0035 gamma = p.getParameter<double>("gamma") * GeV;
0036 amplitude = p.getParameter<double>("amplitude") * millibarn;
0037 suppressionfactor = p.getParameter<double>("reggeSuppression");
0038 hadronlifetime = p.getParameter<double>("hadronLifeTime");
0039 reggemodel = p.getParameter<bool>("reggeModel");
0040 mixing = p.getParameter<double>("mixing");
0041
0042 edm::LogInfo("SimG4CoreCustomPhysics") << "ProcessHelper: Read in physics parameters:"
0043 << "\n Resonant = " << resonant << "\n ResonanceEnergy = " << ek_0 / GeV
0044 << " GeV"
0045 << "\n Gamma = " << gamma / GeV << " GeV"
0046 << "\n Amplitude = " << amplitude / millibarn << " millibarn"
0047 << "ReggeSuppression = " << 100 * suppressionfactor << " %"
0048 << "HadronLifeTime = " << hadronlifetime << " s"
0049 << "ReggeModel = " << reggemodel << "Mixing = " << mixing * 100 << " %";
0050
0051 checkfraction = 0;
0052 n_22 = 0;
0053 n_23 = 0;
0054
0055 while (getline(process_stream, line)) {
0056 std::vector<G4String> tokens;
0057
0058 ReadAndParse(line, tokens, "#");
0059
0060 G4String incident = tokens[0];
0061
0062 G4ParticleDefinition* incidentDef = particleTable->FindParticle(incident);
0063
0064 G4int incidentPDG = incidentDef->GetPDGEncoding();
0065 known_particles[incidentDef] = true;
0066
0067 G4String target = tokens[1];
0068 edm::LogInfo("SimG4CoreCustomPhysics") << "ProcessHelper: Incident " << incident << "; Target " << target;
0069
0070
0071 ReactionProduct prod;
0072 for (size_t i = 2; i != tokens.size(); i++) {
0073 G4String part = tokens[i];
0074 if (particleTable->contains(part)) {
0075 prod.push_back(particleTable->FindParticle(part)->GetPDGEncoding());
0076 } else {
0077 G4Exception("G4ProcessHelper",
0078 "UnkownParticle",
0079 FatalException,
0080 "Initialization: The reaction product list contained an unknown particle");
0081 }
0082 }
0083 if (target == "proton") {
0084 pReactionMap[incidentPDG].push_back(prod);
0085 } else if (target == "neutron") {
0086 nReactionMap[incidentPDG].push_back(prod);
0087 } else {
0088 G4Exception("G4ProcessHelper",
0089 "IllegalTarget",
0090 FatalException,
0091 "Initialization: The reaction product list contained an illegal target particle");
0092 }
0093 }
0094
0095 process_stream.close();
0096
0097 for (auto part : fParticleFactory->getCustomParticles()) {
0098 CustomParticle* particle = dynamic_cast<CustomParticle*>(part);
0099 if (particle) {
0100 edm::LogInfo("SimG4CoreCustomPhysics") << "ProcessHelper: Lifetime of " << part->GetParticleName() << " set to "
0101 << particle->GetPDGLifeTime() / s << " s;"
0102 << " isStable: " << particle->GetPDGStable();
0103 }
0104 }
0105 }
0106
0107 G4ProcessHelper::~G4ProcessHelper() {}
0108
0109 G4bool G4ProcessHelper::ApplicabilityTester(const G4ParticleDefinition& aPart) {
0110 const G4ParticleDefinition* aP = &aPart;
0111 if (known_particles[aP])
0112 return true;
0113 return false;
0114 }
0115
0116 G4double G4ProcessHelper::GetInclusiveCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement) {
0117
0118
0119
0120
0121 G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
0122 double boost = (aParticle->GetKineticEnergy() + aParticle->GetMass()) / aParticle->GetMass();
0123
0124 G4double theXsec = 0;
0125 G4String name = aParticle->GetDefinition()->GetParticleName();
0126 if (!reggemodel) {
0127
0128 if (CustomPDGParser::s_isRGlueball(thePDGCode)) {
0129 theXsec = 24 * millibarn;
0130 } else {
0131 std::vector<G4int> nq = CustomPDGParser::s_containedQuarks(thePDGCode);
0132
0133 for (std::vector<G4int>::iterator it = nq.begin(); it != nq.end(); it++) {
0134
0135 if (*it == 1 || *it == 2)
0136 theXsec += 12 * millibarn;
0137 if (*it == 3)
0138 theXsec += 6 * millibarn;
0139 }
0140 }
0141 } else {
0142 double R = Regge(boost);
0143 double P = Pom(boost);
0144 if (thePDGCode > 0) {
0145 if (CustomPDGParser::s_isMesonino(thePDGCode))
0146 theXsec = (P + R) * millibarn;
0147 if (CustomPDGParser::s_isSbaryon(thePDGCode))
0148 theXsec = 2 * P * millibarn;
0149 if (CustomPDGParser::s_isRMeson(thePDGCode) || CustomPDGParser::s_isRGlueball(thePDGCode))
0150 theXsec = (R + 2 * P) * millibarn;
0151 if (CustomPDGParser::s_isRBaryon(thePDGCode))
0152 theXsec = 3 * P * millibarn;
0153 } else {
0154 if (CustomPDGParser::s_isMesonino(thePDGCode))
0155 theXsec = P * millibarn;
0156 if (CustomPDGParser::s_isSbaryon(thePDGCode))
0157 theXsec = (2 * (P + R) + 30 / sqrt(boost)) * millibarn;
0158 if (CustomPDGParser::s_isRMeson(thePDGCode) || CustomPDGParser::s_isRGlueball(thePDGCode))
0159 theXsec = (R + 2 * P) * millibarn;
0160 if (CustomPDGParser::s_isRBaryon(thePDGCode))
0161 theXsec = 3 * P * millibarn;
0162 }
0163 }
0164
0165 if (resonant) {
0166 double e_0 = ek_0 + aParticle->GetDefinition()->GetPDGMass();
0167
0168 e_0 = sqrt(aParticle->GetDefinition()->GetPDGMass() * aParticle->GetDefinition()->GetPDGMass() +
0169 theProton->GetPDGMass() * theProton->GetPDGMass() + 2. * e_0 * theProton->GetPDGMass());
0170
0171
0172 double sqrts = sqrt(aParticle->GetDefinition()->GetPDGMass() * aParticle->GetDefinition()->GetPDGMass() +
0173 theProton->GetPDGMass() * theProton->GetPDGMass() +
0174 2 * aParticle->GetTotalEnergy() * theProton->GetPDGMass());
0175
0176 double res_result = amplitude * (gamma * gamma / 4.) /
0177 ((sqrts - e_0) * (sqrts - e_0) + (gamma * gamma / 4.));
0178
0179 theXsec += res_result;
0180
0181 }
0182
0183
0184 return theXsec * pow(anElement->GetN(), 0.7) * 1.25;
0185 }
0186
0187 ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4ParticleDefinition*& aTarget) {
0188 const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
0189
0190
0191
0192
0193
0194
0195 G4Material* aMaterial = aTrack.GetMaterial();
0196 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
0197 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
0198
0199 G4double NumberOfProtons = 0;
0200 G4double NumberOfNucleons = 0;
0201
0202 for (size_t elm = 0; elm < aMaterial->GetNumberOfElements(); elm++) {
0203
0204 NumberOfProtons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetZ();
0205
0206 NumberOfNucleons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetN();
0207 }
0208
0209 if (G4UniformRand() < NumberOfProtons / NumberOfNucleons) {
0210 theReactionMap = &pReactionMap;
0211 theTarget = theProton;
0212 } else {
0213 theReactionMap = &nReactionMap;
0214 theTarget = theNeutron;
0215 }
0216 aTarget = theTarget;
0217
0218 G4int theIncidentPDG = aDynamicParticle->GetDefinition()->GetPDGEncoding();
0219
0220 if (reggemodel && CustomPDGParser::s_isMesonino(theIncidentPDG) && G4UniformRand() * mixing > 0.5 &&
0221 aDynamicParticle->GetDefinition()->GetPDGCharge() == 0.) {
0222
0223 theIncidentPDG *= -1;
0224 }
0225
0226 bool baryonise = false;
0227
0228 if (reggemodel && G4UniformRand() > 0.9 &&
0229 ((CustomPDGParser::s_isMesonino(theIncidentPDG) && theIncidentPDG > 0) ||
0230 CustomPDGParser::s_isRMeson(theIncidentPDG)))
0231 baryonise = true;
0232
0233
0234 ReactionProductList* aReactionProductList = &((*theReactionMap)[theIncidentPDG]);
0235
0236
0237
0238
0239
0240
0241
0242 G4int N22 = 0;
0243 G4int N23 = 0;
0244
0245
0246 ReactionProductList theReactionProductList;
0247 std::vector<bool> theChargeChangeList;
0248
0249 for (ReactionProductList::iterator prod_it = aReactionProductList->begin(); prod_it != aReactionProductList->end();
0250 prod_it++) {
0251 G4int secondaries = prod_it->size();
0252
0253 if (ReactionIsPossible(*prod_it, aDynamicParticle) &&
0254 (!reggemodel || (baryonise && ReactionGivesBaryon(*prod_it)) ||
0255 (!baryonise && !ReactionGivesBaryon(*prod_it)) || (CustomPDGParser::s_isSbaryon(theIncidentPDG)) ||
0256 (CustomPDGParser::s_isRBaryon(theIncidentPDG)))) {
0257
0258 theReactionProductList.push_back(*prod_it);
0259 if (secondaries == 2) {
0260 N22++;
0261 } else if (secondaries == 3) {
0262 N23++;
0263 } else {
0264 G4cerr << "ReactionProduct has unsupported number of secondaries: " << secondaries << G4endl;
0265 }
0266 }
0267
0268
0269 }
0270
0271
0272 if (theReactionProductList.empty())
0273 G4Exception("G4ProcessHelper",
0274 "NoProcessPossible",
0275 FatalException,
0276 "GetFinalState: No process could be selected from the given list.");
0277
0278
0279 if (reggemodel) {
0280 int n_rps = theReactionProductList.size();
0281 int select = (int)(G4UniformRand() * n_rps);
0282
0283 return theReactionProductList[select];
0284 }
0285
0286
0287
0288 G4double p22 = 0.15;
0289 G4double p23 = 1 - p22;
0290
0291 std::vector<G4double> Probabilities;
0292 std::vector<G4bool> TwotoThreeFlag;
0293
0294 G4double CumulatedProbability = 0;
0295
0296
0297
0298 for (unsigned int i = 0; i != theReactionProductList.size(); i++) {
0299 if (theReactionProductList[i].size() == 2) {
0300 CumulatedProbability += p22 / N22;
0301 TwotoThreeFlag.push_back(false);
0302 } else {
0303 CumulatedProbability += p23 / N23;
0304 TwotoThreeFlag.push_back(true);
0305 }
0306 Probabilities.push_back(CumulatedProbability);
0307
0308 }
0309
0310
0311
0312 for (std::vector<G4double>::iterator it = Probabilities.begin(); it != Probabilities.end(); it++) {
0313 *it /= CumulatedProbability;
0314
0315 }
0316
0317
0318
0319
0320 G4bool selected = false;
0321 G4int tries = 0;
0322
0323
0324
0325 unsigned int i;
0326 while (!selected && tries < 100) {
0327 i = 0;
0328 G4double dice = G4UniformRand();
0329
0330 while (dice > Probabilities[i] && i < theReactionProductList.size()) {
0331
0332 i++;
0333 }
0334
0335
0336
0337 if (!TwotoThreeFlag[i]) {
0338
0339 selected = true;
0340 } else {
0341
0342 if (PhaseSpace(theReactionProductList[i], aDynamicParticle) > G4UniformRand())
0343 selected = true;
0344
0345 }
0346
0347 if (selected && particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge() !=
0348 aDynamicParticle->GetDefinition()->GetPDGCharge()) {
0349
0350
0351
0352
0353
0354
0355 if (G4UniformRand() < suppressionfactor)
0356 selected = false;
0357 }
0358 tries++;
0359
0360 }
0361 if (tries >= 100)
0362 G4cerr << "Could not select process!!!!" << G4endl;
0363
0364
0365
0366
0367
0368 if (theReactionProductList[i].size() == 2) {
0369 n_22++;
0370 } else {
0371 n_23++;
0372 }
0373
0374 checkfraction = (1.0 * n_22) / (n_22 + n_23);
0375
0376
0377
0378 return theReactionProductList[i];
0379 }
0380
0381 G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,
0382 const G4DynamicParticle* aDynamicParticle) {
0383
0384 G4double E_incident = aDynamicParticle->GetTotalEnergy();
0385
0386
0387 G4double m_1 = aDynamicParticle->GetDefinition()->GetPDGMass();
0388 G4double m_2 = theTarget->GetPDGMass();
0389
0390 G4double sqrts = sqrt(m_1 * m_1 + m_2 * (m_2 + 2 * E_incident));
0391
0392
0393 G4double M_after = 0;
0394 for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it != aReaction.end(); r_it++) {
0395
0396 M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
0397 }
0398
0399 return sqrts - M_after;
0400 }
0401
0402 G4bool G4ProcessHelper::ReactionIsPossible(const ReactionProduct& aReaction,
0403 const G4DynamicParticle* aDynamicParticle) {
0404 if (ReactionProductMass(aReaction, aDynamicParticle) > 0)
0405 return true;
0406 return false;
0407 }
0408
0409 G4bool G4ProcessHelper::ReactionGivesBaryon(const ReactionProduct& aReaction) {
0410 for (ReactionProduct::const_iterator it = aReaction.begin(); it != aReaction.end(); it++)
0411 if (CustomPDGParser::s_isSbaryon(*it) || CustomPDGParser::s_isRBaryon(*it))
0412 return true;
0413 return false;
0414 }
0415
0416 G4double G4ProcessHelper::PhaseSpace(const ReactionProduct& aReaction, const G4DynamicParticle* aDynamicParticle) {
0417 G4double qValue = ReactionProductMass(aReaction, aDynamicParticle);
0418 G4double phi = sqrt(1 + qValue / (2 * 0.139 * GeV)) * pow(qValue / (1.1 * GeV), 3. / 2.);
0419 return (phi / (1 + phi));
0420 }
0421
0422 void G4ProcessHelper::ReadAndParse(const G4String& str, std::vector<G4String>& tokens, const G4String& delimiters) {
0423
0424 G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
0425
0426 G4String::size_type pos = str.find_first_of(delimiters, lastPos);
0427
0428 while (G4String::npos != pos || G4String::npos != lastPos) {
0429
0430 G4String temp = str.substr(lastPos, pos - lastPos);
0431 while (temp.c_str()[0] == ' ')
0432 temp.erase(0, 1);
0433 while (temp[temp.size() - 1] == ' ')
0434 temp.erase(temp.size() - 1, 1);
0435
0436 tokens.push_back(temp);
0437
0438 lastPos = str.find_first_not_of(delimiters, pos);
0439
0440 pos = str.find_first_of(delimiters, lastPos);
0441 }
0442 }
0443
0444 double G4ProcessHelper::Regge(const double boost) {
0445 double a = 2.165635078566177;
0446 double b = 0.1467453738547229;
0447 double c = -0.9607903711871166;
0448 return 1.5 * exp(a + b / boost + c * log(boost));
0449 }
0450
0451 double G4ProcessHelper::Pom(const double boost) {
0452 double a = 4.138224000651535;
0453 double b = 1.50377557581421;
0454 double c = -0.05449742257808247;
0455 double d = 0.0008221235048211401;
0456 return a + b * sqrt(boost) + c * boost + d * pow(boost, 1.5);
0457 }