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