Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/particle_hp/src/G4ParticleHPDiscreteTwoBody.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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