| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France << 32 //
>> 33 // INCL++ revision: v5.0_rc3
33 // 34 //
34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1
35 36
36 #include "globals.hh" 37 #include "globals.hh"
37 38
38 #include "G4INCLDeltaDecayChannel.hh" 39 #include "G4INCLDeltaDecayChannel.hh"
39 #include "G4INCLKinematicsUtils.hh" 40 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh" 41 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh" 42 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh" 43 #include "G4INCLGlobals.hh"
43 44
44 namespace G4INCL { 45 namespace G4INCL {
45 46
46 DeltaDecayChannel::DeltaDecayChannel(Particl << 47 DeltaDecayChannel::DeltaDecayChannel(Nucleus *n, Particle *p, ThreeVector const dir)
47 :theParticle(p), incidentDirection(dir) << 48 :theParticle(p), theNucleus(n), incidentDirection(dir)
48 { } 49 { }
49 50
50 DeltaDecayChannel::~DeltaDecayChannel() {} 51 DeltaDecayChannel::~DeltaDecayChannel() {}
51 52
52 G4double DeltaDecayChannel::computeDecayTime 53 G4double DeltaDecayChannel::computeDecayTime(Particle *p) {
53 const G4double m = p->getMass(); 54 const G4double m = p->getMass();
54 const G4double g0 = 115.0; 55 const G4double g0 = 115.0;
55 G4double gg = g0; 56 G4double gg = g0;
56 if(m > 1500.0) gg = 200.0; 57 if(m > 1500.0) gg = 200.0;
57 const G4double geff = p->getEnergy()/m; 58 const G4double geff = p->getEnergy()/m;
58 const G4double qqq = KinematicsUtils::mome 59 const G4double qqq = KinematicsUtils::momentumInCM(m, ParticleTable::effectiveNucleonMass, ParticleTable::effectivePionMass);
59 const G4double psf = std::pow(qqq, 3)/(std << 60 const G4double psf = std::pow(qqq, 3)/(std::pow(qqq, 3) + 5832000.0);
60 const G4double tdel = -G4INCL::PhysicalCon << 61 const G4double tdel = -G4INCL::hc/(gg*psf)*std::log(Random::shoot())*geff;
61 if( m > 1400) return tdel * 1./(1. + std:: << 62 return tdel;
62 return tdel; // fm <<
63 } 63 }
64 64
65 void DeltaDecayChannel::sampleAngles(G4doubl 65 void DeltaDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
66 const G4double hel = theParticle->getHelic 66 const G4double hel = theParticle->getHelicity();
67 unsigned long loopCounter = 0; <<
68 const unsigned long maxLoopCounter = 10000 <<
69 do { 67 do {
70 (*ctet_par) = -1.0 + 2.0*Random::shoot() 68 (*ctet_par) = -1.0 + 2.0*Random::shoot();
71 if(std::abs(*ctet_par) > 1.0) (*ctet_par 69 if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
72 ++loopCounter; << 70 } while(Random::shoot() > ((1.0 + 3.0 * hel * (*ctet_par) * (*ctet_par))
73 } while(loopCounter<maxLoopCounter && Rand << 71 /(1.0 + 3.0 * hel)));
74 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*c 72 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
75 (*phi_par) = Math::twoPi * Random::shoot() 73 (*phi_par) = Math::twoPi * Random::shoot();
76 } 74 }
77 75
78 void DeltaDecayChannel::fillFinalState(Final << 76 FinalState* DeltaDecayChannel::getFinalState() {
79 // SUBROUTINE DECAY2(P1,P2,P3,WP,ij, 77 // SUBROUTINE DECAY2(P1,P2,P3,WP,ij,
80 // s X1,X2,hel,B1,B2,B3) 78 // s X1,X2,hel,B1,B2,B3)
81 79
82 // This routine describes the anisotropic 80 // This routine describes the anisotropic decay of a particle of mass
83 // xi into 2 particles of masses x1,x2. << 81 // xi G4into 2 particles of masses x1,x2.
84 // The anisotropy is supposed to follow a 82 // The anisotropy is supposed to follow a 1+3*hel*(cos(theta))**2
85 // law with respect to the direction of th 83 // law with respect to the direction of the incoming particle.
86 // In the input, p1,p2,p3 is the momentum 84 // In the input, p1,p2,p3 is the momentum of particle xi.
87 // In the output, p1,p2,p3 is the momentum 85 // In the output, p1,p2,p3 is the momentum of particle x1 , while
88 // q1,q2,q3 is the momentum of particle x2 86 // q1,q2,q3 is the momentum of particle x2.
89 87
90 // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200) 88 // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200),QQ4(200),
91 // s YY1(200),YY2(200),YY3(200) 89 // s YY1(200),YY2(200),YY3(200),YM(200),IPI(200)
92 // common/hazard/ial,IY1,IY2,IY3,IY4,IY5 90 // common/hazard/ial,IY1,IY2,IY3,IY4,IY5,IY6,IY7,IY8,IY9,IY10,
93 // s IY11,IY12,IY13,IY14,IY1 91 // s IY11,IY12,IY13,IY14,IY15,IY16,IY17,IY18,IY19
94 92
95 // DATA IY8,IY9,IY10/82345,92345,45681/ << 93 // DATA IY8,IY9,IY10/82345,92345,45681/
96 // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E* 94 // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E**2-(A-C)**2))/E P-N20800
97 // XI=YM(ij) 95 // XI=YM(ij)
98 96
99 // XE=WP 97 // XE=WP P-N20810
100 // B1=P1/XE 98 // B1=P1/XE P-N20820
101 // B2=P2/XE 99 // B2=P2/XE P-N20830
102 // B3=P3/XE << 100 // B3=P3/XE
103 // XQ=PCM(XI,X1,X2) << 101 // XQ=PCM(XI,X1,X2)
104 102
105 const G4double deltaMass = theParticle->ge 103 const G4double deltaMass = theParticle->getMass();
106 104
107 G4double fi, ctet, stet; 105 G4double fi, ctet, stet;
108 sampleAngles(&ctet, &stet, &fi); 106 sampleAngles(&ctet, &stet, &fi);
109 107
110 G4double cfi = std::cos(fi); 108 G4double cfi = std::cos(fi);
111 G4double sfi = std::sin(fi); 109 G4double sfi = std::sin(fi);
112 G4double beta = incidentDirection.mag(); 110 G4double beta = incidentDirection.mag();
113 111
114 G4double q1, q2, q3; 112 G4double q1, q2, q3;
115 G4double sal=0.0; 113 G4double sal=0.0;
116 if (beta >= 1.0e-10) 114 if (beta >= 1.0e-10)
117 sal = incidentDirection.perp()/beta; 115 sal = incidentDirection.perp()/beta;
118 if (sal >= 1.0e-6) { 116 if (sal >= 1.0e-6) {
119 G4double b1 = incidentDirection.getX(); 117 G4double b1 = incidentDirection.getX();
120 G4double b2 = incidentDirection.getY(); 118 G4double b2 = incidentDirection.getY();
121 G4double b3 = incidentDirection.getZ(); 119 G4double b3 = incidentDirection.getZ();
122 G4double cal = b3/beta; 120 G4double cal = b3/beta;
123 G4double t1 = ctet+cal*stet*sfi/sal; 121 G4double t1 = ctet+cal*stet*sfi/sal;
124 G4double t2 = stet/sal; 122 G4double t2 = stet/sal;
125 q1=(b1*t1+b2*t2*cfi)/beta; 123 q1=(b1*t1+b2*t2*cfi)/beta;
126 q2=(b2*t1-b1*t2*cfi)/beta; 124 q2=(b2*t1-b1*t2*cfi)/beta;
127 q3=(b3*t1/beta-t2*sfi); 125 q3=(b3*t1/beta-t2*sfi);
128 } else { 126 } else {
129 q1 = stet*cfi; 127 q1 = stet*cfi;
130 q2 = stet*sfi; 128 q2 = stet*sfi;
131 q3 = ctet; 129 q3 = ctet;
132 } 130 }
133 theParticle->setHelicity(0.0); 131 theParticle->setHelicity(0.0);
134 132
135 ParticleType pionType; 133 ParticleType pionType;
136 #ifdef INCLXX_IN_GEANT4_MODE <<
137 G4int deltaPDGCode = 0; <<
138 #endif <<
139 switch(theParticle->getType()) { 134 switch(theParticle->getType()) {
140 case DeltaPlusPlus: 135 case DeltaPlusPlus:
141 theParticle->setType(Proton); 136 theParticle->setType(Proton);
142 pionType = PiPlus; 137 pionType = PiPlus;
143 #ifdef INCLXX_IN_GEANT4_MODE <<
144 deltaPDGCode = 2224; <<
145 #endif <<
146 break; 138 break;
147 case DeltaPlus: 139 case DeltaPlus:
148 if(Random::shoot() < 1.0/3.0) { 140 if(Random::shoot() < 1.0/3.0) {
149 theParticle->setType(Neutron); 141 theParticle->setType(Neutron);
150 pionType = PiPlus; 142 pionType = PiPlus;
151 } else { 143 } else {
152 theParticle->setType(Proton); 144 theParticle->setType(Proton);
153 pionType = PiZero; 145 pionType = PiZero;
154 } 146 }
155 #ifdef INCLXX_IN_GEANT4_MODE <<
156 deltaPDGCode = 2214; <<
157 #endif <<
158 break; 147 break;
159 case DeltaZero: 148 case DeltaZero:
160 if(Random::shoot() < 1.0/3.0) { 149 if(Random::shoot() < 1.0/3.0) {
161 theParticle->setType(Proton); 150 theParticle->setType(Proton);
162 pionType = PiMinus; 151 pionType = PiMinus;
163 } else { 152 } else {
164 theParticle->setType(Neutron); 153 theParticle->setType(Neutron);
165 pionType = PiZero; 154 pionType = PiZero;
166 } 155 }
167 #ifdef INCLXX_IN_GEANT4_MODE <<
168 deltaPDGCode = 2114; <<
169 #endif <<
170 break; 156 break;
171 case DeltaMinus: 157 case DeltaMinus:
172 theParticle->setType(Neutron); 158 theParticle->setType(Neutron);
173 pionType = PiMinus; 159 pionType = PiMinus;
174 #ifdef INCLXX_IN_GEANT4_MODE <<
175 deltaPDGCode = 1114; <<
176 #endif <<
177 break; 160 break;
178 default: 161 default:
179 INCL_FATAL("Unrecognized delta type; t << 162 FATAL("Unrecognized delta type; type=" << theParticle->getType() << std::endl);
180 pionType = UnknownParticle; << 163 abort();
181 break; 164 break;
182 } 165 }
183 166
184 G4double xq = KinematicsUtils::momentumInC 167 G4double xq = KinematicsUtils::momentumInCM(deltaMass,
185 theParticle->getMass(), 168 theParticle->getMass(),
186 ParticleTable::getINCLMass(pionType)); << 169 ParticleTable::getMass(pionType));
187 170
188 q1 *= xq; 171 q1 *= xq;
189 q2 *= xq; 172 q2 *= xq;
190 q3 *= xq; 173 q3 *= xq;
191 174
192 ThreeVector pionMomentum(q1, q2, q3); 175 ThreeVector pionMomentum(q1, q2, q3);
193 ThreeVector pionPosition(theParticle->getP 176 ThreeVector pionPosition(theParticle->getPosition());
194 Particle *pion = new Particle(pionType, pi 177 Particle *pion = new Particle(pionType, pionMomentum, pionPosition);
195 theParticle->setMomentum(-pionMomentum); 178 theParticle->setMomentum(-pionMomentum);
196 theParticle->adjustEnergyFromMomentum(); 179 theParticle->adjustEnergyFromMomentum();
197 180
198 #ifdef INCLXX_IN_GEANT4_MODE << 181 FinalState *fs = new FinalState;
199 // Set the information about the parent re <<
200 // (as unique, integer ID, we take the rou <<
201 G4int parentResonanceID = static_cast<G4in <<
202 pion->setParentResonancePDGCode(deltaPDGCo <<
203 pion->setParentResonanceID(parentResonance <<
204 theParticle->setParentResonancePDGCode(del <<
205 theParticle->setParentResonanceID(parentRe <<
206 #endif <<
207 <<
208 fs->addModifiedParticle(theParticle); 182 fs->addModifiedParticle(theParticle);
209 fs->addCreatedParticle(pion); 183 fs->addCreatedParticle(pion);
210 // call loren(q1,q2,q3,b1,b2,b3,wq) 184 // call loren(q1,q2,q3,b1,b2,b3,wq)
211 // call loren(p1,p2,p3,b1,b2,b3,wp) 185 // call loren(p1,p2,p3,b1,b2,b3,wp)
212 // qq1(ij)=q1 186 // qq1(ij)=q1
213 // qq2(ij)=q2 187 // qq2(ij)=q2
214 // qq3(ij)=q3 188 // qq3(ij)=q3
215 // qq4(ij)=wq 189 // qq4(ij)=wq
216 // ym(ij)=xi 190 // ym(ij)=xi
217 // RETURN 191 // RETURN P-N21120
218 // END 192 // END P-N21130
>> 193 return fs;
219 } 194 }
220 } 195 }
221 196