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 10.6.p1)


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