Geant4: G4MuonRadiativeDecayChannelWithSpin Class Reference

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

Definition at line 60 of file G4MuonRadiativeDecayChannelWithSpin.cc.

References G4cout, G4endl, G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::SetBR(), G4VDecayChannel::SetDaughter(), G4VDecayChannel::SetNumberOfDaughters(), and G4VDecayChannel::SetParent().

Referenced by G4MuonRadiativeDecayChannelWithSpin().

00062              : G4VDecayChannel("Radiative Muon Decay",1),
00063                parent_polarization()
00064 {
00065   // set names for daughter particles
00066   if (theParentName == "mu+") {
00067     SetBR(theBR);
00068     SetParent("mu+");
00069     SetNumberOfDaughters(4);
00070     SetDaughter(0, "e+");
00071     SetDaughter(1, "gamma");
00072     SetDaughter(2, "nu_e");
00073     SetDaughter(3, "anti_nu_mu");
00074   } else if (theParentName == "mu-") {
00075     SetBR(theBR);
00076     SetParent("mu-");
00077     SetNumberOfDaughters(4);
00078     SetDaughter(0, "e-");
00079     SetDaughter(1, "gamma");
00080     SetDaughter(2, "anti_nu_e");
00081     SetDaughter(3, "nu_mu");
00082   } else {
00083 #ifdef G4VERBOSE
00084     if (GetVerboseLevel()>0) {
00085       G4cout << "G4RadiativeMuonDecayChannel:: constructor :";
00086       G4cout << " parent particle is not muon but ";
00087       G4cout << theParentName << G4endl;
00088     }
00089 #endif
00090   }
00091 }

G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin ( ) [virtual]

Definition at line 93 of file G4MuonRadiativeDecayChannelWithSpin.cc.

00094 {
00095 }

G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin ( const G4MuonRadiativeDecayChannelWithSpin & ) [protected]

Definition at line 97 of file G4MuonRadiativeDecayChannelWithSpin.cc.

References G4MuonRadiativeDecayChannelWithSpin(), and parent_polarization.

00097                                                                                                                         :
00098   G4VDecayChannel(right)
00099 {
00100   parent_polarization = right.parent_polarization;
00101 }

G4DecayProducts * G4MuonRadiativeDecayChannelWithSpin::DecayIt ( G4double ) [virtual]

Implements G4VDecayChannel.

Definition at line 132 of file G4MuonRadiativeDecayChannelWithSpin.cc.

References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4VDecayChannel::FillDaughters(), G4VDecayChannel::FillParent(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4VDecayChannel::GetParentName(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalMomentum(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, G4DecayProducts::PushProducts(), and G4DynamicParticle::Set4Momentum().

00133 {
00134 
00135 #ifdef G4VERBOSE
00136   if (GetVerboseLevel()>1) 
00137                  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
00138 #endif
00139 
00140   if (parent == 0) FillParent();  
00141   if (daughters == 0) FillDaughters();
00142 
00143   // parent mass
00144   G4double parentmass = parent->GetPDGMass();
00145 
00146   G4double EMMU = parentmass;
00147 
00148   //daughters'mass
00149   G4double daughtermass[4]; 
00150   G4double sumofdaughtermass = 0.0;
00151   for (G4int index=0; index<4; index++){
00152     daughtermass[index] = daughters[index]->GetPDGMass();
00153     sumofdaughtermass += daughtermass[index];
00154   }
00155 
00156   G4double EMASS = daughtermass[0];
00157 
00158   //create parent G4DynamicParticle at rest
00159   G4ThreeVector dummy;
00160   G4DynamicParticle * parentparticle = 
00161                                new G4DynamicParticle( parent, dummy, 0.0);
00162   //create G4Decayproducts
00163   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00164   delete parentparticle;
00165 
00166   G4int i = 0;
00167 
00168   G4double eps = EMASS/EMMU;
00169 
00170   G4double som0, Qsqr, x, y, xx, yy, zz;
00171   G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
00172 
00173   do {
00174 
00175 //     leap1:
00176 
00177      i++;
00178 
00179 //     leap2:
00180 
00181      do {
00182 //
00183 //--------------------------------------------------------------------------
00184 //      Build two vectors of random length and random direction, for the
00185 //      positron and the photon.
00186 //      x/y is the length of the vector, xx, yy and zz the components,
00187 //      phi is the azimutal angle, theta the polar angle.
00188 //--------------------------------------------------------------------------
00189 //
00190 //      For the positron
00191 //
00192         x = G4UniformRand();
00193 
00194         rn3dim(xx,yy,zz,x);
00195 
00196         if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(x*x))>0.001){
00197           G4cout << "Norm of x not correct" << G4endl;
00198         }
00199 
00200         phiE = atan4(xx,yy);
00201         cthetaE = zz/x;
00202         G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x;
00203 //
00204 //      What you get:
00205 //
00206 //      x       = positron energy
00207 //      phiE    = azimutal angle of positron momentum
00208 //      cthetaE = cosine of polar angle of positron momentum
00209 //      sthetaE = sine of polar angle of positron momentum
00210 //
00216 //
00217 //-----------------------------------------------------------------------
00218 //
00219 //      For the photon
00220 //
00221         y = G4UniformRand();
00222 
00223         rn3dim(xx,yy,zz,y);
00224 
00225         if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(y*y))>0.001){
00226           G4cout << " Norm of y not correct " << G4endl;
00227         }
00228 
00229         phiG = atan4(xx,yy);
00230         cthetaG = zz/y;
00231         G4double sthetaG = std::sqrt((xx*xx)+(yy*yy))/y;
00232 //
00233 //      What you get:
00234 //
00235 //      y       = photon energy
00236 //      phiG    = azimutal angle of photon momentum
00237 //      cthetaG = cosine of polar angle of photon momentum
00238 //      sthetaG = sine of polar angle of photon momentum
00239 //
00245 //
00246 //-----------------------------------------------------------------------
00247 //
00248 //      Maybe certain restrictions on the kinematical variables:
00249 //
00254 //
00255 //-----------------------------------------------------------------------
00256 //
00257 //      Calculate the angle between positron and photon (cosine)
00258 //
00259         cthetaGE = cthetaE*cthetaG+sthetaE*sthetaG*std::cos(phiE-phiG);
00260 //
00264 //
00265 //-----------------------------------------------------------------------
00266 //
00267         G4double term0 = eps*eps;
00268         G4double term1 = x*((1.0-eps)*(1.0-eps))+2.0*eps;
00269         G4double beta  = std::sqrt( x*((1.0-eps)*(1.0-eps))*
00270                                    (x*((1.0-eps)*(1.0-eps))+4.0*eps))/term1;
00271         G4double delta = 1.0-beta*cthetaGE;
00272 
00273         G4double term3 = y*(1.0-(eps*eps));
00274         G4double term6 = term1*delta*term3;
00275 
00276         Qsqr  = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps));
00277 //
00278 //-----------------------------------------------------------------------
00279 //
00280 //      Check the kinematics.
00281 //
00282      } while ( Qsqr<0.0 || Qsqr>1.0 );
00283 //
00285 //
00286 //   Do the calculation for -1 muon polarization (i.e. mu+)
00287 //
00288      G4double Pmu = -1.0;
00289      if (GetParentName() == "mu-")Pmu = +1.0;
00290 //
00291 //   and for Fronsdal
00292 //
00293 //-----------------------------------------------------------------------
00294 //
00295      som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE);
00296 //
00304 //
00305 //-----------------------------------------------------------------------
00306 //
00308 //
00309 //----------------------------------------------------------------------
00310 //
00311 //   Sample the decay rate
00312 //
00313 
00314   } while (G4UniformRand()*177.0 > som0);
00315 
00317 //
00318 //-----------------------------------------------------------------------
00319 //
00320   G4double E = EMMU/2.*(x*((1.-eps)*(1.-eps))+2.*eps);
00321   G4double G = EMMU/2.*y*(1.-eps*eps);
00322 //
00323 //-----------------------------------------------------------------------
00324 //
00325 
00326   if(E < EMASS) E = EMASS;
00327 
00328   // calculate daughter momentum
00329   G4double daughtermomentum[4];
00330 
00331   daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
00332 
00333   G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
00334   G4double cphiE = std::cos(phiE);
00335   G4double sphiE = std::sin(phiE);
00336 
00337   //Coordinates of the decay positron with respect to the muon spin
00338 
00339   G4double px = sthetaE*cphiE;
00340   G4double py = sthetaE*sphiE;
00341   G4double pz = cthetaE;
00342 
00343   G4ThreeVector direction0(px,py,pz);
00344 
00345   direction0.rotateUz(parent_polarization);
00346 
00347   G4DynamicParticle * daughterparticle0 
00348     = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
00349 
00350   products->PushProducts(daughterparticle0);
00351 
00352   daughtermomentum[1] = G;
00353 
00354   G4double sthetaG = std::sqrt(1.-cthetaG*cthetaG);
00355   G4double cphiG = std::cos(phiG);
00356   G4double sphiG = std::sin(phiG);
00357 
00358   //Coordinates of the decay gamma with respect to the muon spin
00359 
00360   px = sthetaG*cphiG;
00361   py = sthetaG*sphiG;
00362   pz = cthetaG;
00363 
00364   G4ThreeVector direction1(px,py,pz);
00365 
00366   direction1.rotateUz(parent_polarization);
00367 
00368   G4DynamicParticle * daughterparticle1
00369     = new G4DynamicParticle( daughters[1], daughtermomentum[1]*direction1);
00370 
00371   products->PushProducts(daughterparticle1);
00372 
00373   // daughter 3 ,4 (neutrinos)
00374   // create neutrinos in the C.M frame of two neutrinos
00375 
00376   G4double energy2 = parentmass*(1.0 - (x+y)/2.0);
00377 
00378   G4double vmass   = std::sqrt((energy2-
00379                                 (daughtermomentum[0]+daughtermomentum[1]))*
00380                                (energy2+
00381                                 (daughtermomentum[0]+daughtermomentum[1])));
00382   G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
00383   beta = -1.0 * std::min(beta,0.99);
00384 
00385   G4double costhetan = 2.*G4UniformRand()-1.0;
00386   G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00387   G4double phin  = twopi*G4UniformRand()*rad;
00388   G4double sinphin = std::sin(phin);
00389   G4double cosphin = std::cos(phin);
00390 
00391   G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan);
00392 
00393   G4DynamicParticle * daughterparticle2
00394     = new G4DynamicParticle( daughters[2], direction2*(vmass/2.));
00395   G4DynamicParticle * daughterparticle3
00396     = new G4DynamicParticle( daughters[3], direction2*(-1.0*vmass/2.));
00397 
00398   // boost to the muon rest frame
00399 
00400   G4ThreeVector direction34(direction0.x()+direction1.x(),
00401                             direction0.y()+direction1.y(),
00402                             direction0.z()+direction1.z());
00403   direction34 = direction34.unit();
00404 
00405   G4LorentzVector p4 = daughterparticle2->Get4Momentum();
00406   p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
00407   daughterparticle2->Set4Momentum(p4);
00408 
00409   p4 = daughterparticle3->Get4Momentum();
00410   p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
00411   daughterparticle3->Set4Momentum(p4);
00412 
00413   products->PushProducts(daughterparticle2);
00414   products->PushProducts(daughterparticle3);
00415 
00416   daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
00417   daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
00418 
00419 // output message
00420 #ifdef G4VERBOSE
00421   if (GetVerboseLevel()>1) {
00422     G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
00423     G4cout << "  create decay products in rest frame " <<G4endl;
00424     products->DumpInfo();
00425   }
00426 #endif
00427   return products;
00428 }

const G4ThreeVector & G4MuonRadiativeDecayChannelWithSpin::GetPolarization ( ) const [inline]

Definition at line 114 of file G4MuonRadiativeDecayChannelWithSpin.hh.

00115 {
00116   return parent_polarization;
00117 }

G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator= ( const G4MuonRadiativeDecayChannelWithSpin & ) [protected]

Definition at line 103 of file G4MuonRadiativeDecayChannelWithSpin.cc.

References G4VDecayChannel::ClearDaughtersName(), G4VDecayChannel::daughters_name, G4VDecayChannel::kinematics_name, G4VDecayChannel::numberOfDaughters, G4VDecayChannel::parent_name, parent_polarization, G4VDecayChannel::rbranch, and G4VDecayChannel::verboseLevel.

00104 {
00105   if (this != &right) { 
00106     kinematics_name = right.kinematics_name;
00107     verboseLevel = right.verboseLevel;
00108     rbranch = right.rbranch;
00109 
00110     // copy parent name
00111     parent_name = new G4String(*right.parent_name);
00112 
00113     // clear daughters_name array
00114     ClearDaughtersName();
00115 
00116     // recreate array
00117     numberOfDaughters = right.numberOfDaughters;
00118     if ( numberOfDaughters >0 ) {
00119       if (daughters_name !=0) ClearDaughtersName();
00120       daughters_name = new G4String*[numberOfDaughters];
00121       //copy daughters name
00122       for (G4int index=0; index < numberOfDaughters; index++) {
00123           daughters_name[index] = new G4String(*right.daughters_name[index]);
00124       }
00125     }
00126     parent_polarization = right.parent_polarization;
00127   }
00128   return *this;
00129 }

void G4MuonRadiativeDecayChannelWithSpin::SetPolarization ( G4ThreeVector ) [inline]

Definition at line 108 of file G4MuonRadiativeDecayChannelWithSpin.hh.

00109 {
00110   parent_polarization = polar;
00111 }


Public Member Functions
	G4MuonRadiativeDecayChannelWithSpin (const G4String &theParentName, G4double theBR)
virtual	~G4MuonRadiativeDecayChannelWithSpin ()
virtual G4DecayProducts *	DecayIt (G4double)
void	SetPolarization (G4ThreeVector)
const G4ThreeVector &	GetPolarization () const
Protected Member Functions
	G4MuonRadiativeDecayChannelWithSpin (const G4MuonRadiativeDecayChannelWithSpin &)
G4MuonRadiativeDecayChannelWithSpin &	operator= (const G4MuonRadiativeDecayChannelWithSpin &)

G4MuonRadiativeDecayChannelWithSpin Class Reference

Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation