ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PionRadiativeDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PionRadiativeDecayChannel.cc
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 }