ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutrinoElectronCcModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NeutrinoElectronCcModel.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 // Geant4 Header : G4NeutrinoElectronCcModel
28 //
29 // Author : V.Grichine 26.4.17
30 //
31 
33 #include "G4SystemOfUnits.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4IonTable.hh"
37 #include "Randomize.hh"
38 #include "G4NeutrinoE.hh"
39 #include "G4AntiNeutrinoE.hh"
40 
41 #include "G4NeutrinoMu.hh"
42 #include "G4AntiNeutrinoMu.hh"
43 #include "G4NeutrinoTau.hh"
44 #include "G4AntiNeutrinoTau.hh"
45 #include "G4MuonMinus.hh"
46 #include "G4TauMinus.hh"
47 #include "G4HadronicParameters.hh"
48 
49 using namespace std;
50 using namespace CLHEP;
51 
53  : G4HadronicInteraction(name)
54 {
55  SetMinEnergy( 0.0*GeV );
57  SetMinEnergy(1.e-6*eV);
58 
61 
64 
67 
70 
71  // PDG2016: sin^2 theta Weinberg
72 
73  fSin2tW = 0.23129; // 0.2312;
74 
75  fCutEnergy = 0.; // default value
76 
77 }
78 
79 
81 {}
82 
83 
84 void G4NeutrinoElectronCcModel::ModelDescription(std::ostream& outFile) const
85 {
86 
87  outFile << "G4NeutrinoElectronCcModel is a neutrino-electron (neutral current) elastic scattering\n"
88  << "model which uses the standard model \n"
89  << "transfer parameterization. The model is fully relativistic\n";
90 
91 }
92 
94 
96  G4Nucleus & targetNucleus)
97 {
98  G4bool result = false;
99  G4String pName = aPart.GetDefinition()->GetParticleName();
100  if(pName == "anti_nu_mu" || pName == "anti_nu_tau") return result; // no cc for anti_nu_(mu,tau)
101  G4double minEnergy = 0., energy = aPart.GetTotalEnergy();
102  G4double fmass, emass = electron_mass_c2;
103 
104  if( pName == "nu_mu" ) fmass = theMuonMinus->GetPDGMass();
105  else if( pName == "nu_tau" ) fmass = theTauMinus->GetPDGMass();
106  else fmass = emass;
107 
108  minEnergy = (fmass-emass)*(fmass+emass)/emass;
109  SetMinEnergy( minEnergy );
110 
111  if( ( pName == "nu_mu" || pName == "nu_tau" || pName == "anti_nu_e" ) && energy > minEnergy )
112  {
113  result = true;
114  }
115  G4int Z = targetNucleus.GetZ_asInt();
116  Z *= 1;
117 
118  return result;
119 }
120 
122 //
123 //
124 
126  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
127 {
129 
130  const G4HadProjectile* aParticle = &aTrack;
131  G4double energy = aParticle->GetTotalEnergy();
132 
133  G4String pName = aParticle->GetDefinition()->GetParticleName();
134  G4double minEnergy(0.), fmass(0.), emass = electron_mass_c2;
135 
136  if( pName == "nu_mu" ) fmass = theMuonMinus->GetPDGMass();
137  else if( pName == "nu_tau" ) fmass = theTauMinus->GetPDGMass();
138  else fmass = emass;
139 
140  minEnergy = (fmass-emass)*(fmass+emass)/emass;
141 
142  if( energy <= minEnergy )
143  {
146  return &theParticleChange;
147  }
148  G4double massf(0.), massf2(0.); // , emass = electron_mass_c2;
149  G4double sTot = 2.*energy*emass + emass*emass;
150 
151  G4LorentzVector lvp1 = aParticle->Get4Momentum();
152  G4LorentzVector lvt1(0.,0.,0.,electron_mass_c2);
153  G4LorentzVector lvsum = lvp1+lvt1;
154  G4ThreeVector bst = lvsum.boostVector();
155 
156  // sample and make final state in CMS frame
157 
158  G4double cost = SampleCosCMS( aParticle );
159  G4double sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
161 
162  G4ThreeVector eP( sint*std::cos(phi), sint*std::sin(phi), cost );
163 
164  if( pName == "nu_mu" ) massf = theMuonMinus->GetPDGMass();
165  else if( pName == "nu_tau" ) massf = theTauMinus->GetPDGMass();
166 
167  massf2 = massf*massf;
168 
169  G4double epf = 0.5*(sTot - massf2)/sqrt(sTot);
170  // G4double etf = epf*(sTot + massf2)/(sTot - massf2);
171 
172  eP *= epf;
173  G4LorentzVector lvp2( eP, epf );
174  lvp2.boost(bst); // back to lab frame
175 
176  G4LorentzVector lvt2 = lvsum - lvp2; // ?
177 
178  G4DynamicParticle* aNu = nullptr;
179  G4DynamicParticle* aLept = nullptr;
180 
181  if( pName == "nu_mu" || pName == "nu_tau")
182  {
183  aNu = new G4DynamicParticle( theNeutrinoE, lvp2 );
184  }
185  else if( pName == "anti_nu_e" ) aNu = new G4DynamicParticle( theAntiNeutrinoMu, lvp2 ); // s-channel for mu (tau later)
186 
187  if( pName == "nu_mu" || pName == "anti_nu_e")
188  {
189  aLept = new G4DynamicParticle( theMuonMinus, lvt2 );
190  }
191  else if( pName == "nu_tau" ) // || pName == "anti_nu_tau")
192  {
193  aLept = new G4DynamicParticle( theTauMinus, lvt2 );
194  }
195  if(aNu) { theParticleChange.AddSecondary( aNu ); }
196  if(aLept) { theParticleChange.AddSecondary( aLept ); }
197 
198  G4int Z = targetNucleus.GetZ_asInt();
199  Z *= 1;
200 
201  return &theParticleChange;
202 }
203 
205 //
206 // sample recoil electron energy in lab frame
207 
209 {
210  G4double result = 0., cofL, cofR, cofLR, massf2, sTot, emass = electron_mass_c2, emass2;
211 
212  G4double energy = aParticle->GetTotalEnergy();
213 
214  if( energy == 0.) return result; // vmg: < th?? as in xsc
215 
216  G4String pName = aParticle->GetDefinition()->GetParticleName();
217 
218  if( pName == "nu_mu" || pName == "nu_tau")
219  {
220  return 2.*G4UniformRand()-1.; // uniform scattering cos in CMS
221  }
222  else if( pName == "anti_nu_mu" || pName == "anti_nu_tau")
223  {
224  emass2 = emass*emass;
225  sTot = 2.*energy*emass + emass2;
226 
227  cofL = (sTot-emass2)/(sTot+emass2);
228 
229  if(pName == "anti_nu_mu") massf2 = theMuonMinus->GetPDGMass()*theMuonMinus->GetPDGMass();
230  else massf2 = theTauMinus->GetPDGMass()*theTauMinus->GetPDGMass();
231 
232  cofR = (sTot-massf2)/(sTot+massf2);
233 
234  cofLR = cofL*cofR/3.;
235 
236  // cofs of cos 3rd equation
237 
238  G4double a = cofLR;
239  G4double b = 0.5*(cofR+cofL);
240  G4double c = 1.;
241 
242  G4double d = -G4UniformRand()*2.*(1.+ cofLR);
243  d += c - b + a;
244 
245  // G4cout<<a<<" "<<b<<" "<<c<<" "<<d<<G4endl<<G4endl;
246 
247  // cofs of the incomplete 3rd equation
248 
249  G4double p = c/a;
250  p -= b*b/a/a/3.;
251 
252  G4double q = d/a;
253  q -= b*c/a/a/3.;
254  q += 2*b*b*b/a/a/a/27.;
255 
256 
257  // cofs for the incomplete colutions
258 
259  G4double D = p*p*p/3./3./3.;
260  D += q*q/2./2.;
261 
262  // G4cout<<"D = "<<D<<G4endl;
263  if(D < 0.) D = -D;
264  // G4complex A1 = G4complex(- q/2., std::sqrt(-D) );
265  // G4complex A = std::pow(A1,1./3.);
266 
267  // G4complex B1 = G4complex(- q/2., -std::sqrt(-D) );
268  // G4complex B = std::pow(B1,1./3.);
269 
270  G4double A, B;
271 
272  G4double A1 = - q/2. + std::sqrt(D);
273  if (A1 < 0.) A1 = -A1;
274  A = std::pow(A1,1./3.);
275  if (A1 < 0.) A = -A;
276 
277  G4double B1 = - q/2. - std::sqrt(D);
278  // G4double B = std::pow(-B1,1./3.);
279  if(B1 < 0.) B1 = -B1;
280  B = std::pow(B1,1./3.);
281  if(B1 < 0.) B = -B;
282  // G4cout<<"A1 = "<<A1<<"; A = "<<A<<"; B1 = "<<B1<<"; B = "<<B<<G4endl;
283  // roots of the incomplete 3rd equation
284 
285  G4complex y1 = A + B;
286  // G4complex y2 = -0.5*(A + B) + 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
287  // G4complex y3 = -0.5*(A + B) - 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
288 
289  G4complex x1 = y1 - b/a/3.;
290  // G4complex x2 = y2 - b/a/3.;
291  // G4complex x3 = y3 - b/a/3.;
292  // G4cout<<"re_x1 = "<<real(x1)<<" + i*"<<imag(x1)<<G4endl;
293  // G4cout<<"re_x1 = "<<real(x1)<<"; re_x2 = "<<real(x2)<<"; re_x3 = "<<real(x3)<<G4endl;
294  // G4cout<<"im_x1 = "<<imag(x1)<<"; im_x2 = "<<imag(x2)<<"; im_x3 = "<<imag(x3)<<G4endl<<G4endl;
295 
296  result = real(x1);
297  }
298  else
299  {
300  return result;
301  }
302  return result;
303 }
304 
305 //
306 //