Geant4  10.01.p03
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 67971 2013-03-13 10:13:24Z 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
111  if (G4MT_parent == 0) FillParent();
112  if (G4MT_daughters == 0) FillDaughters();
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  do {
130  x = G4UniformRand()*(xmax-xmin) + xmin;
131  w = G4UniformRand()*wmax;
132  t = std::exp(x);
133  w1 = (1.0-4.0*leptonmass*leptonmass/t);
134  if ( w1 > 0.0) {
135  w2 = ( 1.0 + 2.0*leptonmass*leptonmass/t );
136  w3 = ( 1.0 - t/parentmass/parentmass );
137  w3 = w3 * w3 * w3;
138  ww = w3 * w2 * std::sqrt(w1);
139  } else {
140  ww = 0.0;
141  }
142  } while (w > ww);
143 
144  // calculate gamma momentum
145  G4double Pgamma =
146  G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, std::sqrt(t));
147  G4double costheta = 2.*G4UniformRand()-1.0;
148  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
149  G4double phi = twopi*G4UniformRand()*rad;
150  G4ThreeVector gdirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
151 
152  //create G4DynamicParticle for gamma
153  G4DynamicParticle * gammaparticle
154  = new G4DynamicParticle(G4MT_daughters[idGamma] , gdirection, Pgamma);
155 
156  // calcurate beta of (l+ l-)system
157  G4double beta = Pgamma/(parentmass-Pgamma);
158 
159  // calculate lepton momentum in the rest frame of (l+ l-)system
160  G4double Plepton =
161  G4PhaseSpaceDecayChannel::Pmx(std::sqrt(t),leptonmass, leptonmass);
162  G4double Elepton = std::sqrt(Plepton*Plepton + leptonmass*leptonmass );
163  costheta = 2.*G4UniformRand()-1.0;
164  sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
165  phi = twopi*G4UniformRand()*rad;
166  G4ThreeVector ldirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
167  //create G4DynamicParticle for leptons in the rest frame of (l+ l-)system
168  G4DynamicParticle * leptonparticle
170  ldirection, Elepton-leptonmass );
171  G4DynamicParticle * antileptonparticle
173  -1.0*ldirection, Elepton-leptonmass );
174  //boost leptons in the rest frame of the parent
175  G4LorentzVector p4 = leptonparticle->Get4Momentum();
176  p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
177  leptonparticle->Set4Momentum(p4);
178  p4 = antileptonparticle->Get4Momentum();
179  p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
180  antileptonparticle->Set4Momentum(p4);
181 
182  //create G4Decayproducts
183  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
184  delete parentparticle;
185  products->PushProducts(gammaparticle);
186  products->PushProducts(leptonparticle);
187  products->PushProducts(antileptonparticle);
188 
189 #ifdef G4VERBOSE
190  if (GetVerboseLevel()>1) {
191  G4cout << "G4DalitzDecayChannel::DecayIt ";
192  G4cout << " create decay products in rest frame " <<G4endl;
193  products->DumpInfo();
194  }
195 #endif
196  return products;
197 }
198 
199 
200 
201 
202 
void SetBR(G4double value)
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
virtual G4DecayProducts * DecayIt(G4double)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
G4String kinematics_name
void DumpInfo() const
void SetNumberOfDaughters(G4int value)
G4String * parent_name
G4LorentzVector Get4Momentum() const
static const double rad
Definition: G4SIunits.hh:130
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