HadronicProcessHelper Class Reference

#include <HadronicProcessHelper.h>

Public Types
typedef std::map< G4int, ReactionProductList >	ReactionMap
typedef std::vector< G4int >	ReactionProduct
typedef std::vector< ReactionProduct >	ReactionProductList

Public Member Functions
G4bool	applicabilityTester (const G4ParticleDefinition &particle)
ReactionProduct	finalState (const G4Track &track, G4ParticleDefinition *&target)
ReactionProduct	finalState (const G4DynamicParticle *particle, const G4Material *material, G4ParticleDefinition *&target)
	HadronicProcessHelper (const std::string &fileName)
	HadronicProcessHelper (const HadronicProcessHelper &)
G4double	inclusiveCrossSection (const G4DynamicParticle *particle, const G4Element *element)
HadronicProcessHelper *	instance ()
HadronicProcessHelper &	operator= (const HadronicProcessHelper &)

Private Member Functions
G4double	m_phaseSpace (const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle, G4ParticleDefinition *target)
G4bool	m_reactionIsPossible (const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle, G4ParticleDefinition *target)
G4double	m_reactionProductMass (const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle, G4ParticleDefinition *target)
void	m_readAndParse (const G4String &str, std::vector< G4String > &tokens, const G4String &delimiters=" ")

Private Attributes
G4double	m_checkFraction
std::map< const G4ParticleDefinition *, G4bool >	m_knownParticles
G4int	m_n22
G4int	m_n23
G4ParticleDefinition *	m_neutron
ReactionMap	m_neutronReactionMap
G4ParticleTable *	m_particleTable
G4ParticleDefinition *	m_proton
ReactionMap	m_protonReactionMap

Detailed Description

Definition at line 16 of file HadronicProcessHelper.h.

Member Typedef Documentation

typedef std::map<G4int , ReactionProductList> HadronicProcessHelper::ReactionMap

Definition at line 23 of file HadronicProcessHelper.h.

typedef std::vector<G4int> HadronicProcessHelper::ReactionProduct

Definition at line 21 of file HadronicProcessHelper.h.

typedef std::vector<ReactionProduct > HadronicProcessHelper::ReactionProductList

Definition at line 22 of file HadronicProcessHelper.h.

Constructor & Destructor Documentation

HadronicProcessHelper::HadronicProcessHelper

(

const std::string &

fileName

)

Definition at line 16 of file HadronicProcessHelper.cc.

References mps_fire::i, mps_splice::line, parseEventContent::prod, and edmPickEvents::target.

                                                                      {

  m_particleTable = G4ParticleTable::GetParticleTable();
  m_proton = m_particleTable->FindParticle("proton");
  m_neutron = m_particleTable->FindParticle("neutron");

  G4String line;

  std::ifstream processListStream (fileName.c_str());

  while(getline(processListStream,line)){
    std::vector<G4String> tokens;

    //Getting a line
    m_readAndParse(line,tokens,"#");
    
    //Important info
    G4String incident = tokens[0];
    G4ParticleDefinition* incidentDef = m_particleTable->FindParticle(incident);
    G4int incidentPDG = incidentDef->GetPDGEncoding();
    m_knownParticles[incidentDef]=true;
    //    G4cout<<"Incident: "<<incident<<G4endl;
    G4String target = tokens[1];
    //    G4cout<<"Target: "<<target<<G4endl;
    
    // Making a ReactionProduct
    ReactionProduct prod;
    for (size_t i = 2; i != tokens.size();i++){
      G4String part = tokens[i];
      if (m_particleTable->contains(part))
    {
      prod.push_back(m_particleTable->FindParticle(part)->GetPDGEncoding());
    } else {
      G4Exception("HadronicProcessHelper", "UnkownParticle", FatalException,
              "Initialization: The reaction product list contained an unknown particle");
    }
    }
    if (target == "proton")
      {
    m_protonReactionMap[incidentPDG].push_back(prod);
      } else if (target == "neutron") {
    m_neutronReactionMap[incidentPDG].push_back(prod);
      } else {
    G4Exception("HadronicProcessHelper", "IllegalTarget", FatalException,
            "Initialization: The reaction product list contained an illegal target particle");
      }
   }

  processListStream.close();

  m_checkFraction = 0;
  m_n22 = 0;
  m_n23 = 0;
}

HadronicProcessHelper::HadronicProcessHelper

(

const HadronicProcessHelper &

)

Member Function Documentation

G4bool HadronicProcessHelper::applicabilityTester

(

const G4ParticleDefinition &

particle

)

Definition at line 71 of file HadronicProcessHelper.cc.

Referenced by instance().

                                                                                  {
  const G4ParticleDefinition* aP = &aPart; 
  if (m_knownParticles[aP]) return true;
  return false;
}

ReactionProduct HadronicProcessHelper::finalState ( const G4Track &  track, G4ParticleDefinition *&  target  )

inline

Definition at line 39 of file HadronicProcessHelper.h.

References edmPickEvents::target.

  {
   return finalState(track.GetDynamicParticle(),track.GetMaterial(),  target);
  }

HadronicProcessHelper::ReactionProduct HadronicProcessHelper::finalState	(	const G4DynamicParticle *	particle,
		const G4Material *	material,
		G4ParticleDefinition *&	target
	)

Definition at line 107 of file HadronicProcessHelper.cc.

References ecalTB2006H4_GenSimDigiReco_cfg::G4cerr, mps_fire::i, and findQualityFiles::size.

                                                             {

//  const G4DynamicParticle* incidentDynamicParticle = track.GetDynamicParticle();

  //-----------------------------------------------
  // Choose n / p as target
  // and get ReactionProductList pointer
  //-----------------------------------------------
  ReactionMap* m_reactionMap;
  //G4Material* material = track.GetMaterial();
  const G4ElementVector* elementVector = material->GetElementVector() ;
  const G4double* numberOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();

  G4double numberOfProtons=0;
  G4double numberOfNucleons=0;

  //Loop on elements 
  for ( size_t elm=0 ; elm < material->GetNumberOfElements() ; elm++ )
    {
      //Summing number of protons per unit volume
      numberOfProtons += numberOfAtomsPerVolume[elm]*(*elementVector)[elm]->GetZ();
      //Summing nucleons (not neutrons)
      numberOfNucleons += numberOfAtomsPerVolume[elm]*(*elementVector)[elm]->GetN();
    }
  
  //random decision of the target
  if(G4UniformRand()<numberOfProtons/numberOfNucleons)
    { //target is a proton
      m_reactionMap = &m_protonReactionMap;
      target = m_proton;
    } else { //target is a neutron
      m_reactionMap = &m_neutronReactionMap;
      target = m_neutron;
    }
  
  const G4int incidentPDG = incidentDynamicParticle->GetDefinition()->GetPDGEncoding();

  //Making a pointer directly to the ReactionProductList we are looking at. Makes life easier :-)
  ReactionProductList*  reactionProductList = &((*m_reactionMap)[incidentPDG]);

  //-----------------------------------------------
  // Count processes
  // kinematic check
  // compute number of 2 -> 2 and 2 -> 3 processes
  //-----------------------------------------------

  G4int good22 = 0; //Number of 2 -> 2 processes that are possible
  G4int good23 = 0; //Number of 2 -> 3 processes that are possible
  
  //This is the list to be populated
  ReactionProductList goodReactionProductList;

  for (ReactionProductList::iterator prod_it = reactionProductList->begin();
       prod_it != reactionProductList->end();
       prod_it++){
    G4int secondaries = prod_it->size();
    // If the reaction is not possible we will not consider it
    if(m_reactionIsPossible(*prod_it,incidentDynamicParticle,target)){
      // The reaction is possible. Let's store and count it
      goodReactionProductList.push_back(*prod_it);
      if (secondaries == 2){
    good22++;
      } else if (secondaries ==3) {
    good23++;
      } else {
    G4cerr << "ReactionProduct has unsupported number of secondaries: "<<secondaries<<G4endl;
      }
    } /*else {
      G4cout<<"There was an impossible process"<<G4endl;
      }*/
  }
  //  G4cout<<"The size of the ReactionProductList is: "<<theReactionProductList.size()<<G4endl;

  if (goodReactionProductList.empty()) G4Exception("HadronicProcessHelper", "NoProcessPossible", FatalException,
                            "GetFinalState: No process could be selected from the given list.");
  // Fill a probability map. Remember total probability
  // 2->2 is 0.15*1/n_22 2->3 uses phase space
  G4double prob22 = 0.15;
  G4double prob23 = 1-prob22; // :-)

  std::vector<G4double> probabilities;
  std::vector<G4bool> twoToThreeFlag;
  
  G4double cumulatedProbability = 0;

  // To each ReactionProduct we assign a cumulated probability and a flag
  // discerning between 2 -> 2 and 2 -> 3
  size_t numberOfReactions = goodReactionProductList.size();
  for (size_t i = 0; i != numberOfReactions; i++){
    if (goodReactionProductList[i].size() == 2) {
      cumulatedProbability += prob22/good22;
      twoToThreeFlag.push_back(false);
    } else {
      cumulatedProbability += prob23/good23;
      twoToThreeFlag.push_back(true);
    }
    probabilities.push_back(cumulatedProbability);
  }

  //Normalising probabilities to 1
  for (std::vector<G4double>::iterator it = probabilities.begin();
       it != probabilities.end();
       it++)
    {
      *it /= cumulatedProbability;
    }

  // Choosing ReactionProduct
  G4bool selected = false;
  G4int tries = 0;
  //  ReactionProductList::iterator prod_it;

  //Keep looping over the list until we have a choice, or until we have tried 100 times  
  size_t i;
  while(!selected && tries < 100){
    i=0;
    G4double dice = G4UniformRand();
    //select the process using the dice
    while(dice>probabilities[i] && i<numberOfReactions)  i++;

    if(twoToThreeFlag[i]) {
      // 2 -> 3 processes require a phase space lookup
      if (m_phaseSpace(goodReactionProductList[i],incidentDynamicParticle,target)>G4UniformRand()) selected = true;
    } else {
      // 2 -> 2 processes are chosen immediately
      selected = true;
    }
    tries++;
  }
  if(tries>=100) G4cerr<<"Could not select process!!!!"<<G4endl;

  
  //Debugging stuff
  //Updating checkfraction:
  if (goodReactionProductList[i].size()==2) {
    m_n22++;
  } else {
    m_n23++;
  }
  m_checkFraction = (1.0*m_n22)/(m_n22+m_n23);
  
  //return the selected productList
  return goodReactionProductList[i];
}

G4double HadronicProcessHelper::inclusiveCrossSection	(	const G4DynamicParticle *	particle,
		const G4Element *	element
	)

Definition at line 77 of file HadronicProcessHelper.cc.

References funct::pow(), CustomPDGParser::s_containedQuarks(), and CustomPDGParser::s_isRGlueball().

Referenced by instance().

                                                    {

  //We really do need a dedicated class to handle the cross sections. They might not always be constant

  G4int pdgCode = particle->GetDefinition()->GetPDGEncoding();

  //24mb for gluino-balls
  if(CustomPDGParser::s_isRGlueball(pdgCode)) return 24 * millibarn * element->GetN();
  
  //get quark vector
  std::vector<G4int> quarks=CustomPDGParser::s_containedQuarks(pdgCode);
  
  G4double totalNucleonCrossSection = 0;

  for (std::vector<G4int>::iterator it = quarks.begin();
       it != quarks.end();
       it++)
    {
      // 12mb for each 'up' or 'down'
      if (*it == 1 || *it == 2) totalNucleonCrossSection += 12 * millibarn;
      //  6mb for each 'strange'
      if (*it == 3) totalNucleonCrossSection += 6 * millibarn;
    }
  
  //Convert to xsec per nucleus
  //  return totalNucleonCrossSection * element->GetN();
  return totalNucleonCrossSection * pow(element->GetN(),0.7)*1.25;// * 0.523598775598299;
}

HadronicProcessHelper* HadronicProcessHelper::instance ( )

inline

Definition at line 31 of file HadronicProcessHelper.h.

References applicabilityTester(), and inclusiveCrossSection().

{return this;};

G4double HadronicProcessHelper::m_phaseSpace ( const ReactionProduct &  aReaction, const G4DynamicParticle *  aDynamicParticle, G4ParticleDefinition *  target  )

private

Definition at line 278 of file HadronicProcessHelper.cc.

References GeV, funct::pow(), and mathSSE::sqrt().

                                                                            {
  G4double qValue = m_reactionProductMass(aReaction,aDynamicParticle,target);
  G4double phi = sqrt(1+qValue/(2*0.139*GeV))*pow(qValue/(1.1*GeV),3./2.);
  return (phi/(1+phi));
}

G4bool HadronicProcessHelper::m_reactionIsPossible ( const ReactionProduct &  aReaction, const G4DynamicParticle *  aDynamicParticle, G4ParticleDefinition *  target  )

private

Definition at line 272 of file HadronicProcessHelper.cc.

                                                                            {
  if (m_reactionProductMass(aReaction,aDynamicParticle,target)>0) return true;
  return false;
}

G4double HadronicProcessHelper::m_reactionProductMass ( const ReactionProduct &  aReaction, const G4DynamicParticle *  aDynamicParticle, G4ParticleDefinition *  target  )

private

Definition at line 253 of file HadronicProcessHelper.cc.

References mathSSE::sqrt().

                                                                                {
  // Incident energy:
  G4double incidentEnergy = incidentDynamicParticle->GetTotalEnergy();
  G4double m_1 = incidentDynamicParticle->GetDefinition()->GetPDGMass();
  G4double m_2 = target->GetPDGMass();
  //square root of "s"
  G4double sqrtS = sqrt(m_1*m_1 + m_2*(m_2 + 2 * incidentEnergy));
  // Sum of rest masses after reaction:
  G4double productsMass = 0;
  //Loop on reaction producs
  for (ReactionProduct::const_iterator r_it = reactionProd.begin(); r_it !=reactionProd.end(); r_it++){
    //Sum the masses of the products
    productsMass += m_particleTable->FindParticle(*r_it)->GetPDGMass();
  }
  //the result is square root of "s" minus the masses of the products
  return sqrtS - productsMass;
}

void HadronicProcessHelper::m_readAndParse ( const G4String &  str, std::vector< G4String > &  tokens, const G4String &  delimiters = " "  )

private

Definition at line 285 of file HadronicProcessHelper.cc.

References groupFilesInBlocks::temp.

{
  // Skip delimiters at beginning.
  G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
  // Find first "non-delimiter".
  G4String::size_type pos     = str.find_first_of(delimiters, lastPos);
  
  while (G4String::npos != pos || G4String::npos != lastPos)
    {
      //Skipping leading / trailing whitespaces
      G4String temp = str.substr(lastPos, pos - lastPos);
      while(temp.c_str()[0] == ' ') temp.erase(0,1);
      while(temp[temp.size()-1] == ' ') temp.erase(temp.size()-1,1);
      // Found a token, add it to the vector.
      tokens.push_back(temp);
      // Skip delimiters.  Note the "not_of"
      lastPos = str.find_first_not_of(delimiters, pos);
      // Find next "non-delimiter"
      pos = str.find_first_of(delimiters, lastPos);
    }
}

HadronicProcessHelper& HadronicProcessHelper::operator=

(

const HadronicProcessHelper &

)

Member Data Documentation

G4double HadronicProcessHelper::m_checkFraction

private

Definition at line 78 of file HadronicProcessHelper.h.

std::map<const G4ParticleDefinition*,G4bool> HadronicProcessHelper::m_knownParticles

private

Definition at line 67 of file HadronicProcessHelper.h.

G4int HadronicProcessHelper::m_n22

private

Definition at line 79 of file HadronicProcessHelper.h.

G4int HadronicProcessHelper::m_n23

private

Definition at line 80 of file HadronicProcessHelper.h.

G4ParticleDefinition* HadronicProcessHelper::m_neutron

private

Definition at line 49 of file HadronicProcessHelper.h.

ReactionMap HadronicProcessHelper::m_neutronReactionMap

private

Definition at line 73 of file HadronicProcessHelper.h.

G4ParticleTable* HadronicProcessHelper::m_particleTable

private

Definition at line 75 of file HadronicProcessHelper.h.

G4ParticleDefinition* HadronicProcessHelper::m_proton

private

Definition at line 48 of file HadronicProcessHelper.h.

ReactionMap HadronicProcessHelper::m_protonReactionMap

private

Definition at line 70 of file HadronicProcessHelper.h.