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