Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/particles/management/src/G4TauLeptonicDecayChannel.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 /particles/management/src/G4TauLeptonicDecayChannel.cc (Version 11.3.0) and /particles/management/src/G4TauLeptonicDecayChannel.cc (Version 9.4)


  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 // G4TauLeptonicDecayChannel class implementat << 
 27 //                                                 26 //
 28 // Author: H.Kurashige, 30 May 1997            <<  27 // $Id: G4TauLeptonicDecayChannel.cc,v 1.5 2006/06/29 19:26:18 gunter Exp $
 29 // ------------------------------------------- <<  28 // GEANT4 tag $Name: geant4-09-02 $
 30                                                <<  29 //
 31 #include "G4TauLeptonicDecayChannel.hh"        <<  30 // 
                                                   >>  31 // ------------------------------------------------------------
                                                   >>  32 //      GEANT 4 class header file
                                                   >>  33 //
                                                   >>  34 //      History: first implementation, based on object model of
                                                   >>  35 //      30 May  1997 H.Kurashige
                                                   >>  36 //
                                                   >>  37 //      Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige 
                                                   >>  38 // ------------------------------------------------------------
 32                                                    39 
 33 #include "G4DecayProducts.hh"                  << 
 34 #include "G4LorentzRotation.hh"                << 
 35 #include "G4LorentzVector.hh"                  << 
 36 #include "G4ParticleDefinition.hh"                 40 #include "G4ParticleDefinition.hh"
 37 #include "G4PhysicalConstants.hh"              <<  41 #include "G4DecayProducts.hh"
 38 #include "G4SystemOfUnits.hh"                  << 
 39 #include "G4VDecayChannel.hh"                      42 #include "G4VDecayChannel.hh"
                                                   >>  43 #include "G4TauLeptonicDecayChannel.hh"
 40 #include "Randomize.hh"                            44 #include "Randomize.hh"
                                                   >>  45 #include "G4LorentzVector.hh"
                                                   >>  46 #include "G4LorentzRotation.hh"
 41                                                    47 
 42 G4TauLeptonicDecayChannel::G4TauLeptonicDecayC <<  48 
 43                                                <<  49 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel(
 44   : G4VDecayChannel("Tau Leptonic Decay", 1)   <<  50                  const G4String& theParentName, 
                                                   >>  51                  G4double        theBR,
                                                   >>  52                  const G4String& theLeptonName)
                                                   >>  53                    :G4VDecayChannel("Tau Leptonic Decay",1)
 45 {                                                  54 {
 46   // set names for daughter particles              55   // set names for daughter particles
 47   if (theParentName == "tau+") {                   56   if (theParentName == "tau+") {
 48     SetBR(theBR);                                  57     SetBR(theBR);
 49     SetParent("tau+");                             58     SetParent("tau+");
 50     SetNumberOfDaughters(3);                       59     SetNumberOfDaughters(3);
 51     if ((theLeptonName == "e-" || theLeptonNam <<  60     if ((theLeptonName=="e-"||theLeptonName=="e+")){
 52       SetDaughter(0, "e+");                        61       SetDaughter(0, "e+");
 53       SetDaughter(1, "nu_e");                      62       SetDaughter(1, "nu_e");
 54       SetDaughter(2, "anti_nu_tau");               63       SetDaughter(2, "anti_nu_tau");
 55     }                                          <<  64     } else {
 56     else {                                     << 
 57       SetDaughter(0, "mu+");                       65       SetDaughter(0, "mu+");
 58       SetDaughter(1, "nu_mu");                     66       SetDaughter(1, "nu_mu");
 59       SetDaughter(2, "anti_nu_tau");               67       SetDaughter(2, "anti_nu_tau");
 60     }                                          <<  68     } 
 61   }                                            <<  69   } else if (theParentName == "tau-") {
 62   else if (theParentName == "tau-") {          << 
 63     SetBR(theBR);                                  70     SetBR(theBR);
 64     SetParent("tau-");                             71     SetParent("tau-");
 65     SetNumberOfDaughters(3);                       72     SetNumberOfDaughters(3);
 66     if ((theLeptonName == "e-" || theLeptonNam <<  73     if ((theLeptonName=="e-"||theLeptonName=="e+")){
 67       SetDaughter(0, "e-");                        74       SetDaughter(0, "e-");
 68       SetDaughter(1, "anti_nu_e");                 75       SetDaughter(1, "anti_nu_e");
 69       SetDaughter(2, "nu_tau");                    76       SetDaughter(2, "nu_tau");
 70     }                                          <<  77     } else {
 71     else {                                     << 
 72       SetDaughter(0, "mu-");                       78       SetDaughter(0, "mu-");
 73       SetDaughter(1, "anti_nu_mu");                79       SetDaughter(1, "anti_nu_mu");
 74       SetDaughter(2, "nu_tau");                    80       SetDaughter(2, "nu_tau");
 75     }                                          <<  81     } 
 76   }                                            <<  82   } else {
 77   else {                                       << 
 78 #ifdef G4VERBOSE                                   83 #ifdef G4VERBOSE
 79     if (GetVerboseLevel() > 0) {               <<  84     if (GetVerboseLevel()>0) {
 80       G4cout << "G4TauLeptonicDecayChannel:: c     85       G4cout << "G4TauLeptonicDecayChannel:: constructor :";
 81       G4cout << " parent particle is not tau b     86       G4cout << " parent particle is not tau but ";
 82       G4cout << theParentName << G4endl;           87       G4cout << theParentName << G4endl;
 83     }                                              88     }
 84 #endif                                             89 #endif
 85   }                                                90   }
 86 }                                                  91 }
 87                                                    92 
 88 G4TauLeptonicDecayChannel&                     <<  93 G4TauLeptonicDecayChannel::~G4TauLeptonicDecayChannel()
 89 G4TauLeptonicDecayChannel::operator=(const G4T << 
 90 {                                                  94 {
 91   if (this != &right) {                        << 
 92     kinematics_name = right.kinematics_name;   << 
 93     verboseLevel = right.verboseLevel;         << 
 94     rbranch = right.rbranch;                   << 
 95                                                << 
 96     // copy parent name                        << 
 97     parent_name = new G4String(*right.parent_n << 
 98                                                << 
 99     // clear daughters_name array              << 
100     ClearDaughtersName();                      << 
101                                                << 
102     // recreate array                          << 
103     numberOfDaughters = right.numberOfDaughter << 
104     if (numberOfDaughters > 0) {               << 
105       if (daughters_name != nullptr) ClearDaug << 
106       daughters_name = new G4String*[numberOfD << 
107       // copy daughters name                   << 
108       for (G4int index = 0; index < numberOfDa << 
109         daughters_name[index] = new G4String(* << 
110       }                                        << 
111     }                                          << 
112   }                                            << 
113   return *this;                                << 
114 }                                                  95 }
115                                                    96 
116 G4DecayProducts* G4TauLeptonicDecayChannel::De <<  97 G4DecayProducts *G4TauLeptonicDecayChannel::DecayIt(G4double) 
117 {                                                  98 {
118   // this version neglects muon polarization   <<  99   // this version neglects muon polarization 
119   //              assumes the pure V-A couplin    100   //              assumes the pure V-A coupling
120   //              gives incorrect energy spect    101   //              gives incorrect energy spectrum for neutrinos
121                                                << 
122 #ifdef G4VERBOSE                                  102 #ifdef G4VERBOSE
123   if (GetVerboseLevel() > 1) G4cout << "G4TauL << 103   if (GetVerboseLevel()>1) G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
124 #endif                                            104 #endif
125                                                   105 
126   CheckAndFillParent();                        << 106   if (parent == 0) FillParent();  
127   CheckAndFillDaughters();                     << 107   if (daughters == 0) FillDaughters();
128                                                << 108  
129   // parent mass                                  109   // parent mass
130   G4double parentmass = G4MT_parent->GetPDGMas << 110   G4double parentmass = parent->GetPDGMass();
131                                                   111 
132   // daughters'mass                            << 112   //daughters'mass
133   const G4int N_DAUGHTER = 3;                  << 113   G4double daughtermass[3]; 
134   G4double daughtermass[N_DAUGHTER];           << 114   for (G4int index=0; index<3; index++){
135   for (G4int index = 0; index < N_DAUGHTER; ++ << 115     daughtermass[index] = daughters[index]->GetPDGMass();
136     daughtermass[index] = G4MT_daughters[index << 
137   }                                               116   }
138                                                   117 
139   // create parent G4DynamicParticle at rest   << 118    //create parent G4DynamicParticle at rest
140   G4ThreeVector dummy;                            119   G4ThreeVector dummy;
141   auto parentparticle = new G4DynamicParticle( << 120   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
142   // create G4Decayproducts                    << 121   //create G4Decayproducts
143   auto products = new G4DecayProducts(*parentp << 122   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
144   delete parentparticle;                          123   delete parentparticle;
145                                                   124 
146   // calculate daughter momentum                  125   // calculate daughter momentum
147   G4double daughtermomentum[N_DAUGHTER];       << 126   G4double daughtermomentum[3];
                                                   >> 127   G4double energy;
148                                                   128 
149   // calculate lepton momentum                 << 129   // calcurate lepton momentum
150   G4double pmax = (parentmass * parentmass - d << 130   G4double pmax = (parentmass*parentmass-daughtermass[0]*daughtermass[0])/2./parentmass;
151   G4double p, e;                                  131   G4double p, e;
152   G4double r;                                     132   G4double r;
153   const std::size_t MAX_LOOP = 10000;          << 133   do {
154   for (std::size_t loop_counter = 0; loop_coun << 
155     // determine momentum/energy                  134     // determine momentum/energy
156     r = G4UniformRand();                          135     r = G4UniformRand();
157     p = pmax * G4UniformRand();                << 136     p = pmax*G4UniformRand();
158     e = std::sqrt(p * p + daughtermass[0] * da << 137     e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]);
159     if (r < spectrum(p, e, parentmass, daughte << 138   } while (r > spectrum(p,e,parentmass,daughtermass[0]) );
160   }                                            << 139   energy = e- daughtermass[0];
161                                                   140 
162   // create daughter G4DynamicParticle         << 141   //create daughter G4DynamicParticle 
163   // daughter 0 (lepton)                          142   // daughter 0 (lepton)
164   daughtermomentum[0] = p;                        143   daughtermomentum[0] = p;
165   G4double costheta, sintheta, phi, sinphi, co << 144   G4double costheta, sintheta, phi, sinphi, cosphi; 
166   costheta = 2. * G4UniformRand() - 1.0;       << 145   costheta = 2.*G4UniformRand()-1.0;
167   sintheta = std::sqrt((1.0 - costheta) * (1.0 << 146   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
168   phi = twopi * G4UniformRand() * rad;         << 147   phi  = twopi*G4UniformRand()*rad;
169   sinphi = std::sin(phi);                         148   sinphi = std::sin(phi);
170   cosphi = std::cos(phi);                         149   cosphi = std::cos(phi);
171   G4ThreeVector direction0(sintheta * cosphi,  << 150   G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
172   auto daughterparticle =                      << 151   G4DynamicParticle * daughterparticle 
173     new G4DynamicParticle(G4MT_daughters[0], d << 152          = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
174   products->PushProducts(daughterparticle);       153   products->PushProducts(daughterparticle);
175                                                   154 
176   // daughter 1 ,2 (nutrinos)                     155   // daughter 1 ,2 (nutrinos)
177   // create neutrinos in the C.M frame of two     156   // create neutrinos in the C.M frame of two neutrinos
178   G4double energy2 = parentmass - e;           << 157   G4double energy2 = parentmass-e; 
179   G4double vmass = std::sqrt((energy2 - daught << 158   G4double vmass   = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
180   G4double beta = -1.0 * daughtermomentum[0] / << 159   G4double beta = -1.0*daughtermomentum[0]/energy2;
181   G4double costhetan = 2. * G4UniformRand() -  << 160   G4double costhetan = 2.*G4UniformRand()-1.0;
182   G4double sinthetan = std::sqrt((1.0 - costhe << 161   G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
183   G4double phin = twopi * G4UniformRand() * ra << 162   G4double phin  = twopi*G4UniformRand()*rad;
184   G4double sinphin = std::sin(phin);              163   G4double sinphin = std::sin(phin);
185   G4double cosphin = std::cos(phin);              164   G4double cosphin = std::cos(phin);
186                                                   165 
187   G4ThreeVector direction1(sinthetan * cosphin << 166   G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
188   auto daughterparticle1 = new G4DynamicPartic << 167   G4DynamicParticle * daughterparticle1 
189   auto daughterparticle2 =                     << 168          = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
190     new G4DynamicParticle(G4MT_daughters[2], d << 169   G4DynamicParticle * daughterparticle2
                                                   >> 170          = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
191                                                   171 
192   // boost to the muon rest frame                 172   // boost to the muon rest frame
193   G4LorentzVector p4;                             173   G4LorentzVector p4;
194   p4 = daughterparticle1->Get4Momentum();         174   p4 = daughterparticle1->Get4Momentum();
195   p4.boost(direction0.x() * beta, direction0.y << 175   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
196   daughterparticle1->Set4Momentum(p4);            176   daughterparticle1->Set4Momentum(p4);
197   p4 = daughterparticle2->Get4Momentum();         177   p4 = daughterparticle2->Get4Momentum();
198   p4.boost(direction0.x() * beta, direction0.y << 178   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
199   daughterparticle2->Set4Momentum(p4);            179   daughterparticle2->Set4Momentum(p4);
200   products->PushProducts(daughterparticle1);      180   products->PushProducts(daughterparticle1);
201   products->PushProducts(daughterparticle2);      181   products->PushProducts(daughterparticle2);
202   daughtermomentum[1] = daughterparticle1->Get    182   daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
203   daughtermomentum[2] = daughterparticle2->Get    183   daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
204                                                   184 
205   // output message                            << 185 
                                                   >> 186  // output message
206 #ifdef G4VERBOSE                                  187 #ifdef G4VERBOSE
207   if (GetVerboseLevel() > 1) {                 << 188   if (GetVerboseLevel()>1) {
208     G4cout << "G4TauLeptonicDecayChannel::Deca    189     G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
209     G4cout << "  create decay products in rest << 190     G4cout << "  create decay products in rest frame " <<G4endl;
210     products->DumpInfo();                         191     products->DumpInfo();
211   }                                               192   }
212 #endif                                            193 #endif
213   return products;                                194   return products;
214 }                                                 195 }
215                                                   196 
216 G4double G4TauLeptonicDecayChannel::spectrum(G << 197 
                                                   >> 198 
                                                   >> 199 
                                                   >> 200 G4double G4TauLeptonicDecayChannel::spectrum(G4double p,
                                                   >> 201                G4double e,
                                                   >> 202                G4double mtau,
                                                   >> 203                G4double ml)
217 {                                                 204 {
218   G4double f1;                                    205   G4double f1;
219   f1 = 3.0 * e * (mtau * mtau + ml * ml) - 4.0 << 206   f1 = 3.0*e*(mtau*mtau+ml*ml)-4.0*mtau*e*e-2.0*mtau*ml*ml;
220   return p * (f1) / (mtau * mtau * mtau * mtau << 207   return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6);
221 }                                                 208 }
                                                   >> 209 
                                                   >> 210 
                                                   >> 211 
222                                                   212