Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4DalitzDecayChannel class implementation << 27 // 23 // 28 // Author: H.Kurashige, 30 May 1997 << 24 // $Id: G4DalitzDecayChannel.cc,v 1.5 2001/07/11 10:01:59 gunter Exp $ 29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-05-02-patch-01 $ 30 << 26 // 31 #include "G4DalitzDecayChannel.hh" << 27 // >> 28 // ------------------------------------------------------------ >> 29 // GEANT 4 class header file >> 30 // >> 31 // History: first implementation, based on object model of >> 32 // 30 May 1997 H.Kurashige >> 33 // ------------------------------------------------------------ 32 34 33 #include "G4DecayProducts.hh" << 34 #include "G4LorentzRotation.hh" << 35 #include "G4LorentzVector.hh" << 36 #include "G4ParticleDefinition.hh" 35 #include "G4ParticleDefinition.hh" 37 #include "G4PhaseSpaceDecayChannel.hh" << 36 #include "G4DecayProducts.hh" 38 #include "G4PhysicalConstants.hh" << 39 #include "G4SystemOfUnits.hh" << 40 #include "G4VDecayChannel.hh" 37 #include "G4VDecayChannel.hh" >> 38 #include "G4DalitzDecayChannel.hh" >> 39 #include "G4PhaseSpaceDecayChannel.hh" 41 #include "Randomize.hh" 40 #include "Randomize.hh" >> 41 #include "G4LorentzVector.hh" >> 42 #include "G4LorentzRotation.hh" 42 43 43 G4DalitzDecayChannel::G4DalitzDecayChannel(con << 44 G4DalitzDecayChannel::G4DalitzDecayChannel( 44 con << 45 const G4String& theParentName, 45 con << 46 G4double theBR, 46 : G4VDecayChannel("Dalitz Decay", 1) << 47 const G4String& theLeptonName, >> 48 const G4String& theAntiLeptonName) >> 49 :G4VDecayChannel("Dalitz Decay",1) 47 { 50 { >> 51 //#ifdef G4VERBOSE >> 52 //if (GetVerboseLevel()>1) { >> 53 // G4cout << "G4DalitzDecayChannel:: constructor "; >> 54 // G4cout << "addr[" << this << "]" << G4endl; >> 55 //} >> 56 //#endif 48 // set names for daughter particles 57 // set names for daughter particles 49 SetParent(theParentName); 58 SetParent(theParentName); 50 SetBR(theBR); 59 SetBR(theBR); 51 SetNumberOfDaughters(3); 60 SetNumberOfDaughters(3); 52 G4String gammaName = "gamma"; 61 G4String gammaName = "gamma"; 53 SetDaughter(idGamma, gammaName); 62 SetDaughter(idGamma, gammaName); 54 SetDaughter(idLepton, theLeptonName); 63 SetDaughter(idLepton, theLeptonName); 55 SetDaughter(idAntiLepton, theAntiLeptonName) 64 SetDaughter(idAntiLepton, theAntiLeptonName); >> 65 // >> 66 //#ifdef G4VERBOSE >> 67 //if (GetVerboseLevel()>1) DumpInfo(); >> 68 //#endif 56 } 69 } 57 70 58 G4DalitzDecayChannel& G4DalitzDecayChannel::op << 71 G4DalitzDecayChannel::~G4DalitzDecayChannel() 59 { 72 { 60 if (this != &right) { << 61 kinematics_name = right.kinematics_name; << 62 verboseLevel = right.verboseLevel; << 63 rbranch = right.rbranch; << 64 << 65 // copy parent name << 66 parent_name = new G4String(*right.parent_n << 67 << 68 // clear daughters_name array << 69 ClearDaughtersName(); << 70 << 71 // recreate array << 72 numberOfDaughters = right.numberOfDaughter << 73 if (numberOfDaughters > 0) { << 74 if (daughters_name != nullptr) ClearDaug << 75 daughters_name = new G4String*[numberOfD << 76 // copy daughters name << 77 for (G4int index = 0; index < numberOfDa << 78 daughters_name[index] = new G4String(* << 79 } << 80 } << 81 } << 82 return *this; << 83 } 73 } 84 74 85 G4DecayProducts* G4DalitzDecayChannel::DecayIt << 75 G4DecayProducts *G4DalitzDecayChannel::DecayIt(G4double) 86 { 76 { 87 #ifdef G4VERBOSE 77 #ifdef G4VERBOSE 88 if (GetVerboseLevel() > 1) G4cout << "G4Dali << 78 if (GetVerboseLevel()>1) G4cout << "G4DalitzDecayChannel::DecayIt "; 89 #endif << 79 #endif 90 CheckAndFillParent(); << 80 if (parent == 0) FillParent(); 91 CheckAndFillDaughters(); << 81 if (daughters == 0) FillDaughters(); 92 82 93 // parent mass 83 // parent mass 94 G4double parentmass = G4MT_parent->GetPDGMas << 84 G4double parentmass = parent->GetPDGMass(); 95 << 85 96 // create parent G4DynamicParticle at rest << 86 //create parent G4DynamicParticle at rest 97 G4ThreeVector dummy; 87 G4ThreeVector dummy; 98 auto parentparticle = new G4DynamicParticle( << 88 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 99 << 89 100 // daughters' mass << 90 //daughters'mass 101 G4double leptonmass = G4MT_daughters[idLepto << 91 G4double leptonmass = daughters[idLepton]->GetPDGMass(); 102 << 92 103 // Generate t ( = std::exp(x):mass Square of << 93 // Generate t ( = exp(x):mass Square of (l+ l-) system) 104 G4double xmin = 2.0 * std::log(2.0 * leptonm << 94 G4double xmin = 2.0*log(2.0*leptonmass); 105 G4double xmax = 2.0 * std::log(parentmass); << 95 G4double xmax = 2.0*log(parentmass); 106 G4double wmax = 1.5; 96 G4double wmax = 1.5; 107 G4double x, w, ww, w1, w2, w3, t; 97 G4double x, w, ww, w1, w2, w3, t; 108 const std::size_t MAX_LOOP = 10000; << 98 do { 109 for (std::size_t loop_counter = 0; loop_coun << 99 x = G4UniformRand()*(xmax-xmin) + xmin; 110 x = G4UniformRand() * (xmax - xmin) + xmin << 100 w = G4UniformRand()*wmax; 111 w = G4UniformRand() * wmax; << 101 t = exp(x); 112 t = std::exp(x); << 102 w1 = (1.0-4.0*leptonmass*leptonmass/t); 113 w1 = (1.0 - 4.0 * leptonmass * leptonmass << 103 if ( w1 > 0.0) { 114 if (w1 > 0.0) { << 104 w2 = ( 1.0 + 2.0*leptonmass*leptonmass/t ); 115 w2 = (1.0 + 2.0 * leptonmass * leptonmas << 105 w3 = ( 1.0 - t/parentmass/parentmass ); 116 w3 = (1.0 - t / parentmass / parentmass) << 117 w3 = w3 * w3 * w3; 106 w3 = w3 * w3 * w3; 118 ww = w3 * w2 * std::sqrt(w1); << 107 ww = w3 * w2 * sqrt(w1); 119 } << 108 } else { 120 else { << 121 ww = 0.0; 109 ww = 0.0; 122 } 110 } 123 if (w <= ww) break; << 111 } while (w > ww); 124 } << 112 125 << 126 // calculate gamma momentum 113 // calculate gamma momentum 127 G4double Pgamma = G4PhaseSpaceDecayChannel:: << 114 G4double Pgamma = 128 G4double costheta = 2. * G4UniformRand() - 1 << 115 G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, sqrt(t)); 129 G4double sintheta = std::sqrt((1.0 - costhet << 116 G4double costheta = 2.*G4UniformRand()-1.0; 130 G4double phi = twopi * G4UniformRand() * rad << 117 G4double sintheta = sqrt((1.0 - costheta)*(1.0 + costheta)); 131 G4ThreeVector gdirection(sintheta * std::cos << 118 G4double phi = 2.0*M_PI*G4UniformRand()*rad; 132 << 119 G4ThreeVector gdirection(sintheta*cos(phi),sintheta*sin(phi),costheta); 133 // create G4DynamicParticle for gamma << 120 134 auto gammaparticle = new G4DynamicParticle(G << 121 //create G4DynamicParticle for gamma >> 122 G4DynamicParticle * gammaparticle >> 123 = new G4DynamicParticle(daughters[idGamma] , gdirection, Pgamma); 135 124 136 // calculate beta of (l+ l-)system << 125 // calcurate beta of (l+ l-)system 137 G4double beta = Pgamma / (parentmass - Pgamm << 126 G4double beta = Pgamma/(parentmass-Pgamma); 138 127 139 // calculate lepton momentum in the rest fra 128 // calculate lepton momentum in the rest frame of (l+ l-)system 140 G4double Plepton = G4PhaseSpaceDecayChannel: << 129 G4double Plepton = 141 G4double Elepton = std::sqrt(Plepton * Plept << 130 G4PhaseSpaceDecayChannel::Pmx(sqrt(t),leptonmass, leptonmass); 142 costheta = 2. * G4UniformRand() - 1.0; << 131 G4double Elepton = sqrt(Plepton*Plepton + leptonmass*leptonmass ); 143 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 132 costheta = 2.*G4UniformRand()-1.0; 144 phi = twopi * G4UniformRand() * rad; << 133 sintheta = sqrt((1.0 - costheta)*(1.0 + costheta)); 145 G4ThreeVector ldirection(sintheta * std::cos << 134 phi = 2.0*M_PI*G4UniformRand()*rad; 146 // create G4DynamicParticle for leptons in << 135 G4ThreeVector ldirection(sintheta*cos(phi),sintheta*sin(phi),costheta); 147 auto leptonparticle = << 136 //create G4DynamicParticle for leptons in the rest frame of (l+ l-)system 148 new G4DynamicParticle(G4MT_daughters[idLep << 137 G4DynamicParticle * leptonparticle 149 auto antileptonparticle = << 138 = new G4DynamicParticle(daughters[idLepton] , 150 new G4DynamicParticle(G4MT_daughters[idAnt << 139 ldirection, Elepton-leptonmass ); 151 // boost leptons in the rest frame of the pa << 140 G4DynamicParticle * antileptonparticle >> 141 = new G4DynamicParticle(daughters[idAntiLepton] , >> 142 -1.0*ldirection, Elepton-leptonmass ); >> 143 //boost leptons in the rest frame of the parent 152 G4LorentzVector p4 = leptonparticle->Get4Mom 144 G4LorentzVector p4 = leptonparticle->Get4Momentum(); 153 p4.boost(-1.0 * gdirection.x() * beta, -1.0 << 145 p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta); 154 -1.0 * gdirection.z() * beta); << 155 leptonparticle->Set4Momentum(p4); 146 leptonparticle->Set4Momentum(p4); 156 p4 = antileptonparticle->Get4Momentum(); 147 p4 = antileptonparticle->Get4Momentum(); 157 p4.boost(-1.0 * gdirection.x() * beta, -1.0 << 148 p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta); 158 -1.0 * gdirection.z() * beta); << 159 antileptonparticle->Set4Momentum(p4); 149 antileptonparticle->Set4Momentum(p4); 160 150 161 // create G4Decayproducts << 151 //create G4Decayproducts 162 auto products = new G4DecayProducts(*parentp << 152 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 163 delete parentparticle; 153 delete parentparticle; 164 products->PushProducts(gammaparticle); 154 products->PushProducts(gammaparticle); 165 products->PushProducts(leptonparticle); 155 products->PushProducts(leptonparticle); 166 products->PushProducts(antileptonparticle); 156 products->PushProducts(antileptonparticle); 167 157 168 #ifdef G4VERBOSE 158 #ifdef G4VERBOSE 169 if (GetVerboseLevel() > 1) { << 159 if (GetVerboseLevel()>1) { 170 G4cout << "G4DalitzDecayChannel::DecayIt " << 160 G4cout << "G4DalitzDecayChannel::DecayIt "; 171 G4cout << " create decay products in rest << 161 G4cout << " create decay products in rest frame " <<G4endl; 172 products->DumpInfo(); << 162 products->DumpInfo(); 173 } 163 } 174 #endif 164 #endif 175 return products; 165 return products; 176 } 166 } >> 167 >> 168 >> 169 >> 170 >> 171 177 172