| 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; 66 unsigned long loopCounter = 0;
68 const unsigned long maxLoopCounter = 10000 67 const unsigned long maxLoopCounter = 10000000;
69 do { 68 do {
70 (*ctet_par) = -1.0 + 2.0*Random::shoot() 69 (*ctet_par) = -1.0 + 2.0*Random::shoot();
71 if(std::abs(*ctet_par) > 1.0) (*ctet_par 70 if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
72 ++loopCounter; 71 ++loopCounter;
73 } while(loopCounter<maxLoopCounter && Rand 72 } while(loopCounter<maxLoopCounter && Random::shoot() > ((1.0 + 3.0 * hel * (*ctet_par) * (*ctet_par)) /(1.0 + 3.0 * hel))); /* Loop checking, 10.07.2015, D.Mancusi */
74 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*c 73 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
75 (*phi_par) = Math::twoPi * Random::shoot() 74 (*phi_par) = Math::twoPi * Random::shoot();
76 } 75 }
77 76
78 void DeltaDecayChannel::fillFinalState(Final 77 void DeltaDecayChannel::fillFinalState(FinalState *fs) {
79 // SUBROUTINE DECAY2(P1,P2,P3,WP,ij, 78 // SUBROUTINE DECAY2(P1,P2,P3,WP,ij,
80 // s X1,X2,hel,B1,B2,B3) 79 // s X1,X2,hel,B1,B2,B3)
81 80
82 // This routine describes the anisotropic 81 // This routine describes the anisotropic decay of a particle of mass
83 // xi into 2 particles of masses x1,x2. 82 // xi into 2 particles of masses x1,x2.
84 // The anisotropy is supposed to follow a 83 // The anisotropy is supposed to follow a 1+3*hel*(cos(theta))**2
85 // law with respect to the direction of th 84 // law with respect to the direction of the incoming particle.
86 // In the input, p1,p2,p3 is the momentum 85 // In the input, p1,p2,p3 is the momentum of particle xi.
87 // In the output, p1,p2,p3 is the momentum 86 // In the output, p1,p2,p3 is the momentum of particle x1 , while
88 // q1,q2,q3 is the momentum of particle x2 87 // q1,q2,q3 is the momentum of particle x2.
89 88
90 // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200) 89 // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200),QQ4(200),
91 // s YY1(200),YY2(200),YY3(200) 90 // s YY1(200),YY2(200),YY3(200),YM(200),IPI(200)
92 // common/hazard/ial,IY1,IY2,IY3,IY4,IY5 91 // common/hazard/ial,IY1,IY2,IY3,IY4,IY5,IY6,IY7,IY8,IY9,IY10,
93 // s IY11,IY12,IY13,IY14,IY1 92 // s IY11,IY12,IY13,IY14,IY15,IY16,IY17,IY18,IY19
94 93
95 // DATA IY8,IY9,IY10/82345,92345,45681/ 94 // DATA IY8,IY9,IY10/82345,92345,45681/
96 // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E* 95 // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E**2-(A-C)**2))/E P-N20800
97 // XI=YM(ij) 96 // XI=YM(ij)
98 97
99 // XE=WP 98 // XE=WP P-N20810
100 // B1=P1/XE 99 // B1=P1/XE P-N20820
101 // B2=P2/XE 100 // B2=P2/XE P-N20830
102 // B3=P3/XE 101 // B3=P3/XE
103 // XQ=PCM(XI,X1,X2) 102 // XQ=PCM(XI,X1,X2)
104 103
105 const G4double deltaMass = theParticle->ge 104 const G4double deltaMass = theParticle->getMass();
106 105
107 G4double fi, ctet, stet; 106 G4double fi, ctet, stet;
108 sampleAngles(&ctet, &stet, &fi); 107 sampleAngles(&ctet, &stet, &fi);
109 108
110 G4double cfi = std::cos(fi); 109 G4double cfi = std::cos(fi);
111 G4double sfi = std::sin(fi); 110 G4double sfi = std::sin(fi);
112 G4double beta = incidentDirection.mag(); 111 G4double beta = incidentDirection.mag();
113 112
114 G4double q1, q2, q3; 113 G4double q1, q2, q3;
115 G4double sal=0.0; 114 G4double sal=0.0;
116 if (beta >= 1.0e-10) 115 if (beta >= 1.0e-10)
117 sal = incidentDirection.perp()/beta; 116 sal = incidentDirection.perp()/beta;
118 if (sal >= 1.0e-6) { 117 if (sal >= 1.0e-6) {
119 G4double b1 = incidentDirection.getX(); 118 G4double b1 = incidentDirection.getX();
120 G4double b2 = incidentDirection.getY(); 119 G4double b2 = incidentDirection.getY();
121 G4double b3 = incidentDirection.getZ(); 120 G4double b3 = incidentDirection.getZ();
122 G4double cal = b3/beta; 121 G4double cal = b3/beta;
123 G4double t1 = ctet+cal*stet*sfi/sal; 122 G4double t1 = ctet+cal*stet*sfi/sal;
124 G4double t2 = stet/sal; 123 G4double t2 = stet/sal;
125 q1=(b1*t1+b2*t2*cfi)/beta; 124 q1=(b1*t1+b2*t2*cfi)/beta;
126 q2=(b2*t1-b1*t2*cfi)/beta; 125 q2=(b2*t1-b1*t2*cfi)/beta;
127 q3=(b3*t1/beta-t2*sfi); 126 q3=(b3*t1/beta-t2*sfi);
128 } else { 127 } else {
129 q1 = stet*cfi; 128 q1 = stet*cfi;
130 q2 = stet*sfi; 129 q2 = stet*sfi;
131 q3 = ctet; 130 q3 = ctet;
132 } 131 }
133 theParticle->setHelicity(0.0); 132 theParticle->setHelicity(0.0);
134 133
135 ParticleType pionType; 134 ParticleType pionType;
136 #ifdef INCLXX_IN_GEANT4_MODE <<
137 G4int deltaPDGCode = 0; <<
138 #endif <<
139 switch(theParticle->getType()) { 135 switch(theParticle->getType()) {
140 case DeltaPlusPlus: 136 case DeltaPlusPlus:
141 theParticle->setType(Proton); 137 theParticle->setType(Proton);
142 pionType = PiPlus; 138 pionType = PiPlus;
143 #ifdef INCLXX_IN_GEANT4_MODE <<
144 deltaPDGCode = 2224; <<
145 #endif <<
146 break; 139 break;
147 case DeltaPlus: 140 case DeltaPlus:
148 if(Random::shoot() < 1.0/3.0) { 141 if(Random::shoot() < 1.0/3.0) {
149 theParticle->setType(Neutron); 142 theParticle->setType(Neutron);
150 pionType = PiPlus; 143 pionType = PiPlus;
151 } else { 144 } else {
152 theParticle->setType(Proton); 145 theParticle->setType(Proton);
153 pionType = PiZero; 146 pionType = PiZero;
154 } 147 }
155 #ifdef INCLXX_IN_GEANT4_MODE <<
156 deltaPDGCode = 2214; <<
157 #endif <<
158 break; 148 break;
159 case DeltaZero: 149 case DeltaZero:
160 if(Random::shoot() < 1.0/3.0) { 150 if(Random::shoot() < 1.0/3.0) {
161 theParticle->setType(Proton); 151 theParticle->setType(Proton);
162 pionType = PiMinus; 152 pionType = PiMinus;
163 } else { 153 } else {
164 theParticle->setType(Neutron); 154 theParticle->setType(Neutron);
165 pionType = PiZero; 155 pionType = PiZero;
166 } 156 }
167 #ifdef INCLXX_IN_GEANT4_MODE <<
168 deltaPDGCode = 2114; <<
169 #endif <<
170 break; 157 break;
171 case DeltaMinus: 158 case DeltaMinus:
172 theParticle->setType(Neutron); 159 theParticle->setType(Neutron);
173 pionType = PiMinus; 160 pionType = PiMinus;
174 #ifdef INCLXX_IN_GEANT4_MODE <<
175 deltaPDGCode = 1114; <<
176 #endif <<
177 break; 161 break;
178 default: 162 default:
179 INCL_FATAL("Unrecognized delta type; t 163 INCL_FATAL("Unrecognized delta type; type=" << theParticle->getType() << '\n');
180 pionType = UnknownParticle; 164 pionType = UnknownParticle;
181 break; 165 break;
182 } 166 }
183 167
184 G4double xq = KinematicsUtils::momentumInC 168 G4double xq = KinematicsUtils::momentumInCM(deltaMass,
185 theParticle->getMass(), 169 theParticle->getMass(),
186 ParticleTable::getINCLMass(pionType)); 170 ParticleTable::getINCLMass(pionType));
187 171
188 q1 *= xq; 172 q1 *= xq;
189 q2 *= xq; 173 q2 *= xq;
190 q3 *= xq; 174 q3 *= xq;
191 175
192 ThreeVector pionMomentum(q1, q2, q3); 176 ThreeVector pionMomentum(q1, q2, q3);
193 ThreeVector pionPosition(theParticle->getP 177 ThreeVector pionPosition(theParticle->getPosition());
194 Particle *pion = new Particle(pionType, pi 178 Particle *pion = new Particle(pionType, pionMomentum, pionPosition);
195 theParticle->setMomentum(-pionMomentum); 179 theParticle->setMomentum(-pionMomentum);
196 theParticle->adjustEnergyFromMomentum(); 180 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 181
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
219 } 193 }
220 } 194 }
221 195