Geant4  10.01.p02
G4LMsdGenerator.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 // G4LMsdGenerator
29 //
30 //
31 
32 #include "G4DynamicParticle.hh"
33 #include "G4LMsdGenerator.hh"
35 #include "G4ReactionProduct.hh"
36 #include "G4IonTable.hh"
37 #include "G4NucleiProperties.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4HadFinalState.hh"
40 #include "G4KineticTrack.hh"
41 #include "G4DecayKineticTracks.hh"
42 #include "G4KineticTrackVector.hh"
43 #include "G4Log.hh"
44 
45 
47  : G4HadronicInteraction(name)
48 
49 {
50  fPDGencoding = 0;
51 
52  // theParticleChange = new G4HadFinalState;
53 }
54 
56 {
57  // delete theParticleChange;
58 }
59 
60 void G4LMsdGenerator::ModelDescription(std::ostream& outFile) const
61 {
62  outFile << GetModelName() <<" consists of a "
63  << " string model and a stage to de-excite the excited nuclear fragment."
64  << "\n<p>"
65  << "The string model simulates the interaction of\n"
66  << "an incident hadron with a nucleus, forming \n"
67  << "excited strings, decays these strings into hadrons,\n"
68  << "and leaves an excited nucleus. \n"
69  << "<p>The string model:\n";
70 }
71 
73 //
74 // Particle and kinematical limitation od diffraction dissociation
75 
76 G4bool
78  G4Nucleus& targetNucleus )
79 {
80  G4bool applied = false;
81 
82  if( ( aTrack.GetDefinition() == G4Proton::Proton() ||
83  aTrack.GetDefinition() == G4Neutron::Neutron() ) &&
84  targetNucleus.GetA_asInt() >= 1 &&
85  aTrack.GetKineticEnergy() > 1800*CLHEP::MeV ) // 750*CLHEP::MeV )
86  {
87  applied = true;
88  }
89  else if( ( aTrack.GetDefinition() == G4PionPlus::PionPlus() ||
90  aTrack.GetDefinition() == G4PionMinus::PionMinus() ) &&
91  targetNucleus.GetA_asInt() >= 1 &&
92  aTrack.GetKineticEnergy() > 2340*CLHEP::MeV )
93  {
94  applied = true;
95  }
96  else if( ( aTrack.GetDefinition() == G4KaonPlus::KaonPlus() ||
97  aTrack.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
98  targetNucleus.GetA_asInt() >= 1 &&
99  aTrack.GetKineticEnergy() > 1980*CLHEP::MeV )
100  {
101  applied = true;
102  }
103  return applied;
104 }
105 
107 //
108 // Return dissociated particle products and recoil nucleus
109 
112  G4Nucleus& targetNucleus )
113 {
115 
116  const G4HadProjectile* aParticle = &aTrack;
117  G4double eTkin = aParticle->GetKineticEnergy();
118 
119  if( eTkin <= 1.*CLHEP::GeV )
120  {
122  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
123  return &theParticleChange;
124  }
125 
126  G4int A = targetNucleus.GetA_asInt();
127  G4int Z = targetNucleus.GetZ_asInt();
128 
129  G4double plab = aParticle->GetTotalMomentum();
130  G4double plab2 = plab*plab;
131 
132  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
133  G4double partMass = theParticle->GetPDGMass();
134 
135  G4double oldE = partMass + eTkin;
136 
138  G4double targMass2 = targMass*targMass;
139 
140  G4LorentzVector partLV = aParticle->Get4Momentum();
141 
142  G4double sumE = oldE + targMass;
143  G4double sumE2 = sumE*sumE;
144 
145  G4ThreeVector p1 = partLV.vect();
146  // G4cout<<"p1 = "<<p1<<G4endl;
147  G4ParticleMomentum p1unit = p1.unit();
148 
149  G4double Mx = SampleMx(aParticle); // in GeV
150  G4double t = SampleT( Mx);
151 
152  Mx *= CLHEP::GeV;
153 
154  G4double Mx2 = Mx*Mx;
155 
156  // equation for q|| based on sum-E-P and new invariant mass
157 
158  G4double B = sumE2 + targMass2 - Mx2 - plab2;
159 
160  G4double a = 4*(plab2 - sumE2);
161  G4double b = 4*plab*B;
162  G4double c = B*B - 4*sumE2*targMass2;
163  G4double det2 = b*b - 4*a*c;
164  G4double qLong, det, eRetard; // , x2, x3, e2;
165 
166  if( det2 >= 0.)
167  {
168  det = std::sqrt(det2);
169  qLong = (-b - det)/2./a;
170  eRetard = std::sqrt((plab-qLong)*(plab-qLong)+Mx2);
171  }
172  else
173  {
175  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
176  return &theParticleChange;
177  }
179 
180  plab -= qLong;
181 
182  G4ThreeVector pRetard = plab*p1unit;
183 
184  G4ThreeVector pTarg = p1 - pRetard;
185 
186  G4double eTarg = std::sqrt( targMass2 + pTarg.mag2()); // std::sqrt( targMass*targMass + pTarg.mag2() );
187 
188  G4LorentzVector lvRetard(pRetard, eRetard);
189  G4LorentzVector lvTarg(pTarg, eTarg);
190 
191  lvTarg += lvRetard; // sum LV
192 
193  G4ThreeVector bst = lvTarg.boostVector();
194 
195  lvRetard.boost(-bst); // to CNS
196 
197  G4ThreeVector pCMS = lvRetard.vect();
198  G4double momentumCMS = pCMS.mag();
199  G4double tMax = 4.0*momentumCMS*momentumCMS;
200 
201  if( t > tMax ) t = tMax*G4UniformRand();
202 
203  G4double cost = 1. - 2.0*t/tMax;
204 
205 
206  G4double phi = G4UniformRand()*CLHEP::twopi;
207  G4double sint;
208 
209  if( cost > 1.0 || cost < -1.0 ) //
210  {
211  cost = 1.0;
212  sint = 0.0;
213  }
214  else // normal situation
215  {
216  sint = std::sqrt( (1.0-cost)*(1.0+cost) );
217  }
218  G4ThreeVector v1( sint*std::cos(phi), sint*std::sin(phi), cost);
219 
220  v1 *= momentumCMS;
221 
222  G4LorentzVector lvRes( v1.x(),v1.y(),v1.z(), std::sqrt( momentumCMS*momentumCMS + Mx2));
223 
224  lvRes.boost(bst); // to LS
225 
226  lvTarg -= lvRes;
227 
228  G4double eRecoil = lvTarg.e() - targMass;
229 
230  if( eRecoil > 100.*CLHEP::MeV ) // add recoil nucleus
231  {
232  G4ParticleDefinition * recoilDef = 0;
233 
234  if ( Z == 1 && A == 1 ) { recoilDef = G4Proton::Proton(); }
235  else if ( Z == 1 && A == 2 ) { recoilDef = G4Deuteron::Deuteron(); }
236  else if ( Z == 1 && A == 3 ) { recoilDef = G4Triton::Triton(); }
237  else if ( Z == 2 && A == 3 ) { recoilDef = G4He3::He3(); }
238  else if ( Z == 2 && A == 4 ) { recoilDef = G4Alpha::Alpha(); }
239  else
240  {
241  recoilDef =
243  }
244  G4DynamicParticle * aSec = new G4DynamicParticle( recoilDef, lvTarg);
246  }
247  else if( eRecoil > 0.0 )
248  {
250  }
251 
253  FindParticle(fPDGencoding);
254 
255  // G4cout<<fPDGencoding<<", "<<ddPart->GetParticleName()<<", "<<ddPart->GetPDGMass()<<" MeV; lvRes = "<<lvRes<<G4endl;
256 
257  G4KineticTrack ddkt( ddPart, 0., G4ThreeVector(0.,0.,0.), lvRes);
258  G4KineticTrackVector* ddktv = ddkt.Decay();
259 
260 
261  for( unsigned int i = 0; i < ddktv->size(); i++ ) // add products to partchange
262  {
263  G4DynamicParticle * aNew =
264  new G4DynamicParticle( ddktv->operator[](i)->GetDefinition(),
265  ddktv->operator[](i)->Get4Momentum());
266  // G4cout<<" "<<i<<", "<<aNew->GetDefinition()->GetParticleName()<<", "<<aNew->Get4Momentum()<<G4endl;
267 
269  delete ddktv->operator[](i);
270  }
271  delete ddktv;
272 
273  return &theParticleChange;
274 }
275 
277 //
278 // Sample Mx as Roper resonances, set PDG encoding
279 
281 {
282  G4double Mx=0.;
283  G4int i;
284  G4double rand = G4UniformRand();
285 
286  for( i = 0; i < 60; i++)
287  {
288  if( rand >= fProbMx[i][1] ) break;
289  }
290  if(i <= 0) Mx = fProbMx[0][0];
291  else if(i >= 59) Mx = fProbMx[59][0];
292  else Mx = fProbMx[i][0];
293 
294  if ( Mx <= 1.45 )
295  {
296  if( aParticle->GetDefinition() == G4Proton::Proton() )
297  {
298  Mx = 1.44;
299  fPDGencoding = 12212;
300  }
301  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
302  {
303  Mx = 1.44;
304  fPDGencoding = 12112;
305  }
306  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
307  {
308  Mx = 1.3;
309  fPDGencoding = 100211;
310  }
311  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
312  {
313  Mx = 1.3;
314  fPDGencoding = -100211;
315  }
316  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
317  {
318  Mx = 1.27;
319  fPDGencoding = 10323;
320  }
321  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
322  {
323  Mx = 1.27;
324  fPDGencoding = -10323;
325  }
326  }
327  else if ( Mx <= 1.55 )
328  {
329  if( aParticle->GetDefinition() == G4Proton::Proton() )
330  {
331  Mx = 1.52;
332  fPDGencoding = 2124;
333  }
334  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
335  {
336  Mx = 1.52;
337  fPDGencoding = 1214;
338  }
339  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
340  {
341  Mx = 1.45;
342  fPDGencoding = 10211;
343  }
344  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
345  {
346  Mx = 1.45;
347  fPDGencoding = -10211;
348  }
349  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
350  {
351  Mx = 1.46;
352  fPDGencoding = 100321;
353  }
354  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
355  {
356  Mx = 1.46;
357  fPDGencoding = -100321;
358  }
359  }
360  else
361  {
362  if( aParticle->GetDefinition() == G4Proton::Proton() )
363  {
364  Mx = 1.68;
365  fPDGencoding = 12216;
366  }
367  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
368  {
369  Mx = 1.68;
370  fPDGencoding = 12116;
371  }
372  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
373  {
374  // Mx = 1.67;
375  // fPDGencoding = 10215;
376  Mx = 1.45;
377  fPDGencoding = 10211;
378  }
379  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
380  {
381  // Mx = 1.67; // f0 problems->4pi vmg 20.11.14
382  // fPDGencoding = -10215;
383  Mx = 1.45;
384  fPDGencoding = -10211;
385  }
386  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
387  {
388  Mx = 1.68;
389  fPDGencoding = 30323;
390  }
391  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
392  {
393  Mx = 1.68;
394  fPDGencoding = -30323;
395  }
396  }
397  if(fPDGencoding == 0)
398  {
399  Mx = 1.44;
400  fPDGencoding = 12212;
401  }
402  return Mx;
403 }
404 
406 //
407 // Sample t with kinematic limitations of Mx and Tkin
408 
409 G4double G4LMsdGenerator::SampleT( // const G4HadProjectile* aParticle,
410 G4double Mx)
411 {
412  G4double t=0., b=0.;
413  G4int i;
414 
415  for( i = 0; i < 23; ++i)
416  {
417  if( Mx <= fMxBdata[i][0] ) break;
418  }
419  if( i <= 0 ) b = fMxBdata[0][1];
420  else if( i >= 22 ) b = fMxBdata[22][1];
421  else b = fMxBdata[i][1];
422 
423  G4double rand = G4UniformRand();
424 
425  t = -G4Log(rand)/b;
426 
427  t *= (CLHEP::GeV*CLHEP::GeV); // in G4 internal units
428 
429  return t;
430 }
431 
432 
434 //
435 // Integral spectrum of Mx (GeV)
436 
437 const G4double G4LMsdGenerator::fProbMx[60][2] =
438 {
439  {1.000000e+00, 1.000000e+00},
440  {1.025000e+00, 1.000000e+00},
441  {1.050000e+00, 1.000000e+00},
442  {1.075000e+00, 1.000000e+00},
443  {1.100000e+00, 9.975067e-01},
444  {1.125000e+00, 9.934020e-01},
445  {1.150000e+00, 9.878333e-01},
446  {1.175000e+00, 9.805002e-01},
447  {1.200000e+00, 9.716846e-01},
448  {1.225000e+00, 9.604761e-01},
449  {1.250000e+00, 9.452960e-01},
450  {1.275000e+00, 9.265278e-01},
451  {1.300000e+00, 9.053632e-01},
452  {1.325000e+00, 8.775566e-01},
453  {1.350000e+00, 8.441969e-01},
454  {1.375000e+00, 8.076336e-01},
455  {1.400000e+00, 7.682520e-01},
456  {1.425000e+00, 7.238306e-01},
457  {1.450000e+00, 6.769306e-01},
458  {1.475000e+00, 6.303898e-01},
459  {1.500000e+00, 5.824632e-01},
460  {1.525000e+00, 5.340696e-01},
461  {1.550000e+00, 4.873736e-01},
462  {1.575000e+00, 4.422901e-01},
463  {1.600000e+00, 3.988443e-01},
464  {1.625000e+00, 3.583727e-01},
465  {1.650000e+00, 3.205405e-01},
466  {1.675000e+00, 2.856655e-01},
467  {1.700000e+00, 2.537508e-01},
468  {1.725000e+00, 2.247863e-01},
469  {1.750000e+00, 1.985798e-01},
470  {1.775000e+00, 1.750252e-01},
471  {1.800000e+00, 1.539777e-01},
472  {1.825000e+00, 1.352741e-01},
473  {1.850000e+00, 1.187157e-01},
474  {1.875000e+00, 1.040918e-01},
475  {1.900000e+00, 9.118422e-02},
476  {1.925000e+00, 7.980909e-02},
477  {1.950000e+00, 6.979378e-02},
478  {1.975000e+00, 6.097771e-02},
479  {2.000000e+00, 5.322122e-02},
480  {2.025000e+00, 4.639628e-02},
481  {2.050000e+00, 4.039012e-02},
482  {2.075000e+00, 3.510275e-02},
483  {2.100000e+00, 3.044533e-02},
484  {2.125000e+00, 2.633929e-02},
485  {2.150000e+00, 2.271542e-02},
486  {2.175000e+00, 1.951295e-02},
487  {2.200000e+00, 1.667873e-02},
488  {2.225000e+00, 1.416633e-02},
489  {2.250000e+00, 1.193533e-02},
490  {2.275000e+00, 9.950570e-03},
491  {2.300000e+00, 8.181515e-03},
492  {2.325000e+00, 6.601664e-03},
493  {2.350000e+00, 5.188025e-03},
494  {2.375000e+00, 3.920655e-03},
495  {2.400000e+00, 2.782246e-03},
496  {2.425000e+00, 1.757765e-03},
497  {2.450000e+00, 8.341435e-04},
498  {2.475000e+00, 0.000000e+00}
499 };
500 
502 //
503 // Slope b (1/GeV/GeV) vs Mx (GeV) for t-sampling over exp(-b*t)
504 
505 const G4double G4LMsdGenerator::fMxBdata[23][2] =
506 {
507  {1.09014, 17.8620},
508  {1.12590, 19.2831},
509  {1.18549, 17.6907},
510  {1.21693, 16.4760},
511  {1.25194, 15.3867},
512  {1.26932, 14.4236},
513  {1.29019, 13.2931},
514  {1.30755, 12.2882},
515  {1.31790, 11.4509},
516  {1.33888, 10.6969},
517  {1.34911, 9.44130},
518  {1.37711, 8.56148},
519  {1.39101, 7.76593},
520  {1.42608, 6.88582},
521  {1.48593, 6.13019},
522  {1.53179, 5.87723},
523  {1.58111, 5.37308},
524  {1.64105, 4.95217},
525  {1.69037, 4.44803},
526  {1.81742, 3.89879},
527  {1.88096, 3.68693},
528  {1.95509, 3.43278},
529  {2.02219, 3.30445}
530 };
531 
532 
533 
534 //
535 //
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static const double MeV
Definition: G4SIunits.hh:193
static G4double GetNuclearMass(const G4double A, const G4double Z)
CLHEP::Hep3Vector G4ThreeVector
G4KineticTrackVector * Decay()
G4String name
Definition: TRTMaterials.hh:40
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:463
G4double a
Definition: TRTMaterials.hh:39
static const G4double fProbMx[60][2]
const G4String & GetModelName() const
int G4int
Definition: G4Types.hh:78
G4bool IsApplicable(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
void SetStatusChange(G4HadFinalStateStatus aS)
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:93
const G4ParticleDefinition * GetDefinition() const
G4double SampleMx(const G4HadProjectile *aParticle)
G4double SampleT(G4double Mx)
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
const G4double p1
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static const double GeV
Definition: G4SIunits.hh:196
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
void ModelDescription(std::ostream &outFile) const
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4double G4Log(G4double x)
Definition: G4Log.hh:230
void SetEnergyChange(G4double anEnergy)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4LMsdGenerator(const G4String &name="LMsdGenerator")
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void SetLocalEnergyDeposit(G4double aE)
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
double G4double
Definition: G4Types.hh:76
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4ThreeVector G4ParticleMomentum
static G4He3 * He3()
Definition: G4He3.cc:94
void SetMomentumChange(const G4ThreeVector &aV)
static const G4double fMxBdata[23][2]
G4double GetTotalMomentum() const
CLHEP::HepLorentzVector G4LorentzVector
G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)