Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4HadPhaseSpaceNBodyAsai Class Reference

#include <G4HadPhaseSpaceNBodyAsai.hh>

Inheritance diagram for G4HadPhaseSpaceNBodyAsai:
Collaboration diagram for G4HadPhaseSpaceNBodyAsai:

Public Member Functions

 G4HadPhaseSpaceNBodyAsai (G4int verbose=0)
 
virtual ~G4HadPhaseSpaceNBodyAsai ()
 
- Public Member Functions inherited from G4VHadPhaseSpaceAlgorithm
 G4VHadPhaseSpaceAlgorithm (const G4String &algName, G4int verbose=0)
 
virtual ~G4VHadPhaseSpaceAlgorithm ()
 
- Public Member Functions inherited from G4VHadDecayAlgorithm
 G4VHadDecayAlgorithm (const G4String &algName, G4int verbose=0)
 
virtual ~G4VHadDecayAlgorithm ()
 
void Generate (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
 
virtual void SetVerboseLevel (G4int verbose)
 
G4int GetVerboseLevel () const
 
const G4StringGetName () const
 

Protected Member Functions

virtual void GenerateMultiBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
 
- Protected Member Functions inherited from G4VHadPhaseSpaceAlgorithm
virtual void GenerateTwoBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
 
G4ThreeVector UniformVector (G4double mag=1.) const
 
- Protected Member Functions inherited from G4VHadDecayAlgorithm
virtual G4bool IsDecayAllowed (G4double initialMass, const std::vector< G4double > &masses) const
 
G4double TwoBodyMomentum (G4double M0, G4double M1, G4double M2) const
 
G4double UniformTheta () const
 
G4double UniformPhi () const
 
void PrintVector (const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
 

Detailed Description

Definition at line 38 of file G4HadPhaseSpaceNBodyAsai.hh.

Constructor & Destructor Documentation

G4HadPhaseSpaceNBodyAsai::G4HadPhaseSpaceNBodyAsai ( G4int  verbose = 0)
inline

Definition at line 40 of file G4HadPhaseSpaceNBodyAsai.hh.

41  : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpaceNBodyAsai",verbose) {;}
G4VHadPhaseSpaceAlgorithm(const G4String &algName, G4int verbose=0)
virtual G4HadPhaseSpaceNBodyAsai::~G4HadPhaseSpaceNBodyAsai ( )
inlinevirtual

Definition at line 42 of file G4HadPhaseSpaceNBodyAsai.hh.

42 {;}

Member Function Documentation

void G4HadPhaseSpaceNBodyAsai::GenerateMultiBody ( G4double  initialMass,
const std::vector< G4double > &  masses,
std::vector< G4LorentzVector > &  finalState 
)
protectedvirtual

Implements G4VHadDecayAlgorithm.

Definition at line 52 of file G4HadPhaseSpaceNBodyAsai.cc.

54  {
55  if (GetVerboseLevel()) G4cout << GetName() << "::GenerateMultiBody" << G4endl;
56 
57  finalState.clear();
58 
59  //daughters' mass
60  G4int numberOfDaughters = masses.size();
61  G4double sumofmasses =
62  std::accumulate(masses.begin(), masses.end(), 0.);
63 
64  //Calculate daughter momentum
65  std::vector<G4double> daughtermomentum(numberOfDaughters);
66  std::vector<G4double> sm(numberOfDaughters);
67  G4double tmas;
68  G4double weight = 1.0;
69  G4int numberOfTry = 0;
70  G4int i;
71 
72  std::vector<G4double> rd(numberOfDaughters);
73  do {
74  //Generate random number in descending order
75  rd[0] = 1.0;
76  std::generate(rd.begin()+1, rd.end(), uniformRand);
77  std::sort(rd.begin(), rd.end(), std::greater<G4double>());
78 
79  if (GetVerboseLevel()>1) PrintVector(rd,"rd",G4cout);
80 
81  //calcurate virtual mass
82  tmas = initialMass - sumofmasses;
83  G4double temp = sumofmasses;
84  for(i =0; i < numberOfDaughters; i++) {
85  sm[i] = rd[i]*tmas + temp;
86  temp -= masses[i];
87  if (GetVerboseLevel()>1) {
88  G4cout << i << " random number:" << rd[i]
89  << " virtual mass:" << sm[i]/GeV << " GeV/c2" <<G4endl;
90  }
91  }
92 
93  //Calculate daughter momentum
94  weight = 1.0;
95  i = numberOfDaughters-1;
96  daughtermomentum[i] = TwoBodyMomentum(sm[i-1],masses[i-1],sm[i]);
97  if (GetVerboseLevel()>1) {
98  G4cout << " daughter " << i << ": momentum "
99  << daughtermomentum[i]/GeV << " GeV/c" <<G4endl;
100  }
101  for(i =numberOfDaughters-2; i>=0; i--) {
102  // calculate
103  daughtermomentum[i] = TwoBodyMomentum(sm[i],masses[i],sm[i+1]);
104  if(daughtermomentum[i] < 0.0) {
105  // !!! illegal momentum !!!
106  if (GetVerboseLevel()>0) {
107  G4cout << "G4HadPhaseSpaceNBodyAsai::Generate "
108  << " can not calculate daughter momentum "
109  << "\n initialMass " << initialMass/GeV << " GeV/c2"
110  << "\n daughter " << i << ": mass "
111  << masses[i]/GeV << " GeV/c2; momentum "
112  << daughtermomentum[i]/GeV << " GeV/c" << G4endl;
113  }
114  return; // Error detection
115  }
116 
117  // calculate weight of this events
118  weight *= daughtermomentum[i]/sm[i];
119  if (GetVerboseLevel()>1) {
120  G4cout << " daughter " << i << ": momentum "
121  << daughtermomentum[i]/GeV << " GeV/c" <<G4endl;
122  }
123  }
124  if (GetVerboseLevel()>1) {
125  G4cout << " weight: " << weight <<G4endl;
126  }
127 
128  // exit if number of Try exceeds 100
129  if (numberOfTry++ > 100) {
130  if (GetVerboseLevel()>0) {
131  G4cout << "G4HadPhaseSpaceNBodyAsai::Generate "
132  << " can not determine Decay Kinematics " << G4endl;
133  }
134  return; // Error detection
135  }
136  } while (weight > G4UniformRand()); /* Loop checking, 02.11.2015, A.Ribon */
137 
138  if (GetVerboseLevel()>1) {
139  G4cout << "Start calulation of daughters momentum vector "<<G4endl;
140  }
141 
142  G4double beta;
143 
144  finalState.resize(numberOfDaughters);
145 
146  i = numberOfDaughters-2;
147 
148  G4ThreeVector direction = UniformVector(daughtermomentum[i]);
149 
150  finalState[i].setVectM(direction, masses[i]);
151  finalState[i+1].setVectM(-direction, masses[i+1]);
152 
153  for (i = numberOfDaughters-3; i >= 0; i--) {
154  direction = UniformVector();
155 
156  //create daughter particle
157  finalState[i].setVectM(-daughtermomentum[i]*direction, masses[i]);
158 
159  // boost already created particles
160  beta = daughtermomentum[i];
161  beta /= std::sqrt(beta*beta + sm[i+1]*sm[i+1]);
162  for (G4int j = i+1; j<numberOfDaughters; j++) {
163  finalState[j].boost(beta*direction);
164  }
165  }
166 }
G4int GetVerboseLevel() const
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
const G4String & GetName() const
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ThreeVector UniformVector(G4double mag=1.) const

Here is the call graph for this function:


The documentation for this class was generated from the following files: