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 9.5.p2)


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