Geant4  10.02.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 = 0.; // SampleT( aParticle, Mx); // in GeV
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 
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  G4DynamicParticle * aRes = new G4DynamicParticle( ddPart, lvRes);
258  theParticleChange.AddSecondary(aRes); // simply return resonance
259 
260 
261  /*
262  // Recursive decay using methods of G4KineticTrack
263 
264  G4KineticTrack ddkt( ddPart, 0., G4ThreeVector(0.,0.,0.), lvRes);
265  G4KineticTrackVector* ddktv = ddkt.Decay();
266 
267  G4DecayKineticTracks decay( ddktv );
268 
269  for( unsigned int i = 0; i < ddktv->size(); i++ ) // add products to partchange
270  {
271  G4DynamicParticle * aNew =
272  new G4DynamicParticle( ddktv->operator[](i)->GetDefinition(),
273  ddktv->operator[](i)->Get4Momentum());
274  // G4cout<<" "<<i<<", "<<aNew->GetDefinition()->GetParticleName()<<", "<<aNew->Get4Momentum()<<G4endl;
275 
276  theParticleChange.AddSecondary(aNew);
277  delete ddktv->operator[](i);
278  }
279  delete ddktv;
280  */
281  return &theParticleChange;
282 }
283 
285 //
286 // Sample Mx as Roper resonances, set PDG encoding
287 
289 {
290  G4double Mx = 0.;
291  G4int i;
292  G4double rand = G4UniformRand();
293 
294  for( i = 0; i < 60; i++)
295  {
296  if( rand >= fProbMx[i][1] ) break;
297  }
298  if(i <= 0) Mx = fProbMx[0][0];
299  else if(i >= 59) Mx = fProbMx[59][0];
300  else Mx = fProbMx[i][0];
301 
302  fPDGencoding = 0;
303 
304  if ( Mx <= 1.45 )
305  {
306  if( aParticle->GetDefinition() == G4Proton::Proton() )
307  {
308  Mx = 1.44;
309  fPDGencoding = 12212;
310  }
311  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
312  {
313  Mx = 1.44;
314  fPDGencoding = 12112;
315  }
316  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
317  {
318  // Mx = 1.3;
319  // fPDGencoding = 100211;
320  Mx = 1.26;
321  fPDGencoding = 20213; // a1(1260)+
322  }
323  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
324  {
325  // Mx = 1.3;
326  // fPDGencoding = -100211;
327  Mx = 1.26;
328  fPDGencoding = -20213; // a1(1260)-
329  }
330  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
331  {
332  Mx = 1.27;
333  fPDGencoding = 10323;
334  }
335  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
336  {
337  Mx = 1.27;
338  fPDGencoding = -10323;
339  }
340  }
341  else if ( Mx <= 1.55 )
342  {
343  if( aParticle->GetDefinition() == G4Proton::Proton() )
344  {
345  Mx = 1.52;
346  fPDGencoding = 2124;
347  }
348  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
349  {
350  Mx = 1.52;
351  fPDGencoding = 1214;
352  }
353  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
354  {
355  // Mx = 1.45;
356  // fPDGencoding = 10211;
357  Mx = 1.32;
358  fPDGencoding = 215; // a2(1320)+
359  }
360  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
361  {
362  // Mx = 1.45;
363  // fPDGencoding = -10211;
364  Mx = 1.32;
365  fPDGencoding = -215; // a2(1320)-
366  }
367  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
368  {
369  Mx = 1.46;
370  fPDGencoding = 100321;
371  }
372  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
373  {
374  Mx = 1.46;
375  fPDGencoding = -100321;
376  }
377  }
378  else
379  {
380  if( aParticle->GetDefinition() == G4Proton::Proton() )
381  {
382  Mx = 1.68;
383  fPDGencoding = 12216;
384  }
385  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
386  {
387  Mx = 1.68;
388  fPDGencoding = 12116;
389  }
390  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
391  {
392  Mx = 1.67;
393  fPDGencoding = 10215; // pi2(1670)+
394  // Mx = 1.45;
395  // fPDGencoding = 10211;
396  }
397  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
398  {
399  Mx = 1.67; // f0 problems->4pi vmg 20.11.14
400  fPDGencoding = -10215; // pi2(1670)-
401  // Mx = 1.45;
402  // fPDGencoding = -10211;
403  }
404  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
405  {
406  Mx = 1.68;
407  fPDGencoding = 30323;
408  }
409  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
410  {
411  Mx = 1.68;
412  fPDGencoding = -30323;
413  }
414  }
415  if(fPDGencoding == 0)
416  {
417  Mx = 1.44;
418  fPDGencoding = 12212;
419  }
420  G4ParticleDefinition* myResonance =
422 
423  if ( myResonance ) Mx = myResonance->GetPDGMass();
424 
425  // G4cout<<"PDG-ID = "<<fPDGencoding<<"; with mass = "<<Mx/CLHEP::GeV<<" GeV"<<G4endl;
426 
427  return Mx/CLHEP::GeV;
428 }
429 
431 //
432 // Sample t with kinematic limitations of Mx and Tkin
433 
435 G4double Mx)
436 {
437  G4double t=0., b=0., rTkin = 50.*CLHEP::GeV, eTkin = aParticle->GetKineticEnergy();
438  G4int i;
439 
440  for( i = 0; i < 23; ++i)
441  {
442  if( Mx <= fMxBdata[i][0] ) break;
443  }
444  if( i <= 0 ) b = fMxBdata[0][1];
445  else if( i >= 22 ) b = fMxBdata[22][1];
446  else b = fMxBdata[i][1];
447 
448  if( eTkin > rTkin ) b *= 1. + G4Log(eTkin/rTkin);
449 
450  G4double rand = G4UniformRand();
451 
452  t = -G4Log(rand)/b;
453 
454  t *= (CLHEP::GeV*CLHEP::GeV); // in G4 internal units
455 
456  return t;
457 }
458 
459 
461 //
462 // Integral spectrum of Mx (GeV)
463 
464 const G4double G4LMsdGenerator::fProbMx[60][2] =
465 {
466  {1.000000e+00, 1.000000e+00},
467  {1.025000e+00, 1.000000e+00},
468  {1.050000e+00, 1.000000e+00},
469  {1.075000e+00, 1.000000e+00},
470  {1.100000e+00, 9.975067e-01},
471  {1.125000e+00, 9.934020e-01},
472  {1.150000e+00, 9.878333e-01},
473  {1.175000e+00, 9.805002e-01},
474  {1.200000e+00, 9.716846e-01},
475  {1.225000e+00, 9.604761e-01},
476  {1.250000e+00, 9.452960e-01},
477  {1.275000e+00, 9.265278e-01},
478  {1.300000e+00, 9.053632e-01},
479  {1.325000e+00, 8.775566e-01},
480  {1.350000e+00, 8.441969e-01},
481  {1.375000e+00, 8.076336e-01},
482  {1.400000e+00, 7.682520e-01},
483  {1.425000e+00, 7.238306e-01},
484  {1.450000e+00, 6.769306e-01},
485  {1.475000e+00, 6.303898e-01},
486  {1.500000e+00, 5.824632e-01},
487  {1.525000e+00, 5.340696e-01},
488  {1.550000e+00, 4.873736e-01},
489  {1.575000e+00, 4.422901e-01},
490  {1.600000e+00, 3.988443e-01},
491  {1.625000e+00, 3.583727e-01},
492  {1.650000e+00, 3.205405e-01},
493  {1.675000e+00, 2.856655e-01},
494  {1.700000e+00, 2.537508e-01},
495  {1.725000e+00, 2.247863e-01},
496  {1.750000e+00, 1.985798e-01},
497  {1.775000e+00, 1.750252e-01},
498  {1.800000e+00, 1.539777e-01},
499  {1.825000e+00, 1.352741e-01},
500  {1.850000e+00, 1.187157e-01},
501  {1.875000e+00, 1.040918e-01},
502  {1.900000e+00, 9.118422e-02},
503  {1.925000e+00, 7.980909e-02},
504  {1.950000e+00, 6.979378e-02},
505  {1.975000e+00, 6.097771e-02},
506  {2.000000e+00, 5.322122e-02},
507  {2.025000e+00, 4.639628e-02},
508  {2.050000e+00, 4.039012e-02},
509  {2.075000e+00, 3.510275e-02},
510  {2.100000e+00, 3.044533e-02},
511  {2.125000e+00, 2.633929e-02},
512  {2.150000e+00, 2.271542e-02},
513  {2.175000e+00, 1.951295e-02},
514  {2.200000e+00, 1.667873e-02},
515  {2.225000e+00, 1.416633e-02},
516  {2.250000e+00, 1.193533e-02},
517  {2.275000e+00, 9.950570e-03},
518  {2.300000e+00, 8.181515e-03},
519  {2.325000e+00, 6.601664e-03},
520  {2.350000e+00, 5.188025e-03},
521  {2.375000e+00, 3.920655e-03},
522  {2.400000e+00, 2.782246e-03},
523  {2.425000e+00, 1.757765e-03},
524  {2.450000e+00, 8.341435e-04},
525  {2.475000e+00, 0.000000e+00}
526 };
527 
529 //
530 // Slope b (1/GeV/GeV) vs Mx (GeV) for t-sampling over exp(-b*t)
531 
532 const G4double G4LMsdGenerator::fMxBdata[23][2] =
533 {
534  {1.09014, 17.8620},
535  {1.12590, 19.2831},
536  {1.18549, 17.6907},
537  {1.21693, 16.4760},
538  {1.25194, 15.3867},
539  {1.26932, 14.4236},
540  {1.29019, 13.2931},
541  {1.30755, 12.2882},
542  {1.31790, 11.4509},
543  {1.33888, 10.6969},
544  {1.34911, 9.44130},
545  {1.37711, 8.56148},
546  {1.39101, 7.76593},
547  {1.42608, 6.88582},
548  {1.48593, 6.13019},
549  {1.53179, 5.87723},
550  {1.58111, 5.37308},
551  {1.64105, 4.95217},
552  {1.69037, 4.44803},
553  {1.81742, 3.89879},
554  {1.88096, 3.68693},
555  {1.95509, 3.43278},
556  {2.02219, 3.30445}
557 };
558 
559 
560 
561 //
562 //
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static const double MeV
Definition: G4SIunits.hh:211
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
CLHEP::Hep3Vector G4ThreeVector
G4String name
Definition: TRTMaterials.hh:40
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
G4double a
Definition: TRTMaterials.hh:39
double B(double temperature)
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:97
double A(double temperature)
const G4ParticleDefinition * GetDefinition() const
G4double SampleMx(const G4HadProjectile *aParticle)
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
static const double twopi
Definition: G4SIunits.hh:75
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:214
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
void ModelDescription(std::ostream &outFile) const
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
G4double SampleT(const G4HadProjectile *aParticle, G4double Mx)
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)