Geant4 Cross Reference

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


  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 // G4DalitzDecayChannel class implementation       26 // G4DalitzDecayChannel class implementation
 27 //                                                 27 //
 28 // Author: H.Kurashige, 30 May 1997            <<  28 // Author: H.Kurashige, 30 May 1997 
 29 // -------------------------------------------     29 // --------------------------------------------------------------------
 30                                                    30 
 31 #include "G4DalitzDecayChannel.hh"             << 
 32                                                << 
 33 #include "G4DecayProducts.hh"                  << 
 34 #include "G4LorentzRotation.hh"                << 
 35 #include "G4LorentzVector.hh"                  << 
 36 #include "G4ParticleDefinition.hh"             << 
 37 #include "G4PhaseSpaceDecayChannel.hh"         << 
 38 #include "G4PhysicalConstants.hh"                  31 #include "G4PhysicalConstants.hh"
 39 #include "G4SystemOfUnits.hh"                      32 #include "G4SystemOfUnits.hh"
                                                   >>  33 #include "G4ParticleDefinition.hh"
                                                   >>  34 #include "G4DecayProducts.hh"
 40 #include "G4VDecayChannel.hh"                      35 #include "G4VDecayChannel.hh"
                                                   >>  36 #include "G4DalitzDecayChannel.hh"
                                                   >>  37 #include "G4PhaseSpaceDecayChannel.hh"
 41 #include "Randomize.hh"                            38 #include "Randomize.hh"
                                                   >>  39 #include "G4LorentzVector.hh"
                                                   >>  40 #include "G4LorentzRotation.hh"
 42                                                    41 
 43 G4DalitzDecayChannel::G4DalitzDecayChannel(con <<  42 G4DalitzDecayChannel::G4DalitzDecayChannel()
 44                                            con <<  43   : G4VDecayChannel()
 45                                            con <<  44 {
                                                   >>  45 }
                                                   >>  46 
                                                   >>  47 G4DalitzDecayChannel::
                                                   >>  48 G4DalitzDecayChannel(const G4String& theParentName,
                                                   >>  49                            G4double        theBR,
                                                   >>  50                      const G4String& theLeptonName,
                                                   >>  51                      const G4String& theAntiLeptonName)
 46   : G4VDecayChannel("Dalitz Decay", 1)             52   : G4VDecayChannel("Dalitz Decay", 1)
 47 {                                                  53 {
 48   // set names for daughter particles              54   // set names for daughter particles
 49   SetParent(theParentName);                        55   SetParent(theParentName);
 50   SetBR(theBR);                                    56   SetBR(theBR);
 51   SetNumberOfDaughters(3);                         57   SetNumberOfDaughters(3);
 52   G4String gammaName = "gamma";                    58   G4String gammaName = "gamma";
 53   SetDaughter(idGamma, gammaName);                 59   SetDaughter(idGamma, gammaName);
 54   SetDaughter(idLepton, theLeptonName);            60   SetDaughter(idLepton, theLeptonName);
 55   SetDaughter(idAntiLepton, theAntiLeptonName)     61   SetDaughter(idAntiLepton, theAntiLeptonName);
 56 }                                                  62 }
 57                                                    63 
 58 G4DalitzDecayChannel& G4DalitzDecayChannel::op <<  64 G4DalitzDecayChannel::~G4DalitzDecayChannel()
 59 {                                                  65 {
 60   if (this != &right) {                        <<  66 }
                                                   >>  67 
                                                   >>  68 G4DalitzDecayChannel::G4DalitzDecayChannel(const G4DalitzDecayChannel& right)
                                                   >>  69   : G4VDecayChannel(right)
                                                   >>  70 {
                                                   >>  71 }
                                                   >>  72 
                                                   >>  73 G4DalitzDecayChannel&
                                                   >>  74 G4DalitzDecayChannel::operator=(const G4DalitzDecayChannel& right)
                                                   >>  75 {
                                                   >>  76   if (this != &right)
                                                   >>  77   { 
 61     kinematics_name = right.kinematics_name;       78     kinematics_name = right.kinematics_name;
 62     verboseLevel = right.verboseLevel;             79     verboseLevel = right.verboseLevel;
 63     rbranch = right.rbranch;                       80     rbranch = right.rbranch;
 64                                                    81 
 65     // copy parent name                            82     // copy parent name
 66     parent_name = new G4String(*right.parent_n     83     parent_name = new G4String(*right.parent_name);
 67                                                    84 
 68     // clear daughters_name array                  85     // clear daughters_name array
 69     ClearDaughtersName();                          86     ClearDaughtersName();
 70                                                    87 
 71     // recreate array                              88     // recreate array
 72     numberOfDaughters = right.numberOfDaughter     89     numberOfDaughters = right.numberOfDaughters;
 73     if (numberOfDaughters > 0) {               <<  90     if ( numberOfDaughters >0 )
                                                   >>  91     {
 74       if (daughters_name != nullptr) ClearDaug     92       if (daughters_name != nullptr) ClearDaughtersName();
 75       daughters_name = new G4String*[numberOfD     93       daughters_name = new G4String*[numberOfDaughters];
 76       // copy daughters name                   <<  94       //copy daughters name
 77       for (G4int index = 0; index < numberOfDa <<  95       for (G4int index=0; index < numberOfDaughters; ++index)
                                                   >>  96       {
 78         daughters_name[index] = new G4String(*     97         daughters_name[index] = new G4String(*right.daughters_name[index]);
 79       }                                            98       }
 80     }                                              99     }
 81   }                                               100   }
 82   return *this;                                   101   return *this;
 83 }                                                 102 }
 84                                                   103 
 85 G4DecayProducts* G4DalitzDecayChannel::DecayIt << 104 G4DecayProducts* G4DalitzDecayChannel::DecayIt(G4double) 
 86 {                                                 105 {
 87 #ifdef G4VERBOSE                                  106 #ifdef G4VERBOSE
 88   if (GetVerboseLevel() > 1) G4cout << "G4Dali << 107   if (GetVerboseLevel()>1) G4cout << "G4DalitzDecayChannel::DecayIt ";
 89 #endif                                         << 108 #endif 
 90   CheckAndFillParent();                           109   CheckAndFillParent();
 91   CheckAndFillDaughters();                        110   CheckAndFillDaughters();
 92                                                   111 
 93   // parent mass                                  112   // parent mass
 94   G4double parentmass = G4MT_parent->GetPDGMas    113   G4double parentmass = G4MT_parent->GetPDGMass();
 95                                                << 114  
 96   // create parent G4DynamicParticle at rest      115   // create parent G4DynamicParticle at rest
 97   G4ThreeVector dummy;                            116   G4ThreeVector dummy;
 98   auto parentparticle = new G4DynamicParticle( << 117   G4DynamicParticle* parentparticle
 99                                                << 118     = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
100   // daughters' mass                           << 119  
101   G4double leptonmass = G4MT_daughters[idLepto << 120   //daughters' mass
102                                                << 121   G4double leptonmass = G4MT_daughters[idLepton]->GetPDGMass(); 
103   // Generate t ( = std::exp(x):mass Square of << 122 
104   G4double xmin = 2.0 * std::log(2.0 * leptonm << 123   // Generate t ( = std::exp(x):mass Square of (l+ l-) system) 
105   G4double xmax = 2.0 * std::log(parentmass);  << 124   G4double xmin  = 2.0*std::log(2.0*leptonmass);
                                                   >> 125   G4double xmax  = 2.0*std::log(parentmass);
106   G4double wmax = 1.5;                            126   G4double wmax = 1.5;
107   G4double x, w, ww, w1, w2, w3, t;               127   G4double x, w, ww, w1, w2, w3, t;
108   const std::size_t MAX_LOOP = 10000;             128   const std::size_t MAX_LOOP = 10000;
109   for (std::size_t loop_counter = 0; loop_coun << 129   for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
110     x = G4UniformRand() * (xmax - xmin) + xmin << 130   {
111     w = G4UniformRand() * wmax;                << 131     x = G4UniformRand()*(xmax-xmin) + xmin;
                                                   >> 132     w = G4UniformRand()*wmax;
112     t = std::exp(x);                              133     t = std::exp(x);
113     w1 = (1.0 - 4.0 * leptonmass * leptonmass  << 134     w1 = (1.0-4.0*leptonmass*leptonmass/t);
114     if (w1 > 0.0) {                            << 135     if ( w1 > 0.0)
115       w2 = (1.0 + 2.0 * leptonmass * leptonmas << 136     {
116       w3 = (1.0 - t / parentmass / parentmass) << 137       w2 = ( 1.0 + 2.0*leptonmass*leptonmass/t );
                                                   >> 138       w3 = ( 1.0 - t/parentmass/parentmass );
117       w3 = w3 * w3 * w3;                          139       w3 = w3 * w3 * w3;
118       ww = w3 * w2 * std::sqrt(w1);               140       ww = w3 * w2 * std::sqrt(w1);
119     }                                             141     }
120     else {                                     << 142     else
                                                   >> 143     {
121       ww = 0.0;                                   144       ww = 0.0;
122     }                                             145     }
123     if (w <= ww) break;                        << 146     if (w <= ww) break;    
124   }                                               147   }
125                                                   148 
126   // calculate gamma momentum                     149   // calculate gamma momentum
127   G4double Pgamma = G4PhaseSpaceDecayChannel:: << 150   G4double Pgamma = 
128   G4double costheta = 2. * G4UniformRand() - 1 << 151       G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, std::sqrt(t)); 
129   G4double sintheta = std::sqrt((1.0 - costhet << 152   G4double costheta = 2.*G4UniformRand()-1.0;
130   G4double phi = twopi * G4UniformRand() * rad << 153   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
131   G4ThreeVector gdirection(sintheta * std::cos << 154   G4double phi  = twopi*G4UniformRand()*rad;
132                                                << 155   G4ThreeVector gdirection(sintheta*std::cos(phi),
133   // create G4DynamicParticle for gamma        << 156                            sintheta*std::sin(phi),
134   auto gammaparticle = new G4DynamicParticle(G << 157                            costheta);
                                                   >> 158 
                                                   >> 159   // create G4DynamicParticle for gamma 
                                                   >> 160   G4DynamicParticle* gammaparticle
                                                   >> 161     = new G4DynamicParticle(G4MT_daughters[idGamma] , gdirection, Pgamma);
135                                                   162 
136   // calculate beta of (l+ l-)system              163   // calculate beta of (l+ l-)system
137   G4double beta = Pgamma / (parentmass - Pgamm << 164   G4double beta = Pgamma/(parentmass-Pgamma);
138                                                   165 
139   // calculate lepton momentum in the rest fra    166   // calculate lepton momentum in the rest frame of (l+ l-)system
140   G4double Plepton = G4PhaseSpaceDecayChannel: << 167   G4double Plepton = 
141   G4double Elepton = std::sqrt(Plepton * Plept << 168       G4PhaseSpaceDecayChannel::Pmx(std::sqrt(t),leptonmass, leptonmass); 
142   costheta = 2. * G4UniformRand() - 1.0;       << 169   G4double Elepton = std::sqrt(Plepton*Plepton + leptonmass*leptonmass );
143   sintheta = std::sqrt((1.0 - costheta) * (1.0 << 170   costheta = 2.*G4UniformRand()-1.0;
144   phi = twopi * G4UniformRand() * rad;         << 171   sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
145   G4ThreeVector ldirection(sintheta * std::cos << 172   phi = twopi*G4UniformRand()*rad;
                                                   >> 173   G4ThreeVector ldirection(sintheta*std::cos(phi),
                                                   >> 174                            sintheta*std::sin(phi),
                                                   >> 175                            costheta);
146   // create G4DynamicParticle for leptons  in     176   // create G4DynamicParticle for leptons  in the rest frame of (l+ l-)system
147   auto leptonparticle =                        << 177   G4DynamicParticle* leptonparticle 
148     new G4DynamicParticle(G4MT_daughters[idLep << 178     = new G4DynamicParticle(G4MT_daughters[idLepton] , 
149   auto antileptonparticle =                    << 179                             ldirection, Elepton-leptonmass );
150     new G4DynamicParticle(G4MT_daughters[idAnt << 180   G4DynamicParticle* antileptonparticle 
151   // boost leptons in the rest frame of the pa << 181     = new G4DynamicParticle(G4MT_daughters[idAntiLepton] , 
                                                   >> 182                             -1.0*ldirection, Elepton-leptonmass );
                                                   >> 183   //boost leptons in the rest frame of the parent 
152   G4LorentzVector p4 = leptonparticle->Get4Mom    184   G4LorentzVector p4 = leptonparticle->Get4Momentum();
153   p4.boost(-1.0 * gdirection.x() * beta, -1.0  << 185   p4.boost( -1.0*gdirection.x()*beta,
154            -1.0 * gdirection.z() * beta);      << 186             -1.0*gdirection.y()*beta,
                                                   >> 187             -1.0*gdirection.z()*beta);
155   leptonparticle->Set4Momentum(p4);               188   leptonparticle->Set4Momentum(p4);
156   p4 = antileptonparticle->Get4Momentum();        189   p4 = antileptonparticle->Get4Momentum();
157   p4.boost(-1.0 * gdirection.x() * beta, -1.0  << 190   p4.boost( -1.0*gdirection.x()*beta,
158            -1.0 * gdirection.z() * beta);      << 191             -1.0*gdirection.y()*beta,
                                                   >> 192             -1.0*gdirection.z()*beta);
159   antileptonparticle->Set4Momentum(p4);           193   antileptonparticle->Set4Momentum(p4);
160                                                   194 
161   // create G4Decayproducts                    << 195   //create G4Decayproducts
162   auto products = new G4DecayProducts(*parentp << 196   G4DecayProducts* products = new G4DecayProducts(*parentparticle);
163   delete parentparticle;                          197   delete parentparticle;
164   products->PushProducts(gammaparticle);          198   products->PushProducts(gammaparticle);
165   products->PushProducts(leptonparticle);         199   products->PushProducts(leptonparticle);
166   products->PushProducts(antileptonparticle);     200   products->PushProducts(antileptonparticle);
167                                                   201 
168 #ifdef G4VERBOSE                                  202 #ifdef G4VERBOSE
169   if (GetVerboseLevel() > 1) {                 << 203   if (GetVerboseLevel()>1)
170     G4cout << "G4DalitzDecayChannel::DecayIt " << 204   {
171     G4cout << "  create decay products in rest << 205      G4cout << "G4DalitzDecayChannel::DecayIt ";
172     products->DumpInfo();                      << 206      G4cout << "  create decay products in rest frame " <<G4endl;
                                                   >> 207      products->DumpInfo();
173   }                                               208   }
174 #endif                                            209 #endif
175   return products;                                210   return products;
176 }                                                 211 }
177                                                   212