Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/particles/management/src/G4PionRadiativeDecayChannel.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 /particles/management/src/G4PionRadiativeDecayChannel.cc (Version 11.3.0) and /particles/management/src/G4PionRadiativeDecayChannel.cc (Version 10.7.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 // G4PionRadiativeDecayChannel class implement     26 // G4PionRadiativeDecayChannel class implementation
 27 //      GEANT 4 class header file                  27 //      GEANT 4 class header file
 28 //                                                 28 //
 29 // Author: P.Gumplinger, 30 July 2007          <<  29 // Author: P.Gumplinger, 30 July 2007 
 30 // Reference: M. Blecher, TRIUMF/PIENU Technot     30 // Reference: M. Blecher, TRIUMF/PIENU Technote
 31 //            "Inclusion of pi->enug in the Mo <<  31 //            "Inclusion of pi->enug in the Monte Carlo" 
 32 // -------------------------------------------     32 // --------------------------------------------------------------------
 33                                                    33 
 34 #include "G4PionRadiativeDecayChannel.hh"          34 #include "G4PionRadiativeDecayChannel.hh"
 35                                                    35 
 36 #include "G4DecayProducts.hh"                  << 
 37 #include "G4LorentzVector.hh"                  << 
 38 #include "G4PhysicalConstants.hh"                  36 #include "G4PhysicalConstants.hh"
 39 #include "G4SystemOfUnits.hh"                      37 #include "G4SystemOfUnits.hh"
 40 #include "Randomize.hh"                            38 #include "Randomize.hh"
                                                   >>  39 #include "G4DecayProducts.hh"
                                                   >>  40 #include "G4LorentzVector.hh"
 41                                                    41 
 42 namespace                                          42 namespace
 43 {                                                  43 {
 44 const G4double beta = 3.6612e-03;              <<  44   const G4double beta = 3.6612e-03;
 45 const G4double cib = 1.16141e-03;              <<  45   const G4double cib  = 1.16141e-03;
 46 const G4double csdp = 3.45055e-02;             <<  46   const G4double csdp = 3.45055e-02;
 47 const G4double csdm = 5.14122e-03;             <<  47   const G4double csdm = 5.14122e-03;
 48 const G4double cif = 4.63543e-05;              <<  48   const G4double cif  = 4.63543e-05;
 49 const G4double cig = 1.78928e-05;              <<  49   const G4double cig  = 1.78928e-05;
 50 const G4double xl = 2. * 0.1 * MeV / 139.57 *  <<  50   const G4double xl = 2.*0.1*MeV/139.57*MeV;
 51 const G4double yl = ((1. - xl) + std::sqrt((1  <<  51   const G4double yl = ((1.-xl) + std::sqrt((1-xl)*(1-xl)+4*beta*beta))/2.;
 52                                                <<  52 
 53 const G4double xu = 1. - (yl - std::sqrt(yl *  <<  53   const G4double xu = 1. - (yl - std::sqrt(yl*yl-4.*beta*beta))/2.;
 54 const G4double yu = 1. + beta * beta;          <<  54   const G4double yu = 1. + beta*beta;
 55                                                <<  55 
 56 inline G4double D2W(const G4double x, const G4 <<  56   inline G4double D2W(const G4double x,const G4double y)
 57 {                                              <<  57   {
 58   return cib * (1. - y) * (1. + ((1. - x) * (1 <<  58     return cib*(1.-y)*(1.+((1.-x)*(1.-x)))/((x*x)*(x+y-1.)) +
 59          + csdp * (1. - x) * ((x + y - 1.) * ( <<  59         csdp*(1.-x)*((x+y-1.)*(x+y-1.)) +
 60          + cif * (x - 1.) * (1. - y) / x + cig <<  60         csdm*(1.-x)*((1.-y)*(1.-y)) +
                                                   >>  61         cif*(x-1.)*(1.-y)/x +
                                                   >>  62         cig*(1.-y)*(1.-x+(x*x)/(x+y-1.))/x;
                                                   >>  63   }
                                                   >>  64 
                                                   >>  65   const G4double d2wmax = D2W(xl,yl);
 61 }                                                  66 }
 62                                                    67 
 63 const G4double d2wmax = D2W(xl, yl);           <<  68 // --------------------------------------------------------------------
 64 }  // namespace                                <<  69 G4PionRadiativeDecayChannel::G4PionRadiativeDecayChannel()
                                                   >>  70   : G4VDecayChannel()
                                                   >>  71 {
                                                   >>  72 }
 65                                                    73 
 66 G4PionRadiativeDecayChannel::G4PionRadiativeDe <<  74 // --------------------------------------------------------------------
 67                                                <<  75 G4PionRadiativeDecayChannel::
                                                   >>  76 G4PionRadiativeDecayChannel(const G4String& theParentName,
                                                   >>  77                                   G4double        theBR)
 68   : G4VDecayChannel("Radiative Pion Decay", 1)     78   : G4VDecayChannel("Radiative Pion Decay", 1)
 69 {                                                  79 {
 70   // set names for daughter particles              80   // set names for daughter particles
 71   if (theParentName == "pi+") {                <<  81   if (theParentName == "pi+")
                                                   >>  82   {
 72     SetBR(theBR);                                  83     SetBR(theBR);
 73     SetParent("pi+");                              84     SetParent("pi+");
 74     SetNumberOfDaughters(3);                       85     SetNumberOfDaughters(3);
 75     SetDaughter(0, "e+");                          86     SetDaughter(0, "e+");
 76     SetDaughter(1, "gamma");                       87     SetDaughter(1, "gamma");
 77     SetDaughter(2, "nu_e");                        88     SetDaughter(2, "nu_e");
 78   }                                                89   }
 79   else if (theParentName == "pi-") {           <<  90   else if (theParentName == "pi-")
                                                   >>  91   {
 80     SetBR(theBR);                                  92     SetBR(theBR);
 81     SetParent("pi-");                              93     SetParent("pi-");
 82     SetNumberOfDaughters(3);                       94     SetNumberOfDaughters(3);
 83     SetDaughter(0, "e-");                          95     SetDaughter(0, "e-");
 84     SetDaughter(1, "gamma");                       96     SetDaughter(1, "gamma");
 85     SetDaughter(2, "anti_nu_e");                   97     SetDaughter(2, "anti_nu_e");
 86   }                                                98   }
 87   else {                                       <<  99   else
                                                   >> 100   {
 88 #ifdef G4VERBOSE                                  101 #ifdef G4VERBOSE
 89     if (GetVerboseLevel() > 0) {               << 102     if (GetVerboseLevel()>0)
 90       G4cout << "G4RadiativePionDecayChannel:: << 103     {
                                                   >> 104       G4cout << "G4RadiativePionDecayChannel::G4PionRadiativeDecayChannel()"
                                                   >> 105              << G4endl;
 91       G4cout << "Parent particle is not charge    106       G4cout << "Parent particle is not charged pion: ";
 92       G4cout << theParentName << G4endl;          107       G4cout << theParentName << G4endl;
 93     }                                             108     }
 94 #endif                                            109 #endif
 95   }                                               110   }
 96 }                                                 111 }
 97                                                   112 
                                                   >> 113 // --------------------------------------------------------------------
                                                   >> 114 G4PionRadiativeDecayChannel::~G4PionRadiativeDecayChannel()
                                                   >> 115 {
                                                   >> 116 }
                                                   >> 117 
                                                   >> 118 // --------------------------------------------------------------------
                                                   >> 119 G4PionRadiativeDecayChannel::
                                                   >> 120 G4PionRadiativeDecayChannel(const G4PionRadiativeDecayChannel& right)
                                                   >> 121   : G4VDecayChannel(right)
                                                   >> 122 {
                                                   >> 123 }
                                                   >> 124 
 98 G4PionRadiativeDecayChannel&                      125 G4PionRadiativeDecayChannel&
 99 G4PionRadiativeDecayChannel::operator=(const G    126 G4PionRadiativeDecayChannel::operator=(const G4PionRadiativeDecayChannel& right)
100 {                                                 127 {
101   if (this != &right) {                        << 128   if (this != &right)
                                                   >> 129   { 
102     kinematics_name = right.kinematics_name;      130     kinematics_name = right.kinematics_name;
103     verboseLevel = right.verboseLevel;            131     verboseLevel = right.verboseLevel;
104     rbranch = right.rbranch;                      132     rbranch = right.rbranch;
105                                                   133 
106     // copy parent name                           134     // copy parent name
107     parent_name = new G4String(*right.parent_n    135     parent_name = new G4String(*right.parent_name);
108                                                   136 
109     // clear daughters_name array                 137     // clear daughters_name array
110     ClearDaughtersName();                         138     ClearDaughtersName();
111                                                   139 
112     // recreate array                             140     // recreate array
113     numberOfDaughters = right.numberOfDaughter    141     numberOfDaughters = right.numberOfDaughters;
114     if (numberOfDaughters > 0) {               << 142     if ( numberOfDaughters >0 )
                                                   >> 143     {
115       if (daughters_name != nullptr) ClearDaug    144       if (daughters_name != nullptr) ClearDaughtersName();
116       daughters_name = new G4String*[numberOfD    145       daughters_name = new G4String*[numberOfDaughters];
117       // copy daughters name                   << 146       //copy daughters name
118       for (G4int index = 0; index < numberOfDa << 147       for (G4int index=0; index<numberOfDaughters; ++index)
                                                   >> 148       {
119         daughters_name[index] = new G4String(*    149         daughters_name[index] = new G4String(*right.daughters_name[index]);
120       }                                           150       }
121     }                                             151     }
122   }                                               152   }
123   return *this;                                   153   return *this;
124 }                                                 154 }
125                                                   155 
126 G4DecayProducts* G4PionRadiativeDecayChannel:: << 156 // --------------------------------------------------------------------
                                                   >> 157 G4DecayProducts* G4PionRadiativeDecayChannel::DecayIt(G4double) 
127 {                                                 158 {
                                                   >> 159 
128 #ifdef G4VERBOSE                                  160 #ifdef G4VERBOSE
129   if (GetVerboseLevel() > 1) G4cout << "G4Pion << 161   if (GetVerboseLevel()>1) 
                                                   >> 162      G4cout << "G4PionRadiativeDecayChannel::DecayIt ";
130 #endif                                            163 #endif
131                                                   164 
132   CheckAndFillParent();                           165   CheckAndFillParent();
133   CheckAndFillDaughters();                        166   CheckAndFillDaughters();
134                                                   167 
135   // parent mass                                  168   // parent mass
136   G4double parentmass = G4MT_parent->GetPDGMas    169   G4double parentmass = G4MT_parent->GetPDGMass();
137                                                   170 
138   G4double EMPI = parentmass;                     171   G4double EMPI = parentmass;
139                                                   172 
140   // daughters'mass                               173   // daughters'mass
141   const G4int N_DAUGHTER = 3;                  << 174   const G4int N_DAUGHTER=3;
142   G4double daughtermass[N_DAUGHTER];           << 175   G4double daughtermass[N_DAUGHTER]; 
143   // G4double sumofdaughtermass = 0.0;         << 176   G4double sumofdaughtermass = 0.0;
144   for (G4int index = 0; index < N_DAUGHTER; ++ << 177   for (G4int index=0; index<N_DAUGHTER; ++index)
                                                   >> 178   {
145     daughtermass[index] = G4MT_daughters[index    179     daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
146     // sumofdaughtermass += daughtermass[index << 180     sumofdaughtermass += daughtermass[index];
147   }                                               181   }
148                                                   182 
149   G4double EMASS = daughtermass[0];               183   G4double EMASS = daughtermass[0];
150                                                   184 
151   // create parent G4DynamicParticle at rest      185   // create parent G4DynamicParticle at rest
152   G4ThreeVector dummy;                            186   G4ThreeVector dummy;
153   auto parentparticle = new G4DynamicParticle( << 187   G4DynamicParticle* parentparticle
                                                   >> 188     = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
154   // create G4Decayproducts                       189   // create G4Decayproducts
155   auto products = new G4DecayProducts(*parentp << 190   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
156   delete parentparticle;                          191   delete parentparticle;
157                                                   192 
158   G4double x, y;                                  193   G4double x, y;
159                                                   194 
160   const std::size_t MAX_LOOP = 1000;           << 195   const std::size_t MAX_LOOP=1000;
161                                                   196 
162   for (std::size_t loop_counter1 = 0; loop_cou << 197   for (std::size_t loop_counter1=0; loop_counter1<MAX_LOOP; ++loop_counter1)
163     for (std::size_t loop_counter2 = 0; loop_c << 198   {
164       x = xl + G4UniformRand() * (xu - xl);    << 199     for (std::size_t loop_counter2=0; loop_counter2<MAX_LOOP; ++loop_counter2)
165       y = yl + G4UniformRand() * (yu - yl);    << 200     {
166       if (x + y > 1.) break;                   << 201       x = xl + G4UniformRand()*(xu-xl);
                                                   >> 202       y = yl + G4UniformRand()*(yu-yl);
                                                   >> 203       if (x+y > 1.) break;
167     }                                             204     }
168     G4double d2w = D2W(x, y);                  << 205     G4double d2w = D2W(x,y);
169     if (d2w > G4UniformRand() * d2wmax) break; << 206     if (d2w > G4UniformRand()*d2wmax) break;
170   }                                               207   }
171                                                   208 
172   // Calculate the angle between positron and     209   // Calculate the angle between positron and photon (cosine)
173   //                                              210   //
174   G4double cthetaGE =                          << 211   G4double cthetaGE = (y*(x-2.)+2.*(1.-x+beta*beta)) /
175     (y * (x - 2.) + 2. * (1. - x + beta * beta << 212                       (x*std::sqrt(y*y-4.*beta*beta));
176                                                   213 
177   G4double G = x * EMPI / 2.;                  << 214   G4double G = x * EMPI/2.;
178   G4double E = y * EMPI / 2.;                  << 215   G4double E = y * EMPI/2.;
179                                                   216 
180   if (E < EMASS) E = EMASS;                       217   if (E < EMASS) E = EMASS;
181                                                   218 
182   // calculate daughter momentum                  219   // calculate daughter momentum
183   G4double daughtermomentum[2];                   220   G4double daughtermomentum[2];
184                                                   221 
185   daughtermomentum[0] = std::sqrt(E * E - EMAS << 222   daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
186                                                   223 
187   G4double cthetaE = 2. * G4UniformRand() - 1. << 224   G4double cthetaE = 2.*G4UniformRand()-1.;
188   G4double sthetaE = std::sqrt(1. - cthetaE *  << 225   G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
189                                                   226 
190   G4double phiE = twopi * G4UniformRand() * ra << 227   G4double phiE = twopi*G4UniformRand()*rad;
191   G4double cphiE = std::cos(phiE);                228   G4double cphiE = std::cos(phiE);
192   G4double sphiE = std::sin(phiE);                229   G4double sphiE = std::sin(phiE);
193                                                   230 
194   // Coordinates of the decay positron            231   // Coordinates of the decay positron
195   //                                              232   //
196   G4double px = sthetaE * cphiE;               << 233   G4double px = sthetaE*cphiE;
197   G4double py = sthetaE * sphiE;               << 234   G4double py = sthetaE*sphiE;
198   G4double pz = cthetaE;                          235   G4double pz = cthetaE;
199                                                   236 
200   G4ThreeVector direction0(px, py, pz);        << 237   G4ThreeVector direction0(px,py,pz);
201                                                   238 
202   auto daughterparticle0 =                     << 239   G4DynamicParticle * daughterparticle0 
203     new G4DynamicParticle(G4MT_daughters[0], d << 240     = new G4DynamicParticle(G4MT_daughters[0], daughtermomentum[0]*direction0);
204                                                   241 
205   products->PushProducts(daughterparticle0);      242   products->PushProducts(daughterparticle0);
206                                                   243 
207   daughtermomentum[1] = G;                        244   daughtermomentum[1] = G;
208                                                   245 
209   G4double sthetaGE = std::sqrt(1. - cthetaGE  << 246   G4double sthetaGE = std::sqrt(1.-cthetaGE*cthetaGE);
210                                                   247 
211   G4double phiGE = twopi * G4UniformRand() * r << 248   G4double phiGE = twopi*G4UniformRand()*rad;
212   G4double cphiGE = std::cos(phiGE);              249   G4double cphiGE = std::cos(phiGE);
213   G4double sphiGE = std::sin(phiGE);              250   G4double sphiGE = std::sin(phiGE);
214                                                   251 
215   // Coordinates of the decay gamma with respe    252   // Coordinates of the decay gamma with respect to the decay positron
216   //                                              253   //
217   px = sthetaGE * cphiGE;                      << 254   px = sthetaGE*cphiGE;
218   py = sthetaGE * sphiGE;                      << 255   py = sthetaGE*sphiGE;
219   pz = cthetaGE;                                  256   pz = cthetaGE;
220                                                   257 
221   G4ThreeVector direction1(px, py, pz);        << 258   G4ThreeVector direction1(px,py,pz);
222                                                   259 
223   direction1.rotateUz(direction0);                260   direction1.rotateUz(direction0);
224                                                   261 
225   auto daughterparticle1 =                     << 262   G4DynamicParticle * daughterparticle1
226     new G4DynamicParticle(G4MT_daughters[1], d << 263     = new G4DynamicParticle(G4MT_daughters[1], daughtermomentum[1]*direction1);
227                                                   264 
228   products->PushProducts(daughterparticle1);      265   products->PushProducts(daughterparticle1);
229                                                   266 
230   // output message                               267   // output message
231 #ifdef G4VERBOSE                                  268 #ifdef G4VERBOSE
232   if (GetVerboseLevel() > 1) {                 << 269   if (GetVerboseLevel()>1)
                                                   >> 270   {
233     G4cout << "G4PionRadiativeDecayChannel::De    271     G4cout << "G4PionRadiativeDecayChannel::DecayIt() -";
234     G4cout << " create decay products in rest     272     G4cout << " create decay products in rest frame " << G4endl;
235     products->DumpInfo();                         273     products->DumpInfo();
236   }                                               274   }
237 #endif                                            275 #endif
238                                                   276 
239   return products;                                277   return products;
240 }                                                 278 }
241                                                   279