ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DalitzDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DalitzDecayChannel.cc
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 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 30 May 1997 H.Kurashige
34 // ------------------------------------------------------------
35 
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4DecayProducts.hh"
40 #include "G4VDecayChannel.hh"
41 #include "G4DalitzDecayChannel.hh"
43 #include "Randomize.hh"
44 #include "G4LorentzVector.hh"
45 #include "G4LorentzRotation.hh"
46 
49 {
50 }
51 
53  const G4String& theParentName,
54  G4double theBR,
55  const G4String& theLeptonName,
56  const G4String& theAntiLeptonName)
57  :G4VDecayChannel("Dalitz Decay",1)
58 {
59  // set names for daughter particles
60  SetParent(theParentName);
61  SetBR(theBR);
63  G4String gammaName = "gamma";
64  SetDaughter(idGamma, gammaName);
65  SetDaughter(idLepton, theLeptonName);
66  SetDaughter(idAntiLepton, theAntiLeptonName);
67 }
68 
70 {
71 }
72 
74  G4VDecayChannel(right)
75 {
76 }
77 
79 {
80  if (this != &right) {
82  verboseLevel = right.verboseLevel;
83  rbranch = right.rbranch;
84 
85  // copy parent name
86  parent_name = new G4String(*right.parent_name);
87 
88  // clear daughters_name array
90 
91  // recreate array
93  if ( numberOfDaughters >0 ) {
96  //copy daughters name
97  for (G4int index=0; index < numberOfDaughters; index++) {
98  daughters_name[index] = new G4String(*right.daughters_name[index]);
99  }
100  }
101  }
102  return *this;
103 }
104 
106 {
107 #ifdef G4VERBOSE
108  if (GetVerboseLevel()>1) G4cout << "G4DalitzDecayChannel::DecayIt ";
109 #endif
112 
113  // parent mass
114  G4double parentmass = G4MT_parent->GetPDGMass();
115 
116  //create parent G4DynamicParticle at rest
117  G4ThreeVector dummy;
118  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
119 
120  //daughters'mass
121  G4double leptonmass = G4MT_daughters[idLepton]->GetPDGMass();
122 
123  // Generate t ( = std::exp(x):mass Square of (l+ l-) system)
124  G4double xmin = 2.0*std::log(2.0*leptonmass);
125  G4double xmax = 2.0*std::log(parentmass);
126  G4double wmax = 1.5;
127  G4double x, w, ww, w1, w2, w3, t;
128  const size_t MAX_LOOP = 10000;
129  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
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  if (w <= ww) break;
143  }
144 
145  // calculate gamma momentum
146  G4double Pgamma =
147  G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, std::sqrt(t));
148  G4double costheta = 2.*G4UniformRand()-1.0;
149  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
151  G4ThreeVector gdirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
152 
153  //create G4DynamicParticle for gamma
154  G4DynamicParticle * gammaparticle
155  = new G4DynamicParticle(G4MT_daughters[idGamma] , gdirection, Pgamma);
156 
157  // calcurate beta of (l+ l-)system
158  G4double beta = Pgamma/(parentmass-Pgamma);
159 
160  // calculate lepton momentum in the rest frame of (l+ l-)system
161  G4double Plepton =
162  G4PhaseSpaceDecayChannel::Pmx(std::sqrt(t),leptonmass, leptonmass);
163  G4double Elepton = std::sqrt(Plepton*Plepton + leptonmass*leptonmass );
164  costheta = 2.*G4UniformRand()-1.0;
165  sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
166  phi = twopi*G4UniformRand()*rad;
167  G4ThreeVector ldirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
168  //create G4DynamicParticle for leptons in the rest frame of (l+ l-)system
169  G4DynamicParticle * leptonparticle
171  ldirection, Elepton-leptonmass );
172  G4DynamicParticle * antileptonparticle
174  -1.0*ldirection, Elepton-leptonmass );
175  //boost leptons in the rest frame of the parent
176  G4LorentzVector p4 = leptonparticle->Get4Momentum();
177  p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
178  leptonparticle->Set4Momentum(p4);
179  p4 = antileptonparticle->Get4Momentum();
180  p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
181  antileptonparticle->Set4Momentum(p4);
182 
183  //create G4Decayproducts
184  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
185  delete parentparticle;
186  products->PushProducts(gammaparticle);
187  products->PushProducts(leptonparticle);
188  products->PushProducts(antileptonparticle);
189 
190 #ifdef G4VERBOSE
191  if (GetVerboseLevel()>1) {
192  G4cout << "G4DalitzDecayChannel::DecayIt ";
193  G4cout << " create decay products in rest frame " <<G4endl;
194  products->DumpInfo();
195  }
196 #endif
197  return products;
198 }
199 
200 
201 
202 
203