ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPDiscreteTwoBody.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPDiscreteTwoBody.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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 //080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
31 //080709 Bug fix Sampling Legendre expansion by T. Koi
32 //101110 Bug fix in MF=6, LAW=2 case; contribution from E. Mendoza, D. Cano-Ott (CIEMAT)
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 //
36 #include "G4ParticleHPDiscreteTwoBody.hh"
37 #include "G4Gamma.hh"
38 #include "G4Electron.hh"
39 #include "G4Positron.hh"
40 #include "G4Neutron.hh"
41 #include "G4Proton.hh"
42 #include "G4Deuteron.hh"
43 #include "G4Triton.hh"
44 #include "G4He3.hh"
45 #include "G4Alpha.hh"
46 #include "G4ParticleHPVector.hh"
47 #include "G4ParticleHPLegendreStore.hh"
48 
49 G4ReactionProduct * G4ParticleHPDiscreteTwoBody::Sample(G4double anEnergy, G4double massCode, G4double )
50 { // Interpolation still only for the most used parts; rest to be Done @@@@@
51  G4ReactionProduct * result = new G4ReactionProduct;
52  G4int Z = static_cast<G4int>(massCode/1000);
53  G4int A = static_cast<G4int>(massCode-1000*Z);
54 
55  if(massCode==0)
56  {
57  result->SetDefinition(G4Gamma::Gamma());
58  }
59  else if(A==0)
60  {
61  result->SetDefinition(G4Electron::Electron());
62  if(Z==1) result->SetDefinition(G4Positron::Positron());
63  }
64  else if(A==1)
65  {
66  result->SetDefinition(G4Neutron::Neutron());
67  if(Z==1) result->SetDefinition(G4Proton::Proton());
68  }
69  else if(A==2)
70  {
71  result->SetDefinition(G4Deuteron::Deuteron());
72  }
73  else if(A==3)
74  {
75  result->SetDefinition(G4Triton::Triton());
76  if(Z==2) result->SetDefinition(G4He3::He3());
77  }
78  else if(A==4)
79  {
80  result->SetDefinition(G4Alpha::Alpha());
81  if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
82  }
83  else
84  {
85  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPDiscreteTwoBody: Unknown ion case 2");
86  }
87 
88 // get cosine(theta)
89  G4int i(0), it(0);
90  G4double cosTh(0);
91  for(i=0; i<nEnergy; i++)
92  {
93  it = i;
94  if(theCoeff[i].GetEnergy()>anEnergy) break;
95  }
96  if(it==0||it==nEnergy-1)
97  {
98  if(theCoeff[it].GetRepresentation()==0)
99  {
100 //TK Legendre expansion
101  G4ParticleHPLegendreStore theStore(1);
102  theStore.SetCoeff(0, theCoeff);
103  theStore.SetManager(theManager);
104  //cosTh = theStore.SampleMax(anEnergy);
105  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
106  cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
107  }
108  else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
109  {
110  G4ParticleHPVector theStore;
111  G4InterpolationManager aManager;
112  aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
113  theStore.SetInterpolationManager(aManager);
114  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
115  {
116  //101110
117  //theStore.SetX(i, theCoeff[it].GetCoeff(i));
118  //theStore.SetY(i, theCoeff[it].GetCoeff(i));
119  theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
120  theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
121  }
122  cosTh = theStore.Sample();
123  }
124  else if(theCoeff[it].GetRepresentation()==14) //this is LOGLIN
125  {
126  G4ParticleHPVector theStore;
127  G4InterpolationManager aManager;
128  aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
129  theStore.SetInterpolationManager(aManager);
130  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
131  {
132  //101110
133  //theStore.SetX(i, theCoeff[it].GetCoeff(i));
134  //theStore.SetY(i, theCoeff[it].GetCoeff(i));
135  theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
136  theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
137  }
138  cosTh = theStore.Sample();
139  }
140  else
141  {
142  throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering");
143  }
144  }
145  else
146  {
147  if(!bCheckDiffCoeffRepr || theCoeff[it].GetRepresentation() == theCoeff[it-1].GetRepresentation())
148  {
149  if(theCoeff[it].GetRepresentation()==0)
150  {
151 //TK Legendre expansion
152  G4ParticleHPLegendreStore theStore(2);
153  theStore.SetCoeff(0, &(theCoeff[it-1]));
154  theStore.SetCoeff(1, &(theCoeff[it]));
155  G4InterpolationManager aManager;
156  aManager.Init(theManager.GetScheme(it), 2);
157  theStore.SetManager(aManager);
158  //cosTh = theStore.SampleMax(anEnergy);
159 //080709 TKDB
160  cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
161  }
162  else if(theCoeff[it].GetRepresentation()==12) // LINLIN
163  {
164  G4ParticleHPVector theBuff1;
165  G4InterpolationManager aManager1;
166  aManager1.Init(LINLIN, theCoeff[it-1].GetNumberOfPoly()/2);
167  theBuff1.SetInterpolationManager(aManager1);
168  for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i+=2)
169  {
170  //101110
171  //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
172  //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
173  theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
174  theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
175  }
176  G4ParticleHPVector theBuff2;
177  G4InterpolationManager aManager2;
178  aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
179  theBuff2.SetInterpolationManager(aManager2);
180  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
181  {
182  //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
183  //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
184  theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
185  theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
186  }
187 
188  G4double x1 = theCoeff[it-1].GetEnergy();
189  G4double x2 = theCoeff[it].GetEnergy();
190  G4double x = anEnergy;
191  G4double y1, y2, y, mu;
192 
193  G4ParticleHPVector theStore1;
194  theStore1.SetInterpolationManager(aManager1);
195  G4ParticleHPVector theStore2;
196  theStore2.SetInterpolationManager(aManager2);
197  G4ParticleHPVector theStore;
198 
199  // for fixed mu get p1, p2 and interpolate according to x
200  for(i=0; i<theBuff1.GetVectorLength(); i++)
201  {
202  mu = theBuff1.GetX(i);
203  y1 = theBuff1.GetY(i);
204  y2 = theBuff2.GetY(mu);
205  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
206  theStore1.SetData(i, mu, y);
207  }
208  for(i=0; i<theBuff2.GetVectorLength(); i++)
209  {
210  mu = theBuff2.GetX(i);
211  y1 = theBuff2.GetY(i);
212  y2 = theBuff1.GetY(mu);
213  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
214  theStore2.SetData(i, mu, y);
215  }
216  theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
217  cosTh = theStore.Sample();
218  }
219  else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
220  {
221  G4ParticleHPVector theBuff1;
222  G4InterpolationManager aManager1;
223  aManager1.Init(LOGLIN, theCoeff[it-1].GetNumberOfPoly()/2);
224  theBuff1.SetInterpolationManager(aManager1);
225  for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i+=2)
226  {
227  //101110
228  //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
229  //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
230  theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
231  theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
232  }
233 
234  G4ParticleHPVector theBuff2;
235  G4InterpolationManager aManager2;
236  aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
237  theBuff2.SetInterpolationManager(aManager2);
238  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
239  {
240  //101110
241  //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
242  //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
243  theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
244  theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
245  }
246 
247  G4double x1 = theCoeff[it-1].GetEnergy();
248  G4double x2 = theCoeff[it].GetEnergy();
249  G4double x = anEnergy;
250  G4double y1, y2, y, mu;
251 
252  G4ParticleHPVector theStore1;
253  theStore1.SetInterpolationManager(aManager1);
254  G4ParticleHPVector theStore2;
255  theStore2.SetInterpolationManager(aManager2);
256  G4ParticleHPVector theStore;
257 
258  // for fixed mu get p1, p2 and interpolate according to x
259  for(i=0; i<theBuff1.GetVectorLength(); i++)
260  {
261  mu = theBuff1.GetX(i);
262  y1 = theBuff1.GetY(i);
263  y2 = theBuff2.GetY(mu);
264  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
265  theStore1.SetData(i, mu, y);
266  }
267  for(i=0; i<theBuff2.GetVectorLength(); i++)
268  {
269  mu = theBuff2.GetX(i);
270  y1 = theBuff2.GetY(i);
271  y2 = theBuff1.GetY(mu);
272  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
273  theStore2.SetData(i, mu, y);
274  }
275  theStore.Merge(&theStore1, &theStore2);
276  cosTh = theStore.Sample();
277  }
278  else
279  {
280  throw G4HadronicException(__FILE__, __LINE__, "Two neighbouring distributions with different interpolation");
281  }
282  }
283  else
284  {
285  G4cout << " theCoeff[it].GetRepresent MEM " << it << " " << &theCoeff[it] << " " << &theCoeff[it-1] << G4endl;
286  G4cout << " theCoeff[it].GetRepresent " << it << " " << theCoeff[it].GetRepresentation() << " != " << theCoeff[it-1].GetRepresentation() << G4endl;
287 
288  throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering, case 2");
289  }
290  }
291 
292 // now get the energy from kinematics and Q-value.
293 
294  //G4double restEnergy = anEnergy+GetQValue();
295 
296 // assumed to be in CMS @@@@@@@@@@@@@@@@@
297 
298  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
299  //G4double residualMass = GetTarget()->GetMass() + GetNeutron()->GetMass()
300  // - result->GetMass() - GetQValue();
301  //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
302  G4double A1 = GetTarget()->GetMass()/GetProjectileRP()->GetMass();
303  G4double A1prim = result->GetMass()/GetProjectileRP()->GetMass();
304  //G4double E1 = (A1+1)*(A1+1)/A1/A1*anEnergy;
305  //Bug fix Bugzilla #1815
306  G4double E1 = anEnergy;
307  G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
308 
309  result->SetKineticEnergy(kinE); // non relativistic @@
310  G4double phi = CLHEP::twopi*G4UniformRand();
311  G4double theta = std::acos(cosTh);
312  G4double sinth = std::sin(theta);
313  G4double mtot = result->GetTotalMomentum();
314  G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
315  result->SetMomentum(tempVector);
316 
317 // some garbage collection
318 
319 // return the result
320  return result;
321 }