Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4AngularDistributionPP.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 #include "G4AngularDistributionPP.hh"
 27 #include "G4PhysicalConstants.hh"
 28 #include "G4SystemOfUnits.hh"
 29 #include "Randomize.hh"
 30 #include "G4ios.hh"
 31 
 32 // Initialization of static data arrays:
 33 #include "G4AngularDistributionPPData.hh"
 34 #include "Randomize.hh"
 35 
 36 
 37 G4double G4AngularDistributionPP::CosTheta(G4double S, G4double m_1, G4double m_2) const
 38 {
 39     G4int verboseLevel=1;
 40 
 41     G4double ek= ((S - sqr(m_1) -sqr(m_2) )/(2*m_1) - m_1 )/GeV   ;    // kinetic energy in GeV
 42 
 43     // Find energy bin
 44 
 45     G4int je1 = 0;
 46     G4int je2 = NENERGY - 1;
 47    G4int iterationsLeft=2*NENERGY +1;
 48     do {
 49       G4int midBin = (je1 + je2)/2;
 50       if (ek < elab[midBin])
 51         je2 = midBin;
 52       else
 53         je1 = midBin;
 54     } while ( (je2 - je1) > 1 && --iterationsLeft > 0 );    /* Loop checking, 30-Oct-2015, G.Folger */
 55    if ( iterationsLeft <= 0 ) {
 56        G4Exception("G4AngularDistributionPP", "im_r_matrix020", FatalException,
 57              "Problem with energy bin (elab) data");
 58 
 59    }
 60     //    G4int j;
 61     //std::abs(ek-elab[je1]) < std::abs(ek-elab[je2]) ? j = je1 : j = je2;
 62     G4double delab = elab[je2] - elab[je1];
 63 
 64     // Sample the angle
 65 
 66     G4double sample = G4UniformRand();
 67     G4int ke1 = 0;
 68     G4int ke2 = NANGLE - 1;
 69     G4double dsig = sig[je2][0] - sig[je1][0];
 70     G4double rc = dsig/delab;
 71     G4double b = sig[je1][0] - rc*elab[je1];
 72     G4double sigint1 = rc*ek + b;
 73     G4double sigint2 = 0.;
 74 
 75     if (verboseLevel > 1) G4cout << "sample=" << sample << G4endl
 76                                  << ek << " " << ke1 << " " << ke2 << " "
 77                                  << sigint1 << " " << sigint2 << G4endl;
 78 
 79     iterationsLeft= 2*NANGLE +1;
 80      do {
 81       G4int midBin = (ke1 + ke2)/2;
 82       dsig = sig[je2][midBin] - sig[je1][midBin];
 83       rc = dsig/delab;
 84       b = sig[je1][midBin] - rc*elab[je1];
 85       G4double sigint = rc*ek + b;
 86       if (sample < sigint) {
 87         ke2 = midBin;
 88         sigint2 = sigint;
 89       }
 90       else {
 91         ke1 = midBin;
 92         sigint1 = sigint;
 93       }
 94       if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " "
 95                                   << sigint1 << " " << sigint2 << G4endl;
 96     } while ( (ke2 - ke1) > 1 && --iterationsLeft > 0);    /* Loop checking, 30-Oct-2015, G.Folger */
 97    if ( iterationsLeft <= 0 ) {
 98        G4Exception("G4AngularDistributionPP", "im_r_matrix021", FatalException,
 99              "Problem with angular distribution (sig) data");
100    }
101 
102     // sigint1 and sigint2 should be recoverable from above loop
103 
104     //    G4double dsig = sig[je2][ke1] - sig[je1][ke1];
105     //    G4double rc = dsig/delab;
106     //    G4double b = sig[je1][ke1] - rc*elab[je1];
107     //    G4double sigint1 = rc*ek + b;
108 
109     //    G4double dsig = sig[je2][ke2] - sig[je1][ke2];
110     //    G4double rc = dsig/delab;
111     //    G4double b = sig[je1][ke2] - rc*elab[je1];
112     //    G4double sigint2 = rc*ek + b;
113 
114     dsig = sigint2 - sigint1;
115     rc = 1./dsig;
116     b = ke1 - rc*sigint1;
117     G4double kint = rc*sample + b;
118     G4double theta = (0.5 + kint)*pi/180.;
119 
120     //    G4int k;
121     //std::abs(sample-sig[j][ke1]) < std::abs(sample-sig[j][ke2]) ? k = ke1 : k = ke2;
122     //    G4double theta = (0.5 + k)*pi/180.;
123 
124     if (verboseLevel > 1) {
125       G4cout << "   energy bin " << je1 << " energy=" << elab[je1] << G4endl;
126       G4cout << "   angle bin " << kint << " angle=" << theta/degree << G4endl;
127     }
128     G4double costh= std::cos(theta);
129 //    G4cout << "scattering angle hpw "<<costh<<" "<<s<<" "<<m1<<" "<<m2<<G4endl;
130     return costh;
131 }
132 
133 G4double G4AngularDistributionPP::Phi() const
134 {
135     return twopi * G4UniformRand();
136 }
137