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