Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // Geant4 Header : G4NeutrinoElectronNcModel 28 // 29 // Author : V.Grichine 6.4.17 30 // 31 32 #include "G4NeutrinoElectronNcModel.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4ParticleTable.hh" 35 #include "G4ParticleDefinition.hh" 36 #include "G4IonTable.hh" 37 #include "Randomize.hh" 38 #include "G4Electron.hh" 39 #include "G4HadronicParameters.hh" 40 #include "G4PhysicsModelCatalog.hh" 41 42 using namespace std; 43 using namespace CLHEP; 44 45 G4NeutrinoElectronNcModel::G4NeutrinoElectronN 46 : G4HadronElastic(name) 47 { 48 secID = G4PhysicsModelCatalog::GetModelID( " 49 50 SetMinEnergy( 0.0*GeV ); 51 SetMaxEnergy( G4HadronicParameters::Instance 52 SetLowestEnergyLimit(1.e-6*eV); 53 54 theElectron = G4Electron::Electron(); 55 // PDG2016: sin^2 theta Weinberg 56 57 fSin2tW = 0.23129; // 0.2312; 58 59 fCutEnergy = 0.; // default value 60 } 61 62 63 G4NeutrinoElectronNcModel::~G4NeutrinoElectron 64 {} 65 66 void G4NeutrinoElectronNcModel::ModelDescripti 67 { 68 outFile << "G4NeutrinoElectronNcModel is a n 69 << "model which uses the standard model \n 70 << "transfer parameterization. The model 71 } 72 73 ////////////////////////////////////////////// 74 75 G4bool G4NeutrinoElectronNcModel::IsApplicable 76 { 77 G4bool result = false; 78 G4String pName = aTrack.GetDefinition()->Get 79 G4double minEnergy = 0.; 80 G4double energy = aTrack.GetTotalEnergy(); 81 82 if( fCutEnergy > 0. ) // min detected recoil 83 { 84 minEnergy = 0.5*(fCutEnergy+sqrt(fCutEnerg 85 } 86 if( ( pName == "nu_e" || pName == "anti_nu 87 pName == "nu_mu" || pName == "anti_nu 88 pName == "nu_tau" || pName == "anti_nu 89 energy > minEnergy 90 { 91 result = true; 92 } 93 return result; 94 } 95 96 ////////////////////////////////////////////// 97 // 98 // 99 100 G4HadFinalState* G4NeutrinoElectronNcModel::Ap 101 const G4HadProjectile& aTrack, G4Nucleus& 102 { 103 theParticleChange.Clear(); 104 105 const G4HadProjectile* aParticle = &aTrack; 106 G4double nuTkin = aParticle->GetKineticEnerg 107 108 if( nuTkin <= LowestEnergyLimit() ) 109 { 110 theParticleChange.SetEnergyChange(nuTkin); 111 theParticleChange.SetMomentumChange(aTrack 112 return &theParticleChange; 113 } 114 // sample and make final state in lab frame 115 116 G4double eTkin = SampleElectronTkin( aPartic 117 118 if( eTkin > fCutEnergy ) 119 { 120 G4double ePlab = sqrt( eTkin*(eTkin + 2.*e 121 122 G4double cost2 = eTkin*(nuTkin + electron 123 cost2 /= nuTkin*nuTkin*(eTkin + 2 124 125 if( cost2 > 1. ) cost2 = 1.; 126 if( cost2 < 0. ) cost2 = 0.; 127 128 G4double cost = sqrt(cost2); 129 G4double sint = std::sqrt( (1.0 - cost)*(1 130 G4double phi = G4UniformRand()*CLHEP::two 131 132 G4ThreeVector eP( sint*std::cos(phi), sint 133 eP *= ePlab; 134 G4LorentzVector lvt2( eP, eTkin + electron 135 G4DynamicParticle * aSec = new G4DynamicPa 136 theParticleChange.AddSecondary( aSec, secI 137 138 G4LorentzVector lvp1 = aParticle->Get4Mome 139 G4LorentzVector lvt1(0.,0.,0.,electron_mas 140 G4LorentzVector lvsum = lvp1+lvt1; 141 142 G4LorentzVector lvp2 = lvsum-lvt2; 143 G4double nuTkin2 = lvp2.e()-aParticle->Get 144 theParticleChange.SetEnergyChange(nuTkin2) 145 theParticleChange.SetMomentumChange(lvp2.v 146 } 147 else if( eTkin > 0.0 ) 148 { 149 theParticleChange.SetLocalEnergyDeposit( e 150 nuTkin -= eTkin; 151 152 if( nuTkin > 0. ) 153 { 154 theParticleChange.SetEnergyChange( nuTki 155 theParticleChange.SetMomentumChange( aTr 156 } 157 } 158 else 159 { 160 theParticleChange.SetEnergyChange( nuTkin 161 theParticleChange.SetMomentumChange( aTrac 162 } 163 return &theParticleChange; 164 } 165 166 ////////////////////////////////////////////// 167 // 168 // sample recoil electron energy in lab frame 169 170 G4double G4NeutrinoElectronNcModel::SampleElec 171 { 172 G4double result = 0., xi, cofL, cofR, cofL2, 173 174 G4double energy = aParticle->GetTotalEnergy( 175 if( energy == 0.) return result; // vmg: < t 176 177 G4String pName = aParticle->GetDefinition() 178 179 if( pName == "nu_e") 180 { 181 cofL = 0.5 + fSin2tW; 182 cofR = fSin2tW; 183 } 184 else if( pName == "anti_nu_e") 185 { 186 cofL = fSin2tW; 187 cofR = 0.5 + fSin2tW; 188 } 189 else if( pName == "nu_mu") 190 { 191 cofL = -0.5 + fSin2tW; 192 cofR = fSin2tW; 193 } 194 else if( pName == "anti_nu_mu") 195 { 196 cofL = fSin2tW; 197 cofR = -0.5 + fSin2tW; 198 } 199 else if( pName == "nu_tau") // vmg: nu_tau a 200 { 201 cofL = -0.5 + fSin2tW; 202 cofR = fSin2tW; 203 } 204 else if( pName == "anti_nu_tau") 205 { 206 cofL = fSin2tW; 207 cofR = -0.5 + fSin2tW; 208 } 209 else 210 { 211 return result; 212 } 213 xi = 0.5*electron_mass_c2/energy; 214 215 cofL2 = cofL*cofL; 216 cofR2 = cofR*cofR; 217 cofLR = cofL*cofR; 218 219 // cofs of Tkin/Enu 3rd equation 220 221 G4double a = cofR2/3.; 222 G4double b = -(cofR2+cofLR*xi); 223 G4double c = cofL2+cofR2; 224 225 G4double xMax = 1./(1. + xi); 226 G4double xMax2 = xMax*xMax; 227 G4double xMax3 = xMax*xMax2; 228 229 G4double d = -( a*xMax3 + b*xMax2 + c*xMax 230 d *= G4UniformRand(); 231 232 // G4cout<<a<<" "<<b<<" "<<c<<" "<<d<< 233 234 // cofs of the incomplete 3rd equation 235 236 G4double p = c/a; 237 p -= b*b/a/a/3.; 238 G4double q = d/a; 239 q -= b*c/a/a/3.; 240 q += 2*b*b*b/a/a/a/27.; 241 242 243 // cofs for the incomplete colutions 244 245 G4double D = p*p*p/3./3./3.; 246 D += q*q/2./2.; 247 248 // G4cout<<"D = "<<D<<G4endl; 249 // D = -D; 250 // G4complex A1 = G4complex(- q/2., std:: 251 // G4complex A = std::pow(A1,1./3.); 252 253 // G4complex B1 = G4complex(- q/2., -std: 254 // G4complex B = std::pow(B1,1./3.); 255 256 G4double A1 = - q/2. + std::sqrt(D); 257 G4double A = std::pow(A1,1./3.); 258 259 G4double B1 = - q/2. - std::sqrt(D); 260 G4double B = std::pow(-B1,1./3.); 261 B = -B; 262 263 // roots of the incomplete 3rd equation 264 265 G4complex y1 = A + B; 266 // G4complex y2 = -0.5*(A + B) + 0.5*std::sq 267 // G4complex y3 = -0.5*(A + B) - 0.5*std::sq 268 269 G4complex x1 = y1 - b/a/3.; 270 // G4complex x2 = y2 - b/a/3.; 271 // G4complex x3 = y3 - b/a/3.; 272 273 // G4cout<<"re_x1 = "<<real(x1)<<"; re_x2 = 274 // G4cout<<"im_x1 = "<<imag(x1)<<"; im_x2 = 275 276 result = real(x1)*energy; 277 278 return result; 279 } 280 281 // 282 // 283 /////////////////////////// 284