ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuonDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuonDecayChannel.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 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 30 May 1997 H.Kurashige
34 //
35 // - Fix bug in calculation of electron energy in DecayIt()
36 // 28 Feb 01 H.Kurashige
37 // - Adding V-A fluxes for neutrinos using a new algorithm, 2005
38 // M. Melissas ( melissas AT cppm.in2p3.fr)
39 // J. Brunner ( brunner AT cppm.in2p3.fr)
40 // ------------------------------------------------------------
41 
42 #include "G4ParticleDefinition.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4DecayProducts.hh"
46 #include "G4VDecayChannel.hh"
47 #include "G4MuonDecayChannel.hh"
48 #include "Randomize.hh"
49 #include "G4LorentzVector.hh"
50 #include "G4LorentzRotation.hh"
51 #include "G4RotationMatrix.hh"
52 
55 {
56 }
57 
59  G4double theBR)
60  :G4VDecayChannel("Muon Decay",1)
61 {
62  // set names for daughter particles
63  if (theParentName == "mu+") {
64  SetBR(theBR);
65  SetParent("mu+");
67  SetDaughter(0, "e+");
68  SetDaughter(1, "nu_e");
69  SetDaughter(2, "anti_nu_mu");
70  } else if (theParentName == "mu-") {
71  SetBR(theBR);
72  SetParent("mu-");
74  SetDaughter(0, "e-");
75  SetDaughter(1, "anti_nu_e");
76  SetDaughter(2, "nu_mu");
77  } else {
78 #ifdef G4VERBOSE
79  if (GetVerboseLevel()>0) {
80  G4cout << "G4MuonDecayChannel:: constructor :";
81  G4cout << " parent particle is not muon but ";
82  G4cout << theParentName << G4endl;
83  }
84 #endif
85  }
86 }
87 
89  G4VDecayChannel(right)
90 {
91 }
92 
94 {
95 }
96 
98 {
99  if (this != &right) {
101  verboseLevel = right.verboseLevel;
102  rbranch = right.rbranch;
103 
104  // copy parent name
105  parent_name = new G4String(*right.parent_name);
106 
107  // clear daughters_name array
109 
110  // recreate array
112  if ( numberOfDaughters >0 ) {
115  //copy daughters name
116  for (G4int index=0; index < numberOfDaughters; index++) {
117  daughters_name[index] = new G4String(*right.daughters_name[index]);
118  }
119  }
120  }
121  return *this;
122 }
123 
125 {
126  // this version neglects muon polarization,and electron mass
127  // assumes the pure V-A coupling
128  // the Neutrinos are correctly V-A.
129 #ifdef G4VERBOSE
130  if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt ";
131 #endif
132 
135 
136  // parent mass
137  G4double parentmass = G4MT_parent->GetPDGMass();
138  const int N_DAUGHTER=3;
139 
140  //daughters'mass
141  G4double daughtermass[N_DAUGHTER];
142  G4double sumofdaughtermass = 0.0;
143  for (G4int index=0; index<N_DAUGHTER; 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[N_DAUGHTER];
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  const size_t MAX_LOOP=1000;
167  //Generating Random Energy
168  for (size_t loop1=0; loop1 <MAX_LOOP; ++loop1){
169  Ee=G4UniformRand();
170  for (size_t loop2 =0; loop2<MAX_LOOP; ++loop2){
171  x=xmax*G4UniformRand();
172  gam=G4UniformRand();
173  if (gam <= x*(1.-x)) break;
174  x = xmax;
175  }
176  Ene=x;
177  if ( Ene >= (1.-Ee)) break;
178  Ene = 1.-Ee;
179  }
180  G4double Enm=(2.-Ee-Ene);
181 
182 
183  //initialisation of rotation parameters
184 
185  G4double costheta,sintheta,rphi,rtheta,rpsi;
186  costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
187  sintheta=std::sqrt(1.-costheta*costheta);
188 
189 
190  rphi=twopi*G4UniformRand()*rad;
191  rtheta=(std::acos(2.*G4UniformRand()-1.));
192  rpsi=twopi*G4UniformRand()*rad;
193 
194  G4RotationMatrix rot;
195  rot.set(rphi,rtheta,rpsi);
196 
197  //electron 0
198  daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
199  G4ThreeVector direction0(0.0,0.0,1.0);
200 
201  direction0 *= rot;
202 
203  G4DynamicParticle * daughterparticle = new G4DynamicParticle ( G4MT_daughters[0], direction0 * daughtermomentum[0]);
204 
205  products->PushProducts(daughterparticle);
206 
207  //electronic neutrino 1
208 
209  daughtermomentum[1]=std::sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]);
210  G4ThreeVector direction1(sintheta,0.0,costheta);
211 
212  direction1 *= rot;
213 
214  G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( G4MT_daughters[1], direction1 * daughtermomentum[1]);
215  products->PushProducts(daughterparticle1);
216 
217  //muonnic neutrino 2
218 
219  daughtermomentum[2]=std::sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]);
220  G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
221 
222  direction2 *= rot;
223 
224  G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( G4MT_daughters[2],
225  direction2 * daughtermomentum[2]);
226  products->PushProducts(daughterparticle2);
227 
228  // output message
229 #ifdef G4VERBOSE
230  if (GetVerboseLevel()>1) {
231  G4cout << "G4MuonDecayChannel::DecayIt ";
232  G4cout << " create decay products in rest frame " <<G4endl;
233  products->DumpInfo();
234  }
235 #endif
236  return products;
237 }