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