Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/highenergy/src/G4AnnihiToMuPair.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/highenergy/src/G4AnnihiToMuPair.cc (Version 11.3.0) and /processes/electromagnetic/highenergy/src/G4AnnihiToMuPair.cc (Version 9.3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
                                                   >>  27 // $Id: G4AnnihiToMuPair.cc,v 1.6 2009/11/09 18:24:07 vnivanch Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-03 $
 27 //                                                 29 //
 28 //         ------------ G4AnnihiToMuPair physi     30 //         ------------ G4AnnihiToMuPair physics process ------
 29 //         by H.Burkhardt, S. Kelner and R. Ko     31 //         by H.Burkhardt, S. Kelner and R. Kokoulin, November 2002
 30 // -------------------------------------------     32 // -----------------------------------------------------------------------------
 31 //                                                 33 //
 32 //....oooOO0OOooo........oooOO0OOooo........oo     34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......//
 33 //                                                 35 //
 34 // 27.01.03 : first implementation (hbu)           36 // 27.01.03 : first implementation (hbu)
 35 // 04.02.03 : cosmetic simplifications (mma)       37 // 04.02.03 : cosmetic simplifications (mma)
 36 // 25.10.04 : migrade to new interfaces of Par     38 // 25.10.04 : migrade to new interfaces of ParticleChange (vi)
 37 // 28.02.18 : cross section now including SSS  << 
 38 //                                                 39 //
 39 //....oooOO0OOooo........oooOO0OOooo........oo     40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 40                                                    41 
 41 #include "G4AnnihiToMuPair.hh"                     42 #include "G4AnnihiToMuPair.hh"
 42                                                    43 
 43 #include "G4Exp.hh"                            << 
 44 #include "G4LossTableManager.hh"               << 
 45 #include "G4Material.hh"                       << 
 46 #include "G4MuonMinus.hh"                      << 
 47 #include "G4MuonPlus.hh"                       << 
 48 #include "G4PhysicalConstants.hh"              << 
 49 #include "G4Positron.hh"                       << 
 50 #include "G4Step.hh"                           << 
 51 #include "G4SystemOfUnits.hh"                  << 
 52 #include "G4TauMinus.hh"                       << 
 53 #include "G4TauPlus.hh"                        << 
 54 #include "G4ios.hh"                                44 #include "G4ios.hh"
 55 #include "Randomize.hh"                            45 #include "Randomize.hh"
 56                                                    46 
                                                   >>  47 #include "G4Positron.hh"
                                                   >>  48 #include "G4MuonPlus.hh"
                                                   >>  49 #include "G4MuonMinus.hh"
                                                   >>  50 #include "G4Material.hh"
                                                   >>  51 #include "G4Step.hh"
                                                   >>  52 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 58                                                    54 
                                                   >>  55 using namespace std;
                                                   >>  56 
 59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4Str     57 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4String& processName,
 60     G4ProcessType type):G4VDiscreteProcess (pr     58     G4ProcessType type):G4VDiscreteProcess (processName, type)
 61 {                                                  59 {
 62   //e+ Energy threshold                        <<  60  //e+ Energy threshold
 63   if(processName == "AnnihiToTauPair") {       <<  61  const G4double Mu_massc2 = G4MuonPlus::MuonPlus()->GetPDGMass();
 64     SetProcessSubType(fAnnihilationToTauTau);  <<  62  LowestEnergyLimit  = 2*Mu_massc2*Mu_massc2/electron_mass_c2 - electron_mass_c2;
 65     part1 = G4TauPlus::TauPlus();              <<  63  
 66     part2 = G4TauMinus::TauMinus();            <<  64  //modele ok up to 1000 TeV due to neglected Z-interference
 67     fInfo = "e+e->tau+tau-";                   <<  65  HighestEnergyLimit = 1000*TeV;
 68   } else {                                     <<  66  
 69     SetProcessSubType(fAnnihilationToMuMu);    <<  67  CurrentSigma = 0.0;
 70     part1 = G4MuonPlus::MuonPlus();            <<  68  CrossSecFactor = 1.;
 71     part2 = G4MuonMinus::MuonMinus();          <<  69  SetProcessSubType(6);
 72   }                                            << 
 73   fMass = part1->GetPDGMass();                 << 
 74   fLowEnergyLimit = 2. * fMass * fMass / CLHEP << 
 75                                                << 
 76   // model is ok up to 1000 TeV due to neglect << 
 77   fHighEnergyLimit = 1000. * TeV;              << 
 78                                                    70 
 79   fCurrentSigma = 0.0;                         << 
 80   fCrossSecFactor = 1.;                        << 
 81   fManager = G4LossTableManager::Instance();   << 
 82   fManager->Register(this);                    << 
 83 }                                                  71 }
 84                                                    72 
 85 //....oooOO0OOooo........oooOO0OOooo........oo     73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 86                                                    74 
 87 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empt     75 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor
 88 {                                              <<  76 { }
 89   fManager->DeRegister(this);                  << 
 90 }                                              << 
 91                                                    77 
 92 //....oooOO0OOooo........oooOO0OOooo........oo     78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 93                                                    79 
 94 G4bool G4AnnihiToMuPair::IsApplicable(const G4     80 G4bool G4AnnihiToMuPair::IsApplicable(const G4ParticleDefinition& particle)
 95 {                                                  81 {
 96   return ( &particle == G4Positron::Positron()     82   return ( &particle == G4Positron::Positron() );
 97 }                                                  83 }
 98                                                    84 
 99 //....oooOO0OOooo........oooOO0OOooo........oo     85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100                                                    86 
101 void G4AnnihiToMuPair::BuildPhysicsTable(const     87 void G4AnnihiToMuPair::BuildPhysicsTable(const G4ParticleDefinition&)
                                                   >>  88 // Build cross section and mean free path tables
                                                   >>  89 //here no tables, just calling PrintInfoDefinition
102 {                                                  90 {
                                                   >>  91   CurrentSigma = 0.0;
103   PrintInfoDefinition();                           92   PrintInfoDefinition();
104 }                                                  93 }
105                                                    94 
106 //....oooOO0OOooo........oooOO0OOooo........oo     95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107                                                    96 
108 void G4AnnihiToMuPair::SetCrossSecFactor(G4dou     97 void G4AnnihiToMuPair::SetCrossSecFactor(G4double fac)
109 // Set the factor to artificially increase the     98 // Set the factor to artificially increase the cross section
110 {                                                  99 { 
111   fCrossSecFactor = fac;                       << 100   CrossSecFactor = fac;
112   //G4cout << "The cross section for AnnihiToM << 101   G4cout << "The cross section for AnnihiToMuPair is artificially "
113   //       << "increased by the CrossSecFactor << 102          << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
114 }                                                 103 }
115                                                   104 
116 //....oooOO0OOooo........oooOO0OOooo........oo    105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117                                                   106 
118 G4double G4AnnihiToMuPair::ComputeCrossSection << 107 G4double G4AnnihiToMuPair::ComputeCrossSectionPerAtom(G4double Epos, G4double Z)
119 // Calculates the microscopic cross section in    108 // Calculates the microscopic cross section in GEANT4 internal units.
120 // It gives a good description from threshold     109 // It gives a good description from threshold to 1000 GeV
121 {                                                 110 {
122   G4double rmuon = CLHEP::elm_coupling/fMass;  << 111   static const G4double Mmuon = G4MuonPlus::MuonPlus()->GetPDGMass();
123   G4double sig0 = CLHEP::pi*rmuon*rmuon/3.;    << 112   static const G4double Rmuon = elm_coupling/Mmuon; //classical particle radius
124   const G4double pial = CLHEP::pi*CLHEP::fine_ << 113   static const G4double Sig0  = pi*Rmuon*Rmuon/3.;  //constant in crossSection
125                                                   114 
126   if (e <= fLowEnergyLimit) return 0.0;        << 115   G4double CrossSection = 0.;
                                                   >> 116   if (Epos < LowestEnergyLimit) return CrossSection;
127                                                   117    
128   const G4double xi = fLowEnergyLimit/e;       << 118   G4double xi = LowestEnergyLimit/Epos;
129   const G4double piaxi = pial * std::sqrt(xi); << 119   G4double SigmaEl = Sig0*xi*(1.+xi/2.)*sqrt(1.-xi); // per electron
130   G4double sigma = sig0 * xi * (1. + xi*0.5);  << 120   CrossSection = SigmaEl*Z;         // number of electrons per atom
131   //G4cout << "### xi= " << xi << " piaxi=" << << 121   return CrossSection;
132                                                << 
133   // argument of the exponent below 0.1 or abo << 
134   // Sigma per electron * number of electrons  << 
135   if(xi <= 1.0 - 100*piaxi*piaxi) {            << 
136     sigma *= std::sqrt(1.0 - xi);              << 
137   }                                            << 
138   else if (xi >= 1.0 - 0.01 * piaxi * piaxi) { << 
139     sigma *= piaxi;                            << 
140   }                                            << 
141   else {                                       << 
142     sigma *= piaxi / (1. - G4Exp(-piaxi / std: << 
143   }                                            << 
144   // G4cout << "### sigma= " << sigma << G4end << 
145   return sigma;                                << 
146 }                                                 122 }
147                                                   123 
148 //....oooOO0OOooo........oooOO0OOooo........oo    124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149                                                   125 
150 G4double G4AnnihiToMuPair::ComputeCrossSection << 126 G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double PositronEnergy, 
151                                                << 127              const G4Material* aMaterial)
152 {                                                 128 {
153   return ComputeCrossSectionPerElectron(energy << 129   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
154 }                                              << 130   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
155                                                   131 
156 //....oooOO0OOooo........oooOO0OOooo........oo << 132   G4double SIGMA = 0.0;
157                                                   133 
158 G4double G4AnnihiToMuPair::CrossSectionPerVolu << 134   for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; ++i )
159              const G4Material* aMaterial)      << 135   {
160 {                                              << 136     G4double AtomicZ = (*theElementVector)[i]->GetZ();
161   return ComputeCrossSectionPerElectron(energy << 137     SIGMA += NbOfAtomsPerVolume[i] *
                                                   >> 138       ComputeCrossSectionPerAtom(PositronEnergy,AtomicZ);
                                                   >> 139   }
                                                   >> 140   return SIGMA;
162 }                                                 141 }
163                                                   142 
164 //....oooOO0OOooo........oooOO0OOooo........oo    143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165                                                   144 
166 G4double G4AnnihiToMuPair::GetMeanFreePath(con    145 G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack,
167                                            G4d    146                                            G4double, G4ForceCondition*)
                                                   >> 147 
168 // returns the positron mean free path in GEAN    148 // returns the positron mean free path in GEANT4 internal units
                                                   >> 149 
169 {                                                 150 {
170   const G4DynamicParticle* aDynamicPositron =     151   const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
171   G4double energy = aDynamicPositron->GetTotal << 152   G4double PositronEnergy = aDynamicPositron->GetKineticEnergy()
172   const G4Material* aMaterial = aTrack.GetMate << 153                                               +electron_mass_c2;
173                                                << 154   G4Material* aMaterial = aTrack.GetMaterial();
174   // cross section before step                 << 155   CurrentSigma = CrossSectionPerVolume(PositronEnergy, aMaterial);
175   fCurrentSigma = CrossSectionPerVolume(energy << 
176                                                   156 
177   // increase the CrossSection by CrossSecFact    157   // increase the CrossSection by CrossSecFactor (default 1)
178   return (fCurrentSigma > 0.0) ? 1.0/(fCurrent << 158   G4double mfp = DBL_MAX;
                                                   >> 159   if(CurrentSigma > DBL_MIN) mfp = 1.0/(CurrentSigma*CrossSecFactor);
                                                   >> 160 
                                                   >> 161   return mfp;
179 }                                                 162 }
180                                                   163 
181 //....oooOO0OOooo........oooOO0OOooo........oo    164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182                                                   165 
183 G4VParticleChange* G4AnnihiToMuPair::PostStepD    166 G4VParticleChange* G4AnnihiToMuPair::PostStepDoIt(const G4Track& aTrack,
184                                                   167                                                   const G4Step&  aStep)
185 //                                                168 //
186 // generation of e+e- -> mu+mu-                   169 // generation of e+e- -> mu+mu-
187 //                                                170 //
188 {                                                 171 {
                                                   >> 172 
189   aParticleChange.Initialize(aTrack);             173   aParticleChange.Initialize(aTrack);
                                                   >> 174   static const G4double Mele=electron_mass_c2;
                                                   >> 175   static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
190                                                   176 
191   // current Positron energy and direction, re    177   // current Positron energy and direction, return if energy too low
192   const G4DynamicParticle *aDynamicPositron =     178   const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle();
193   const G4double Mele = CLHEP::electron_mass_c << 179   G4double Epos = aDynamicPositron->GetKineticEnergy() + Mele; 
194   G4double Epos = aDynamicPositron->GetTotalEn << 
195   G4double xs = CrossSectionPerVolume(Epos, aT << 
196                                                   180 
197   // test of cross section                        181   // test of cross section
198   if(xs > 0.0 && fCurrentSigma*G4UniformRand() << 182   if(CurrentSigma*G4UniformRand() > 
199     return G4VDiscreteProcess::PostStepDoIt(aT << 183      CrossSectionPerVolume(Epos, aTrack.GetMaterial())) 
                                                   >> 184     {
                                                   >> 185       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
                                                   >> 186     }
                                                   >> 187 
                                                   >> 188   if (Epos < LowestEnergyLimit) {
                                                   >> 189      return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
200   }                                               190   }
201                                                   191 
202   const G4ThreeVector PosiDirection = aDynamic << 192   G4ParticleMomentum PositronDirection = 
203   G4double xi = fLowEnergyLimit/Epos; // xi is << 193                                        aDynamicPositron->GetMomentumDirection();
204                                       // goes  << 194   G4double xi = LowestEnergyLimit/Epos; // xi is always less than 1,
                                                   >> 195                                         // goes to 0 at high Epos
205                                                   196 
206   // generate cost; probability function 1+cos << 197   // generate cost
207   //                                              198   //
208   G4double cost;                                  199   G4double cost;
209   do { cost = 2.*G4UniformRand()-1.; }         << 200   do cost = 2.*G4UniformRand()-1.;
210   // Loop checking, 07-Aug-2015, Vladimir Ivan << 
211   while (2.*G4UniformRand() > 1.+xi+cost*cost*    201   while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) ); 
212   G4double sint = std::sqrt(1.-cost*cost);     << 202                                                        //1+cost**2 at high Epos
                                                   >> 203   G4double sint = sqrt(1.-cost*cost);
213                                                   204 
214   // generate phi                                 205   // generate phi
215   //                                              206   //
216   G4double phi = 2.*CLHEP::pi*G4UniformRand(); << 207   G4double phi=2.*pi*G4UniformRand();
217                                                   208 
218   G4double Ecm   = std::sqrt(0.5*Mele*(Epos+Me << 209   G4double Ecm   = sqrt(0.5*Mele*(Epos+Mele));
219   G4double Pcm   = std::sqrt(Ecm*Ecm - fMass*f << 210   G4double Pcm   = sqrt(Ecm*Ecm-Mmuon*Mmuon);
220   G4double beta  = std::sqrt((Epos-Mele)/(Epos << 211   G4double beta  = sqrt((Epos-Mele)/(Epos+Mele));
221   G4double gamma = Ecm/Mele;                   << 212   G4double gamma = Ecm/Mele;                    // =sqrt((Epos+Mele)/(2.*Mele));
222   G4double Pt    = Pcm*sint;                      213   G4double Pt    = Pcm*sint;
223                                                   214   
224   // energy and momentum of the muons in the L    215   // energy and momentum of the muons in the Lab
225   //                                              216   //
226   G4double EmuPlus   = gamma*(Ecm + cost*beta* << 217   G4double EmuPlus   = gamma*(     Ecm+cost*beta*Pcm);
227   G4double EmuMinus  = gamma*(Ecm - cost*beta* << 218   G4double EmuMinus  = gamma*(     Ecm-cost*beta*Pcm);
228   G4double PmuPlusZ  = gamma*(beta*Ecm + cost* << 219   G4double PmuPlusZ  = gamma*(beta*Ecm+cost*     Pcm);
229   G4double PmuMinusZ = gamma*(beta*Ecm - cost* << 220   G4double PmuMinusZ = gamma*(beta*Ecm-cost*     Pcm);
230   G4double PmuPlusX  = Pt*std::cos(phi);       << 221   G4double PmuPlusX  = Pt*cos(phi);
231   G4double PmuPlusY  = Pt*std::sin(phi);       << 222   G4double PmuPlusY  = Pt*sin(phi);
232   G4double PmuMinusX =-PmuPlusX;               << 223   G4double PmuMinusX =-Pt*cos(phi);
233   G4double PmuMinusY =-PmuPlusY;               << 224   G4double PmuMinusY =-Pt*sin(phi);
234   // absolute momenta                             225   // absolute momenta
235   G4double PmuPlus  = std::sqrt(Pt*Pt+PmuPlusZ << 226   G4double PmuPlus  = sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ );
236   G4double PmuMinus = std::sqrt(Pt*Pt+PmuMinus << 227   G4double PmuMinus = sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ);
237                                                   228 
238   // mu+ mu- directions for Positron in z-dire    229   // mu+ mu- directions for Positron in z-direction
239   //                                              230   //
240   G4ThreeVector MuPlusDirection(PmuPlusX / Pmu << 231   G4ThreeVector
241   G4ThreeVector MuMinusDirection(PmuMinusX / P << 232     MuPlusDirection ( PmuPlusX/PmuPlus, PmuPlusY/PmuPlus,  PmuPlusZ/PmuPlus  );
                                                   >> 233   G4ThreeVector
                                                   >> 234     MuMinusDirection(PmuMinusX/PmuMinus,PmuMinusY/PmuMinus,PmuMinusZ/PmuMinus);
242                                                   235 
243   // rotate to actual Positron direction          236   // rotate to actual Positron direction
244   //                                              237   //
245   MuPlusDirection.rotateUz(PosiDirection);     << 238   MuPlusDirection.rotateUz(PositronDirection);
246   MuMinusDirection.rotateUz(PosiDirection);    << 239   MuMinusDirection.rotateUz(PositronDirection);
247                                                   240 
248   aParticleChange.SetNumberOfSecondaries(2);      241   aParticleChange.SetNumberOfSecondaries(2);
249                                                << 
250   // create G4DynamicParticle object for the p    242   // create G4DynamicParticle object for the particle1
251   auto aParticle1 = new G4DynamicParticle(part << 243   G4DynamicParticle* aParticle1= new G4DynamicParticle(
                                                   >> 244                          G4MuonPlus::MuonPlus(),MuPlusDirection,EmuPlus-Mmuon);
252   aParticleChange.AddSecondary(aParticle1);       245   aParticleChange.AddSecondary(aParticle1);
253   // create G4DynamicParticle object for the p    246   // create G4DynamicParticle object for the particle2
254   auto aParticle2 = new G4DynamicParticle(part << 247   G4DynamicParticle* aParticle2= new G4DynamicParticle(
                                                   >> 248                      G4MuonMinus::MuonMinus(),MuMinusDirection,EmuMinus-Mmuon);
255   aParticleChange.AddSecondary(aParticle2);       249   aParticleChange.AddSecondary(aParticle2);
256                                                   250 
257   // Kill the incident positron                   251   // Kill the incident positron 
258   //                                              252   //
259   aParticleChange.ProposeEnergy(0.);              253   aParticleChange.ProposeEnergy(0.); 
260   aParticleChange.ProposeTrackStatus(fStopAndK    254   aParticleChange.ProposeTrackStatus(fStopAndKill);
261                                                   255 
262   return &aParticleChange;                        256   return &aParticleChange;
263 }                                                 257 }
264                                                   258 
265 //....oooOO0OOooo........oooOO0OOooo........oo    259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266                                                   260 
267 void G4AnnihiToMuPair::PrintInfoDefinition()      261 void G4AnnihiToMuPair::PrintInfoDefinition()
268 {                                                 262 {
269   G4String comments = fInfo + " annihilation,  << 263   G4String comments ="e+e->mu+mu- annihilation, atomic e- at rest, SubType=.";
270   G4cout << G4endl << GetProcessName() << ":   << 264   G4cout << G4endl << GetProcessName() << ":  " << comments 
271   G4cout << "        threshold at " << fLowEne << 265    << GetProcessSubType() << G4endl;
272          << " good description up to " << fHig << 266   G4cout << "        threshold at " << LowestEnergyLimit/GeV << " GeV"
273          << G4endl;                            << 267          << " good description up to "
                                                   >> 268          << HighestEnergyLimit/TeV << " TeV for all Z." << G4endl;
274 }                                                 269 }
275                                                   270 
276 //....oooOO0OOooo........oooOO0OOooo........oo    271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277                                                   272