Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4AngularDistributionNP.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 "G4AngularDistributionNP.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 "G4AngularDistributionNPData.hh"
 34 #include "Randomize.hh"
 35 
 36 
 37 G4double G4AngularDistributionNP::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("G4AngularDistributionNP", "im_r_matrix010", 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     iterationsLeft= 2*NANGLE +1;
 79     do {
 80       G4int midBin = (ke1 + ke2)/2;
 81       dsig = sig[je2][midBin] - sig[je1][midBin];
 82       rc = dsig/delab;
 83       b = sig[je1][midBin] - rc*elab[je1];
 84       G4double sigint = rc*ek + b;
 85       if (sample < sigint) {
 86         ke2 = midBin;
 87         sigint2 = sigint;
 88       }
 89       else {
 90         ke1 = midBin;
 91         sigint1 = sigint;
 92       }
 93       if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " "
 94                                   << sigint1 << " " << sigint2 << G4endl;
 95     } while ( (ke2 - ke1) > 1 && --iterationsLeft > 0);    /* Loop checking, 30-Oct-2015, G.Folger */
 96    if ( iterationsLeft <= 0 ) {
 97        G4Exception("G4AngularDistributionNP", "im_r_matrix011", FatalException,
 98              "Problem with angular distribution (sig) data");
 99    }
100 
101     // sigint1 and sigint2 should be recoverable from above loop
102 
103     //    G4double dsig = sig[je2][ke1] - sig[je1][ke1];
104     //    G4double rc = dsig/delab;
105     //    G4double b = sig[je1][ke1] - rc*elab[je1];
106     //    G4double sigint1 = rc*ek + b;
107 
108     //    G4double dsig = sig[je2][ke2] - sig[je1][ke2];
109     //    G4double rc = dsig/delab;
110     //    G4double b = sig[je1][ke2] - rc*elab[je1];
111     //    G4double sigint2 = rc*ek + b;
112 
113     dsig = sigint2 - sigint1;
114     rc = 1./dsig;
115     b = ke1 - rc*sigint1;
116     G4double kint = rc*sample + b;
117     G4double theta = (0.5 + kint)*pi/180.;
118 
119     //    G4int k;
120     //std::abs(sample-sig[j][ke1]) < std::abs(sample-sig[j][ke2]) ? k = ke1 : k = ke2;
121     //    G4double theta = (0.5 + k)*pi/180.;
122 
123     if (verboseLevel > 1) {
124       G4cout << "   energy bin " << je1 << " energy=" << elab[je1] << G4endl;
125       G4cout << "   angle bin " << kint << " angle=" << theta/degree << G4endl;
126     }
127     G4double costh= std::cos(theta);
128     return costh;
129 }
130 
131 G4double G4AngularDistributionNP::Phi() const
132 {
133     return twopi * G4UniformRand();
134 }
135