Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLDeltaProductionChannel.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 // INCL++ intra-nuclear cascade model
27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Davide Mancusi, CEA
29 // Alain Boudard, CEA
30 // Sylvie Leray, CEA
31 // Joseph Cugnon, University of Liege
32 //
33 #define INCLXX_IN_GEANT4_MODE 1
34 
35 #include "globals.hh"
36 
37 #include "G4INCLDeltaProductionChannel.hh"
38 #include "G4INCLKinematicsUtils.hh"
39 #include "G4INCLBinaryCollisionAvatar.hh"
40 #include "G4INCLRandom.hh"
41 #include "G4INCLGlobals.hh"
42 #include "G4INCLLogger.hh"
43 
44 namespace G4INCL {
45 
46  DeltaProductionChannel::DeltaProductionChannel(Particle *p1,
47  Particle *p2,
48  Nucleus *n)
49  :theNucleus(n), particle1(p1), particle2(p2)
50  {}
51 
52  DeltaProductionChannel::~DeltaProductionChannel() {}
53 
54  G4double DeltaProductionChannel::sampleDeltaMass(G4double ecm) {
55  const G4double ramass = 0.0;
56  const G4int maxTries = 100000;
57  G4int nTries = 0;
58  deltaProd101: G4double rndm = Random::shoot();
59  nTries++;
60  G4double y = std::tan(Math::pi*(rndm-0.5));
61  G4double x = 1232.+0.5*130.*y+ramass;
62  if (x < ParticleTable::effectiveDeltaDecayThreshold && (nTries < maxTries))
63  goto deltaProd101;
64  if (ecm < x + ParticleTable::effectiveNucleonMass + 1.0 && (nTries < maxTries)) goto deltaProd101;
65 
66  // generation of the delta mass with the penetration factor
67  // (see prc56(1997)2431)
68  y=ecm*ecm;
69  G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
70  G4double q3=std::pow(std::sqrt(q2), 3.);
71  G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3
72  y=x*x;
73  q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
74  q3=std::pow(std::sqrt(q2), 3.);
75  G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3
76  rndm = Random::shoot();
77  if (rndm > f3/f3max && (nTries < maxTries)) goto deltaProd101;
78  if(nTries >= maxTries) {
79  WARN("DeltaProductionChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Delta mass " << x << " MeV with CM energy " << ecm << " MeV may be unphysical." << std::endl);
80  }
81  return x;
82  }
83 
84  FinalState* DeltaProductionChannel::getFinalState() {
93  // 100 IF (K4.NE.1) GO TO 101 // ThA K4 = 2 by default
94  // ParticleType p1TypeOld = particle1->getType();
95  // ParticleType p2TypeOld = particle2->getType();
96  G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2);
97 
98  const G4int isospin = ParticleTable::getIsospin(particle1->getType()) +
99  ParticleTable::getIsospin(particle2->getType());
100 
101  // Calculate the outcome of the channel:
102  G4double pin = particle1->getMomentum().mag();
103  G4double rndm = 0.0, b = 0.0;
104 
105  G4double xmdel = sampleDeltaMass(ecm);
106  // deltaProduction103: // This label is not used
107  G4double pnorm = KinematicsUtils::momentumInCM(ecm, ParticleTable::effectiveNucleonMass, xmdel);
108  if (pnorm <= 0.0) pnorm=0.000001;
109  G4int index=0;
110  G4int index2=0;
111  rndm = Random::shoot();
112  if (rndm < 0.5) index=1;
113  if (isospin == 0) { // pn case
114  rndm = Random::shoot();
115  if (rndm < 0.5) index2=1;
116  }
117 
118  // G4double x=0.001*0.5*ecm*std::sqrt(ecm*ecm-4.*ParticleTable::effectiveNucleonMass2)
119  // / ParticleTable::effectiveNucleonMass;
120  G4double x = 0.001 * KinematicsUtils::momentumInLab(ecm*ecm, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
121  if(x < 1.4) {
122  b=(5.287/(1.+std::exp((1.3-x)/0.05)))*1.e-6;
123  } else {
124  b=(4.65+0.706*(x-1.4))*1.e-6;
125  }
126  G4double xkh = 2.*b*pin*pnorm;
127  rndm = Random::shoot();
128  G4double ctet=1.0+std::log(1.-rndm*(1.-std::exp(-2.*xkh)))/xkh;
129  if(std::abs(ctet) > 1.0) ctet = Math::sign(ctet);
130  G4double stet = std::sqrt(1.-ctet*ctet);
131 
132  rndm = Random::shoot();
133  G4double fi = Math::twoPi*rndm;
134  G4double cfi = std::cos(fi);
135  G4double sfi = std::sin(fi);
136  // delta production: correction of the angular distribution 02/09/02
137 
138  G4double xx = particle1->getMomentum().perp2();
139  G4double zz = std::pow(particle1->getMomentum().getZ(), 2);
140  G4double xp1, xp2, xp3;
141  if (xx >= zz*1.e-8) {
142  G4double yn = std::sqrt(xx);
143  G4double zn = yn*pin;
144  G4double ex[3], ey[3], ez[3];
145  G4double p1 = particle1->getMomentum().getX();
146  G4double p2 = particle1->getMomentum().getY();
147  G4double p3 = particle1->getMomentum().getZ();
148  ez[0] = p1/pin;
149  ez[1] = p2/pin;
150  ez[2] = p3/pin;
151  ex[0] = p2/yn;
152  ex[1] = -p1/yn;
153  ex[2] = 0.0;
154  ey[0] = p1*p3/zn;
155  ey[1] = p2*p3/zn;
156  ey[2] = -xx/zn;
157  xp1 = (ex[0]*cfi*stet+ey[0]*sfi*stet+ez[0]*ctet)*pnorm;
158  xp2 = (ex[1]*cfi*stet+ey[1]*sfi*stet+ez[1]*ctet)*pnorm;
159  xp3 = (ex[2]*cfi*stet+ey[2]*sfi*stet+ez[2]*ctet)*pnorm;
160  }else {
161  xp1=pnorm*stet*cfi;
162  xp2=pnorm*stet*sfi;
163  xp3=pnorm*ctet;
164  }
165  // end of correction angular distribution of delta production
166  G4double e3 = std::sqrt(xp1*xp1+xp2*xp2+xp3*xp3
167  +ParticleTable::effectiveNucleonMass2);
168  // if(k4.ne.0) go to 161
169 
170  // long-lived delta
171  G4int m1 = 0;
172  G4int m2 = 0;
173  if (index != 1) {
174  ThreeVector mom(xp1, xp2, xp3);
175  particle1->setMomentum(mom);
176  // e1=ecm-eout1
177  m1=1;
178  } else {
179  ThreeVector mom(-xp1, -xp2, -xp3);
180  particle1->setMomentum(mom);
181  // e1=ecm-eout1
182  m1=1;
183  }
184 
185  particle1->setEnergy(ecm - e3);
186  particle2->setEnergy(e3);
187  particle2->setMomentum(-particle1->getMomentum());
188 
189  // SYMMETRIZATION OF CHARGES IN pn -> N DELTA
190  // THE TEST ON "INDEX" ABOVE SYMETRIZES THE EXCITATION OF ONE
191  // OF THE NUCLEONS WITH RESPECT TO THE DELTA EXCITATION
192  // (SEE NOTE 16/10/97)
193  G4int is1 = ParticleTable::getIsospin(particle1->getType());
194  G4int is2 = ParticleTable::getIsospin(particle2->getType());
195  if (isospin == 0) {
196  if(index2 == 1) {
197  G4int isi=is1;
198  is1=is2;
199  is2=isi;
200  }
201  particle1->setHelicity(0.0);
202  } else {
203  rndm = Random::shoot();
204  if (rndm >= 0.25) {
205  is1=3*is1*m1-(1-m1)*is1;
206  is2=3*is2*m2-(1-m2)*is2;
207  }
208  particle1->setHelicity(ctet*ctet);
209  }
210 
211  if(is1 == ParticleTable::getIsospin(Proton) && m1 == 0) {
212  particle1->setType(Proton);
213  } else if(is1 == ParticleTable::getIsospin(Neutron) && m1 == 0) {
214  particle1->setType(Neutron);
215  } else if(is1 == ParticleTable::getIsospin(DeltaMinus) && m1 == 1) {
216  particle1->setType(DeltaMinus);
217  } else if(is1 == ParticleTable::getIsospin(DeltaZero) && m1 == 1) {
218  particle1->setType(DeltaZero);
219  } else if(is1 == ParticleTable::getIsospin(DeltaPlus) && m1 == 1) {
220  particle1->setType(DeltaPlus);
221  } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus) && m1 == 1) {
222  particle1->setType(DeltaPlusPlus);
223  }
224 
225  if(is2 == ParticleTable::getIsospin(Proton) && m2 == 0) {
226  particle2->setType(Proton);
227  } else if(is2 == ParticleTable::getIsospin(Neutron) && m2 == 0) {
228  particle2->setType(Neutron);
229  } else if(is2 == ParticleTable::getIsospin(DeltaMinus) && m2 == 1) {
230  particle2->setType(DeltaMinus);
231  } else if(is2 == ParticleTable::getIsospin(DeltaZero) && m2 == 1) {
232  particle2->setType(DeltaZero);
233  } else if(is2 == ParticleTable::getIsospin(DeltaPlus) && m2 == 1) {
234  particle2->setType(DeltaPlus);
235  } else if(is2 == ParticleTable::getIsospin(DeltaPlusPlus) && m2 == 1) {
236  particle2->setType(DeltaPlusPlus);
237  }
238 
239  if(particle1->isDelta()) particle1->setMass(xmdel);
240  if(particle2->isDelta()) particle2->setMass(xmdel);
241 
242  FinalState *fs = new FinalState;
243  fs->addModifiedParticle(particle1);
244  fs->addModifiedParticle(particle2);
245  return fs;
246  }
247 }
Generated on Sat May 25 2013 14:33:59 for Geant4 by   doxygen 1.8.4