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 G4DecayProductsDecayIt (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 G4StringGetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinitionGetParent ()
 
G4ParticleDefinitionGetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4StringGetParentName () const
 
const G4StringGetDaughterName (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 G4ThreeVectorGetPolarization () const
 

Protected Member Functions

 G4MuonRadiativeDecayChannelWithSpin (const G4MuonRadiativeDecayChannelWithSpin &)
 
G4MuonRadiativeDecayChannelWithSpinoperator= (const G4MuonRadiativeDecayChannelWithSpin &)
 
- Protected Member Functions inherited from G4VDecayChannel
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChanneloperator= (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
 
G4Stringparent_name
 
G4String ** daughters_name
 
G4double rangeMass
 
G4ThreeVector parent_polarization
 
G4ParticleTableparticletable
 
G4int verboseLevel
 
G4ParticleDefinitionG4MT_parent
 
G4ParticleDefinition ** G4MT_daughters
 
G4double G4MT_parent_mass
 
G4doubleG4MT_daughters_mass
 
G4doubleG4MT_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: