Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/util/src/G4HadPhaseSpaceGenbod.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 GENBOD (CERNLIB W515) method.
 28 //
 29 // Author:  Michael Kelsey (SLAC) <kelsey@slac.stanford.edu>
 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 passing to std::fill
 45   G4double uniformRand() { return G4UniformRand(); }
 46 }
 47 
 48 
 49 // Constructor initializes everything to zero
 50 
 51 G4HadPhaseSpaceGenbod::G4HadPhaseSpaceGenbod(G4int verbose)
 52   : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpaceGenbod",verbose),
 53     nFinal(0), totalMass(0.), massExcess(0.), weightMax(0.), nTrials(0) {;}
 54 
 55 
 56 // C++ re-implementation of GENBOD.F (Raubold-Lynch method)
 57 
 58 void G4HadPhaseSpaceGenbod::
 59 GenerateMultiBody(G4double initialMass,
 60       const std::vector<G4double>& masses,
 61       std::vector<G4LorentzVector>& finalState) {
 62   if (GetVerboseLevel()) G4cout << GetName() << "::GenerateMultiBody" << G4endl;
 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 distribution
 71     ++nTrials;
 72     FillRandomBuffer();
 73     FillEnergySteps(initialMass, masses);
 74   } while ( (!AcceptEvent()) && nTrials < maxNumberOfLoops );  /* Loop checking, 02.11.2015, A.Ribon */
 75   if ( nTrials >= maxNumberOfLoops ) {
 76     G4ExceptionDescription ed;
 77     ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
 78     G4Exception( " G4HadPhaseSpaceGenbod::GenerateMultiBody ", "HAD_GENBOD_001", FatalException, ed );
 79   }
 80   GenerateMomenta(masses, finalState);
 81 }
 82 
 83 void G4HadPhaseSpaceGenbod::
 84 Initialize(G4double initialMass, const std::vector<G4double>& masses) {
 85   if (GetVerboseLevel()>1) G4cout << GetName() << "::Initialize" << G4endl;
 86 
 87   nFinal = masses.size();
 88   msum.resize(nFinal, 0.);    // Initialize buffers for filling
 89   msq.resize(nFinal, 0.);
 90 
 91   std::partial_sum(masses.begin(), masses.end(), msum.begin());
 92   std::transform(masses.begin(), masses.end(), masses.begin(), msq.begin(),
 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 << " massExcess " << massExcess
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() << "::FillRandomBuffer" << G4endl;
112 
113   rndm.resize(nFinal-2,0.); // Final states generated in sorted order
114   std::generate(rndm.begin(), rndm.end(), uniformRand);
115   std::sort(rndm.begin(), rndm.end());
116   if (GetVerboseLevel()>2) PrintVector(rndm, "rndm", G4cout);
117 }
118 
119 
120 // Final state effective masses, min to max
121 
122 void 
123 G4HadPhaseSpaceGenbod::FillEnergySteps(G4double initialMass,
124                const std::vector<G4double>& masses) {
125   if (GetVerboseLevel()>1) G4cout << GetName() << "::FillEnergySteps" << G4endl;
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[i]);
133     pd.push_back(TwoBodyMomentum(meff[i], meff[i-1], masses[i]));
134   }
135   meff.push_back(initialMass);
136   pd.push_back(TwoBodyMomentum(meff[nFinal-1], meff[nFinal-2], masses[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 (used with accept/reject)
146 
147 void 
148 G4HadPhaseSpaceGenbod::ComputeWeightScale(const std::vector<G4double>& masses) {
149   if (GetVerboseLevel()>1) 
150     G4cout << GetName() << "::ComputeWeightScale" << G4endl;
151 
152   weightMax = 1.;
153   for (size_t i=1; i<nFinal; i++) {
154     weightMax *= TwoBodyMomentum(massExcess+msum[i], msum[i-1], masses[i]);
155   }
156 
157   if (GetVerboseLevel()>2) G4cout << " weightMax = " << weightMax << G4endl;
158 }
159 
160 
161 // Event weight computed as either constant or Fermi-dependent cross-section
162 
163 G4double G4HadPhaseSpaceGenbod::ComputeWeight() const {
164   if (GetVerboseLevel()>1) G4cout << GetName() << "::ComputeWeight" << G4endl;
165 
166   return (std::accumulate(pd.begin(), pd.end(), 1./weightMax,
167         std::multiplies<G4double>()));
168 }
169 
170 G4bool G4HadPhaseSpaceGenbod::AcceptEvent() const {
171   if (GetVerboseLevel()>1) 
172     G4cout << GetName() << "::AcceptEvent? " << nTrials << G4endl;
173 
174   return (G4UniformRand() <= ComputeWeight());
175 }
176 
177 
178 // Final state momentum vectors in CMS system, using Raubold-Lynch method
179 
180 void G4HadPhaseSpaceGenbod::
181 GenerateMomenta(const std::vector<G4double>& masses,
182     std::vector<G4LorentzVector>& finalState) {
183   if (GetVerboseLevel()>1) G4cout << GetName() << "::GenerateMomenta" << G4endl;
184 
185   finalState.resize(nFinal);  // Preallocate vectors for convenience below
186 
187   for (size_t i=0; i<nFinal; i++) {
188     AccumulateFinalState(i, masses, finalState);
189     if (GetVerboseLevel()>2)
190       G4cout << " finalState[" << i << "] " << finalState[i] << G4endl;
191   }
192 }
193 
194 // Process final state daughters up to current index
195 
196 void G4HadPhaseSpaceGenbod::
197 AccumulateFinalState(size_t i, 
198          const std::vector<G4double>& masses,
199          std::vector<G4LorentzVector>& finalState) {
200   if (GetVerboseLevel()>2)
201     G4cout << GetName() << "::AccumulateFinalState " << i << G4endl;
202 
203   if (i==0) {     // First final state particle left alone
204     finalState[i].setVectM(G4ThreeVector(0.,pd[i],0.),masses[i]);
205     return;
206   }
207   
208   finalState[i].setVectM(G4ThreeVector(0.,-pd[i-1],0.),masses[i]);
209   G4double phi = G4UniformRand() * twopi;
210   G4double theta = std::acos(2.*G4UniformRand() - 1.);
211 
212   if (GetVerboseLevel() > 2) {
213     G4cout << " initialized Py " << -pd[i-1] << " phi " << phi
214      << " theta " << theta << G4endl;
215   }
216 
217   G4double esys=0.,beta=0.,gamma=1.;
218   if (i < nFinal-1) {     // Do not boost final particle
219     esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[i]);
220     beta = pd[i] / esys;
221     gamma = esys / meff[i];
222 
223     if (GetVerboseLevel()>2)
224       G4cout << " esys " << esys << " beta " << beta << " gamma " << gamma
225        << G4endl;
226   }
227 
228   for (size_t j=0; j<=i; j++) {   // Accumulate rotations
229     finalState[j].rotateZ(theta).rotateY(phi);
230     finalState[j].setY(gamma*(finalState[j].y() + beta*finalState[j].e()));
231     if (GetVerboseLevel()>2) G4cout << " j " << j << " " << finalState[j] << G4endl;
232   }
233 }
234