Geant4 Cross Reference |
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 // 27 // Multibody "phase space" generator using Makoto Asai's NBody method. 28 // 29 // Author: Michael Kelsey (SLAC) <kelsey@slac.stanford.edu> 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 function 46 G4double uniformRand() { return G4UniformRand(); } 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() << "::GenerateMultiBody" << G4endl; 55 56 finalState.clear(); 57 58 //daughters' mass 59 G4int numberOfDaughters = (G4int)masses.size(); 60 G4double sumofmasses = 61 std::accumulate(masses.begin(), masses.end(), 0.); 62 63 //Calculate daughter momentum 64 std::vector<G4double> daughtermomentum(numberOfDaughters); 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 order 74 rd[0] = 1.0; 75 std::generate(rd.begin()+1, rd.end(), uniformRand); 76 std::sort(rd.begin(), rd.end(), std::greater<G4double>()); 77 78 if (GetVerboseLevel()>1) PrintVector(rd,"rd",G4cout); 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[i] 88 << " virtual mass:" << sm[i]/GeV << " GeV/c2" <<G4endl; 89 } 90 } 91 92 //Calculate daughter momentum 93 weight = 1.0; 94 i = numberOfDaughters-1; 95 daughtermomentum[i] = TwoBodyMomentum(sm[i-1],masses[i-1],sm[i]); 96 if (GetVerboseLevel()>1) { 97 G4cout << " daughter " << i << ": momentum " 98 << daughtermomentum[i]/GeV << " GeV/c" <<G4endl; 99 } 100 for(i =numberOfDaughters-2; i>=0; i--) { 101 // calculate 102 daughtermomentum[i] = TwoBodyMomentum(sm[i],masses[i],sm[i+1]); 103 if(daughtermomentum[i] < 0.0) { 104 // !!! illegal momentum !!! 105 if (GetVerboseLevel()>0) { 106 G4cout << "G4HadPhaseSpaceNBodyAsai::Generate " 107 << " can not calculate daughter momentum " 108 << "\n initialMass " << initialMass/GeV << " GeV/c2" 109 << "\n daughter " << i << ": mass " 110 << masses[i]/GeV << " GeV/c2; momentum " 111 << daughtermomentum[i]/GeV << " GeV/c" << G4endl; 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" <<G4endl; 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::Generate " 131 << " can not determine Decay Kinematics " << G4endl; 132 } 133 return; // Error detection 134 } 135 } while (weight > G4UniformRand()); /* Loop checking, 02.11.2015, A.Ribon */ 136 137 if (GetVerboseLevel()>1) { 138 G4cout << "Start calculation of daughters momentum vector "<<G4endl; 139 } 140 141 G4double beta; 142 143 finalState.resize(numberOfDaughters); 144 145 i = numberOfDaughters-2; 146 147 G4ThreeVector direction = UniformVector(daughtermomentum[i]); 148 149 finalState[i].setVectM(direction, masses[i]); 150 finalState[i+1].setVectM(-direction, masses[i+1]); 151 152 for (i = numberOfDaughters-3; i >= 0; i--) { 153 direction = UniformVector(); 154 155 //create daughter particle 156 finalState[i].setVectM(-daughtermomentum[i]*direction, masses[i]); 157 158 // boost already created particles 159 beta = daughtermomentum[i]; 160 beta /= std::sqrt(beta*beta + sm[i+1]*sm[i+1]); 161 for (G4int j = i+1; j<numberOfDaughters; j++) { 162 finalState[j].boost(beta*direction); 163 } 164 } 165 } 166