Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4MuonRadiativeDecayChannelWithSpin Class Reference

#include <G4MuonRadiativeDecayChannelWithSpin.hh>

Inheritance diagram for G4MuonRadiativeDecayChannelWithSpin:
Collaboration diagram for G4MuonRadiativeDecayChannelWithSpin:

Public Member Functions

 G4MuonRadiativeDecayChannelWithSpin (const G4String &theParentName, G4double theBR)
 
virtual ~G4MuonRadiativeDecayChannelWithSpin ()
 
virtual G4DecayProducts * DecayIt (G4double)
 
- Public Member Functions inherited from G4VDecayChannel
 G4VDecayChannel (const G4String &aName, G4int Verbose=1)
 
 G4VDecayChannel (const G4String &aName, const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="", const G4String &theDaughterName4="")
 
virtual ~G4VDecayChannel ()
 
G4int operator== (const G4VDecayChannel &right) const
 
G4int operator!= (const G4VDecayChannel &right) const
 
G4int operator< (const G4VDecayChannel &right) const
 
const G4String & GetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinition * GetParent ()
 
G4ParticleDefinition * GetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4String & GetParentName () const
 
const G4String & GetDaughterName (G4int anIndex) const
 
G4double GetParentMass () const
 
G4double GetDaughterMass (G4int anIndex) const
 
void SetParent (const G4ParticleDefinition *particle_type)
 
void SetParent (const G4String &particle_name)
 
void SetBR (G4double value)
 
void SetNumberOfDaughters (G4int value)
 
void SetDaughter (G4int anIndex, const G4ParticleDefinition *particle_type)
 
void SetDaughter (G4int anIndex, const G4String &particle_name)
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
void DumpInfo ()
 
G4double GetRangeMass () const
 
void SetRangeMass (G4double val)
 
virtual G4bool IsOKWithParentMass (G4double parentMass)
 
void SetPolarization (const G4ThreeVector &)
 
const G4ThreeVector & GetPolarization () const
 

Protected Member Functions

 G4MuonRadiativeDecayChannelWithSpin (const G4MuonRadiativeDecayChannelWithSpin &)
 
G4MuonRadiativeDecayChannelWithSpin & operator= (const G4MuonRadiativeDecayChannelWithSpin &)
 
- Protected Member Functions inherited from G4VDecayChannel
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChannel & operator= (const G4VDecayChannel &)
 
G4double DynamicalMass (G4double massPDG, G4double width, G4double maxDev=+1.) const
 

Additional Inherited Members

- Protected Attributes inherited from G4VDecayChannel
G4String kinematics_name
 
G4double rbranch
 
G4int numberOfDaughters
 
G4String * parent_name
 
G4String ** daughters_name
 
G4double rangeMass
 
G4ThreeVector parent_polarization
 
G4ParticleTable * particletable
 
G4int verboseLevel
 
G4ParticleDefinition * G4MT_parent
 
G4ParticleDefinition ** G4MT_daughters
 
G4double G4MT_parent_mass
 
G4double * G4MT_daughters_mass
 
G4double * G4MT_daughters_width
 
G4Mutex daughtersMutex
 
G4Mutex parentMutex
 
- Static Protected Attributes inherited from G4VDecayChannel
static const G4String noName = " "
 

Detailed Description

Definition at line 61 of file G4MuonRadiativeDecayChannelWithSpin.hh.

Constructor & Destructor Documentation

G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin ( const G4String &  theParentName,
G4double  theBR 
)

Definition at line 59 of file G4MuonRadiativeDecayChannelWithSpin.cc.

61  : G4VDecayChannel("Radiative 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, "gamma");
70  SetDaughter(2, "nu_e");
71  SetDaughter(3, "anti_nu_mu");
72  } else if (theParentName == "mu-") {
73  SetBR(theBR);
74  SetParent("mu-");
76  SetDaughter(0, "e-");
77  SetDaughter(1, "gamma");
78  SetDaughter(2, "anti_nu_e");
79  SetDaughter(3, "nu_mu");
80  } else {
81 #ifdef G4VERBOSE
82  if (GetVerboseLevel()>0) {
83  G4cout << "G4RadiativeMuonDecayChannel:: constructor :";
84  G4cout << " parent particle is not muon but ";
85  G4cout << theParentName << G4endl;
86  }
87 #endif
88  }
89 }
void SetBR(G4double value)
G4GLOB_DLL std::ostream G4cout
void SetNumberOfDaughters(G4int value)
G4int GetVerboseLevel() const
void SetParent(const G4ParticleDefinition *particle_type)
void SetDaughter(G4int anIndex, const G4ParticleDefinition *particle_type)
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin ( )
virtual

Definition at line 91 of file G4MuonRadiativeDecayChannelWithSpin.cc.

92 {
93 }
G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin ( const G4MuonRadiativeDecayChannelWithSpin &  right)
protected

Definition at line 95 of file G4MuonRadiativeDecayChannelWithSpin.cc.

95  :
96  G4VDecayChannel(right)
97 {
98 }

Here is the call graph for this function:

Member Function Documentation

G4DecayProducts * G4MuonRadiativeDecayChannelWithSpin::DecayIt ( G4double  )
virtual

if(i<10000000)goto leap1:

Implements G4VDecayChannel.

Definition at line 129 of file G4MuonRadiativeDecayChannelWithSpin.cc.

130 {
131 
132 #ifdef G4VERBOSE
133  if (GetVerboseLevel()>1)
134  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
135 #endif
136 
139 
140  // parent mass
141  G4double parentmass = G4MT_parent->GetPDGMass();
142 
143  G4double EMMU = parentmass;
144 
145  //daughters'mass
146  G4double daughtermass[4];
147  G4double sumofdaughtermass = 0.0;
148  for (G4int index=0; index<4; index++){
149  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
150  sumofdaughtermass += daughtermass[index];
151  }
152 
153  G4double EMASS = daughtermass[0];
154 
155  //create parent G4DynamicParticle at rest
156  G4ThreeVector dummy;
157  G4DynamicParticle * parentparticle =
158  new G4DynamicParticle( G4MT_parent, dummy, 0.0);
159  //create G4Decayproducts
160  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
161  delete parentparticle;
162 
163  G4int i = 0;
164 
165  G4double eps = EMASS/EMMU;
166 
167  G4double som0, x, y, xx, yy, zz;
168  G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
169  const size_t MAX_LOOP=10000;
170 
171  for (size_t loop_counter1=0; loop_counter1 <MAX_LOOP; ++loop_counter1){
172 
173 // leap1:
174 
175  i++;
176 
177 // leap2:
178 
179  for (size_t loop_counter2=0; loop_counter2 <MAX_LOOP; ++loop_counter2){
180 //
181 //--------------------------------------------------------------------------
182 // Build two vectors of random length and random direction, for the
183 // positron and the photon.
184 // x/y is the length of the vector, xx, yy and zz the components,
185 // phi is the azimutal angle, theta the polar angle.
186 //--------------------------------------------------------------------------
187 //
188 // For the positron
189 //
190  x = G4UniformRand();
191 
192  rn3dim(xx,yy,zz,x);
193 
194  if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(x*x))>0.001){
195  G4cout << "Norm of x not correct" << G4endl;
196  }
197 
198  phiE = atan4(xx,yy);
199  cthetaE = zz/x;
200  G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x;
201 //
202 // What you get:
203 //
204 // x = positron energy
205 // phiE = azimutal angle of positron momentum
206 // cthetaE = cosine of polar angle of positron momentum
207 // sthetaE = sine of polar angle of positron momentum
208 //
214 //
215 //-----------------------------------------------------------------------
216 //
217 // For the photon
218 //
219  y = G4UniformRand();
220 
221  rn3dim(xx,yy,zz,y);
222 
223  if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(y*y))>0.001){
224  G4cout << " Norm of y not correct " << G4endl;
225  }
226 
227  phiG = atan4(xx,yy);
228  cthetaG = zz/y;
229  G4double sthetaG = std::sqrt((xx*xx)+(yy*yy))/y;
230 //
231 // What you get:
232 //
233 // y = photon energy
234 // phiG = azimutal angle of photon momentum
235 // cthetaG = cosine of polar angle of photon momentum
236 // sthetaG = sine of polar angle of photon momentum
237 //
243 //
244 //-----------------------------------------------------------------------
245 //
246 // Maybe certain restrictions on the kinematical variables:
247 //
252 //
253 //-----------------------------------------------------------------------
254 //
255 // Calculate the angle between positron and photon (cosine)
256 //
257  cthetaGE = cthetaE*cthetaG+sthetaE*sthetaG*std::cos(phiE-phiG);
258 //
262 //
263 //-----------------------------------------------------------------------
264 //
265  G4double term0 = eps*eps;
266  G4double term1 = x*((1.0-eps)*(1.0-eps))+2.0*eps;
267  G4double beta = std::sqrt( x*((1.0-eps)*(1.0-eps))*
268  (x*((1.0-eps)*(1.0-eps))+4.0*eps))/term1;
269  G4double delta = 1.0-beta*cthetaGE;
270 
271  G4double term3 = y*(1.0-(eps*eps));
272  G4double term6 = term1*delta*term3;
273 
274  G4double Qsqr = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps));
275 //
276 //-----------------------------------------------------------------------
277 //
278 // Check the kinematics.
279 //
280  if ( Qsqr>=0.0 && Qsqr<=1.0 ) break;
281  }
282 //
284 //
285 // Do the calculation for -1 muon polarization (i.e. mu+)
286 //
287  G4double Pmu = -1.0;
288  if (GetParentName() == "mu-")Pmu = +1.0;
289 //
290 // and for Fronsdal
291 //
292 //-----------------------------------------------------------------------
293 //
294  som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE);
295 //
303 //
304 //-----------------------------------------------------------------------
305 //
307 //
308 //----------------------------------------------------------------------
309 //
310 // Sample the decay rate
311 //
312 
313  if (G4UniformRand()*177.0 <= som0) break;
314  }
315 
317 //
318 //-----------------------------------------------------------------------
319 //
320  G4double E = EMMU/2.*(x*((1.-eps)*(1.-eps))+2.*eps);
321  G4double G = EMMU/2.*y*(1.-eps*eps);
322 //
323 //-----------------------------------------------------------------------
324 //
325 
326  if(E < EMASS) E = EMASS;
327 
328  // calculate daughter momentum
329  G4double daughtermomentum[4];
330 
331  daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
332 
333  G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
334  G4double cphiE = std::cos(phiE);
335  G4double sphiE = std::sin(phiE);
336 
337  //Coordinates of the decay positron with respect to the muon spin
338 
339  G4double px = sthetaE*cphiE;
340  G4double py = sthetaE*sphiE;
341  G4double pz = cthetaE;
342 
343  G4ThreeVector direction0(px,py,pz);
344 
345  direction0.rotateUz(parent_polarization);
346 
347  G4DynamicParticle * daughterparticle0
348  = new G4DynamicParticle( G4MT_daughters[0], daughtermomentum[0]*direction0);
349 
350  products->PushProducts(daughterparticle0);
351 
352  daughtermomentum[1] = G;
353 
354  G4double sthetaG = std::sqrt(1.-cthetaG*cthetaG);
355  G4double cphiG = std::cos(phiG);
356  G4double sphiG = std::sin(phiG);
357 
358  //Coordinates of the decay gamma with respect to the muon spin
359 
360  px = sthetaG*cphiG;
361  py = sthetaG*sphiG;
362  pz = cthetaG;
363 
364  G4ThreeVector direction1(px,py,pz);
365 
366  direction1.rotateUz(parent_polarization);
367 
368  G4DynamicParticle * daughterparticle1
369  = new G4DynamicParticle( G4MT_daughters[1], daughtermomentum[1]*direction1);
370 
371  products->PushProducts(daughterparticle1);
372 
373  // daughter 3 ,4 (neutrinos)
374  // create neutrinos in the C.M frame of two neutrinos
375 
376  G4double energy2 = parentmass*(1.0 - (x+y)/2.0);
377 
378  G4double vmass2 = energy2*energy2 -
379  (daughtermomentum[0]*direction0+daughtermomentum[1]*direction1)*
380  (daughtermomentum[0]*direction0+daughtermomentum[1]*direction1);
381  G4double vmass = std::sqrt(vmass2);
382 
383  G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
384  beta = -1.0 * std::min(beta,0.99);
385 
386  G4double costhetan = 2.*G4UniformRand()-1.0;
387  G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
388  G4double phin = twopi*G4UniformRand()*rad;
389  G4double sinphin = std::sin(phin);
390  G4double cosphin = std::cos(phin);
391 
392  G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan);
393 
394  G4DynamicParticle * daughterparticle2
395  = new G4DynamicParticle( G4MT_daughters[2], direction2*(vmass/2.));
396  G4DynamicParticle * daughterparticle3
397  = new G4DynamicParticle( G4MT_daughters[3], direction2*(-1.0*vmass/2.));
398 
399  // boost to the muon rest frame
400 
401  G4ThreeVector direction34(direction0.x()+direction1.x(),
402  direction0.y()+direction1.y(),
403  direction0.z()+direction1.z());
404  direction34 = direction34.unit();
405 
406  G4LorentzVector p4 = daughterparticle2->Get4Momentum();
407  p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
408  daughterparticle2->Set4Momentum(p4);
409 
410  p4 = daughterparticle3->Get4Momentum();
411  p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
412  daughterparticle3->Set4Momentum(p4);
413 
414  products->PushProducts(daughterparticle2);
415  products->PushProducts(daughterparticle3);
416 
417  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
418  daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
419 
420 // output message
421 #ifdef G4VERBOSE
422  if (GetVerboseLevel()>1) {
423  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
424  G4cout << " create decay products in rest frame " <<G4endl;
425  products->DumpInfo();
426  }
427 #endif
428  return products;
429 }
void CheckAndFillDaughters()
const G4String & GetParentName() const
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
static const G4double eps
static constexpr double rad
Definition: G4SIunits.hh:149
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void DumpInfo() const
HepLorentzVector & boost(double, double, double)
G4ThreeVector parent_polarization
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4int GetVerboseLevel() const
G4double GetPDGMass() const
Hep3Vector unit() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator= ( const G4MuonRadiativeDecayChannelWithSpin &  right)
protected

Definition at line 100 of file G4MuonRadiativeDecayChannelWithSpin.cc.

101 {
102  if (this != &right) {
104  verboseLevel = right.verboseLevel;
105  rbranch = right.rbranch;
106 
107  // copy parent name
108  parent_name = new G4String(*right.parent_name);
109 
110  // clear daughters_name array
112 
113  // recreate array
115  if ( numberOfDaughters >0 ) {
118  //copy daughters name
119  for (G4int index=0; index < numberOfDaughters; index++) {
120  daughters_name[index] = new G4String(*right.daughters_name[index]);
121  }
122  }
124  }
125  return *this;
126 }
int G4int
Definition: G4Types.hh:78
G4String kinematics_name
G4String * parent_name
G4ThreeVector parent_polarization
G4String ** daughters_name

Here is the call graph for this function:


The documentation for this class was generated from the following files: