| 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 <<
33 // 32 //
34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1
35 34
36 #include "globals.hh" 35 #include "globals.hh"
37 36
38 #include "G4INCLDeltaDecayChannel.hh" 37 #include "G4INCLDeltaDecayChannel.hh"
39 #include "G4INCLKinematicsUtils.hh" 38 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh" 39 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh" 40 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh" 41 #include "G4INCLGlobals.hh"
43 42
44 namespace G4INCL { 43 namespace G4INCL {
45 44
46 DeltaDecayChannel::DeltaDecayChannel(Particl << 45 DeltaDecayChannel::DeltaDecayChannel(Particle *p, ThreeVector const dir)
47 :theParticle(p), incidentDirection(dir) 46 :theParticle(p), incidentDirection(dir)
48 { } 47 { }
49 48
50 DeltaDecayChannel::~DeltaDecayChannel() {} 49 DeltaDecayChannel::~DeltaDecayChannel() {}
51 50
52 G4double DeltaDecayChannel::computeDecayTime 51 G4double DeltaDecayChannel::computeDecayTime(Particle *p) {
53 const G4double m = p->getMass(); 52 const G4double m = p->getMass();
54 const G4double g0 = 115.0; 53 const G4double g0 = 115.0;
55 G4double gg = g0; 54 G4double gg = g0;
56 if(m > 1500.0) gg = 200.0; 55 if(m > 1500.0) gg = 200.0;
57 const G4double geff = p->getEnergy()/m; 56 const G4double geff = p->getEnergy()/m;
58 const G4double qqq = KinematicsUtils::mome 57 const G4double qqq = KinematicsUtils::momentumInCM(m, ParticleTable::effectiveNucleonMass, ParticleTable::effectivePionMass);
59 const G4double psf = std::pow(qqq, 3)/(std << 58 const G4double psf = std::pow(qqq, 3)/(std::pow(qqq, 3) + 5832000.0);
60 const G4double tdel = -G4INCL::PhysicalCon << 59 const G4double tdel = -G4INCL::PhysicalConstants::hc/(gg*psf)*std::log(Random::shoot())*geff;
61 if( m > 1400) return tdel * 1./(1. + std:: << 60 return tdel;
62 return tdel; // fm <<
63 } 61 }
64 62
65 void DeltaDecayChannel::sampleAngles(G4doubl 63 void DeltaDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
66 const G4double hel = theParticle->getHelic 64 const G4double hel = theParticle->getHelicity();
67 unsigned long loopCounter = 0; <<
68 const unsigned long maxLoopCounter = 10000 <<
69 do { 65 do {
70 (*ctet_par) = -1.0 + 2.0*Random::shoot() 66 (*ctet_par) = -1.0 + 2.0*Random::shoot();
71 if(std::abs(*ctet_par) > 1.0) (*ctet_par 67 if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
72 ++loopCounter; << 68 } while(Random::shoot() > ((1.0 + 3.0 * hel * (*ctet_par) * (*ctet_par))
73 } while(loopCounter<maxLoopCounter && Rand << 69 /(1.0 + 3.0 * hel)));
74 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*c 70 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
75 (*phi_par) = Math::twoPi * Random::shoot() 71 (*phi_par) = Math::twoPi * Random::shoot();
76 } 72 }
77 73
78 void DeltaDecayChannel::fillFinalState(Final << 74 FinalState* DeltaDecayChannel::getFinalState() {
79 // SUBROUTINE DECAY2(P1,P2,P3,WP,ij, 75 // SUBROUTINE DECAY2(P1,P2,P3,WP,ij,
80 // s X1,X2,hel,B1,B2,B3) 76 // s X1,X2,hel,B1,B2,B3)
81 77
82 // This routine describes the anisotropic 78 // This routine describes the anisotropic decay of a particle of mass
83 // xi into 2 particles of masses x1,x2. 79 // xi into 2 particles of masses x1,x2.
84 // The anisotropy is supposed to follow a 80 // The anisotropy is supposed to follow a 1+3*hel*(cos(theta))**2
85 // law with respect to the direction of th 81 // law with respect to the direction of the incoming particle.
86 // In the input, p1,p2,p3 is the momentum 82 // In the input, p1,p2,p3 is the momentum of particle xi.
87 // In the output, p1,p2,p3 is the momentum 83 // In the output, p1,p2,p3 is the momentum of particle x1 , while
88 // q1,q2,q3 is the momentum of particle x2 84 // q1,q2,q3 is the momentum of particle x2.
89 85
90 // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200) 86 // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200),QQ4(200),
91 // s YY1(200),YY2(200),YY3(200) 87 // s YY1(200),YY2(200),YY3(200),YM(200),IPI(200)
92 // common/hazard/ial,IY1,IY2,IY3,IY4,IY5 88 // common/hazard/ial,IY1,IY2,IY3,IY4,IY5,IY6,IY7,IY8,IY9,IY10,
93 // s IY11,IY12,IY13,IY14,IY1 89 // s IY11,IY12,IY13,IY14,IY15,IY16,IY17,IY18,IY19
94 90
95 // DATA IY8,IY9,IY10/82345,92345,45681/ 91 // DATA IY8,IY9,IY10/82345,92345,45681/
96 // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E* 92 // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E**2-(A-C)**2))/E P-N20800
97 // XI=YM(ij) 93 // XI=YM(ij)
98 94
99 // XE=WP 95 // XE=WP P-N20810
100 // B1=P1/XE 96 // B1=P1/XE P-N20820
101 // B2=P2/XE 97 // B2=P2/XE P-N20830
102 // B3=P3/XE 98 // B3=P3/XE
103 // XQ=PCM(XI,X1,X2) 99 // XQ=PCM(XI,X1,X2)
104 100
105 const G4double deltaMass = theParticle->ge 101 const G4double deltaMass = theParticle->getMass();
106 102
107 G4double fi, ctet, stet; 103 G4double fi, ctet, stet;
108 sampleAngles(&ctet, &stet, &fi); 104 sampleAngles(&ctet, &stet, &fi);
109 105
110 G4double cfi = std::cos(fi); 106 G4double cfi = std::cos(fi);
111 G4double sfi = std::sin(fi); 107 G4double sfi = std::sin(fi);
112 G4double beta = incidentDirection.mag(); 108 G4double beta = incidentDirection.mag();
113 109
114 G4double q1, q2, q3; 110 G4double q1, q2, q3;
115 G4double sal=0.0; 111 G4double sal=0.0;
116 if (beta >= 1.0e-10) 112 if (beta >= 1.0e-10)
117 sal = incidentDirection.perp()/beta; 113 sal = incidentDirection.perp()/beta;
118 if (sal >= 1.0e-6) { 114 if (sal >= 1.0e-6) {
119 G4double b1 = incidentDirection.getX(); 115 G4double b1 = incidentDirection.getX();
120 G4double b2 = incidentDirection.getY(); 116 G4double b2 = incidentDirection.getY();
121 G4double b3 = incidentDirection.getZ(); 117 G4double b3 = incidentDirection.getZ();
122 G4double cal = b3/beta; 118 G4double cal = b3/beta;
123 G4double t1 = ctet+cal*stet*sfi/sal; 119 G4double t1 = ctet+cal*stet*sfi/sal;
124 G4double t2 = stet/sal; 120 G4double t2 = stet/sal;
125 q1=(b1*t1+b2*t2*cfi)/beta; 121 q1=(b1*t1+b2*t2*cfi)/beta;
126 q2=(b2*t1-b1*t2*cfi)/beta; 122 q2=(b2*t1-b1*t2*cfi)/beta;
127 q3=(b3*t1/beta-t2*sfi); 123 q3=(b3*t1/beta-t2*sfi);
128 } else { 124 } else {
129 q1 = stet*cfi; 125 q1 = stet*cfi;
130 q2 = stet*sfi; 126 q2 = stet*sfi;
131 q3 = ctet; 127 q3 = ctet;
132 } 128 }
133 theParticle->setHelicity(0.0); 129 theParticle->setHelicity(0.0);
134 130
135 ParticleType pionType; 131 ParticleType pionType;
136 #ifdef INCLXX_IN_GEANT4_MODE <<
137 G4int deltaPDGCode = 0; <<
138 #endif <<
139 switch(theParticle->getType()) { 132 switch(theParticle->getType()) {
140 case DeltaPlusPlus: 133 case DeltaPlusPlus:
141 theParticle->setType(Proton); 134 theParticle->setType(Proton);
142 pionType = PiPlus; 135 pionType = PiPlus;
143 #ifdef INCLXX_IN_GEANT4_MODE <<
144 deltaPDGCode = 2224; <<
145 #endif <<
146 break; 136 break;
147 case DeltaPlus: 137 case DeltaPlus:
148 if(Random::shoot() < 1.0/3.0) { 138 if(Random::shoot() < 1.0/3.0) {
149 theParticle->setType(Neutron); 139 theParticle->setType(Neutron);
150 pionType = PiPlus; 140 pionType = PiPlus;
151 } else { 141 } else {
152 theParticle->setType(Proton); 142 theParticle->setType(Proton);
153 pionType = PiZero; 143 pionType = PiZero;
154 } 144 }
155 #ifdef INCLXX_IN_GEANT4_MODE <<
156 deltaPDGCode = 2214; <<
157 #endif <<
158 break; 145 break;
159 case DeltaZero: 146 case DeltaZero:
160 if(Random::shoot() < 1.0/3.0) { 147 if(Random::shoot() < 1.0/3.0) {
161 theParticle->setType(Proton); 148 theParticle->setType(Proton);
162 pionType = PiMinus; 149 pionType = PiMinus;
163 } else { 150 } else {
164 theParticle->setType(Neutron); 151 theParticle->setType(Neutron);
165 pionType = PiZero; 152 pionType = PiZero;
166 } 153 }
167 #ifdef INCLXX_IN_GEANT4_MODE <<
168 deltaPDGCode = 2114; <<
169 #endif <<
170 break; 154 break;
171 case DeltaMinus: 155 case DeltaMinus:
172 theParticle->setType(Neutron); 156 theParticle->setType(Neutron);
173 pionType = PiMinus; 157 pionType = PiMinus;
174 #ifdef INCLXX_IN_GEANT4_MODE <<
175 deltaPDGCode = 1114; <<
176 #endif <<
177 break; 158 break;
178 default: 159 default:
179 INCL_FATAL("Unrecognized delta type; t << 160 INCL_FATAL("Unrecognized delta type; type=" << theParticle->getType() << std::endl);
180 pionType = UnknownParticle; 161 pionType = UnknownParticle;
181 break; 162 break;
182 } 163 }
183 164
184 G4double xq = KinematicsUtils::momentumInC 165 G4double xq = KinematicsUtils::momentumInCM(deltaMass,
185 theParticle->getMass(), 166 theParticle->getMass(),
186 ParticleTable::getINCLMass(pionType)); 167 ParticleTable::getINCLMass(pionType));
187 168
188 q1 *= xq; 169 q1 *= xq;
189 q2 *= xq; 170 q2 *= xq;
190 q3 *= xq; 171 q3 *= xq;
191 172
192 ThreeVector pionMomentum(q1, q2, q3); 173 ThreeVector pionMomentum(q1, q2, q3);
193 ThreeVector pionPosition(theParticle->getP 174 ThreeVector pionPosition(theParticle->getPosition());
194 Particle *pion = new Particle(pionType, pi 175 Particle *pion = new Particle(pionType, pionMomentum, pionPosition);
195 theParticle->setMomentum(-pionMomentum); 176 theParticle->setMomentum(-pionMomentum);
196 theParticle->adjustEnergyFromMomentum(); 177 theParticle->adjustEnergyFromMomentum();
197 178
198 #ifdef INCLXX_IN_GEANT4_MODE << 179 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); 180 fs->addModifiedParticle(theParticle);
209 fs->addCreatedParticle(pion); 181 fs->addCreatedParticle(pion);
210 // call loren(q1,q2,q3,b1,b2,b3,wq) 182 // call loren(q1,q2,q3,b1,b2,b3,wq)
211 // call loren(p1,p2,p3,b1,b2,b3,wp) 183 // call loren(p1,p2,p3,b1,b2,b3,wp)
212 // qq1(ij)=q1 184 // qq1(ij)=q1
213 // qq2(ij)=q2 185 // qq2(ij)=q2
214 // qq3(ij)=q3 186 // qq3(ij)=q3
215 // qq4(ij)=wq 187 // qq4(ij)=wq
216 // ym(ij)=xi 188 // ym(ij)=xi
217 // RETURN 189 // RETURN P-N21120
218 // END 190 // END P-N21130
>> 191 return fs;
219 } 192 }
220 } 193 }
221 194