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