CustomPhysics/src/G4ProcessHelper.cc

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     //Getting a line
0058     ReadAndParse(line, tokens, "#");
0059     //Important info
0060     G4String incident = tokens[0];
0061
0062     G4ParticleDefinition* incidentDef = particleTable->FindParticle(incident);
0063     //particleTable->DumpTable();
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     // Making a ReactionProduct
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   //We really do need a dedicated class to handle the cross sections. They might not always be constant
0118
0119   //Disassemble the PDG-code
0120
0121   G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
0122   double boost = (aParticle->GetKineticEnergy() + aParticle->GetMass()) / aParticle->GetMass();
0123   //  G4cout<<"thePDGCode: "<<thePDGCode<<G4endl;
0124   G4double theXsec = 0;
0125   G4String name = aParticle->GetDefinition()->GetParticleName();
0126   if (!reggemodel) {
0127     //Flat cross section
0128     if (CustomPDGParser::s_isRGlueball(thePDGCode)) {
0129       theXsec = 24 * millibarn;
0130     } else {
0131       std::vector<G4int> nq = CustomPDGParser::s_containedQuarks(thePDGCode);
0132       //    edm::LogInfo("SimG4CoreCustomPhysics")<<"Number of quarks: "<<nq.size()<<G4endl;
0133       for (std::vector<G4int>::iterator it = nq.begin(); it != nq.end(); it++) {
0134         //    edm::LogInfo("SimG4CoreCustomPhysics")<<"Quarkvector: "<<*it<<G4endl;
0135         if (*it == 1 || *it == 2)
0136           theXsec += 12 * millibarn;
0137         if (*it == 3)
0138           theXsec += 6 * millibarn;
0139       }
0140     }
0141   } else {  //reggemodel
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   //Adding resonance
0165   if (resonant) {
0166     double e_0 = ek_0 + aParticle->GetDefinition()->GetPDGMass();  //Now total energy
0167
0168     e_0 = sqrt(aParticle->GetDefinition()->GetPDGMass() * aParticle->GetDefinition()->GetPDGMass() +
0169                theProton->GetPDGMass() * theProton->GetPDGMass() + 2. * e_0 * theProton->GetPDGMass());
0170     //      edm::LogInfo("SimG4CoreCustomPhysics")<<e_0/GeV<<G4endl;
0171     //      edm::LogInfo("SimG4CoreCustomPhysics")<<ek_0/GeV<<" "<<aParticle->GetDefinition()->GetPDGMass()/GeV<<" "<<theProton->GetPDGMass()/GeV<<G4endl;
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.));  //Non-relativistic Breit Wigner
0178
0179     theXsec += res_result;
0180     //      if(fabs(aParticle->GetKineticEnergy()/GeV-200)<10)  std::cout<<sqrts/GeV<<" "<<theXsec/millibarn<<std::endl;
0181   }
0182
0183   //  std::cout<<"Xsec/nucleon: "<<theXsec/millibarn<<"millibarn, Total Xsec: "<<theXsec * anElement->GetN()/millibarn<<" millibarn"<<std::endl;
0184   return theXsec * pow(anElement->GetN(), 0.7) * 1.25;  // * 0.523598775598299;
0185 }
0186
0187 ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4ParticleDefinition*& aTarget) {
0188   const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
0189
0190   //-----------------------------------------------
0191   // Choose n / p as target
0192   // and get ReactionProductList pointer
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     //Summing number of protons per unit volume
0204     NumberOfProtons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetZ();
0205     //Summing nucleons (not neutrons)
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     //      G4cout<<"Oscillating..."<<G4endl;
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   //Making a pointer directly to the ReactionProductList we are looking at. Makes life easier :-)
0234   ReactionProductList* aReactionProductList = &((*theReactionMap)[theIncidentPDG]);
0235
0236   //-----------------------------------------------
0237   // Count processes
0238   // kinematic check
0239   // compute number of 2 -> 2 and 2 -> 3 processes
0240   //-----------------------------------------------
0241
0242   G4int N22 = 0;  //Number of 2 -> 2 processes
0243   G4int N23 = 0;  //Number of 2 -> 3 processes. Couldn't think of more informative names
0244
0245   //This is the list to be populated
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     // If the reaction is not possible we will not consider it
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       // The reaction is possible. Let's store and count it
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     } /*else {
0267       edm::LogInfo("SimG4CoreCustomPhysics")<<"There was an impossible process"<<G4endl;
0268       }*/
0269   }
0270   //  edm::LogInfo("SimG4CoreCustomPhysics")<<"The size of the ReactionProductList is: "<<theReactionProductList.size()<<G4endl;
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   // For the Regge model no phase space considerations. We pick a process at random
0279   if (reggemodel) {
0280     int n_rps = theReactionProductList.size();
0281     int select = (int)(G4UniformRand() * n_rps);
0282     //      G4cout<<"Possible: "<<n_rps<<", chosen: "<<select<<G4endl;
0283     return theReactionProductList[select];
0284   }
0285
0286   // Fill a probability map. Remember total probability
0287   // 2->2 is 0.15*1/n_22 2->3 uses phase space
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   // To each ReactionProduct we assign a cumulated probability and a flag
0297   // discerning between 2 -> 2 and 2 -> 3
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     //    edm::LogInfo("SimG4CoreCustomPhysics")<<"Pushing back cumulated probability: "<<CumulatedProbability<<G4endl;
0308   }
0309
0310   //Renormalising probabilities
0311   //  edm::LogInfo("SimG4CoreCustomPhysics")<<"Probs: ";
0312   for (std::vector<G4double>::iterator it = Probabilities.begin(); it != Probabilities.end(); it++) {
0313     *it /= CumulatedProbability;
0314     //      edm::LogInfo("SimG4CoreCustomPhysics")<<*it<<" ";
0315   }
0316   //  edm::LogInfo("SimG4CoreCustomPhysics")<<G4endl;
0317
0318   // Choosing ReactionProduct
0319
0320   G4bool selected = false;
0321   G4int tries = 0;
0322   //  ReactionProductList::iterator prod_it;
0323
0324   //Keep looping over the list until we have a choice, or until we have tried 100 times
0325   unsigned int i;
0326   while (!selected && tries < 100) {
0327     i = 0;
0328     G4double dice = G4UniformRand();
0329     // edm::LogInfo("SimG4CoreCustomPhysics")<<"What's the dice?"<<dice<<G4endl;
0330     while (dice > Probabilities[i] && i < theReactionProductList.size()) {
0331       //      edm::LogInfo("SimG4CoreCustomPhysics")<<"i: "<<i<<G4endl;
0332       i++;
0333     }
0334
0335     //    edm::LogInfo("SimG4CoreCustomPhysics")<<"Chosen i: "<<i<<G4endl;
0336
0337     if (!TwotoThreeFlag[i]) {
0338       // 2 -> 2 processes are chosen immediately
0339       selected = true;
0340     } else {
0341       // 2 -> 3 processes require a phase space lookup
0342       if (PhaseSpace(theReactionProductList[i], aDynamicParticle) > G4UniformRand())
0343         selected = true;
0344       //selected = true;
0345     }
0346     //    double suppressionfactor=0.5;
0347     if (selected && particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge() !=
0348                         aDynamicParticle->GetDefinition()->GetPDGCharge()) {
0349       /*
0350     edm::LogInfo("SimG4CoreCustomPhysics")<<"Incoming particle "<<aDynamicParticle->GetDefinition()->GetParticleName()
0351           <<" has charge "<<aDynamicParticle->GetDefinition()->GetPDGCharge()<<G4endl;
0352     edm::LogInfo("SimG4CoreCustomPhysics")<<"Suggested particle "<<particleTable->FindParticle(theReactionProductList[i][0])->GetParticleName()
0353           <<" has charge "<<particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge()<<G4endl;
0354     */
0355       if (G4UniformRand() < suppressionfactor)
0356         selected = false;
0357     }
0358     tries++;
0359     //    edm::LogInfo("SimG4CoreCustomPhysics")<<"Tries: "<<tries<<G4endl;
0360   }
0361   if (tries >= 100)
0362     G4cerr << "Could not select process!!!!" << G4endl;
0363
0364   //  edm::LogInfo("SimG4CoreCustomPhysics")<<"So far so good"<<G4endl;
0365   //  edm::LogInfo("SimG4CoreCustomPhysics")<<"Sec's: "<<theReactionProductList[i].size()<<G4endl;
0366
0367   //Updating checkfraction:
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   //  edm::LogInfo("SimG4CoreCustomPhysics")<<"n_22: "<<n_22<<" n_23: "<<n_23<<" Checkfraction: "<<checkfraction<<G4endl;
0376   //  edm::LogInfo("SimG4CoreCustomPhysics") <<"Biig number: "<<n_22+n_23<<G4endl;
0377   //Return the chosen ReactionProduct
0378   return theReactionProductList[i];
0379 }
0380
0381 G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,
0382                                               const G4DynamicParticle* aDynamicParticle) {
0383   // Incident energy:
0384   G4double E_incident = aDynamicParticle->GetTotalEnergy();
0385   //edm::LogInfo("SimG4CoreCustomPhysics")<<"Total energy: "<<E_incident<<" Kinetic: "<<aDynamicParticle->GetKineticEnergy()<<G4endl;
0386   // sqrt(s)= sqrt(m_1^2 + m_2^2 + 2 E_1 m_2)
0387   G4double m_1 = aDynamicParticle->GetDefinition()->GetPDGMass();
0388   G4double m_2 = theTarget->GetPDGMass();
0389   //edm::LogInfo("SimG4CoreCustomPhysics")<<"M_R: "<<m_1/GeV<<" GeV, M_np: "<<m_2/GeV<<" GeV"<<G4endl;
0390   G4double sqrts = sqrt(m_1 * m_1 + m_2 * (m_2 + 2 * E_incident));
0391   //edm::LogInfo("SimG4CoreCustomPhysics")<<"sqrt(s) = "<<sqrts/GeV<<" GeV"<<G4endl;
0392   // Sum of rest masses after reaction:
0393   G4double M_after = 0;
0394   for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it != aReaction.end(); r_it++) {
0395     //edm::LogInfo("SimG4CoreCustomPhysics")<<"Mass contrib: "<<(particleTable->FindParticle(*r_it)->GetPDGMass())/MeV<<" MeV"<<G4endl;
0396     M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
0397   }
0398   //edm::LogInfo("SimG4CoreCustomPhysics")<<"Intending to return this ReactionProductMass: "<<(sqrts - M_after)/MeV<<" MeV"<<G4endl;
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   // Skip delimiters at beginning.
0424   G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
0425   // Find first "non-delimiter".
0426   G4String::size_type pos = str.find_first_of(delimiters, lastPos);
0427
0428   while (G4String::npos != pos || G4String::npos != lastPos) {
0429     //Skipping leading / trailing whitespaces
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     // Found a token, add it to the vector.
0436     tokens.push_back(temp);
0437     // Skip delimiters.  Note the "not_of"
0438     lastPos = str.find_first_not_of(delimiters, pos);
0439     // Find next "non-delimiter"
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 }