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 Mak 28 // 29 // Author: Michael Kelsey (SLAC) <kelsey@slac 30 31 #include "G4HadPhaseSpaceNBodyAsai.hh" 32 #include "G4LorentzVector.hh" 33 #include "G4PhysicalConstants.hh" 34 #include "G4SystemOfUnits.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 // This wraps the existing #define in a true 46 G4double uniformRand() { return G4UniformRan 47 } 48 49 50 void G4HadPhaseSpaceNBodyAsai:: 51 GenerateMultiBody(G4double initialMass, 52 const std::vector<G4double>& masses, 53 std::vector<G4LorentzVector>& finalState 54 if (GetVerboseLevel()) G4cout << GetName() < 55 56 finalState.clear(); 57 58 //daughters' mass 59 G4int numberOfDaughters = (G4int)masses.size 60 G4double sumofmasses = 61 std::accumulate(masses.begin(), masses.end 62 63 //Calculate daughter momentum 64 std::vector<G4double> daughtermomentum(numbe 65 std::vector<G4double> sm(numberOfDaughters); 66 G4double tmas; 67 G4double weight = 1.0; 68 G4int numberOfTry = 0; 69 G4int i; 70 71 std::vector<G4double> rd(numberOfDaughters); 72 do { 73 //Generate random number in descending ord 74 rd[0] = 1.0; 75 std::generate(rd.begin()+1, rd.end(), unif 76 std::sort(rd.begin(), rd.end(), std::great 77 78 if (GetVerboseLevel()>1) PrintVector(rd,"r 79 80 //calcurate virtual mass 81 tmas = initialMass - sumofmasses; 82 G4double temp = sumofmasses; 83 for(i =0; i < numberOfDaughters; i++) { 84 sm[i] = rd[i]*tmas + temp; 85 temp -= masses[i]; 86 if (GetVerboseLevel()>1) { 87 G4cout << i << " random number:" << rd 88 << " virtual mass:" << sm[i]/GeV << " 89 } 90 } 91 92 //Calculate daughter momentum 93 weight = 1.0; 94 i = numberOfDaughters-1; 95 daughtermomentum[i] = TwoBodyMomentum(sm[i 96 if (GetVerboseLevel()>1) { 97 G4cout << " daughter " << i << ": moment 98 << daughtermomentum[i]/GeV << " GeV/c" 99 } 100 for(i =numberOfDaughters-2; i>=0; i--) { 101 // calculate 102 daughtermomentum[i] = TwoBodyMomentum(sm 103 if(daughtermomentum[i] < 0.0) { 104 // !!! illegal momentum !!! 105 if (GetVerboseLevel()>0) { 106 G4cout << "G4HadPhaseSpaceNBodyAsai: 107 << " can not calculate daughter momentum 108 << "\n initialMass " << initialMass/GeV < 109 << "\n daughter " << i << ": mass " 110 << masses[i]/GeV << " GeV/c2; momentum " 111 << daughtermomentum[i]/GeV << " GeV/c" << 112 } 113 return; // Error detection 114 } 115 116 // calculate weight of this events 117 weight *= daughtermomentum[i]/sm[i]; 118 if (GetVerboseLevel()>1) { 119 G4cout << " daughter " << i << ": momentum " 120 << daughtermomentum[i]/GeV << " GeV/c 121 } 122 } 123 if (GetVerboseLevel()>1) { 124 G4cout << " weight: " << weight <<G4endl 125 } 126 127 // exit if number of Try exceeds 100 128 if (numberOfTry++ > 100) { 129 if (GetVerboseLevel()>0) { 130 G4cout << "G4HadPhaseSpaceNBodyAsai::G 131 << " can not determine Decay Kinemati 132 } 133 return; // Error detection 134 } 135 } while (weight > G4UniformRand()); /* Loop 136 137 if (GetVerboseLevel()>1) { 138 G4cout << "Start calculation of daughter 139 } 140 141 G4double beta; 142 143 finalState.resize(numberOfDaughters); 144 145 i = numberOfDaughters-2; 146 147 G4ThreeVector direction = UniformVector(daug 148 149 finalState[i].setVectM(direction, masses[i]) 150 finalState[i+1].setVectM(-direction, masses[ 151 152 for (i = numberOfDaughters-3; i >= 0; i--) 153 direction = UniformVector(); 154 155 //create daughter particle 156 finalState[i].setVectM(-daughtermomentum[i 157 158 // boost already created particles 159 beta = daughtermomentum[i]; 160 beta /= std::sqrt(beta*beta + sm[i+1]*sm[i 161 for (G4int j = i+1; j<numberOfDaughters; j 162 finalState[j].boost(beta*direction); 163 } 164 } 165 } 166