Geant4  10.03
G4DalitzDecayChannel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4DalitzDecayChannel.cc 95906 2016-03-02 10:56:50Z gcosmo $
28 //
29 //
30 // ------------------------------------------------------------
31 // GEANT 4 class header file
32 //
33 // History: first implementation, based on object model of
34 // 30 May 1997 H.Kurashige
35 // ------------------------------------------------------------
36 
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4DecayProducts.hh"
41 #include "G4VDecayChannel.hh"
42 #include "G4DalitzDecayChannel.hh"
44 #include "Randomize.hh"
45 #include "G4LorentzVector.hh"
46 #include "G4LorentzRotation.hh"
47 
50 {
51 }
52 
54  const G4String& theParentName,
55  G4double theBR,
56  const G4String& theLeptonName,
57  const G4String& theAntiLeptonName)
58  :G4VDecayChannel("Dalitz Decay",1)
59 {
60  // set names for daughter particles
61  SetParent(theParentName);
62  SetBR(theBR);
64  G4String gammaName = "gamma";
65  SetDaughter(idGamma, gammaName);
66  SetDaughter(idLepton, theLeptonName);
67  SetDaughter(idAntiLepton, theAntiLeptonName);
68 }
69 
71 {
72 }
73 
75  G4VDecayChannel(right)
76 {
77 }
78 
80 {
81  if (this != &right) {
83  verboseLevel = right.verboseLevel;
84  rbranch = right.rbranch;
85 
86  // copy parent name
87  parent_name = new G4String(*right.parent_name);
88 
89  // clear daughters_name array
91 
92  // recreate array
94  if ( numberOfDaughters >0 ) {
97  //copy daughters name
98  for (G4int index=0; index < numberOfDaughters; index++) {
99  daughters_name[index] = new G4String(*right.daughters_name[index]);
100  }
101  }
102  }
103  return *this;
104 }
105 
107 {
108 #ifdef G4VERBOSE
109  if (GetVerboseLevel()>1) G4cout << "G4DalitzDecayChannel::DecayIt ";
110 #endif
113 
114  // parent mass
115  G4double parentmass = G4MT_parent->GetPDGMass();
116 
117  //create parent G4DynamicParticle at rest
118  G4ThreeVector dummy;
119  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
120 
121  //daughters'mass
122  G4double leptonmass = G4MT_daughters[idLepton]->GetPDGMass();
123 
124  // Generate t ( = std::exp(x):mass Square of (l+ l-) system)
125  G4double xmin = 2.0*std::log(2.0*leptonmass);
126  G4double xmax = 2.0*std::log(parentmass);
127  G4double wmax = 1.5;
128  G4double x, w, ww, w1, w2, w3, t;
129  const size_t MAX_LOOP = 10000;
130  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
131  x = G4UniformRand()*(xmax-xmin) + xmin;
132  w = G4UniformRand()*wmax;
133  t = std::exp(x);
134  w1 = (1.0-4.0*leptonmass*leptonmass/t);
135  if ( w1 > 0.0) {
136  w2 = ( 1.0 + 2.0*leptonmass*leptonmass/t );
137  w3 = ( 1.0 - t/parentmass/parentmass );
138  w3 = w3 * w3 * w3;
139  ww = w3 * w2 * std::sqrt(w1);
140  } else {
141  ww = 0.0;
142  }
143  if (w <= ww) break;
144  }
145 
146  // calculate gamma momentum
147  G4double Pgamma =
148  G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, std::sqrt(t));
149  G4double costheta = 2.*G4UniformRand()-1.0;
150  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
152  G4ThreeVector gdirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
153 
154  //create G4DynamicParticle for gamma
155  G4DynamicParticle * gammaparticle
156  = new G4DynamicParticle(G4MT_daughters[idGamma] , gdirection, Pgamma);
157 
158  // calcurate beta of (l+ l-)system
159  G4double beta = Pgamma/(parentmass-Pgamma);
160 
161  // calculate lepton momentum in the rest frame of (l+ l-)system
162  G4double Plepton =
163  G4PhaseSpaceDecayChannel::Pmx(std::sqrt(t),leptonmass, leptonmass);
164  G4double Elepton = std::sqrt(Plepton*Plepton + leptonmass*leptonmass );
165  costheta = 2.*G4UniformRand()-1.0;
166  sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
167  phi = twopi*G4UniformRand()*rad;
168  G4ThreeVector ldirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
169  //create G4DynamicParticle for leptons in the rest frame of (l+ l-)system
170  G4DynamicParticle * leptonparticle
172  ldirection, Elepton-leptonmass );
173  G4DynamicParticle * antileptonparticle
175  -1.0*ldirection, Elepton-leptonmass );
176  //boost leptons in the rest frame of the parent
177  G4LorentzVector p4 = leptonparticle->Get4Momentum();
178  p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
179  leptonparticle->Set4Momentum(p4);
180  p4 = antileptonparticle->Get4Momentum();
181  p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
182  antileptonparticle->Set4Momentum(p4);
183 
184  //create G4Decayproducts
185  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
186  delete parentparticle;
187  products->PushProducts(gammaparticle);
188  products->PushProducts(leptonparticle);
189  products->PushProducts(antileptonparticle);
190 
191 #ifdef G4VERBOSE
192  if (GetVerboseLevel()>1) {
193  G4cout << "G4DalitzDecayChannel::DecayIt ";
194  G4cout << " create decay products in rest frame " <<G4endl;
195  products->DumpInfo();
196  }
197 #endif
198  return products;
199 }
200 
201 
202 
203 
204 
void CheckAndFillDaughters()
void SetBR(G4double value)
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
virtual G4DecayProducts * DecayIt(G4double)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
static constexpr double rad
Definition: G4SIunits.hh:149
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4String kinematics_name
void DumpInfo() const
void SetNumberOfDaughters(G4int value)
G4String * parent_name
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4int GetVerboseLevel() const
G4double GetPDGMass() const
void SetParent(const G4ParticleDefinition *particle_type)
void SetDaughter(G4int anIndex, const G4ParticleDefinition *particle_type)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4DalitzDecayChannel & operator=(const G4DalitzDecayChannel &)
G4String ** daughters_name
static G4double Pmx(G4double e, G4double p1, G4double p2)
CLHEP::HepLorentzVector G4LorentzVector