G4MuonDecayChannelWithSpin Class Reference

#include <G4MuonDecayChannelWithSpin.hh>

Inheritance diagram for G4MuonDecayChannelWithSpin:

G4MuonDecayChannel G4VDecayChannel

Public Member Functions

 G4MuonDecayChannelWithSpin (const G4String &theParentName, G4double theBR)
virtual ~G4MuonDecayChannelWithSpin ()
virtual G4DecayProducts * DecayIt (G4double)
void SetPolarization (G4ThreeVector)
const G4ThreeVector & GetPolarization () const

Protected Member Functions

 G4MuonDecayChannelWithSpin (const G4MuonDecayChannelWithSpin &)
G4MuonDecayChannelWithSpin & operator= (const G4MuonDecayChannelWithSpin &)

Detailed Description

Definition at line 50 of file G4MuonDecayChannelWithSpin.hh.

Constructor & Destructor Documentation

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

Definition at line 59 of file G4MuonDecayChannelWithSpin.cc.

References G4MuonDecayChannelWithSpin().

Referenced by G4MuonDecayChannelWithSpin().

00061   : G4MuonDecayChannel(theParentName,theBR),
00062     parent_polarization(),
00063     EMMU( 0.*MeV),
00064     EMASS( 0.*MeV) 
00065 {
00066 }

G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin (  )  [virtual]

Definition at line 68 of file G4MuonDecayChannelWithSpin.cc.

00069 {
00070 }

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

Definition at line 72 of file G4MuonDecayChannelWithSpin.cc.

References EMASS, EMMU, G4MuonDecayChannelWithSpin(), and parent_polarization.

00072                                                                                              :
00073   G4MuonDecayChannel(right)
00074 {
00075   parent_polarization = right.parent_polarization;
00076   EMMU  = right.EMMU;
00077   EMASS = right.EMASS;
00078 }

Member Function Documentation

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

Reimplemented from G4MuonDecayChannel.

Definition at line 111 of file G4MuonDecayChannelWithSpin.cc.

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

00112 {
00113   // This version assumes V-A coupling with 1st order radiative correctons,
00114   //              the standard model Michel parameter values, but
00115   //              gives incorrect energy spectrum for neutrinos
00116 
00117 #ifdef G4VERBOSE
00118   if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
00119 #endif
00120 
00121   if (parent == 0) FillParent();  
00122   if (daughters == 0) FillDaughters();
00123 
00124   // parent mass
00125   G4double parentmass = parent->GetPDGMass();
00126 
00127   EMMU = parentmass;
00128 
00129   //daughters'mass
00130   G4double daughtermass[3]; 
00131   G4double sumofdaughtermass = 0.0;
00132   for (G4int index=0; index<3; index++){
00133     daughtermass[index] = daughters[index]->GetPDGMass();
00134     sumofdaughtermass += daughtermass[index];
00135   }
00136 
00137   EMASS = daughtermass[0];
00138 
00139   //create parent G4DynamicParticle at rest
00140   G4ThreeVector dummy;
00141   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00142   //create G4Decayproducts
00143   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00144   delete parentparticle;
00145 
00146   // calcurate electron energy
00147 
00148   G4double michel_rho   = 0.75; //Standard Model Michel rho
00149   G4double michel_delta = 0.75; //Standard Model Michel delta
00150   G4double michel_xsi   = 1.00; //Standard Model Michel xsi
00151   G4double michel_eta   = 0.00; //Standard Model eta
00152 
00153   G4double rndm, x, ctheta;
00154 
00155   G4double FG; 
00156   G4double FG_max = 2.00;
00157 
00158   G4double W_mue  = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
00159   G4double x0     =           EMASS/W_mue;
00160 
00161   G4double x0_squared = x0*x0;
00162 
00163   // ***************************************************
00164   //     x0 <= x <= 1.   and   -1 <= y <= 1
00165   //
00166   //     F(x,y) = f(x)*g(x,y);   g(x,y) = 1.+g(x)*y
00167   // ***************************************************
00168 
00169   // ***** sampling F(x,y) directly (brute force) *****
00170 
00171   do{
00172 
00173     // Sample the positron energy by sampling from F
00174 
00175     rndm = G4UniformRand();
00176 
00177     x = x0 + rndm*(1.-x0);
00178 
00179     G4double x_squared = x*x;
00180 
00181     G4double F_IS, F_AS, G_IS, G_AS;
00182 
00183     F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
00184     F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
00185 
00186     G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
00187     G_IS = G_IS + michel_eta*(1.-x)*x0;
00188 
00189     G_AS = 3.*(michel_xsi-1.)*(1.-x);
00190     G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
00191     G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
00192 
00193     F_IS = F_IS + G_IS;
00194     F_AS = F_AS + G_AS;
00195 
00196     // *** Radiative Corrections ***
00197 
00198     G4double R_IS = F_c(x,x0);
00199 
00200     G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
00201 
00202     // *** Radiative Corrections ***
00203 
00204     G4double R_AS = F_theta(x,x0);
00205 
00206     rndm = G4UniformRand();
00207 
00208     ctheta = 2.*rndm-1.;
00209 
00210     G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
00211 
00212     FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
00213 
00214     if(FG>FG_max){
00215       G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
00216       FG_max = FG;
00217     }
00218 
00219     rndm = G4UniformRand();
00220 
00221   }while(FG<rndm*FG_max);
00222 
00223   G4double energy = x * W_mue;
00224 
00225   rndm = G4UniformRand();
00226 
00227   G4double phi = twopi * rndm;
00228 
00229   if(energy < EMASS) energy = EMASS;
00230 
00231   // calculate daughter momentum
00232   G4double daughtermomentum[3];
00233 
00234   daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
00235 
00236   G4double stheta = std::sqrt(1.-ctheta*ctheta);
00237   G4double cphi = std::cos(phi);
00238   G4double sphi = std::sin(phi);
00239 
00240   //Coordinates of the decay positron with respect to the muon spin
00241 
00242   G4double px = stheta*cphi;
00243   G4double py = stheta*sphi;
00244   G4double pz = ctheta;
00245 
00246   G4ThreeVector direction0(px,py,pz);
00247 
00248   direction0.rotateUz(parent_polarization);
00249 
00250   G4DynamicParticle * daughterparticle0 
00251     = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
00252 
00253   products->PushProducts(daughterparticle0);
00254 
00255 
00256   // daughter 1 ,2 (neutrinos)
00257   // create neutrinos in the C.M frame of two neutrinos
00258   G4double energy2 = parentmass*(1.0 - x/2.0); 
00259   G4double vmass   = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
00260   G4double beta = -1.0*daughtermomentum[0]/energy2;
00261   G4double costhetan = 2.*G4UniformRand()-1.0;
00262   G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00263   G4double phin  = twopi*G4UniformRand()*rad;
00264   G4double sinphin = std::sin(phin);
00265   G4double cosphin = std::cos(phin);
00266 
00267   G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
00268   G4DynamicParticle * daughterparticle1 
00269     = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
00270   G4DynamicParticle * daughterparticle2
00271     = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
00272 
00273   // boost to the muon rest frame
00274   G4LorentzVector p4;
00275   p4 = daughterparticle1->Get4Momentum();
00276   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
00277   daughterparticle1->Set4Momentum(p4);
00278   p4 = daughterparticle2->Get4Momentum();
00279   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
00280   daughterparticle2->Set4Momentum(p4);
00281   products->PushProducts(daughterparticle1);
00282   products->PushProducts(daughterparticle2);
00283   daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
00284   daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
00285 
00286   // output message
00287 #ifdef G4VERBOSE
00288   if (GetVerboseLevel()>1) {
00289     G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
00290     G4cout << "  create decay products in rest frame " <<G4endl;
00291     products->DumpInfo();
00292   }
00293 #endif
00294   return products;
00295 }

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

Definition at line 101 of file G4MuonDecayChannelWithSpin.hh.

00102 {
00103   return parent_polarization;
00104 }

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

Definition at line 80 of file G4MuonDecayChannelWithSpin.cc.

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

00081 {
00082   if (this != &right) { 
00083     kinematics_name = right.kinematics_name;
00084     verboseLevel = right.verboseLevel;
00085     rbranch = right.rbranch;
00086 
00087     // copy parent name
00088     parent_name = new G4String(*right.parent_name);
00089 
00090     // clear daughters_name array
00091     ClearDaughtersName();
00092 
00093     // recreate array
00094     numberOfDaughters = right.numberOfDaughters;
00095     if ( numberOfDaughters >0 ) {
00096       if (daughters_name !=0) ClearDaughtersName();
00097       daughters_name = new G4String*[numberOfDaughters];
00098       //copy daughters name
00099       for (G4int index=0; index < numberOfDaughters; index++) {
00100           daughters_name[index] = new G4String(*right.daughters_name[index]);
00101       }
00102     }
00103     parent_polarization = right.parent_polarization;
00104     EMMU  = right.EMMU;
00105     EMASS = right.EMASS;
00106   }
00107   return *this;
00108 }

void G4MuonDecayChannelWithSpin::SetPolarization ( G4ThreeVector   )  [inline]

Definition at line 96 of file G4MuonDecayChannelWithSpin.hh.

Referenced by G4DecayWithSpin::DecayIt().

00097 {
00098   parent_polarization = polar;
00099 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:32 2013 for Geant4 by  doxygen 1.4.7