| 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 // $Id: G4Pythia6Decayer.cc 98498 2016-07-20 15:39:59Z ihrivnac $
26 // 27 //
27 /// \file eventgenerator/pythia/decayer6/src/G 28 /// \file eventgenerator/pythia/decayer6/src/G4Pythia6Decayer.cc
28 /// \brief Implementation of the G4Pythia6Deca 29 /// \brief Implementation of the G4Pythia6Decayer class
29 30
30 // ------------------------------------------- 31 // ----------------------------------------------------------------------------
31 // According to TPythia6Decayer class in Root: 32 // According to TPythia6Decayer class in Root:
32 // http://root.cern.ch/ 33 // http://root.cern.ch/
33 // see http://root.cern.ch/root/License.html 34 // see http://root.cern.ch/root/License.html
34 // ------------------------------------------- 35 // ----------------------------------------------------------------------------
35 36
36 #include "G4Pythia6Decayer.hh" 37 #include "G4Pythia6Decayer.hh"
37 <<
38 #include "Pythia6.hh" 38 #include "Pythia6.hh"
39 39
>> 40 #include "G4DynamicParticle.hh"
40 #include "G4DecayProducts.hh" 41 #include "G4DecayProducts.hh"
41 #include "G4DecayTable.hh" 42 #include "G4DecayTable.hh"
42 #include "G4DynamicParticle.hh" <<
43 #include "G4ParticleTable.hh" 43 #include "G4ParticleTable.hh"
44 #include "G4SystemOfUnits.hh" <<
45 #include "G4Track.hh" 44 #include "G4Track.hh"
>> 45 #include "G4SystemOfUnits.hh"
46 46
47 #include <CLHEP/Vector/LorentzVector.h> 47 #include <CLHEP/Vector/LorentzVector.h>
>> 48
48 #include <cmath> 49 #include <cmath>
49 50
50 const EDecayType G4Pythia6Decayer::fgkDefaultD << 51 const EDecayType G4Pythia6Decayer::fgkDefaultDecayType = kAll;
51 52
52 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 54
54 G4Pythia6Decayer::G4Pythia6Decayer() 55 G4Pythia6Decayer::G4Pythia6Decayer()
55 : G4VExtDecayer("G4Pythia6Decayer"), 56 : G4VExtDecayer("G4Pythia6Decayer"),
56 fMessenger(this), 57 fMessenger(this),
57 fVerboseLevel(0), 58 fVerboseLevel(0),
58 fDecayType(fgkDefaultDecayType), 59 fDecayType(fgkDefaultDecayType),
59 fDecayProductsArray(0) 60 fDecayProductsArray(0)
60 { 61 {
61 /// Standard constructor << 62 /// Standard constructor
62 63
63 fDecayProductsArray = new ParticleVector(); 64 fDecayProductsArray = new ParticleVector();
64 << 65
65 ForceDecay(fDecayType); 66 ForceDecay(fDecayType);
66 } 67 }
67 68
68 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 70
70 G4Pythia6Decayer::~G4Pythia6Decayer() << 71 G4Pythia6Decayer::~G4Pythia6Decayer()
71 { 72 {
72 /// Destructor << 73 /// Destructor
73 74
74 delete fDecayProductsArray; 75 delete fDecayProductsArray;
75 } 76 }
76 77
77 // 78 //
78 // private methods 79 // private methods
79 // 80 //
80 81
81 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 83
83 G4ParticleDefinition* G4Pythia6Decayer::GetPar << 84 G4ParticleDefinition* G4Pythia6Decayer::
84 << 85 GetParticleDefinition(const Pythia6Particle* particle, G4bool warn) const
85 { 86 {
86 /// Return G4 particle definition for given << 87 /// Return G4 particle definition for given TParticle
87 88
88 // get particle definition from G4ParticleTa 89 // get particle definition from G4ParticleTable
89 G4int pdgEncoding = particle->fKF; 90 G4int pdgEncoding = particle->fKF;
90 G4ParticleTable* particleTable = G4ParticleT << 91 G4ParticleTable* particleTable
91 G4ParticleDefinition* particleDefinition = 0 << 92 = G4ParticleTable::GetParticleTable();
92 if (pdgEncoding != 0) particleDefinition = p << 93 G4ParticleDefinition* particleDefinition = 0;
93 << 94 if (pdgEncoding != 0)
94 if (particleDefinition == 0 && warn) { << 95 particleDefinition = particleTable->FindParticle(pdgEncoding);
95 G4cerr << "G4Pythia6Decayer: GetParticleDe << 96
96 << "G4ParticleTable::FindParticle() << 97 if ( particleDefinition == 0 && warn) {
97 << " failed." << std::endl; << 98 G4cerr
98 } << 99 << "G4Pythia6Decayer: GetParticleDefinition: " << std::endl
99 << 100 << "G4ParticleTable::FindParticle() for particle with PDG = "
>> 101 << pdgEncoding
>> 102 << " failed." << std::endl;
>> 103 }
>> 104
100 return particleDefinition; 105 return particleDefinition;
101 } 106 }
102 107
103 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 109
105 G4DynamicParticle* G4Pythia6Decayer::CreateDyn << 110 G4DynamicParticle*
106 { << 111 G4Pythia6Decayer::CreateDynamicParticle(const Pythia6Particle* particle) const
107 /// Create G4DynamicParticle. << 112 {
>> 113 /// Create G4DynamicParticle.
108 114
109 // get particle properties 115 // get particle properties
110 const G4ParticleDefinition* particleDefiniti << 116 const G4ParticleDefinition* particleDefinition
111 if (!particleDefinition) return 0; << 117 = GetParticleDefinition(particle);
112 << 118 if ( ! particleDefinition ) return 0;
>> 119
113 G4ThreeVector momentum = GetParticleMomentum 120 G4ThreeVector momentum = GetParticleMomentum(particle);
114 121
115 // create G4DynamicParticle 122 // create G4DynamicParticle
116 G4DynamicParticle* dynamicParticle = new G4D << 123 G4DynamicParticle* dynamicParticle
117 << 124 = new G4DynamicParticle(particleDefinition, momentum);
>> 125
118 return dynamicParticle; 126 return dynamicParticle;
119 } 127 }
120 128
121 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 130
123 G4ThreeVector G4Pythia6Decayer::GetParticlePos << 131 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(
>> 132 const Pythia6Particle* particle) const
124 { 133 {
125 /// Return particle vertex position. << 134 /// Return particle vertex position.
126 135
127 G4ThreeVector position = << 136 G4ThreeVector position
128 G4ThreeVector(particle->fVx * cm, particle << 137 = G4ThreeVector(particle->fVx * cm,
>> 138 particle->fVy * cm,
>> 139 particle->fVz * cm);
129 return position; 140 return position;
130 } << 141 }
131 << 142
132 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 144
134 G4ThreeVector G4Pythia6Decayer::GetParticleMom << 145 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(
>> 146 const Pythia6Particle* particle) const
135 { 147 {
136 /// Return particle momentum. << 148 /// Return particle momentum.
137 149
138 G4ThreeVector momentum = << 150 G4ThreeVector momentum
139 G4ThreeVector(particle->fPx * GeV, particl << 151 = G4ThreeVector(particle->fPx * GeV,
>> 152 particle->fPy * GeV,
>> 153 particle->fPz * GeV);
140 return momentum; 154 return momentum;
141 } 155 }
142 156
143 //....oooOO0OOooo........oooOO0OOooo........oo 157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 158
145 G4int G4Pythia6Decayer::CountProducts(G4int ch 159 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
146 { 160 {
147 /// Count number of decay products << 161 /// Count number of decay products
148 162
149 G4int np = 0; << 163 G4int np = 0;
150 for (G4int i = 1; i <= 5; i++) << 164 for ( G4int i=1; i<=5; i++ )
151 if (std::abs(Pythia6::Instance()->GetKFDP( << 165 if ( std::abs(Pythia6::Instance()->GetKFDP(channel,i) ) == particle )
152 return np; << 166 np++;
>> 167 return np;
153 } 168 }
154 169
155 //....oooOO0OOooo........oooOO0OOooo........oo 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156 171
157 void G4Pythia6Decayer::ForceParticleDecay(G4in << 172 void
>> 173 G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
158 { 174 {
159 /// Force decay of particle into products wi << 175 /// Force decay of particle into products with multiplicity mult
160 176
161 Pythia6* pythia6 = Pythia6::Instance(); << 177 Pythia6* pythia6 = Pythia6::Instance();
162 178
163 G4int kc = pythia6->Pycomp(particle); << 179 G4int kc = pythia6->Pycomp(particle);
164 pythia6->SetMDCY(kc, 1, 1); << 180 pythia6->SetMDCY(kc,1,1);
165 181
166 G4int ifirst = pythia6->GetMDCY(kc, 2); << 182 G4int ifirst = pythia6->GetMDCY(kc,2);
167 G4int ilast = ifirst + pythia6->GetMDCY(kc, << 183 G4int ilast = ifirst + pythia6->GetMDCY(kc,3)-1;
168 184
169 // << 185 //
170 // Loop over decay channels << 186 // Loop over decay channels
171 for (G4int channel = ifirst; channel <= ilas << 187 for (G4int channel= ifirst; channel <= ilast; channel++) {
172 if (CountProducts(channel, product) >= mul << 188 if (CountProducts(channel,product) >= mult) {
173 pythia6->SetMDME(channel, 1, 1); << 189 pythia6->SetMDME(channel,1,1);
174 } << 190 } else {
175 else { << 191 pythia6->SetMDME(channel,1,0);
176 pythia6->SetMDME(channel, 1, 0); << 192 }
177 } << 193 }
178 } <<
179 } 194 }
180 195
181 //....oooOO0OOooo........oooOO0OOooo........oo 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 197
183 void G4Pythia6Decayer::ForceParticleDecay(G4in << 198 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products,
>> 199 G4int* mult, G4int npart)
184 { 200 {
185 /// Force decay of particle into products wi << 201 /// Force decay of particle into products with multiplicity mult
186 202
187 Pythia6* pythia6 = Pythia6::Instance(); << 203 Pythia6* pythia6 = Pythia6::Instance();
188 204
189 G4int kc = pythia6->Pycomp(particle); << 205 G4int kc = pythia6->Pycomp(particle);
190 pythia6->SetMDCY(kc, 1, 1); << 206 pythia6->SetMDCY(kc,1,1);
191 G4int ifirst = pythia6->GetMDCY(kc, 2); << 207 G4int ifirst = pythia6->GetMDCY(kc,2);
192 G4int ilast = ifirst + pythia6->GetMDCY(kc, << 208 G4int ilast = ifirst+pythia6->GetMDCY(kc,3)-1;
193 // << 209 //
194 // Loop over decay channels << 210 // Loop over decay channels
195 for (G4int channel = ifirst; channel <= ilas << 211 for (G4int channel = ifirst; channel <= ilast; channel++) {
196 G4int nprod = 0; << 212 G4int nprod = 0;
197 for (G4int i = 0; i < npart; i++) << 213 for (G4int i = 0; i < npart; i++)
198 nprod += (CountProducts(channel, product << 214 nprod += (CountProducts(channel, products[i]) >= mult[i]);
199 if (nprod) << 215 if (nprod)
200 pythia6->SetMDME(channel, 1, 1); << 216 pythia6->SetMDME(channel,1,1);
201 else { << 217 else {
202 pythia6->SetMDME(channel, 1, 0); << 218 pythia6->SetMDME(channel,1,0);
203 } << 219 }
204 } << 220 }
205 } 221 }
206 222
207 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208 224
209 void G4Pythia6Decayer::ForceHadronicD() 225 void G4Pythia6Decayer::ForceHadronicD()
210 { 226 {
211 /// Force golden D decay modes << 227 /// Force golden D decay modes
212 228
213 const G4int kNHadrons = 4; << 229 const G4int kNHadrons = 4;
214 G4int channel; << 230 G4int channel;
215 G4int hadron[kNHadrons] = {411, 421, 431, 41 << 231 G4int hadron[kNHadrons] = {411, 421, 431, 4112};
216 << 232
217 // for D+ -> K0* (-> K- pi+) pi+ << 233 // for D+ -> K0* (-> K- pi+) pi+
218 G4int iKstar0 = 313; << 234 G4int iKstar0 = 313;
219 G4int iKstarbar0 = -313; << 235 G4int iKstarbar0 = -313;
220 G4int iKPlus = 321; << 236 G4int iKPlus = 321;
221 G4int iKMinus = -321; << 237 G4int iKMinus = -321;
222 G4int iPiPlus = 211; << 238 G4int iPiPlus = 211;
223 G4int iPiMinus = -211; << 239 G4int iPiMinus = -211;
224 << 240
225 G4int products[2] = {iKPlus, iPiMinus}, mult << 241 G4int products[2] = {iKPlus, iPiMinus}, mult[2] = {1, 1};
226 ForceParticleDecay(iKstar0, products, mult, << 242 ForceParticleDecay(iKstar0, products, mult, 2);
227 << 243
228 // for Ds -> Phi pi+ << 244 // for Ds -> Phi pi+
229 G4int iPhi = 333; << 245 G4int iPhi = 333;
230 ForceParticleDecay(iPhi, iKPlus, 2); // Phi << 246 ForceParticleDecay(iPhi,iKPlus,2); // Phi->K+K-
231 << 247
232 G4int decayP1[kNHadrons][3] = { << 248 G4int decayP1[kNHadrons][3] = {
233 {iKMinus, iPiPlus, iPiPlus}, {iKMinus, iPi << 249 {iKMinus, iPiPlus, iPiPlus},
234 G4int decayP2[kNHadrons][3] = { << 250 {iKMinus, iPiPlus, 0 },
235 {iKstarbar0, iPiPlus, 0}, {-1, -1, -1}, {i << 251 {iKPlus , iKstarbar0, 0 },
236 << 252 {-1 , -1 , -1 }
237 Pythia6* pythia6 = Pythia6::Instance(); << 253 };
238 for (G4int ihadron = 0; ihadron < kNHadrons; << 254 G4int decayP2[kNHadrons][3] = {
239 G4int kc = pythia6->Pycomp(hadron[ihadron] << 255 {iKstarbar0, iPiPlus, 0 },
240 pythia6->SetMDCY(kc, 1, 1); << 256 {-1 , -1 , -1 },
241 G4int ifirst = pythia6->GetMDCY(kc, 2); << 257 {iPhi , iPiPlus, 0 },
242 G4int ilast = ifirst + pythia6->GetMDCY(kc << 258 {-1 , -1 , -1 }
243 << 259 };
244 for (channel = ifirst; channel <= ilast; c << 260
245 if ((pythia6->GetKFDP(channel, 1) == dec << 261 Pythia6* pythia6 = Pythia6::Instance();
246 && pythia6->GetKFDP(channel, 2) == << 262 for ( G4int ihadron = 0; ihadron < kNHadrons; ihadron++ ) {
247 && pythia6->GetKFDP(channel, 3) == << 263 G4int kc = pythia6->Pycomp(hadron[ihadron]);
248 && pythia6->GetKFDP(channel, 4) == << 264 pythia6->SetMDCY(kc,1,1);
249 || (pythia6->GetKFDP(channel, 1) == << 265 G4int ifirst = pythia6->GetMDCY(kc,2);
250 && pythia6->GetKFDP(channel, 2) << 266 G4int ilast = ifirst + pythia6->GetMDCY(kc,3)-1;
251 && pythia6->GetKFDP(channel, 3) << 267
252 && pythia6->GetKFDP(channel, 4) << 268 for (channel = ifirst; channel <= ilast; channel++) {
253 { << 269 if ((pythia6->GetKFDP(channel,1) == decayP1[ihadron][0] &&
254 pythia6->SetMDME(channel, 1, 1); << 270 pythia6->GetKFDP(channel,2) == decayP1[ihadron][1] &&
255 } << 271 pythia6->GetKFDP(channel,3) == decayP1[ihadron][2] &&
256 else { << 272 pythia6->GetKFDP(channel,4) == 0) ||
257 pythia6->SetMDME(channel, 1, 0); << 273 (pythia6->GetKFDP(channel,1) == decayP2[ihadron][0] &&
258 } // selected channel ? << 274 pythia6->GetKFDP(channel,2) == decayP2[ihadron][1] &&
259 } // decay channels << 275 pythia6->GetKFDP(channel,3) == decayP2[ihadron][2] &&
260 } // hadrons << 276 pythia6->GetKFDP(channel,4) == 0)) {
>> 277 pythia6->SetMDME(channel,1,1);
>> 278 } else {
>> 279 pythia6->SetMDME(channel,1,0);
>> 280 } // selected channel ?
>> 281 } // decay channels
>> 282 } // hadrons
261 } 283 }
262 284
263 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264 286
265 void G4Pythia6Decayer::ForceOmega() 287 void G4Pythia6Decayer::ForceOmega()
266 { 288 {
267 /// Force Omega -> Lambda K- Decay << 289 /// Force Omega -> Lambda K- Decay
268 290
269 Pythia6* pythia6 = Pythia6::Instance(); << 291 Pythia6* pythia6 = Pythia6::Instance();
270 292
271 G4int iLambda0 = 3122; << 293 G4int iLambda0 = 3122;
272 G4int iKMinus = -321; << 294 G4int iKMinus = -321;
273 295
274 G4int kc = pythia6->Pycomp(3334); << 296 G4int kc = pythia6->Pycomp(3334);
275 pythia6->SetMDCY(kc, 1, 1); << 297 pythia6->SetMDCY(kc,1,1);
276 G4int ifirst = pythia6->GetMDCY(kc, 2); << 298 G4int ifirst = pythia6->GetMDCY(kc,2);
277 G4int ilast = ifirst + pythia6->GetMDCY(kc, << 299 G4int ilast = ifirst + pythia6->GetMDCY(kc,3)-1;
278 << 300
279 for (G4int channel = ifirst; channel <= ilas << 301 for (G4int channel = ifirst; channel <= ilast; channel++) {
280 if (pythia6->GetKFDP(channel, 1) == iLambd << 302 if (pythia6->GetKFDP(channel,1) == iLambda0 &&
281 && pythia6->GetKFDP(channel, 3) == 0) << 303 pythia6->GetKFDP(channel,2) == iKMinus &&
282 pythia6->SetMDME(channel, 1, 1); << 304 pythia6->GetKFDP(channel,3) == 0)
283 else << 305 pythia6->SetMDME(channel,1,1);
284 pythia6->SetMDME(channel, 1, 0); << 306 else
285 // selected channel ? << 307 pythia6->SetMDME(channel,1,0);
286 } // decay channels << 308 // selected channel ?
>> 309 } // decay channels
287 } 310 }
288 311
289 //....oooOO0OOooo........oooOO0OOooo........oo 312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
290 313
291 void G4Pythia6Decayer::ForceDecay(EDecayType d 314 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
292 { 315 {
293 /// Force a particle decay mode << 316 /// Force a particle decay mode
294 <<
295 Pythia6::Instance()->SetMSTJ(21, 2); <<
296 317
297 if (fDecayType == kNoDecayHeavy) return; << 318 Pythia6::Instance()->SetMSTJ(21,2);
298 319
299 // << 320 if ( fDecayType == kNoDecayHeavy ) return;
300 // select mode <<
301 G4int products[3]; <<
302 G4int mult[3]; <<
303 321
304 switch (decayType) { << 322 //
305 case kHardMuons: << 323 // select mode
306 products[0] = 13; << 324 G4int products[3];
307 products[1] = 443; << 325 G4int mult[3];
>> 326
>> 327 switch ( decayType ) {
>> 328
>> 329 case kHardMuons:
>> 330 products[0] = 13;
>> 331 products[1] = 443;
308 products[2] = 100443; 332 products[2] = 100443;
309 mult[0] = 1; 333 mult[0] = 1;
310 mult[1] = 1; 334 mult[1] = 1;
311 mult[2] = 1; 335 mult[2] = 1;
312 ForceParticleDecay(511, products, mult, << 336 ForceParticleDecay( 511, products, mult, 3);
313 ForceParticleDecay(521, products, mult, << 337 ForceParticleDecay( 521, products, mult, 3);
314 ForceParticleDecay(531, products, mult, << 338 ForceParticleDecay( 531, products, mult, 3);
315 ForceParticleDecay(5122, products, mult, << 339 ForceParticleDecay( 5122, products, mult, 3);
316 ForceParticleDecay(5132, products, mult, << 340 ForceParticleDecay( 5132, products, mult, 3);
317 ForceParticleDecay(5232, products, mult, << 341 ForceParticleDecay( 5232, products, mult, 3);
318 ForceParticleDecay(5332, products, mult, << 342 ForceParticleDecay( 5332, products, mult, 3);
319 ForceParticleDecay(100443, 443, 1); // << 343 ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
320 ForceParticleDecay(443, 13, 2); // J/Ps << 344 ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
321 << 345
322 ForceParticleDecay(411, 13, 1); // D+/- << 346 ForceParticleDecay( 411,13,1); // D+/-
323 ForceParticleDecay(421, 13, 1); // D0 << 347 ForceParticleDecay( 421,13,1); // D0
324 ForceParticleDecay(431, 13, 1); // D_s << 348 ForceParticleDecay( 431,13,1); // D_s
325 ForceParticleDecay(4122, 13, 1); // Lam << 349 ForceParticleDecay( 4122,13,1); // Lambda_c
326 ForceParticleDecay(4132, 13, 1); // Xsi << 350 ForceParticleDecay( 4132,13,1); // Xsi_c
327 ForceParticleDecay(4232, 13, 1); // Sig << 351 ForceParticleDecay( 4232,13,1); // Sigma_c
328 ForceParticleDecay(4332, 13, 1); // Ome << 352 ForceParticleDecay( 4332,13,1); // Omega_c
329 break; << 353 break;
330 << 354
331 case kSemiMuonic: << 355 case kSemiMuonic:
332 ForceParticleDecay(411, 13, 1); // D+/- << 356 ForceParticleDecay( 411,13,1); // D+/-
333 ForceParticleDecay(421, 13, 1); // D0 << 357 ForceParticleDecay( 421,13,1); // D0
334 ForceParticleDecay(431, 13, 1); // D_s << 358 ForceParticleDecay( 431,13,1); // D_s
335 ForceParticleDecay(4122, 13, 1); // Lam << 359 ForceParticleDecay( 4122,13,1); // Lambda_c
336 ForceParticleDecay(4132, 13, 1); // Xsi << 360 ForceParticleDecay( 4132,13,1); // Xsi_c
337 ForceParticleDecay(4232, 13, 1); // Sig << 361 ForceParticleDecay( 4232,13,1); // Sigma_c
338 ForceParticleDecay(4332, 13, 1); // Ome << 362 ForceParticleDecay( 4332,13,1); // Omega_c
339 ForceParticleDecay(511, 13, 1); // B0 << 363 ForceParticleDecay( 511,13,1); // B0
340 ForceParticleDecay(521, 13, 1); // B+/- << 364 ForceParticleDecay( 521,13,1); // B+/-
341 ForceParticleDecay(531, 13, 1); // B_s << 365 ForceParticleDecay( 531,13,1); // B_s
342 ForceParticleDecay(5122, 13, 1); // Lam << 366 ForceParticleDecay( 5122,13,1); // Lambda_b
343 ForceParticleDecay(5132, 13, 1); // Xsi << 367 ForceParticleDecay( 5132,13,1); // Xsi_b
344 ForceParticleDecay(5232, 13, 1); // Sig << 368 ForceParticleDecay( 5232,13,1); // Sigma_b
345 ForceParticleDecay(5332, 13, 1); // Ome << 369 ForceParticleDecay( 5332,13,1); // Omega_b
346 break; << 370 break;
347 << 371
348 case kDiMuon: << 372 case kDiMuon:
349 ForceParticleDecay(113, 13, 2); // rho << 373 ForceParticleDecay( 113,13,2); // rho
350 ForceParticleDecay(221, 13, 2); // eta << 374 ForceParticleDecay( 221,13,2); // eta
351 ForceParticleDecay(223, 13, 2); // omeg << 375 ForceParticleDecay( 223,13,2); // omega
352 ForceParticleDecay(333, 13, 2); // phi << 376 ForceParticleDecay( 333,13,2); // phi
353 ForceParticleDecay(443, 13, 2); // J/Ps << 377 ForceParticleDecay( 443,13,2); // J/Psi
354 ForceParticleDecay(100443, 13, 2); // P << 378 ForceParticleDecay(100443,13,2);// Psi'
355 ForceParticleDecay(553, 13, 2); // Upsi << 379 ForceParticleDecay( 553,13,2); // Upsilon
356 ForceParticleDecay(100553, 13, 2); // U << 380 ForceParticleDecay(100553,13,2);// Upsilon'
357 ForceParticleDecay(200553, 13, 2); // U << 381 ForceParticleDecay(200553,13,2);// Upsilon''
358 break; << 382 break;
359 << 383
360 case kSemiElectronic: << 384 case kSemiElectronic:
361 ForceParticleDecay(411, 11, 1); // D+/- << 385 ForceParticleDecay( 411,11,1); // D+/-
362 ForceParticleDecay(421, 11, 1); // D0 << 386 ForceParticleDecay( 421,11,1); // D0
363 ForceParticleDecay(431, 11, 1); // D_s << 387 ForceParticleDecay( 431,11,1); // D_s
364 ForceParticleDecay(4122, 11, 1); // Lam << 388 ForceParticleDecay( 4122,11,1); // Lambda_c
365 ForceParticleDecay(4132, 11, 1); // Xsi << 389 ForceParticleDecay( 4132,11,1); // Xsi_c
366 ForceParticleDecay(4232, 11, 1); // Sig << 390 ForceParticleDecay( 4232,11,1); // Sigma_c
367 ForceParticleDecay(4332, 11, 1); // Ome << 391 ForceParticleDecay( 4332,11,1); // Omega_c
368 ForceParticleDecay(511, 11, 1); // B0 << 392 ForceParticleDecay( 511,11,1); // B0
369 ForceParticleDecay(521, 11, 1); // B+/- << 393 ForceParticleDecay( 521,11,1); // B+/-
370 ForceParticleDecay(531, 11, 1); // B_s << 394 ForceParticleDecay( 531,11,1); // B_s
371 ForceParticleDecay(5122, 11, 1); // Lam << 395 ForceParticleDecay( 5122,11,1); // Lambda_b
372 ForceParticleDecay(5132, 11, 1); // Xsi << 396 ForceParticleDecay( 5132,11,1); // Xsi_b
373 ForceParticleDecay(5232, 11, 1); // Sig << 397 ForceParticleDecay( 5232,11,1); // Sigma_b
374 ForceParticleDecay(5332, 11, 1); // Ome << 398 ForceParticleDecay( 5332,11,1); // Omega_b
375 break; << 399 break;
376 << 400
377 case kDiElectron: << 401 case kDiElectron:
378 ForceParticleDecay(113, 11, 2); // rho << 402 ForceParticleDecay( 113,11,2); // rho
379 ForceParticleDecay(333, 11, 2); // phi << 403 ForceParticleDecay( 333,11,2); // phi
380 ForceParticleDecay(221, 11, 2); // eta << 404 ForceParticleDecay( 221,11,2); // eta
381 ForceParticleDecay(223, 11, 2); // omeg << 405 ForceParticleDecay( 223,11,2); // omega
382 ForceParticleDecay(443, 11, 2); // J/Ps << 406 ForceParticleDecay( 443,11,2); // J/Psi
383 ForceParticleDecay(100443, 11, 2); // P << 407 ForceParticleDecay(100443,11,2);// Psi'
384 ForceParticleDecay(553, 11, 2); // Upsi << 408 ForceParticleDecay( 553,11,2); // Upsilon
385 ForceParticleDecay(100553, 11, 2); // U << 409 ForceParticleDecay(100553,11,2);// Upsilon'
386 ForceParticleDecay(200553, 11, 2); // U << 410 ForceParticleDecay(200553,11,2);// Upsilon''
387 break; 411 break;
388 412
389 case kBJpsiDiMuon: << 413 case kBJpsiDiMuon:
390 414
391 products[0] = 443; << 415 products[0] = 443;
392 products[1] = 100443; 416 products[1] = 100443;
393 mult[0] = 1; 417 mult[0] = 1;
394 mult[1] = 1; 418 mult[1] = 1;
395 419
396 ForceParticleDecay(511, products, mult, << 420 ForceParticleDecay( 511, products, mult, 2); // B0 -> J/Psi (Psi') X
397 ForceParticleDecay(521, products, mult, << 421 ForceParticleDecay( 521, products, mult, 2); // B+/- -> J/Psi (Psi') X
398 ForceParticleDecay(531, products, mult, << 422 ForceParticleDecay( 531, products, mult, 2); // B_s -> J/Psi (Psi') X
399 ForceParticleDecay(5122, products, mult, << 423 ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi')X
400 ForceParticleDecay(100443, 443, 1); // << 424 ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
401 ForceParticleDecay(443, 13, 2); // J/Ps << 425 ForceParticleDecay( 443,13,2); // J/Psi -> mu+ mu-
402 break; 426 break;
403 427
404 case kBPsiPrimeDiMuon: << 428 case kBPsiPrimeDiMuon:
405 ForceParticleDecay(511, 100443, 1); // << 429 ForceParticleDecay( 511,100443,1); // B0
406 ForceParticleDecay(521, 100443, 1); // << 430 ForceParticleDecay( 521,100443,1); // B+/-
407 ForceParticleDecay(531, 100443, 1); // << 431 ForceParticleDecay( 531,100443,1); // B_s
408 ForceParticleDecay(5122, 100443, 1); // << 432 ForceParticleDecay( 5122,100443,1); // Lambda_b
409 ForceParticleDecay(100443, 13, 2); // P << 433 ForceParticleDecay(100443,13,2); // Psi'
410 break; 434 break;
411 435
412 case kBJpsiDiElectron: << 436 case kBJpsiDiElectron:
413 ForceParticleDecay(511, 443, 1); // B0 << 437 ForceParticleDecay( 511,443,1); // B0
414 ForceParticleDecay(521, 443, 1); // B+/ << 438 ForceParticleDecay( 521,443,1); // B+/-
415 ForceParticleDecay(531, 443, 1); // B_s << 439 ForceParticleDecay( 531,443,1); // B_s
416 ForceParticleDecay(5122, 443, 1); // La << 440 ForceParticleDecay( 5122,443,1); // Lambda_b
417 ForceParticleDecay(443, 11, 2); // J/Ps << 441 ForceParticleDecay( 443,11,2); // J/Psi
418 break; 442 break;
419 443
420 case kBJpsi: << 444 case kBJpsi:
421 ForceParticleDecay(511, 443, 1); // B0 << 445 ForceParticleDecay( 511,443,1); // B0
422 ForceParticleDecay(521, 443, 1); // B+/ << 446 ForceParticleDecay( 521,443,1); // B+/-
423 ForceParticleDecay(531, 443, 1); // B_s << 447 ForceParticleDecay( 531,443,1); // B_s
424 ForceParticleDecay(5122, 443, 1); // La << 448 ForceParticleDecay( 5122,443,1); // Lambda_b
425 break; 449 break;
426 450
427 case kBPsiPrimeDiElectron: << 451 case kBPsiPrimeDiElectron:
428 ForceParticleDecay(511, 100443, 1); // << 452 ForceParticleDecay( 511,100443,1); // B0
429 ForceParticleDecay(521, 100443, 1); // << 453 ForceParticleDecay( 521,100443,1); // B+/-
430 ForceParticleDecay(531, 100443, 1); // << 454 ForceParticleDecay( 531,100443,1); // B_s
431 ForceParticleDecay(5122, 100443, 1); // << 455 ForceParticleDecay( 5122,100443,1); // Lambda_b
432 ForceParticleDecay(100443, 11, 2); // P << 456 ForceParticleDecay(100443,11,2); // Psi'
433 break; 457 break;
434 458
435 case kPiToMu: << 459 case kPiToMu:
436 ForceParticleDecay(211, 13, 1); // pi-> << 460 ForceParticleDecay(211,13,1); // pi->mu
437 break; 461 break;
438 462
439 case kKaToMu: << 463 case kKaToMu:
440 ForceParticleDecay(321, 13, 1); // K->m << 464 ForceParticleDecay(321,13,1); // K->mu
441 break; 465 break;
442 466
443 case kWToMuon: << 467 case kWToMuon:
444 ForceParticleDecay(24, 13, 1); // W -> << 468 ForceParticleDecay( 24, 13,1); // W -> mu
445 break; 469 break;
446 470
447 case kWToCharm: << 471 case kWToCharm:
448 ForceParticleDecay(24, 4, 1); // W -> c << 472 ForceParticleDecay( 24, 4,1); // W -> c
449 break; 473 break;
450 474
451 case kWToCharmToMuon: << 475 case kWToCharmToMuon:
452 ForceParticleDecay(24, 4, 1); // W -> c << 476 ForceParticleDecay( 24, 4,1); // W -> c
453 ForceParticleDecay(411, 13, 1); // D+/- << 477 ForceParticleDecay( 411,13,1); // D+/- -> mu
454 ForceParticleDecay(421, 13, 1); // D0 << 478 ForceParticleDecay( 421,13,1); // D0 -> mu
455 ForceParticleDecay(431, 13, 1); // D_s << 479 ForceParticleDecay( 431,13,1); // D_s -> mu
456 ForceParticleDecay(4122, 13, 1); // Lam << 480 ForceParticleDecay( 4122,13,1); // Lambda_c
457 ForceParticleDecay(4132, 13, 1); // Xsi << 481 ForceParticleDecay( 4132,13,1); // Xsi_c
458 ForceParticleDecay(4232, 13, 1); // Sig << 482 ForceParticleDecay( 4232,13,1); // Sigma_c
459 ForceParticleDecay(4332, 13, 1); // Ome << 483 ForceParticleDecay( 4332,13,1); // Omega_c
460 break; 484 break;
461 485
462 case kZDiMuon: << 486 case kZDiMuon:
463 ForceParticleDecay(23, 13, 2); // Z -> << 487 ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
464 break; 488 break;
465 489
466 case kHadronicD: << 490 case kHadronicD:
467 ForceHadronicD(); 491 ForceHadronicD();
468 break; 492 break;
469 493
470 case kPhiKK: << 494 case kPhiKK:
471 ForceParticleDecay(333, 321, 2); // Phi << 495 ForceParticleDecay(333,321,2); // Phi->K+K-
472 break; 496 break;
473 497
474 case kOmega: << 498 case kOmega:
475 ForceOmega(); 499 ForceOmega();
476 500
477 case kAll: << 501 case kAll:
478 break; 502 break;
479 503
480 case kNoDecay: << 504 case kNoDecay:
481 Pythia6::Instance()->SetMSTJ(21, 0); << 505 Pythia6::Instance()->SetMSTJ(21,0);
482 break; 506 break;
483 507
484 case kNoDecayHeavy: << 508 case kNoDecayHeavy: break;
485 break; <<
486 509
487 case kMaxDecay: << 510 case kMaxDecay: break;
488 break; << 511 }
489 } <<
490 } 512 }
491 513
492 //....oooOO0OOooo........oooOO0OOooo........oo 514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
493 515
494 void G4Pythia6Decayer::Decay(G4int pdg, const 516 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
495 { 517 {
496 /// Decay a particle of type IDPART (PDG cod << 518 /// Decay a particle of type IDPART (PDG code) and momentum P.
497 519
498 Pythia6::Instance()->Py1ent(0, pdg, p.e(), p << 520 Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
499 } 521 }
500 522
501 //....oooOO0OOooo........oooOO0OOooo........oo 523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502 524
503 G4int G4Pythia6Decayer::ImportParticles(Partic 525 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
504 { 526 {
505 /// Get the decay products into the passed P << 527 /// Get the decay products into the passed PARTICLES vector
506 528
507 return Pythia6::Instance()->ImportParticles( << 529 return Pythia6::Instance()->ImportParticles(particles,"All");
508 } 530 }
509 531
510 // 532 //
511 // public methods 533 // public methods
512 // 534 //
513 535
514 //....oooOO0OOooo........oooOO0OOooo........oo 536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
515 537
516 G4DecayProducts* G4Pythia6Decayer::ImportDecay 538 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
517 { 539 {
518 /// Import decay products << 540 /// Import decay products
519 541
520 // get particle momentum 542 // get particle momentum
521 G4ThreeVector momentum = track.GetMomentum() << 543 G4ThreeVector momentum = track.GetMomentum();
522 G4double etot = track.GetDynamicParticle()-> << 544 G4double etot = track.GetDynamicParticle()->GetTotalEnergy();;
523 ; << 545 CLHEP::HepLorentzVector p;
524 CLHEP::HepLorentzVector p; <<
525 p[0] = momentum.x() / GeV; 546 p[0] = momentum.x() / GeV;
526 p[1] = momentum.y() / GeV; 547 p[1] = momentum.y() / GeV;
527 p[2] = momentum.z() / GeV; 548 p[2] = momentum.z() / GeV;
528 p[3] = etot / GeV; << 549 p[3] = etot / GeV;
529 << 550
530 // get particle PDG 551 // get particle PDG
531 // ask G4Pythia6Decayer to get PDG encoding << 552 // ask G4Pythia6Decayer to get PDG encoding
532 // (in order to get PDG from extended TDatab 553 // (in order to get PDG from extended TDatabasePDG
533 // in case the standard PDG code is not defi 554 // in case the standard PDG code is not defined)
534 G4ParticleDefinition* particleDef = track.Ge 555 G4ParticleDefinition* particleDef = track.GetDefinition();
535 G4int pdgEncoding = particleDef->GetPDGEncod 556 G4int pdgEncoding = particleDef->GetPDGEncoding();
536 557
537 // let Pythia6Decayer decay the particle 558 // let Pythia6Decayer decay the particle
538 // and import the decay products 559 // and import the decay products
539 Decay(pdgEncoding, p); 560 Decay(pdgEncoding, p);
540 G4int nofParticles = ImportParticles(fDecayP 561 G4int nofParticles = ImportParticles(fDecayProductsArray);
541 << 562
542 if (fVerboseLevel > 0) { << 563 if ( fVerboseLevel > 0 ) {
543 G4cout << "nofParticles: " << nofParticles << 564 G4cout << "nofParticles: " << nofParticles << G4endl;
544 } << 565 }
545 << 566
546 // convert decay products Pythia6Particle ty << 567 // convert decay products Pythia6Particle type
547 // to G4DecayProducts << 568 // to G4DecayProducts
548 G4DecayProducts* decayProducts = new G4Decay << 569 G4DecayProducts* decayProducts
>> 570 = new G4DecayProducts(*(track.GetDynamicParticle()));
549 571
550 G4int counter = 0; 572 G4int counter = 0;
551 for (G4int i = 0; i < nofParticles; i++) { << 573 for (G4int i=0; i<nofParticles; i++) {
>> 574
552 // get particle from ParticleVector 575 // get particle from ParticleVector
553 Pythia6Particle* particle = (*fDecayProduc 576 Pythia6Particle* particle = (*fDecayProductsArray)[i];
554 << 577
555 G4int status = particle->fKS; 578 G4int status = particle->fKS;
556 G4int pdg = particle->fKF; 579 G4int pdg = particle->fKF;
557 if (status > 0 && status < 11 && std::abs( << 580 if ( status>0 && status<11 &&
558 && std::abs(pdg) != 16) << 581 std::abs(pdg)!=12 && std::abs(pdg)!=14 && std::abs(pdg)!=16 ) {
559 { <<
560 // pass to tracking final particles only 582 // pass to tracking final particles only;
561 // skip neutrinos 583 // skip neutrinos
562 584
563 if (fVerboseLevel > 0) { << 585 if ( fVerboseLevel > 0 ) {
564 G4cout << " " << i << "th particle PD 586 G4cout << " " << i << "th particle PDG: " << pdg << " ";
565 } << 587 }
566 << 588
567 // create G4DynamicParticle << 589 // create G4DynamicParticle
568 G4DynamicParticle* dynamicParticle = Cre << 590 G4DynamicParticle* dynamicParticle
>> 591 = CreateDynamicParticle(particle);
569 592
570 if (dynamicParticle) { 593 if (dynamicParticle) {
571 if (fVerboseLevel > 0) { << 594
572 G4cout << " G4 particle name: " << << 595 if ( fVerboseLevel > 0 ) {
>> 596 G4cout << " G4 particle name: "
>> 597 << dynamicParticle->GetDefinition()->GetParticleName()
573 << G4endl; 598 << G4endl;
574 } << 599 }
575 600
576 // add dynamicParticle to decayProduct 601 // add dynamicParticle to decayProducts
577 decayProducts->PushProducts(dynamicPar 602 decayProducts->PushProducts(dynamicParticle);
578 << 603
579 counter++; 604 counter++;
580 } 605 }
581 } << 606 }
582 } << 607 }
583 if (fVerboseLevel > 0) { << 608 if ( fVerboseLevel > 0 ) {
584 G4cout << "nofParticles for tracking: " << << 609 G4cout << "nofParticles for tracking: " << counter << G4endl;
585 } << 610 }
586 << 611
587 return decayProducts; 612 return decayProducts;
588 } 613 }
589 << 614
590 //....oooOO0OOooo........oooOO0OOooo........oo 615 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591 616
592 void G4Pythia6Decayer::ForceDecayType(EDecayTy 617 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
593 { << 618 {
594 /// Force a given decay type << 619 /// Force a given decay type
595 620
596 // Do nothing if the decay type is not diffe 621 // Do nothing if the decay type is not different from current one
597 if (decayType == fDecayType) return; << 622 if ( decayType == fDecayType ) return;
598 << 623
599 fDecayType = decayType; << 624 fDecayType = decayType;
600 ForceDecay(fDecayType); 625 ForceDecay(fDecayType);
601 } 626 }
602 627