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