Geant4  10.02.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  const G4int maxNumberOfLoops = 10000;
70  nTrials = 0;
71  do { // Apply accept/reject to get distribution
72  ++nTrials;
74  FillEnergySteps(initialMass, masses);
75  } while ( (!AcceptEvent()) && nTrials < maxNumberOfLoops ); /* Loop checking, 02.11.2015, A.Ribon */
76  if ( nTrials >= maxNumberOfLoops ) {
78  ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
79  G4Exception( " G4HadPhaseSpaceGenbod::GenerateMultiBody ", "HAD_GENBOD_001", FatalException, ed );
80  }
81  GenerateMomenta(masses, finalState);
82 }
83 
85 Initialize(G4double initialMass, const std::vector<G4double>& masses) {
86  if (GetVerboseLevel()>1) G4cout << GetName() << "::Initialize" << G4endl;
87 
88  nFinal = masses.size();
89  msum.resize(nFinal, 0.); // Initialize buffers for filling
90  msq.resize(nFinal, 0.);
91 
92  std::partial_sum(masses.begin(), masses.end(), msum.begin());
93  std::transform(masses.begin(), masses.end(), masses.begin(), msq.begin(),
94  std::multiplies<G4double>());
95  totalMass = msum.back();
96  massExcess = initialMass - totalMass;
97 
98  if (GetVerboseLevel()>2) {
99  PrintVector(msum, "msum", G4cout);
100  PrintVector(msq, "msq", G4cout);
101  G4cout << " totalMass " << totalMass << " massExcess " << massExcess
102  << G4endl;
103  }
104 
105  ComputeWeightScale(masses);
106 }
107 
108 
109 // Generate ordered list of random numbers
110 
112  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillRandomBuffer" << G4endl;
113 
114  rndm.resize(nFinal-2,0.); // Final states generated in sorted order
115  std::generate(rndm.begin(), rndm.end(), uniformRand);
116  std::sort(rndm.begin(), rndm.end());
117  if (GetVerboseLevel()>2) PrintVector(rndm, "rndm", G4cout);
118 }
119 
120 
121 // Final state effective masses, min to max
122 
123 void
125  const std::vector<G4double>& masses) {
126  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillEnergySteps" << G4endl;
127 
128  meff.clear();
129  pd.clear();
130 
131  meff.push_back(masses[0]);
132  for (size_t i=1; i<nFinal-1; i++) {
133  meff.push_back(rndm[i-1]*massExcess + msum[i]);
134  pd.push_back(TwoBodyMomentum(meff[i], meff[i-1], masses[i]));
135  }
136  meff.push_back(initialMass);
137  pd.push_back(TwoBodyMomentum(meff[nFinal-1], meff[nFinal-2], masses[nFinal-1]));
138 
139  if (GetVerboseLevel()>2) {
140  PrintVector(meff,"meff",G4cout);
141  PrintVector(pd,"pd",G4cout);
142  }
143 }
144 
145 
146 // Maximum possible weight for final state (used with accept/reject)
147 
148 void
149 G4HadPhaseSpaceGenbod::ComputeWeightScale(const std::vector<G4double>& masses) {
150  if (GetVerboseLevel()>1)
151  G4cout << GetName() << "::ComputeWeightScale" << G4endl;
152 
153  weightMax = 1.;
154  for (size_t i=1; i<nFinal; i++) {
155  weightMax *= TwoBodyMomentum(massExcess+msum[i], msum[i-1], masses[i]);
156  }
157 
158  if (GetVerboseLevel()>2) G4cout << " weightMax = " << weightMax << G4endl;
159 }
160 
161 
162 // Event weight computed as either constant or Fermi-dependent cross-section
163 
165  if (GetVerboseLevel()>1) G4cout << GetName() << "::ComputeWeight" << G4endl;
166 
167  return (std::accumulate(pd.begin(), pd.end(), 1./weightMax,
168  std::multiplies<G4double>()));
169 }
170 
172  if (GetVerboseLevel()>1)
173  G4cout << GetName() << "::AcceptEvent? " << nTrials << G4endl;
174 
175  return (G4UniformRand() <= ComputeWeight());
176 }
177 
178 
179 // Final state momentum vectors in CMS system, using Raubold-Lynch method
180 
182 GenerateMomenta(const std::vector<G4double>& masses,
183  std::vector<G4LorentzVector>& finalState) {
184  if (GetVerboseLevel()>1) G4cout << GetName() << "::GenerateMomenta" << G4endl;
185 
186  finalState.resize(nFinal); // Preallocate vectors for convenience below
187 
188  for (size_t i=0; i<nFinal; i++) {
189  AccumulateFinalState(i, masses, finalState);
190  if (GetVerboseLevel()>2)
191  G4cout << " finalState[" << i << "] " << finalState[i] << G4endl;
192  }
193 }
194 
195 // Process final state daughters up to current index
196 
199  const std::vector<G4double>& masses,
200  std::vector<G4LorentzVector>& finalState) {
201  if (GetVerboseLevel()>2)
202  G4cout << GetName() << "::AccumulateFinalState " << i << G4endl;
203 
204  if (i==0) { // First final state particle left alone
205  finalState[i].setVectM(G4ThreeVector(0.,pd[i],0.),masses[i]);
206  return;
207  }
208 
209  finalState[i].setVectM(G4ThreeVector(0.,-pd[i-1],0.),masses[i]);
210  G4double phi = G4UniformRand() * twopi;
211  G4double theta = std::acos(2.*G4UniformRand() - 1.);
212 
213  if (GetVerboseLevel() > 2) {
214  G4cout << " initialized Py " << -pd[i-1] << " phi " << phi
215  << " theta " << theta << G4endl;
216  }
217 
218  G4double esys=0.,beta=0.,gamma=1.;
219  if (i < nFinal-1) { // Do not boost final particle
220  esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[i]);
221  beta = pd[i] / esys;
222  gamma = esys / meff[i];
223 
224  if (GetVerboseLevel()>2)
225  G4cout << " esys " << esys << " beta " << beta << " gamma " << gamma
226  << G4endl;
227  }
228 
229  for (size_t j=0; j<=i; j++) { // Accumulate rotations
230  finalState[j].rotateZ(theta).rotateY(phi);
231  finalState[j].setY(gamma*(finalState[j].y() + beta*finalState[j].e()));
232  if (GetVerboseLevel()>2) G4cout << " j " << j << " " << finalState[j] << G4endl;
233  }
234 }
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
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:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
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