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.p4)


  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: not supported by cvs2svn $
 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];
148                                                   127 
149   // calculate lepton momentum                 << 128   // calcurate lepton momentum
150   G4double pmax = (parentmass * parentmass - d << 129   G4double pmax = (parentmass*parentmass-daughtermass[0]*daughtermass[0])/2./parentmass;
151   G4double p, e;                                  130   G4double p, e;
152   G4double r;                                     131   G4double r;
153   const std::size_t MAX_LOOP = 10000;          << 132   do {
154   for (std::size_t loop_counter = 0; loop_coun << 
155     // determine momentum/energy                  133     // determine momentum/energy
156     r = G4UniformRand();                          134     r = G4UniformRand();
157     p = pmax * G4UniformRand();                << 135     p = pmax*G4UniformRand();
158     e = std::sqrt(p * p + daughtermass[0] * da << 136     e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]);
159     if (r < spectrum(p, e, parentmass, daughte << 137   } while (r > spectrum(p,e,parentmass,daughtermass[0]) );
160   }                                            << 
161                                                   138 
162   // create daughter G4DynamicParticle         << 139   //create daughter G4DynamicParticle 
163   // daughter 0 (lepton)                          140   // daughter 0 (lepton)
164   daughtermomentum[0] = p;                        141   daughtermomentum[0] = p;
165   G4double costheta, sintheta, phi, sinphi, co << 142   G4double costheta, sintheta, phi, sinphi, cosphi; 
166   costheta = 2. * G4UniformRand() - 1.0;       << 143   costheta = 2.*G4UniformRand()-1.0;
167   sintheta = std::sqrt((1.0 - costheta) * (1.0 << 144   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
168   phi = twopi * G4UniformRand() * rad;         << 145   phi  = twopi*G4UniformRand()*rad;
169   sinphi = std::sin(phi);                         146   sinphi = std::sin(phi);
170   cosphi = std::cos(phi);                         147   cosphi = std::cos(phi);
171   G4ThreeVector direction0(sintheta * cosphi,  << 148   G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
172   auto daughterparticle =                      << 149   G4DynamicParticle * daughterparticle 
173     new G4DynamicParticle(G4MT_daughters[0], d << 150          = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
174   products->PushProducts(daughterparticle);       151   products->PushProducts(daughterparticle);
175                                                   152 
176   // daughter 1 ,2 (nutrinos)                     153   // daughter 1 ,2 (nutrinos)
177   // create neutrinos in the C.M frame of two     154   // create neutrinos in the C.M frame of two neutrinos
178   G4double energy2 = parentmass - e;           << 155   G4double energy2 = parentmass-e; 
179   G4double vmass = std::sqrt((energy2 - daught << 156   G4double vmass   = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
180   G4double beta = -1.0 * daughtermomentum[0] / << 157   G4double beta = -1.0*daughtermomentum[0]/energy2;
181   G4double costhetan = 2. * G4UniformRand() -  << 158   G4double costhetan = 2.*G4UniformRand()-1.0;
182   G4double sinthetan = std::sqrt((1.0 - costhe << 159   G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
183   G4double phin = twopi * G4UniformRand() * ra << 160   G4double phin  = twopi*G4UniformRand()*rad;
184   G4double sinphin = std::sin(phin);              161   G4double sinphin = std::sin(phin);
185   G4double cosphin = std::cos(phin);              162   G4double cosphin = std::cos(phin);
186                                                   163 
187   G4ThreeVector direction1(sinthetan * cosphin << 164   G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
188   auto daughterparticle1 = new G4DynamicPartic << 165   G4DynamicParticle * daughterparticle1 
189   auto daughterparticle2 =                     << 166          = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
190     new G4DynamicParticle(G4MT_daughters[2], d << 167   G4DynamicParticle * daughterparticle2
                                                   >> 168          = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
191                                                   169 
192   // boost to the muon rest frame                 170   // boost to the muon rest frame
193   G4LorentzVector p4;                             171   G4LorentzVector p4;
194   p4 = daughterparticle1->Get4Momentum();         172   p4 = daughterparticle1->Get4Momentum();
195   p4.boost(direction0.x() * beta, direction0.y << 173   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
196   daughterparticle1->Set4Momentum(p4);            174   daughterparticle1->Set4Momentum(p4);
197   p4 = daughterparticle2->Get4Momentum();         175   p4 = daughterparticle2->Get4Momentum();
198   p4.boost(direction0.x() * beta, direction0.y << 176   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
199   daughterparticle2->Set4Momentum(p4);            177   daughterparticle2->Set4Momentum(p4);
200   products->PushProducts(daughterparticle1);      178   products->PushProducts(daughterparticle1);
201   products->PushProducts(daughterparticle2);      179   products->PushProducts(daughterparticle2);
202   daughtermomentum[1] = daughterparticle1->Get    180   daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
203   daughtermomentum[2] = daughterparticle2->Get    181   daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
204                                                   182 
205   // output message                            << 183 
                                                   >> 184  // output message
206 #ifdef G4VERBOSE                                  185 #ifdef G4VERBOSE
207   if (GetVerboseLevel() > 1) {                 << 186   if (GetVerboseLevel()>1) {
208     G4cout << "G4TauLeptonicDecayChannel::Deca    187     G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
209     G4cout << "  create decay products in rest << 188     G4cout << "  create decay products in rest frame " <<G4endl;
210     products->DumpInfo();                         189     products->DumpInfo();
211   }                                               190   }
212 #endif                                            191 #endif
213   return products;                                192   return products;
214 }                                                 193 }
215                                                   194 
216 G4double G4TauLeptonicDecayChannel::spectrum(G << 195 
                                                   >> 196 
                                                   >> 197 
                                                   >> 198 G4double G4TauLeptonicDecayChannel::spectrum(G4double p,
                                                   >> 199                G4double e,
                                                   >> 200                G4double mtau,
                                                   >> 201                G4double ml)
217 {                                                 202 {
218   G4double f1;                                    203   G4double f1;
219   f1 = 3.0 * e * (mtau * mtau + ml * ml) - 4.0 << 204   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 << 205   return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6);
221 }                                                 206 }
                                                   >> 207 
                                                   >> 208 
                                                   >> 209 
222                                                   210