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 ]

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