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