Geant4  10.02.p01
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 
49  : G4VDecayChannel(),
50  beta(0.),cib(0.),csdp(0.),csdm(0.),cif(0.),cig(0.),
51  xl(0.), yl(0.), xu(0.), yu(0.), d2wmax(0.)
52 {
53 }
54 
57  G4double theBR)
58  : G4VDecayChannel("Radiative Pion Decay",1)
59 {
60  // set names for daughter particles
61  if (theParentName == "pi+") {
62  SetBR(theBR);
63  SetParent("pi+");
65  SetDaughter(0, "e+");
66  SetDaughter(1, "gamma");
67  SetDaughter(2, "nu_e");
68  } else if (theParentName == "pi-") {
69  SetBR(theBR);
70  SetParent("pi-");
72  SetDaughter(0, "e-");
73  SetDaughter(1, "gamma");
74  SetDaughter(2, "anti_nu_e");
75  } else {
76 #ifdef G4VERBOSE
77  if (GetVerboseLevel()>0) {
78  G4cout << "G4RadiativePionDecayChannel:: constructor :";
79  G4cout << " parent particle is not charged pion but ";
80  G4cout << theParentName << G4endl;
81  }
82 #endif
83  }
84 
85  beta = 3.6612e-03;
86 
87  cib = 1.16141e-03;
88  csdp = 3.45055e-02;
89  csdm = 5.14122e-03;
90  cif = 4.63543e-05;
91  cig = 1.78928e-05;
92 
93  xl = 2.*0.1*MeV/139.57*MeV;
94  yl = ((1.-xl) + std::sqrt((1-xl)*(1-xl)+4*beta*beta))/2.;
95 
96  xu = 1. - (yl - std::sqrt(yl*yl-4.*beta*beta))/2.;
97  yu = 1. + beta*beta;
98 
99  d2wmax = D2W(xl,yl);
100 
101 }
102 
104 {
105 }
107  :G4VDecayChannel(right),
108  beta(right.beta),cib(right.cib),csdp(right.csdp),
109  csdm(right.csdm),cif(right.cif),cig(right.cig),
110  xl(right.xl), yl(right.yl), xu(right.xu), yu(right.yu),
111  d2wmax(right.d2wmax)
112 {
113 }
114 
116 {
117  if (this != &right) {
119  verboseLevel = right.verboseLevel;
120  rbranch = right.rbranch;
121 
122  // copy parent name
123  parent_name = new G4String(*right.parent_name);
124 
125  // clear daughters_name array
127 
128  // recreate array
130  if ( numberOfDaughters >0 ) {
133  //copy daughters name
134  for (G4int index=0; index < numberOfDaughters; index++) {
135  daughters_name[index] = new G4String(*right.daughters_name[index]);
136  }
137  }
138  beta = right.beta;
139  cib = right.cib;
140  csdp = right.csdp;
141  csdm = right.csdm;
142  cif = right.cif;
143  cig = right.cig;
144  xl = right.xl;
145  yl = right.yl;
146  xu = right.xu;
147  yu = right.yu;
148  d2wmax = right.d2wmax;
149  }
150  return *this;
151 }
152 
154 {
155 
156 #ifdef G4VERBOSE
157  if (GetVerboseLevel()>1)
158  G4cout << "G4PionRadiativeDecayChannel::DecayIt ";
159 #endif
160 
161  if (G4MT_parent == 0) FillParent();
162  if (G4MT_daughters == 0) FillDaughters();
163 
164  // parent mass
165  G4double parentmass = G4MT_parent->GetPDGMass();
166 
167  G4double EMPI = parentmass;
168 
169  //daughters'mass
170  const G4int N_DAUGHTER=3;
171  G4double daughtermass[N_DAUGHTER];
172  G4double sumofdaughtermass = 0.0;
173  for (G4int index=0; index<N_DAUGHTER; index++){
174  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
175  sumofdaughtermass += daughtermass[index];
176  }
177 
178  G4double EMASS = daughtermass[0];
179 
180  //create parent G4DynamicParticle at rest
181  G4ThreeVector dummy;
182  G4DynamicParticle * parentparticle =
183  new G4DynamicParticle( G4MT_parent, dummy, 0.0);
184  //create G4Decayproducts
185  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
186  delete parentparticle;
187 
188  G4double x, y;
189 
190  const size_t MAX_LOOP=1000;
191 
192  for (size_t loop_counter1=0; loop_counter1<MAX_LOOP; ++loop_counter1){
193  for (size_t loop_counter2=0; loop_counter2<MAX_LOOP; ++loop_counter2){
194  x = xl + G4UniformRand()*(xu-xl);
195  y = yl + G4UniformRand()*(yu-yl);
196  if (x+y > 1.) break;
197  }
198  G4double d2w = D2W(x,y);
199  if (d2w > G4UniformRand()*d2wmax) break;
200  }
201 
202 //-----------------------------------------------------------------------
203 //
204 // Calculate the angle between positron and photon (cosine)
205 //
206  G4double cthetaGE = (y*(x-2.)+2.*(1.-x+beta*beta)) /
207  (x*std::sqrt(y*y-4.*beta*beta));
208 
209 //
210 //-----------------------------------------------------------------------
211 //
212  G4double G = x * EMPI/2.;
213  G4double E = y * EMPI/2.;
214 //
215 //-----------------------------------------------------------------------
216 //
217 
218  if (E < EMASS) E = EMASS;
219 
220  // calculate daughter momentum
221  G4double daughtermomentum[2];
222 
223  daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
224 
225  G4double cthetaE = 2.*G4UniformRand()-1.;
226  G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
227 
228  G4double phiE = twopi*G4UniformRand()*rad;
229  G4double cphiE = std::cos(phiE);
230  G4double sphiE = std::sin(phiE);
231 
232  //Coordinates of the decay positron
233 
234  G4double px = sthetaE*cphiE;
235  G4double py = sthetaE*sphiE;
236  G4double pz = cthetaE;
237 
238  G4ThreeVector direction0(px,py,pz);
239 
240  G4DynamicParticle * daughterparticle0
241  = new G4DynamicParticle( G4MT_daughters[0], daughtermomentum[0]*direction0);
242 
243  products->PushProducts(daughterparticle0);
244 
245  daughtermomentum[1] = G;
246 
247  G4double sthetaGE = std::sqrt(1.-cthetaGE*cthetaGE);
248 
249  G4double phiGE = twopi*G4UniformRand()*rad;
250  G4double cphiGE = std::cos(phiGE);
251  G4double sphiGE = std::sin(phiGE);
252 
253  //Coordinates of the decay gamma with respect to the decay positron
254 
255  px = sthetaGE*cphiGE;
256  py = sthetaGE*sphiGE;
257  pz = cthetaGE;
258 
259  G4ThreeVector direction1(px,py,pz);
260 
261  direction1.rotateUz(direction0);
262 
263  G4DynamicParticle * daughterparticle1
264  = new G4DynamicParticle( G4MT_daughters[1], daughtermomentum[1]*direction1);
265 
266  products->PushProducts(daughterparticle1);
267 
268 // output message
269 #ifdef G4VERBOSE
270  if (GetVerboseLevel()>1) {
271  G4cout << "G4PionRadiativeDecayChannel::DecayIt ";
272  G4cout << " create decay products in rest frame " <<G4endl;
273  products->DumpInfo();
274  }
275 #endif
276 
277  return products;
278 
279 }
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 &)
G4double D2W(const G4double x, const G4double y)
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