Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:44

0001 #include "Randomize.hh"
0002 #include "G4ParticleTable.hh"
0003 
0004 #include <iostream>
0005 #include <fstream>
0006 #include <string>
0007 
0008 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
0009 #include "SimG4Core/CustomPhysics/interface/HadronicProcessHelper.h"
0010 
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/FileInPath.h"
0013 
0014 using namespace CLHEP;
0015 
0016 HadronicProcessHelper::HadronicProcessHelper(const std::string& fileName) {
0017   m_particleTable = G4ParticleTable::GetParticleTable();
0018   m_proton = m_particleTable->FindParticle("proton");
0019   m_neutron = m_particleTable->FindParticle("neutron");
0020 
0021   G4String line;
0022 
0023   std::ifstream processListStream(fileName.c_str());
0024 
0025   while (getline(processListStream, line)) {
0026     std::vector<G4String> tokens;
0027 
0028     //Getting a line
0029     m_readAndParse(line, tokens, "#");
0030 
0031     //Important info
0032     G4String incident = tokens[0];
0033     G4ParticleDefinition* incidentDef = m_particleTable->FindParticle(incident);
0034     G4int incidentPDG = incidentDef->GetPDGEncoding();
0035     m_knownParticles[incidentDef] = true;
0036     //    G4cout<<"Incident: "<<incident<<G4endl;
0037     G4String target = tokens[1];
0038     //    G4cout<<"Target: "<<target<<G4endl;
0039 
0040     // Making a ReactionProduct
0041     ReactionProduct prod;
0042     for (size_t i = 2; i != tokens.size(); i++) {
0043       G4String part = tokens[i];
0044       if (m_particleTable->contains(part)) {
0045         prod.push_back(m_particleTable->FindParticle(part)->GetPDGEncoding());
0046       } else {
0047         G4Exception("HadronicProcessHelper",
0048                     "UnkownParticle",
0049                     FatalException,
0050                     "Initialization: The reaction product list contained an unknown particle");
0051       }
0052     }
0053     if (target == "proton") {
0054       m_protonReactionMap[incidentPDG].push_back(prod);
0055     } else if (target == "neutron") {
0056       m_neutronReactionMap[incidentPDG].push_back(prod);
0057     } else {
0058       G4Exception("HadronicProcessHelper",
0059                   "IllegalTarget",
0060                   FatalException,
0061                   "Initialization: The reaction product list contained an illegal target particle");
0062     }
0063   }
0064 
0065   processListStream.close();
0066 
0067   m_checkFraction = 0;
0068   m_n22 = 0;
0069   m_n23 = 0;
0070 }
0071 
0072 G4bool HadronicProcessHelper::applicabilityTester(const G4ParticleDefinition& aPart) {
0073   const G4ParticleDefinition* aP = &aPart;
0074   if (m_knownParticles[aP])
0075     return true;
0076   return false;
0077 }
0078 
0079 G4double HadronicProcessHelper::inclusiveCrossSection(const G4DynamicParticle* particle, const G4Element* element) {
0080   //We really do need a dedicated class to handle the cross sections. They might not always be constant
0081 
0082   G4int pdgCode = particle->GetDefinition()->GetPDGEncoding();
0083 
0084   //24mb for gluino-balls
0085   if (CustomPDGParser::s_isRGlueball(pdgCode))
0086     return 24 * millibarn * element->GetN();
0087 
0088   //get quark vector
0089   std::vector<G4int> quarks = CustomPDGParser::s_containedQuarks(pdgCode);
0090 
0091   G4double totalNucleonCrossSection = 0;
0092 
0093   for (std::vector<G4int>::iterator it = quarks.begin(); it != quarks.end(); it++) {
0094     // 12mb for each 'up' or 'down'
0095     if (*it == 1 || *it == 2)
0096       totalNucleonCrossSection += 12 * millibarn;
0097     //  6mb for each 'strange'
0098     if (*it == 3)
0099       totalNucleonCrossSection += 6 * millibarn;
0100   }
0101 
0102   //Convert to xsec per nucleus
0103   //  return totalNucleonCrossSection * element->GetN();
0104   return totalNucleonCrossSection * pow(element->GetN(), 0.7) * 1.25;  // * 0.523598775598299;
0105 }
0106 
0107 HadronicProcessHelper::ReactionProduct HadronicProcessHelper::finalState(
0108     const G4DynamicParticle* incidentDynamicParticle, const G4Material* material, G4ParticleDefinition*& target) {
0109   //  const G4DynamicParticle* incidentDynamicParticle = track.GetDynamicParticle();
0110 
0111   //-----------------------------------------------
0112   // Choose n / p as target
0113   // and get ReactionProductList pointer
0114   //-----------------------------------------------
0115   ReactionMap* m_reactionMap;
0116   //G4Material* material = track.GetMaterial();
0117   const G4ElementVector* elementVector = material->GetElementVector();
0118   const G4double* numberOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
0119 
0120   G4double numberOfProtons = 0;
0121   G4double numberOfNucleons = 0;
0122 
0123   //Loop on elements
0124   for (size_t elm = 0; elm < material->GetNumberOfElements(); elm++) {
0125     //Summing number of protons per unit volume
0126     numberOfProtons += numberOfAtomsPerVolume[elm] * (*elementVector)[elm]->GetZ();
0127     //Summing nucleons (not neutrons)
0128     numberOfNucleons += numberOfAtomsPerVolume[elm] * (*elementVector)[elm]->GetN();
0129   }
0130 
0131   //random decision of the target
0132   if (G4UniformRand() < numberOfProtons / numberOfNucleons) {  //target is a proton
0133     m_reactionMap = &m_protonReactionMap;
0134     target = m_proton;
0135   } else {  //target is a neutron
0136     m_reactionMap = &m_neutronReactionMap;
0137     target = m_neutron;
0138   }
0139 
0140   const G4int incidentPDG = incidentDynamicParticle->GetDefinition()->GetPDGEncoding();
0141 
0142   //Making a pointer directly to the ReactionProductList we are looking at. Makes life easier :-)
0143   ReactionProductList* reactionProductList = &((*m_reactionMap)[incidentPDG]);
0144 
0145   //-----------------------------------------------
0146   // Count processes
0147   // kinematic check
0148   // compute number of 2 -> 2 and 2 -> 3 processes
0149   //-----------------------------------------------
0150 
0151   G4int good22 = 0;  //Number of 2 -> 2 processes that are possible
0152   G4int good23 = 0;  //Number of 2 -> 3 processes that are possible
0153 
0154   //This is the list to be populated
0155   ReactionProductList goodReactionProductList;
0156 
0157   for (ReactionProductList::iterator prod_it = reactionProductList->begin(); prod_it != reactionProductList->end();
0158        prod_it++) {
0159     G4int secondaries = prod_it->size();
0160     // If the reaction is not possible we will not consider it
0161     if (m_reactionIsPossible(*prod_it, incidentDynamicParticle, target)) {
0162       // The reaction is possible. Let's store and count it
0163       goodReactionProductList.push_back(*prod_it);
0164       if (secondaries == 2) {
0165         good22++;
0166       } else if (secondaries == 3) {
0167         good23++;
0168       } else {
0169         G4cerr << "ReactionProduct has unsupported number of secondaries: " << secondaries << G4endl;
0170       }
0171     } /*else {
0172       G4cout<<"There was an impossible process"<<G4endl;
0173       }*/
0174   }
0175   //  G4cout<<"The size of the ReactionProductList is: "<<theReactionProductList.size()<<G4endl;
0176 
0177   if (goodReactionProductList.empty())
0178     G4Exception("HadronicProcessHelper",
0179                 "NoProcessPossible",
0180                 FatalException,
0181                 "GetFinalState: No process could be selected from the given list.");
0182   // Fill a probability map. Remember total probability
0183   // 2->2 is 0.15*1/n_22 2->3 uses phase space
0184   G4double prob22 = 0.15;
0185   G4double prob23 = 1 - prob22;  // :-)
0186 
0187   std::vector<G4double> probabilities;
0188   std::vector<G4bool> twoToThreeFlag;
0189 
0190   G4double cumulatedProbability = 0;
0191 
0192   // To each ReactionProduct we assign a cumulated probability and a flag
0193   // discerning between 2 -> 2 and 2 -> 3
0194   size_t numberOfReactions = goodReactionProductList.size();
0195   for (size_t i = 0; i != numberOfReactions; i++) {
0196     if (goodReactionProductList[i].size() == 2) {
0197       cumulatedProbability += prob22 / good22;
0198       twoToThreeFlag.push_back(false);
0199     } else {
0200       cumulatedProbability += prob23 / good23;
0201       twoToThreeFlag.push_back(true);
0202     }
0203     probabilities.push_back(cumulatedProbability);
0204   }
0205 
0206   //Normalising probabilities to 1
0207   for (std::vector<G4double>::iterator it = probabilities.begin(); it != probabilities.end(); it++) {
0208     *it /= cumulatedProbability;
0209   }
0210 
0211   // Choosing ReactionProduct
0212   G4bool selected = false;
0213   G4int tries = 0;
0214   //  ReactionProductList::iterator prod_it;
0215 
0216   //Keep looping over the list until we have a choice, or until we have tried 100 times
0217   size_t i;
0218   while (!selected && tries < 100) {
0219     i = 0;
0220     G4double dice = G4UniformRand();
0221     //select the process using the dice
0222     while (dice > probabilities[i] && i < numberOfReactions)
0223       i++;
0224 
0225     if (twoToThreeFlag[i]) {
0226       // 2 -> 3 processes require a phase space lookup
0227       if (m_phaseSpace(goodReactionProductList[i], incidentDynamicParticle, target) > G4UniformRand())
0228         selected = true;
0229     } else {
0230       // 2 -> 2 processes are chosen immediately
0231       selected = true;
0232     }
0233     tries++;
0234   }
0235   if (tries >= 100)
0236     G4cerr << "Could not select process!!!!" << G4endl;
0237 
0238   //Debugging stuff
0239   //Updating checkfraction:
0240   if (goodReactionProductList[i].size() == 2) {
0241     m_n22++;
0242   } else {
0243     m_n23++;
0244   }
0245   m_checkFraction = (1.0 * m_n22) / (m_n22 + m_n23);
0246 
0247   //return the selected productList
0248   return goodReactionProductList[i];
0249 }
0250 
0251 G4double HadronicProcessHelper::m_reactionProductMass(const ReactionProduct& reactionProd,
0252                                                       const G4DynamicParticle* incidentDynamicParticle,
0253                                                       G4ParticleDefinition* target) {
0254   // Incident energy:
0255   G4double incidentEnergy = incidentDynamicParticle->GetTotalEnergy();
0256   G4double m_1 = incidentDynamicParticle->GetDefinition()->GetPDGMass();
0257   G4double m_2 = target->GetPDGMass();
0258   //square root of "s"
0259   G4double sqrtS = sqrt(m_1 * m_1 + m_2 * (m_2 + 2 * incidentEnergy));
0260   // Sum of rest masses after reaction:
0261   G4double productsMass = 0;
0262   //Loop on reaction producs
0263   for (ReactionProduct::const_iterator r_it = reactionProd.begin(); r_it != reactionProd.end(); r_it++) {
0264     //Sum the masses of the products
0265     productsMass += m_particleTable->FindParticle(*r_it)->GetPDGMass();
0266   }
0267   //the result is square root of "s" minus the masses of the products
0268   return sqrtS - productsMass;
0269 }
0270 
0271 G4bool HadronicProcessHelper::m_reactionIsPossible(const ReactionProduct& aReaction,
0272                                                    const G4DynamicParticle* aDynamicParticle,
0273                                                    G4ParticleDefinition* target) {
0274   if (m_reactionProductMass(aReaction, aDynamicParticle, target) > 0)
0275     return true;
0276   return false;
0277 }
0278 
0279 G4double HadronicProcessHelper::m_phaseSpace(const ReactionProduct& aReaction,
0280                                              const G4DynamicParticle* aDynamicParticle,
0281                                              G4ParticleDefinition* target) {
0282   G4double qValue = m_reactionProductMass(aReaction, aDynamicParticle, target);
0283   G4double phi = sqrt(1 + qValue / (2 * 0.139 * GeV)) * pow(qValue / (1.1 * GeV), 3. / 2.);
0284   return (phi / (1 + phi));
0285 }
0286 
0287 void HadronicProcessHelper::m_readAndParse(const G4String& str,
0288                                            std::vector<G4String>& tokens,
0289                                            const G4String& delimiters) {
0290   // Skip delimiters at beginning.
0291   G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
0292   // Find first "non-delimiter".
0293   G4String::size_type pos = str.find_first_of(delimiters, lastPos);
0294 
0295   while (G4String::npos != pos || G4String::npos != lastPos) {
0296     //Skipping leading / trailing whitespaces
0297     G4String temp = str.substr(lastPos, pos - lastPos);
0298     while (temp.c_str()[0] == ' ')
0299       temp.erase(0, 1);
0300     while (temp[temp.size() - 1] == ' ')
0301       temp.erase(temp.size() - 1, 1);
0302     // Found a token, add it to the vector.
0303     tokens.push_back(temp);
0304     // Skip delimiters.  Note the "not_of"
0305     lastPos = str.find_first_not_of(delimiters, pos);
0306     // Find next "non-delimiter"
0307     pos = str.find_first_of(delimiters, lastPos);
0308   }
0309 }