Geant4  10.01.p01
G4HadPhaseSpaceGenbod.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id$
27 //
28 // Multibody "phase space" generator using GENBOD (CERNLIB W515) method.
29 //
30 // Author: Michael Kelsey (SLAC) <kelsey@slac.stanford.edu>
31 
32 #include "G4HadPhaseSpaceGenbod.hh"
33 #include "G4LorentzVector.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4ThreeVector.hh"
36 #include "Randomize.hh"
37 #include <algorithm>
38 #include <functional>
39 #include <iterator>
40 #include <numeric>
41 #include <vector>
42 
43 
44 namespace {
45  // Wrap #define in a true function, for passing to std::fill
46  G4double uniformRand() { return G4UniformRand(); }
47 }
48 
49 
50 // Constructor initializes everything to zero
51 
53  : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpaceGenbod",verbose),
54  nFinal(0), totalMass(0.), massExcess(0.), weightMax(0.), nTrials(0) {;}
55 
56 
57 // C++ re-implementation of GENBOD.F (Raubold-Lynch method)
58 
61  const std::vector<G4double>& masses,
62  std::vector<G4LorentzVector>& finalState) {
63  if (GetVerboseLevel()) G4cout << GetName() << "::GenerateMultiBody" << G4endl;
64 
65  finalState.clear();
66 
67  Initialize(initialMass, masses);
68 
69  nTrials = 0;
70  do { // Apply accept/reject to get distribution
71  ++nTrials;
73  FillEnergySteps(initialMass, masses);
74  } while (!AcceptEvent()); // FIXME: Do we need a limit on nTrials?
75 
76  GenerateMomenta(masses, finalState);
77 }
78 
80 Initialize(G4double initialMass, const std::vector<G4double>& masses) {
81  if (GetVerboseLevel()>1) G4cout << GetName() << "::Initialize" << G4endl;
82 
83  nFinal = masses.size();
84  msum.resize(nFinal, 0.); // Initialize buffers for filling
85  msq.resize(nFinal, 0.);
86 
87  std::partial_sum(masses.begin(), masses.end(), msum.begin());
88  std::transform(masses.begin(), masses.end(), masses.begin(), msq.begin(),
89  std::multiplies<G4double>());
90  totalMass = msum.back();
91  massExcess = initialMass - totalMass;
92 
93  if (GetVerboseLevel()>2) {
94  PrintVector(msum, "msum", G4cout);
95  PrintVector(msq, "msq", G4cout);
96  G4cout << " totalMass " << totalMass << " massExcess " << massExcess
97  << G4endl;
98  }
99 
100  ComputeWeightScale(masses);
101 }
102 
103 
104 // Generate ordered list of random numbers
105 
107  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillRandomBuffer" << G4endl;
108 
109  rndm.resize(nFinal-2,0.); // Final states generated in sorted order
110  std::generate(rndm.begin(), rndm.end(), uniformRand);
111  std::sort(rndm.begin(), rndm.end());
112  if (GetVerboseLevel()>2) PrintVector(rndm, "rndm", G4cout);
113 }
114 
115 
116 // Final state effective masses, min to max
117 
118 void
120  const std::vector<G4double>& masses) {
121  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillEnergySteps" << G4endl;
122 
123  meff.clear();
124  pd.clear();
125 
126  meff.push_back(masses[0]);
127  for (size_t i=1; i<nFinal-1; i++) {
128  meff.push_back(rndm[i-1]*massExcess + msum[i]);
129  pd.push_back(TwoBodyMomentum(meff[i], meff[i-1], masses[i]));
130  }
131  meff.push_back(initialMass);
132  pd.push_back(TwoBodyMomentum(meff[nFinal-1], meff[nFinal-2], masses[nFinal-1]));
133 
134  if (GetVerboseLevel()>2) {
135  PrintVector(meff,"meff",G4cout);
136  PrintVector(pd,"pd",G4cout);
137  }
138 }
139 
140 
141 // Maximum possible weight for final state (used with accept/reject)
142 
143 void
144 G4HadPhaseSpaceGenbod::ComputeWeightScale(const std::vector<G4double>& masses) {
145  if (GetVerboseLevel()>1)
146  G4cout << GetName() << "::ComputeWeightScale" << G4endl;
147 
148  weightMax = 1.;
149  for (size_t i=1; i<nFinal; i++) {
150  weightMax *= TwoBodyMomentum(massExcess+msum[i], msum[i-1], masses[i]);
151  }
152 
153  if (GetVerboseLevel()>2) G4cout << " weightMax = " << weightMax << G4endl;
154 }
155 
156 
157 // Event weight computed as either constant or Fermi-dependent cross-section
158 
160  if (GetVerboseLevel()>1) G4cout << GetName() << "::ComputeWeight" << G4endl;
161 
162  return (std::accumulate(pd.begin(), pd.end(), 1./weightMax,
163  std::multiplies<G4double>()));
164 }
165 
167  if (GetVerboseLevel()>1)
168  G4cout << GetName() << "::AcceptEvent? " << nTrials << G4endl;
169 
170  return (G4UniformRand() <= ComputeWeight());
171 }
172 
173 
174 // Final state momentum vectors in CMS system, using Raubold-Lynch method
175 
177 GenerateMomenta(const std::vector<G4double>& masses,
178  std::vector<G4LorentzVector>& finalState) {
179  if (GetVerboseLevel()>1) G4cout << GetName() << "::GenerateMomenta" << G4endl;
180 
181  finalState.resize(nFinal); // Preallocate vectors for convenience below
182 
183  for (size_t i=0; i<nFinal; i++) {
184  AccumulateFinalState(i, masses, finalState);
185  if (GetVerboseLevel()>2)
186  G4cout << " finalState[" << i << "] " << finalState[i] << G4endl;
187  }
188 }
189 
190 // Process final state daughters up to current index
191 
194  const std::vector<G4double>& masses,
195  std::vector<G4LorentzVector>& finalState) {
196  if (GetVerboseLevel()>2)
197  G4cout << GetName() << "::AccumulateFinalState " << i << G4endl;
198 
199  if (i==0) { // First final state particle left alone
200  finalState[i].setVectM(G4ThreeVector(0.,pd[i],0.),masses[i]);
201  return;
202  }
203 
204  finalState[i].setVectM(G4ThreeVector(0.,-pd[i-1],0.),masses[i]);
205  G4double phi = G4UniformRand() * twopi;
206  G4double theta = std::acos(2.*G4UniformRand() - 1.);
207 
208  if (GetVerboseLevel() > 2) {
209  G4cout << " initialized Py " << -pd[i-1] << " phi " << phi
210  << " theta " << theta << G4endl;
211  }
212 
213  G4double esys=0.,beta=0.,gamma=1.;
214  if (i < nFinal-1) { // Do not boost final particle
215  esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[i]);
216  beta = pd[i] / esys;
217  gamma = esys / meff[i];
218 
219  if (GetVerboseLevel()>2)
220  G4cout << " esys " << esys << " beta " << beta << " gamma " << gamma
221  << G4endl;
222  }
223 
224  for (size_t j=0; j<=i; j++) { // Accumulate rotations
225  finalState[j].rotateZ(theta).rotateY(phi);
226  finalState[j].setY(gamma*(finalState[j].y() + beta*finalState[j].e()));
227  if (GetVerboseLevel()>2) G4cout << " j " << j << " " << finalState[j] << G4endl;
228  }
229 }
CLHEP::Hep3Vector G4ThreeVector
void FillEnergySteps(G4double initialMass, const std::vector< G4double > &masses)
void GenerateMomenta(const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
void ComputeWeightScale(const std::vector< G4double > &masses)
void AccumulateFinalState(size_t i, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4int GetVerboseLevel() const
G4double ComputeWeight() const
virtual void GenerateMultiBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
int G4int
Definition: G4Types.hh:78
G4HadPhaseSpaceGenbod(G4int verbose=0)
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
std::vector< G4double > msq
const G4String & GetName() const
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
void Initialize(G4double initialMass, const std::vector< G4double > &masses)
std::vector< G4double > pd
std::vector< G4double > rndm
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
std::vector< G4double > msum
#define G4endl
Definition: G4ios.hh:61
std::vector< G4double > meff
double G4double
Definition: G4Types.hh:76