Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLEtaNElasticChannel.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 ]

Diff markup

Differences between /processes/hadronic/models/inclxx/incl_physics/src/G4INCLEtaNElasticChannel.cc (Version 11.3.0) and /processes/hadronic/models/inclxx/incl_physics/src/G4INCLEtaNElasticChannel.cc (Version 10.3.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // INCL++ intra-nuclear cascade model              26 // INCL++ intra-nuclear cascade model
 27 // Alain Boudard, CEA-Saclay, France               27 // Alain Boudard, CEA-Saclay, France
 28 // Joseph Cugnon, University of Liege, Belgium     28 // Joseph Cugnon, University of Liege, Belgium
 29 // Jean-Christophe David, CEA-Saclay, France       29 // Jean-Christophe David, CEA-Saclay, France
 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H     30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
 31 // Sylvie Leray, CEA-Saclay, France                31 // Sylvie Leray, CEA-Saclay, France
 32 // Davide Mancusi, CEA-Saclay, France              32 // Davide Mancusi, CEA-Saclay, France
 33 //                                                 33 //
 34 #define INCLXX_IN_GEANT4_MODE 1                    34 #define INCLXX_IN_GEANT4_MODE 1
 35                                                    35 
 36 #include "globals.hh"                              36 #include "globals.hh"
 37                                                    37 
 38 #include "G4INCLEtaNElasticChannel.hh"             38 #include "G4INCLEtaNElasticChannel.hh"
 39 #include "G4INCLKinematicsUtils.hh"                39 #include "G4INCLKinematicsUtils.hh"
 40 #include "G4INCLBinaryCollisionAvatar.hh"          40 #include "G4INCLBinaryCollisionAvatar.hh"
 41 #include "G4INCLRandom.hh"                         41 #include "G4INCLRandom.hh"
 42 #include "G4INCLGlobals.hh"                        42 #include "G4INCLGlobals.hh"
 43 #include "G4INCLLogger.hh"                         43 #include "G4INCLLogger.hh"
 44                                                    44 
 45 namespace G4INCL {                                 45 namespace G4INCL {
 46                                                    46 
 47     EtaNElasticChannel::EtaNElasticChannel(Par     47     EtaNElasticChannel::EtaNElasticChannel(Particle *p1, Particle *p2)
 48     : particle1(p1), particle2(p2)                 48     : particle1(p1), particle2(p2)
 49     {                                              49     {
 50                                                    50 
 51     }                                              51     }
 52                                                    52 
 53     EtaNElasticChannel::~EtaNElasticChannel(){     53     EtaNElasticChannel::~EtaNElasticChannel(){
 54                                                    54 
 55     }                                              55     }
 56                                                    56 
 57     void EtaNElasticChannel::fillFinalState(Fi     57     void EtaNElasticChannel::fillFinalState(FinalState *fs) {
 58         Particle * nucleon;                        58         Particle * nucleon;
 59         Particle * eta;                            59         Particle * eta;
 60         if(particle1->isNucleon()) {               60         if(particle1->isNucleon()) {
 61             nucleon = particle1;                   61             nucleon = particle1;
 62             eta = particle2;                       62             eta = particle2;
 63         } else {                                   63         } else {
 64             nucleon = particle2;                   64             nucleon = particle2;
 65             eta = particle1;                       65             eta = particle1;
 66         }                                          66         }
 67                                                    67 
 68           G4double plab=KinematicsUtils::momen     68           G4double plab=KinematicsUtils::momentumInLab(particle1, particle2);
 69                                                    69 
 70           G4double sh=nucleon->getEnergy()+eta     70           G4double sh=nucleon->getEnergy()+eta->getEnergy();
 71           G4double mn=nucleon->getMass();          71           G4double mn=nucleon->getMass();
 72           G4double me=eta->getMass();              72           G4double me=eta->getMass();
 73           G4double en=(sh*sh+mn*mn-me*me)/(2*s     73           G4double en=(sh*sh+mn*mn-me*me)/(2*sh);
 74           nucleon->setEnergy(en);                  74           nucleon->setEnergy(en);
 75           G4double ee=std::sqrt(en*en-mn*mn+me     75           G4double ee=std::sqrt(en*en-mn*mn+me*me);
 76           eta->setEnergy(ee);                      76           eta->setEnergy(ee);
 77           G4double pn=std::sqrt(en*en-mn*mn);      77           G4double pn=std::sqrt(en*en-mn*mn);         
 78                                                    78 
 79           ThreeVector mom_nucleon;                 79           ThreeVector mom_nucleon;
 80                                                    80 
 81           if (plab < 250.) {                       81           if (plab < 250.) {
 82 // Isotropy                                        82 // Isotropy
 83           mom_nucleon = Random::normVector(pn)     83           mom_nucleon = Random::normVector(pn);
 84           }                                        84           }         
 85                                                    85 
 86 // From Kamano                                     86 // From Kamano
 87      else {                                        87      else {
 88                                                    88 
 89             const G4double pi=std::acos(-1.0); <<  89             const G4double pi=std::acos(-1.0);    
 90             G4double x1;                           90             G4double x1;
 91             G4double u1;                           91             G4double u1;
 92             G4double fteta;                        92             G4double fteta;
 93             G4double teta;                         93             G4double teta;
 94             G4double fi;                           94             G4double fi;
 95                                                    95 
 96             G4double a0;                           96             G4double a0;
 97             G4double a1;                           97             G4double a1;
 98             G4double a2;                           98             G4double a2;
 99             G4double a3;                           99             G4double a3;
100             G4double a4;                          100             G4double a4;
101             G4double a5;                          101             G4double a5;
102             G4double a6;                          102             G4double a6;
103                                                   103             
104             if (plab > 1400.) plab=1400.; // n    104             if (plab > 1400.) plab=1400.; // no information on angular distributions above plab=1400 MeV
105       G4double p6=std::pow(plab, 6);              105       G4double p6=std::pow(plab, 6); 
106             G4double p5=std::pow(plab, 5);        106             G4double p5=std::pow(plab, 5); 
107             G4double p4=std::pow(plab, 4);        107             G4double p4=std::pow(plab, 4); 
108             G4double p3=std::pow(plab, 3);        108             G4double p3=std::pow(plab, 3); 
109             G4double p2=std::pow(plab, 2);        109             G4double p2=std::pow(plab, 2); 
110             G4double p1=plab;                     110             G4double p1=plab; 
111                                                   111             
112 // a6                                             112 // a6
113             if (plab < 300.) {                    113             if (plab < 300.) {
114        a6=-8.384000E-08*p1 - 1.15452E-04;         114        a6=-8.384000E-08*p1 - 1.15452E-04;
115             }                                     115             }
116       else if (plab < 500.){                      116       else if (plab < 500.){
117               a6=1.593966E-13*p4 - 2.619560E-1    117               a6=1.593966E-13*p4 - 2.619560E-10*p3 + 1.564701E-07*p2 - 3.986627E-05*p1 + 3.622575E-03;
118       }                                           118       }
119             else {                                119             else {
120               a6=6.143615E-20*p6 - 3.157181E-1    120               a6=6.143615E-20*p6 - 3.157181E-16*p5 + 6.348289E-13*p4 - 6.117961E-10*p3 + 2.764542E-07*p2 - 4.391048E-05*p1 - 1.443857E-03;
121             }                                     121             }
122 // a5                                             122 // a5
123             if (plab < 650.) {                    123             if (plab < 650.) {
124               a5=-9.021076E-18*p6 + 2.176771E-    124               a5=-9.021076E-18*p6 + 2.176771E-14*p5 - 2.136095E-11*p4 + 1.100580E-08*p3 - 3.150857E-06*p2 + 4.761016E-04*p1 - 2.969608E-02;
125             }                                     125             }
126       else if (plab < 950.){                      126       else if (plab < 950.){
127               a5=4.424756E-18*p6 - 1.756295E-1    127               a5=4.424756E-18*p6 - 1.756295E-14*p5 + 2.625428E-11*p4 - 1.678272E-08*p3 + 2.227237E-06*p2 + 2.146666E-03*p1 - 7.065712E-01;
128       }                                           128       }
129             else {                                129             else {
130               a5=2.209585E-19*p6 - 1.546647E-1    130               a5=2.209585E-19*p6 - 1.546647E-15*p5 + 4.578142E-12*p4 - 7.303856E-09*p3 + 6.604074E-06*p2 - 3.205628E-03*p1 + 6.534893E-01;
131             }                                     131             }
132 // a4                                             132 // a4
133             if (plab < 700.) {                    133             if (plab < 700.) {
134               a4=4.826684E-17*p6 - 1.534471E-1    134               a4=4.826684E-17*p6 - 1.534471E-13*p5 + 1.907868E-10*p4 - 1.192317E-07*p3 + 3.988902E-05*p2 - 6.822100E-03*p1 + 4.684685E-01;
135             }                                     135             }
136             else {                                136             else {
137               a4=-3.245143E-18*p6 + 2.174395E-    137               a4=-3.245143E-18*p6 + 2.174395E-14*p5 - 6.012288E-11*p4 + 8.772790E-08*p3 - 7.113554E-05*p2 + 3.029285E-02*p1 - 5.237677E+00;
138             }                                     138             }
139 // a3                                             139 // a3
140             if (plab < 650.) {                    140             if (plab < 650.) {
141               a3=3.783071E-17*p6 - 1.151454E-1    141               a3=3.783071E-17*p6 - 1.151454E-13*p5 + 1.357165E-10*p4 - 8.036891E-08*p3 + 2.572396E-05*p2 - 4.245566E-03*p1 + 2.832772E-01;
142             }                                     142             }
143             else {                                143             else {
144               a3=-5.063316E-18*p6 + 3.223757E-    144               a3=-5.063316E-18*p6 + 3.223757E-14*p5 - 8.435635E-11*p4 + 1.159487E-07*p3 - 8.812510E-05*p2 + 3.500692E-02*p1 - 5.624556E+00;
145             }                                     145             }
146 // a2                                             146 // a2
147             if (plab < 500.) {                    147             if (plab < 500.) {
148               a2=-6.085067E-14*p5 + 1.354078E-    148               a2=-6.085067E-14*p5 + 1.354078E-10*p4 - 1.124158E-07*p3 + 4.292106E-05*p2 - 7.218145E-03*p1 + 4.584962E-01;
149             }                                     149             }
150             else if (plab < 750.) {               150             else if (plab < 750.) {
151               a2= 9.512730E-11*p4 - 2.362724E-    151               a2= 9.512730E-11*p4 - 2.362724E-07*p3 + 2.171883E-04*p2 - 8.742722E-02*p1 + 1.309433E+01;
152             }                                     152             }
153             else {                                153             else {
154               a2=-4.228889E-18*p6 + 2.798222E-    154               a2=-4.228889E-18*p6 + 2.798222E-14*p5 - 7.640831E-11*p4 + 1.100124E-07*p3 - 8.778573E-05*p2 + 3.652772E-02*p1 - 6.025497E+00;
155             }                                     155             }
156 // a1                                             156 // a1
157             if (plab < 500.) {                    157             if (plab < 500.) {
158               a1=-1.524408E-14*p5 + 3.007021E-    158               a1=-1.524408E-14*p5 + 3.007021E-11*p4 - 2.129570E-08*p3 + 5.607250E-06*p2 - 3.001598E-04*p1 + 8.701280E-04;
159             }                                     159             }
160             else if (plab < 750.) {               160             else if (plab < 750.) {
161               a1=-3.255396E-11*p4 + 8.168681E-    161               a1=-3.255396E-11*p4 + 8.168681E-08*p3 - 7.447474E-05*p2 + 2.917630E-02*p1 - 4.152037E+00;
162             }                                     162             }
163             else {                                163             else {
164               a1=9.964504E-19*p6 - 6.380168E-1    164               a1=9.964504E-19*p6 - 6.380168E-15*p5 + 1.638691E-11*p4 - 2.107063E-08*p3 + 1.347462E-05*p2 - 3.318304E-03*p1 - 5.030932E-02;
165             }                                     165             }
166 // a0                                             166 // a0
167               a0=-3.220143E-17*p6 + 1.789654E-    167               a0=-3.220143E-17*p6 + 1.789654E-13*p5 - 3.912863E-10*p4 + 4.181510E-07*p3 - 2.147259E-04*p2 + 3.856266E-02*p1 + 2.609971E+00;
168                                                   168             
169             G4double interg1=2.*(a6/7. + a4/5.    169             G4double interg1=2.*(a6/7. + a4/5. + a2/3. + a0); // (integral to normalize)
170             G4double f1=(a6+a5+a4+a3+a2+a1+a0)    170             G4double f1=(a6+a5+a4+a3+a2+a1+a0)/interg1; // (Max normalized)   
171                                                   171             
172             G4int passe1=0;                       172             G4int passe1=0;
173             while (passe1==0) {                   173             while (passe1==0) {
174               // Sample x from -1 to 1            174               // Sample x from -1 to 1
175               x1=Random::shoot();                 175               x1=Random::shoot();
176               if (Random::shoot() > 0.5) x1=-x    176               if (Random::shoot() > 0.5) x1=-x1;
177                                                   177               
178               // Sample u from 0 to 1             178               // Sample u from 0 to 1
179               u1=Random::shoot();                 179               u1=Random::shoot();
180               fteta=(a6*x1*x1*x1*x1*x1*x1+a5*x    180               fteta=(a6*x1*x1*x1*x1*x1*x1+a5*x1*x1*x1*x1*x1+a4*x1*x1*x1*x1+a3*x1*x1*x1+a2*x1*x1+a1*x1+a0)/interg1;
181               // The condition                    181               // The condition
182               if (u1*f1 < fteta) {                182               if (u1*f1 < fteta) {
183                 teta=std::acos(x1);               183                 teta=std::acos(x1);
184                 //        std::cout << x1  <<     184                 //        std::cout << x1  << " " << fteta << " "<< f1/interg1 << " " << u1 << " " << interg1 << std::endl;
185                 passe1=1;                         185                 passe1=1;
186               }                                   186               }
187             }                                     187             }
188                                                   188             
189             fi=(2.0*pi)*Random::shoot();          189             fi=(2.0*pi)*Random::shoot();    
190                                                   190             
191             ThreeVector mom_nucleon1(             191             ThreeVector mom_nucleon1(
192                pn*std::sin(teta)*std::cos(fi), << 192                                                             pn*std::sin(teta)*std::cos(fi),
193                pn*std::sin(teta)*std::sin(fi), << 193                                                             pn*std::sin(teta)*std::sin(fi),
194                pn*std::cos(teta)               << 194                                                             pn*std::cos(teta)
195             );                                 << 195                                                             );
196                                                   196             
197             mom_nucleon = -mom_nucleon1 ;         197             mom_nucleon = -mom_nucleon1 ;
198                                                   198           
199           }                                       199           }
200                                                   200 
201           nucleon->setMomentum(mom_nucleon);      201           nucleon->setMomentum(mom_nucleon);
202           eta->setMomentum(-mom_nucleon);         202           eta->setMomentum(-mom_nucleon);
203                                                   203 
204          fs->addModifiedParticle(nucleon);        204          fs->addModifiedParticle(nucleon);
205           fs->addModifiedParticle(eta);           205           fs->addModifiedParticle(eta);
206                                                   206     
207         }                                         207         }
208 }                                                 208 }
209                                                   209