Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
86  aTrack.GetKineticEnergy() > 300*CLHEP::MeV
87  ) // 750*CLHEP::MeV )
88  {
89  applied = true;
90  }
91  else if( ( aTrack.GetDefinition() == G4PionPlus::PionPlus() ||
92  aTrack.GetDefinition() == G4PionMinus::PionMinus() ) &&
93  targetNucleus.GetA_asInt() >= 1 &&
94  aTrack.GetKineticEnergy() > 2340*CLHEP::MeV )
95  {
96  applied = true;
97  }
98  else if( ( aTrack.GetDefinition() == G4KaonPlus::KaonPlus() ||
99  aTrack.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
100  targetNucleus.GetA_asInt() >= 1 &&
101  aTrack.GetKineticEnergy() > 1980*CLHEP::MeV )
102  {
103  applied = true;
104  }
105  return applied;
106 }
107 
109 //
110 // Return dissociated particle products and recoil nucleus
111 
114  G4Nucleus& targetNucleus )
115 {
117 
118  const G4HadProjectile* aParticle = &aTrack;
119  G4double eTkin = aParticle->GetKineticEnergy();
120 
121  if( eTkin <= 1.*CLHEP::GeV && aTrack.GetDefinition() != G4Proton::Proton())
122  {
125  return &theParticleChange;
126  }
127 
128  G4int A = targetNucleus.GetA_asInt();
129  G4int Z = targetNucleus.GetZ_asInt();
130 
131  G4double plab = aParticle->GetTotalMomentum();
132  G4double plab2 = plab*plab;
133 
134  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
135  G4double partMass = theParticle->GetPDGMass();
136 
137  G4double oldE = partMass + eTkin;
138 
140  G4double targMass2 = targMass*targMass;
141 
142  G4LorentzVector partLV = aParticle->Get4Momentum();
143 
144  G4double sumE = oldE + targMass;
145  G4double sumE2 = sumE*sumE;
146 
147  G4ThreeVector p1 = partLV.vect();
148  // G4cout<<"p1 = "<<p1<<G4endl;
149  G4ParticleMomentum p1unit = p1.unit();
150 
151  G4double Mx = SampleMx(aParticle); // in GeV
152  G4double t = SampleT( aParticle, Mx); // in GeV
153 
154  Mx *= CLHEP::GeV;
155 
156  G4double Mx2 = Mx*Mx;
157 
158  // equation for q|| based on sum-E-P and new invariant mass
159 
160  G4double B = sumE2 + targMass2 - Mx2 - plab2;
161 
162  G4double a = 4*(plab2 - sumE2);
163  G4double b = 4*plab*B;
164  G4double c = B*B - 4*sumE2*targMass2;
165  G4double det2 = b*b - 4*a*c;
166  G4double qLong, det, eRetard; // , x2, x3, e2;
167 
168  if( det2 >= 0.)
169  {
170  det = std::sqrt(det2);
171  qLong = (-b - det)/2./a;
172  eRetard = std::sqrt((plab-qLong)*(plab-qLong)+Mx2);
173  }
174  else
175  {
178  return &theParticleChange;
179  }
181 
182  plab -= qLong;
183 
184  G4ThreeVector pRetard = plab*p1unit;
185 
186  G4ThreeVector pTarg = p1 - pRetard;
187 
188  G4double eTarg = std::sqrt( targMass2 + pTarg.mag2()); // std::sqrt( targMass*targMass + pTarg.mag2() );
189 
190  G4LorentzVector lvRetard(pRetard, eRetard);
191  G4LorentzVector lvTarg(pTarg, eTarg);
192 
193  lvTarg += lvRetard; // sum LV
194 
195  G4ThreeVector bst = lvTarg.boostVector();
196 
197  lvRetard.boost(-bst); // to CNS
198 
199  G4ThreeVector pCMS = lvRetard.vect();
200  G4double momentumCMS = pCMS.mag();
201  G4double tMax = 4.0*momentumCMS*momentumCMS;
202 
203  if( t > tMax ) t = tMax*G4UniformRand();
204 
205  G4double cost = 1. - 2.0*t/tMax;
206 
207 
209  G4double sint;
210 
211  if( cost > 1.0 || cost < -1.0 ) //
212  {
213  cost = 1.0;
214  sint = 0.0;
215  }
216  else // normal situation
217  {
218  sint = std::sqrt( (1.0-cost)*(1.0+cost) );
219  }
220  G4ThreeVector v1( sint*std::cos(phi), sint*std::sin(phi), cost);
221 
222  v1 *= momentumCMS;
223 
224  G4LorentzVector lvRes( v1.x(),v1.y(),v1.z(), std::sqrt( momentumCMS*momentumCMS + Mx2));
225 
226  lvRes.boost(bst); // to LS
227 
228  lvTarg -= lvRes;
229 
230  G4double eRecoil = lvTarg.e() - targMass;
231 
232  if( eRecoil > 100.*CLHEP::MeV ) // add recoil nucleus
233  {
234  G4ParticleDefinition * recoilDef = 0;
235 
236  if ( Z == 1 && A == 1 ) { recoilDef = G4Proton::Proton(); }
237  else if ( Z == 1 && A == 2 ) { recoilDef = G4Deuteron::Deuteron(); }
238  else if ( Z == 1 && A == 3 ) { recoilDef = G4Triton::Triton(); }
239  else if ( Z == 2 && A == 3 ) { recoilDef = G4He3::He3(); }
240  else if ( Z == 2 && A == 4 ) { recoilDef = G4Alpha::Alpha(); }
241  else
242  {
243  recoilDef =
245  }
246  G4DynamicParticle * aSec = new G4DynamicParticle( recoilDef, lvTarg);
248  }
249  else if( eRecoil > 0.0 )
250  {
252  }
253 
255  FindParticle(fPDGencoding);
256 
257  // G4cout<<fPDGencoding<<", "<<ddPart->GetParticleName()<<", "<<ddPart->GetPDGMass()<<" MeV; lvRes = "<<lvRes<<G4endl;
258 
259  // G4DynamicParticle * aRes = new G4DynamicParticle( ddPart, lvRes);
260  // theParticleChange.AddSecondary(aRes); // simply return resonance
261 
262 
263 
264  // Recursive decay using methods of G4KineticTrack
265 
266  G4KineticTrack ddkt( ddPart, 0., G4ThreeVector(0.,0.,0.), lvRes);
267  G4KineticTrackVector* ddktv = ddkt.Decay();
268 
269  G4DecayKineticTracks decay( ddktv );
270 
271  for( unsigned int i = 0; i < ddktv->size(); i++ ) // add products to partchange
272  {
273  G4DynamicParticle * aNew =
274  new G4DynamicParticle( ddktv->operator[](i)->GetDefinition(),
275  ddktv->operator[](i)->Get4Momentum());
276 
277  // G4cout<<" "<<i<<", "<<aNew->GetDefinition()->GetParticleName()<<", "<<aNew->Get4Momentum()<<G4endl;
278 
280  delete ddktv->operator[](i);
281  }
282  delete ddktv;
283 
284  return &theParticleChange;
285 }
286 
288 //
289 // Sample Mx as Roper resonances, set PDG encoding
290 
292 {
293  G4double Mx = 0.;
294  G4int i;
295  G4double rand = G4UniformRand();
296 
297  for( i = 0; i < 60; i++)
298  {
299  if( rand >= fProbMx[i][1] ) break;
300  }
301  if(i <= 0) Mx = fProbMx[0][0];
302  else if(i >= 59) Mx = fProbMx[59][0];
303  else Mx = fProbMx[i][0];
304 
305  fPDGencoding = 0;
306 
307  if ( Mx <= 1.45 )
308  {
309  if( aParticle->GetDefinition() == G4Proton::Proton() )
310  {
311  Mx = 1.44;
312  // fPDGencoding = 12212;
313  fPDGencoding = 2214;
314  }
315  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
316  {
317  Mx = 1.44;
318  fPDGencoding = 12112;
319  }
320  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
321  {
322  // Mx = 1.3;
323  // fPDGencoding = 100211;
324  Mx = 1.26;
325  fPDGencoding = 20213; // a1(1260)+
326  }
327  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
328  {
329  // Mx = 1.3;
330  // fPDGencoding = -100211;
331  Mx = 1.26;
332  fPDGencoding = -20213; // a1(1260)-
333  }
334  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
335  {
336  Mx = 1.27;
337  fPDGencoding = 10323;
338  }
339  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
340  {
341  Mx = 1.27;
342  fPDGencoding = -10323;
343  }
344  }
345  else if ( Mx <= 1.55 )
346  {
347  if( aParticle->GetDefinition() == G4Proton::Proton() )
348  {
349  Mx = 1.52;
350  // fPDGencoding = 2124;
351  fPDGencoding = 2214;
352  }
353  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
354  {
355  Mx = 1.52;
356  fPDGencoding = 1214;
357  }
358  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
359  {
360  // Mx = 1.45;
361  // fPDGencoding = 10211;
362  Mx = 1.32;
363  fPDGencoding = 215; // a2(1320)+
364  }
365  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
366  {
367  // Mx = 1.45;
368  // fPDGencoding = -10211;
369  Mx = 1.32;
370  fPDGencoding = -215; // a2(1320)-
371  }
372  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
373  {
374  Mx = 1.46;
375  fPDGencoding = 100321;
376  }
377  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
378  {
379  Mx = 1.46;
380  fPDGencoding = -100321;
381  }
382  }
383  else
384  {
385  if( aParticle->GetDefinition() == G4Proton::Proton() )
386  {
387  Mx = 1.68;
388  // fPDGencoding = 12216;
389  fPDGencoding = 2214;
390  }
391  else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
392  {
393  Mx = 1.68;
394  fPDGencoding = 12116;
395  }
396  else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
397  {
398  Mx = 1.67;
399  fPDGencoding = 10215; // pi2(1670)+
400  // Mx = 1.45;
401  // fPDGencoding = 10211;
402  }
403  else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
404  {
405  Mx = 1.67; // f0 problems->4pi vmg 20.11.14
406  fPDGencoding = -10215; // pi2(1670)-
407  // Mx = 1.45;
408  // fPDGencoding = -10211;
409  }
410  else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
411  {
412  Mx = 1.68;
413  fPDGencoding = 30323;
414  }
415  else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
416  {
417  Mx = 1.68;
418  fPDGencoding = -30323;
419  }
420  }
421  if(fPDGencoding == 0)
422  {
423  Mx = 1.44;
424  // fPDGencoding = 12212;
425  fPDGencoding = 2214;
426  }
427  G4ParticleDefinition* myResonance =
429 
430  if ( myResonance ) Mx = myResonance->GetPDGMass();
431 
432  // G4cout<<"PDG-ID = "<<fPDGencoding<<"; with mass = "<<Mx/CLHEP::GeV<<" GeV"<<G4endl;
433 
434  return Mx/CLHEP::GeV;
435 }
436 
438 //
439 // Sample t with kinematic limitations of Mx and Tkin
440 
442 G4double Mx)
443 {
444  G4double t=0., b=0., rTkin = 50.*CLHEP::GeV, eTkin = aParticle->GetKineticEnergy();
445  G4int i;
446 
447  for( i = 0; i < 23; ++i)
448  {
449  if( Mx <= fMxBdata[i][0] ) break;
450  }
451  if( i <= 0 ) b = fMxBdata[0][1];
452  else if( i >= 22 ) b = fMxBdata[22][1];
453  else b = fMxBdata[i][1];
454 
455  if( eTkin > rTkin ) b *= 1. + G4Log(eTkin/rTkin);
456 
457  G4double rand = G4UniformRand();
458 
459  t = -G4Log(rand)/b;
460 
461  t *= (CLHEP::GeV*CLHEP::GeV); // in G4 internal units
462 
463  return t;
464 }
465 
466 
468 //
469 // Integral spectrum of Mx (GeV)
470 
471 const G4double G4LMsdGenerator::fProbMx[60][2] =
472 {
473  {1.000000e+00, 1.000000e+00},
474  {1.025000e+00, 1.000000e+00},
475  {1.050000e+00, 1.000000e+00},
476  {1.075000e+00, 1.000000e+00},
477  {1.100000e+00, 9.975067e-01},
478  {1.125000e+00, 9.934020e-01},
479  {1.150000e+00, 9.878333e-01},
480  {1.175000e+00, 9.805002e-01},
481  {1.200000e+00, 9.716846e-01},
482  {1.225000e+00, 9.604761e-01},
483  {1.250000e+00, 9.452960e-01},
484  {1.275000e+00, 9.265278e-01},
485  {1.300000e+00, 9.053632e-01},
486  {1.325000e+00, 8.775566e-01},
487  {1.350000e+00, 8.441969e-01},
488  {1.375000e+00, 8.076336e-01},
489  {1.400000e+00, 7.682520e-01},
490  {1.425000e+00, 7.238306e-01},
491  {1.450000e+00, 6.769306e-01},
492  {1.475000e+00, 6.303898e-01},
493  {1.500000e+00, 5.824632e-01},
494  {1.525000e+00, 5.340696e-01},
495  {1.550000e+00, 4.873736e-01},
496  {1.575000e+00, 4.422901e-01},
497  {1.600000e+00, 3.988443e-01},
498  {1.625000e+00, 3.583727e-01},
499  {1.650000e+00, 3.205405e-01},
500  {1.675000e+00, 2.856655e-01},
501  {1.700000e+00, 2.537508e-01},
502  {1.725000e+00, 2.247863e-01},
503  {1.750000e+00, 1.985798e-01},
504  {1.775000e+00, 1.750252e-01},
505  {1.800000e+00, 1.539777e-01},
506  {1.825000e+00, 1.352741e-01},
507  {1.850000e+00, 1.187157e-01},
508  {1.875000e+00, 1.040918e-01},
509  {1.900000e+00, 9.118422e-02},
510  {1.925000e+00, 7.980909e-02},
511  {1.950000e+00, 6.979378e-02},
512  {1.975000e+00, 6.097771e-02},
513  {2.000000e+00, 5.322122e-02},
514  {2.025000e+00, 4.639628e-02},
515  {2.050000e+00, 4.039012e-02},
516  {2.075000e+00, 3.510275e-02},
517  {2.100000e+00, 3.044533e-02},
518  {2.125000e+00, 2.633929e-02},
519  {2.150000e+00, 2.271542e-02},
520  {2.175000e+00, 1.951295e-02},
521  {2.200000e+00, 1.667873e-02},
522  {2.225000e+00, 1.416633e-02},
523  {2.250000e+00, 1.193533e-02},
524  {2.275000e+00, 9.950570e-03},
525  {2.300000e+00, 8.181515e-03},
526  {2.325000e+00, 6.601664e-03},
527  {2.350000e+00, 5.188025e-03},
528  {2.375000e+00, 3.920655e-03},
529  {2.400000e+00, 2.782246e-03},
530  {2.425000e+00, 1.757765e-03},
531  {2.450000e+00, 8.341435e-04},
532  {2.475000e+00, 0.000000e+00}
533 };
534 
536 //
537 // Slope b (1/GeV/GeV) vs Mx (GeV) for t-sampling over exp(-b*t)
538 
539 const G4double G4LMsdGenerator::fMxBdata[23][2] =
540 {
541  {1.09014, 17.8620},
542  {1.12590, 19.2831},
543  {1.18549, 17.6907},
544  {1.21693, 16.4760},
545  {1.25194, 15.3867},
546  {1.26932, 14.4236},
547  {1.29019, 13.2931},
548  {1.30755, 12.2882},
549  {1.31790, 11.4509},
550  {1.33888, 10.6969},
551  {1.34911, 9.44130},
552  {1.37711, 8.56148},
553  {1.39101, 7.76593},
554  {1.42608, 6.88582},
555  {1.48593, 6.13019},
556  {1.53179, 5.87723},
557  {1.58111, 5.37308},
558  {1.64105, 4.95217},
559  {1.69037, 4.44803},
560  {1.81742, 3.89879},
561  {1.88096, 3.68693},
562  {1.95509, 3.43278},
563  {2.02219, 3.30445}
564 };
565 
566 
567 
568 //
569 //
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
const XML_Char * name
Definition: expat.h:151
Hep3Vector boostVector() const
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
CLHEP::Hep3Vector G4ThreeVector
double x() const
G4KineticTrackVector * Decay()
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
double B(double temperature)
const G4String & GetModelName() const
int G4int
Definition: G4Types.hh:78
G4bool IsApplicable(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
double z() const
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
void SetStatusChange(G4HadFinalStateStatus aS)
tuple b
Definition: test.py:12
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
const G4ParticleDefinition * GetDefinition() const
G4double SampleMx(const G4HadProjectile *aParticle)
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
G4double GetKineticEnergy() const
static constexpr double MeV
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
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
static constexpr double GeV
void SetEnergyChange(G4double anEnergy)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4LMsdGenerator(const G4String &name="LMsdGenerator")
Hep3Vector unit() const
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
double y() const
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4double SampleT(const G4HadProjectile *aParticle, G4double Mx)
double mag2() const
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
tuple c
Definition: test.py:13
static G4He3 * He3()
Definition: G4He3.cc:94
double mag() const
void SetMomentumChange(const G4ThreeVector &aV)
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4double GetTotalMomentum() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)