Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // Multibody "phase space" generator using GEN 28 // 29 // Author: Michael Kelsey (SLAC) <kelsey@slac 30 31 #include "G4HadPhaseSpaceGenbod.hh" 32 #include "G4LorentzVector.hh" 33 #include "G4PhysicalConstants.hh" 34 #include "G4ThreeVector.hh" 35 #include "Randomize.hh" 36 #include <algorithm> 37 #include <functional> 38 #include <iterator> 39 #include <numeric> 40 #include <vector> 41 42 43 namespace { 44 // Wrap #define in a true function, for pass 45 G4double uniformRand() { return G4UniformRan 46 } 47 48 49 // Constructor initializes everything to zero 50 51 G4HadPhaseSpaceGenbod::G4HadPhaseSpaceGenbod(G 52 : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpace 53 nFinal(0), totalMass(0.), massExcess(0.), 54 55 56 // C++ re-implementation of GENBOD.F (Raubold- 57 58 void G4HadPhaseSpaceGenbod:: 59 GenerateMultiBody(G4double initialMass, 60 const std::vector<G4double>& masses, 61 std::vector<G4LorentzVector>& finalState 62 if (GetVerboseLevel()) G4cout << GetName() < 63 64 finalState.clear(); 65 66 Initialize(initialMass, masses); 67 68 const G4int maxNumberOfLoops = 10000; 69 nTrials = 0; 70 do { // Apply accept/reject to get di 71 ++nTrials; 72 FillRandomBuffer(); 73 FillEnergySteps(initialMass, masses); 74 } while ( (!AcceptEvent()) && nTrials < maxN 75 if ( nTrials >= maxNumberOfLoops ) { 76 G4ExceptionDescription ed; 77 ed << " Failed sampling after maxNumberOfL 78 G4Exception( " G4HadPhaseSpaceGenbod::Gene 79 } 80 GenerateMomenta(masses, finalState); 81 } 82 83 void G4HadPhaseSpaceGenbod:: 84 Initialize(G4double initialMass, const std::ve 85 if (GetVerboseLevel()>1) G4cout << GetName() 86 87 nFinal = masses.size(); 88 msum.resize(nFinal, 0.); // Initialize bu 89 msq.resize(nFinal, 0.); 90 91 std::partial_sum(masses.begin(), masses.end( 92 std::transform(masses.begin(), masses.end(), 93 std::multiplies<G4double>()); 94 totalMass = msum.back(); 95 massExcess = initialMass - totalMass; 96 97 if (GetVerboseLevel()>2) { 98 PrintVector(msum, "msum", G4cout); 99 PrintVector(msq, "msq", G4cout); 100 G4cout << " totalMass " << totalMass << " 101 << G4endl; 102 } 103 104 ComputeWeightScale(masses); 105 } 106 107 108 // Generate ordered list of random numbers 109 110 void G4HadPhaseSpaceGenbod::FillRandomBuffer() 111 if (GetVerboseLevel()>1) G4cout << GetName() 112 113 rndm.resize(nFinal-2,0.); // Final states ge 114 std::generate(rndm.begin(), rndm.end(), unif 115 std::sort(rndm.begin(), rndm.end()); 116 if (GetVerboseLevel()>2) PrintVector(rndm, " 117 } 118 119 120 // Final state effective masses, min to max 121 122 void 123 G4HadPhaseSpaceGenbod::FillEnergySteps(G4doubl 124 const std::vector<G4double>& ma 125 if (GetVerboseLevel()>1) G4cout << GetName() 126 127 meff.clear(); 128 pd.clear(); 129 130 meff.push_back(masses[0]); 131 for (size_t i=1; i<nFinal-1; i++) { 132 meff.push_back(rndm[i-1]*massExcess + msum 133 pd.push_back(TwoBodyMomentum(meff[i], meff 134 } 135 meff.push_back(initialMass); 136 pd.push_back(TwoBodyMomentum(meff[nFinal-1], 137 138 if (GetVerboseLevel()>2) { 139 PrintVector(meff,"meff",G4cout); 140 PrintVector(pd,"pd",G4cout); 141 } 142 } 143 144 145 // Maximum possible weight for final state (us 146 147 void 148 G4HadPhaseSpaceGenbod::ComputeWeightScale(cons 149 if (GetVerboseLevel()>1) 150 G4cout << GetName() << "::ComputeWeightSca 151 152 weightMax = 1.; 153 for (size_t i=1; i<nFinal; i++) { 154 weightMax *= TwoBodyMomentum(massExcess+ms 155 } 156 157 if (GetVerboseLevel()>2) G4cout << " weightM 158 } 159 160 161 // Event weight computed as either constant or 162 163 G4double G4HadPhaseSpaceGenbod::ComputeWeight( 164 if (GetVerboseLevel()>1) G4cout << GetName() 165 166 return (std::accumulate(pd.begin(), pd.end() 167 std::multiplies<G4double>())); 168 } 169 170 G4bool G4HadPhaseSpaceGenbod::AcceptEvent() co 171 if (GetVerboseLevel()>1) 172 G4cout << GetName() << "::AcceptEvent? " < 173 174 return (G4UniformRand() <= ComputeWeight()); 175 } 176 177 178 // Final state momentum vectors in CMS system, 179 180 void G4HadPhaseSpaceGenbod:: 181 GenerateMomenta(const std::vector<G4double>& m 182 std::vector<G4LorentzVector>& finalState) 183 if (GetVerboseLevel()>1) G4cout << GetName() 184 185 finalState.resize(nFinal); // Preallocate v 186 187 for (size_t i=0; i<nFinal; i++) { 188 AccumulateFinalState(i, masses, finalState 189 if (GetVerboseLevel()>2) 190 G4cout << " finalState[" << i << "] " << 191 } 192 } 193 194 // Process final state daughters up to current 195 196 void G4HadPhaseSpaceGenbod:: 197 AccumulateFinalState(size_t i, 198 const std::vector<G4double>& masses, 199 std::vector<G4LorentzVector>& finalSt 200 if (GetVerboseLevel()>2) 201 G4cout << GetName() << "::AccumulateFinalS 202 203 if (i==0) { // First final state particl 204 finalState[i].setVectM(G4ThreeVector(0.,pd 205 return; 206 } 207 208 finalState[i].setVectM(G4ThreeVector(0.,-pd[ 209 G4double phi = G4UniformRand() * twopi; 210 G4double theta = std::acos(2.*G4UniformRand( 211 212 if (GetVerboseLevel() > 2) { 213 G4cout << " initialized Py " << -pd[i-1] < 214 << " theta " << theta << G4endl; 215 } 216 217 G4double esys=0.,beta=0.,gamma=1.; 218 if (i < nFinal-1) { // Do not boost fina 219 esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[ 220 beta = pd[i] / esys; 221 gamma = esys / meff[i]; 222 223 if (GetVerboseLevel()>2) 224 G4cout << " esys " << esys << " beta " < 225 << G4endl; 226 } 227 228 for (size_t j=0; j<=i; j++) { // Accumulat 229 finalState[j].rotateZ(theta).rotateY(phi); 230 finalState[j].setY(gamma*(finalState[j].y( 231 if (GetVerboseLevel()>2) G4cout << " j " < 232 } 233 } 234