Geant4  10.01.p02
G4MuonDecayChannel.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 // $Id: G4MuonDecayChannel.cc 67971 2013-03-13 10:13:24Z gcosmo $
28 //
29 //
30 // ------------------------------------------------------------
31 // GEANT 4 class header file
32 //
33 // History: first implementation, based on object model of
34 // 30 May 1997 H.Kurashige
35 //
36 // Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige
37 //2005
38 // M. Melissas ( melissas AT cppm.in2p3.fr)
39 // J. Brunner ( brunner AT cppm.in2p3.fr)
40 // Adding V-A fluxes for neutrinos using a new algortithm :
41 // ------------------------------------------------------------
42 
43 #include "G4ParticleDefinition.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4DecayProducts.hh"
47 #include "G4VDecayChannel.hh"
48 #include "G4MuonDecayChannel.hh"
49 #include "Randomize.hh"
50 #include "G4LorentzVector.hh"
51 #include "G4LorentzRotation.hh"
52 #include "G4RotationMatrix.hh"
53 
56 {
57 }
58 
60  G4double theBR)
61  :G4VDecayChannel("Muon Decay",1)
62 {
63  // set names for daughter particles
64  if (theParentName == "mu+") {
65  SetBR(theBR);
66  SetParent("mu+");
68  SetDaughter(0, "e+");
69  SetDaughter(1, "nu_e");
70  SetDaughter(2, "anti_nu_mu");
71  } else if (theParentName == "mu-") {
72  SetBR(theBR);
73  SetParent("mu-");
75  SetDaughter(0, "e-");
76  SetDaughter(1, "anti_nu_e");
77  SetDaughter(2, "nu_mu");
78  } else {
79 #ifdef G4VERBOSE
80  if (GetVerboseLevel()>0) {
81  G4cout << "G4MuonDecayChannel:: constructor :";
82  G4cout << " parent particle is not muon but ";
83  G4cout << theParentName << G4endl;
84  }
85 #endif
86  }
87 }
88 
90  G4VDecayChannel(right)
91 {
92 }
93 
95 {
96 }
97 
99 {
100  if (this != &right) {
102  verboseLevel = right.verboseLevel;
103  rbranch = right.rbranch;
104 
105  // copy parent name
106  parent_name = new G4String(*right.parent_name);
107 
108  // clear daughters_name array
110 
111  // recreate array
113  if ( numberOfDaughters >0 ) {
116  //copy daughters name
117  for (G4int index=0; index < numberOfDaughters; index++) {
118  daughters_name[index] = new G4String(*right.daughters_name[index]);
119  }
120  }
121  }
122  return *this;
123 }
124 
126 {
127  // this version neglects muon polarization,and electron mass
128  // assumes the pure V-A coupling
129  // the Neutrinos are correctly V-A.
130 #ifdef G4VERBOSE
131  if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt ";
132 #endif
133 
134  if (G4MT_parent == 0) FillParent();
135  if (G4MT_daughters == 0) FillDaughters();
136 
137  // parent mass
138  G4double parentmass = G4MT_parent->GetPDGMass();
139 
140  //daughters'mass
141  G4double daughtermass[3];
142  G4double sumofdaughtermass = 0.0;
143  for (G4int index=0; index<3; index++){
144  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
145  sumofdaughtermass += daughtermass[index];
146  }
147 
148  //create parent G4DynamicParticle at rest
149  G4ThreeVector dummy;
150  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
151  //create G4Decayproducts
152  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
153  delete parentparticle;
154 
155  // calculate daughter momentum
156  G4double daughtermomentum[3];
157  // calcurate electron energy
158  G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
159  G4double x;
160 
161  G4double Ee,Ene;
162 
163  G4double gam;
164  G4double EMax=parentmass/2-daughtermass[0];
165 
166 
167  //Generating Random Energy
168 do {
169  Ee=G4UniformRand();
170  do{
171  x=xmax*G4UniformRand();
172  gam=G4UniformRand();
173  }while (gam >x*(1.-x));
174  Ene=x;
175  } while ( Ene < (1.-Ee));
176  G4double Enm=(2.-Ee-Ene);
177 
178 
179  //initialisation of rotation parameters
180 
181  G4double costheta,sintheta,rphi,rtheta,rpsi;
182  costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
183  sintheta=std::sqrt(1.-costheta*costheta);
184 
185 
186  rphi=twopi*G4UniformRand()*rad;
187  rtheta=(std::acos(2.*G4UniformRand()-1.));
188  rpsi=twopi*G4UniformRand()*rad;
189 
190  G4RotationMatrix rot;
191  rot.set(rphi,rtheta,rpsi);
192 
193  //electron 0
194  daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
195  G4ThreeVector direction0(0.0,0.0,1.0);
196 
197  direction0 *= rot;
198 
199  G4DynamicParticle * daughterparticle = new G4DynamicParticle ( G4MT_daughters[0], direction0 * daughtermomentum[0]);
200 
201  products->PushProducts(daughterparticle);
202 
203  //electronic neutrino 1
204 
205  daughtermomentum[1]=std::sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]);
206  G4ThreeVector direction1(sintheta,0.0,costheta);
207 
208  direction1 *= rot;
209 
210  G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( G4MT_daughters[1], direction1 * daughtermomentum[1]);
211  products->PushProducts(daughterparticle1);
212 
213  //muonnic neutrino 2
214 
215  daughtermomentum[2]=std::sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]);
216  G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
217 
218  direction2 *= rot;
219 
220  G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( G4MT_daughters[2],
221  direction2 * daughtermomentum[2]);
222  products->PushProducts(daughterparticle2);
223 
224 
225 
226 
227  // output message
228 #ifdef G4VERBOSE
229  if (GetVerboseLevel()>1) {
230  G4cout << "G4MuonDecayChannel::DecayIt ";
231  G4cout << " create decay products in rest frame " <<G4endl;
232  products->DumpInfo();
233  }
234 #endif
235  return products;
236 }
237 
238 
239 
240 
241 
242 
void SetBR(G4double value)
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
G4String kinematics_name
void DumpInfo() const
void SetNumberOfDaughters(G4int value)
G4String * parent_name
static const double rad
Definition: G4SIunits.hh:130
G4int GetVerboseLevel() const
G4double GetPDGMass() const
virtual G4DecayProducts * DecayIt(G4double)
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
G4MuonDecayChannel & operator=(const G4MuonDecayChannel &)
G4String ** daughters_name