source: trunk/source/processes/hadronic/models/high_energy/src/G4HEVector.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 33.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4HEVector.cc,v 1.20 2006/11/22 23:52:32 dennis Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29//
30//
31
32#include "globals.hh"
33#include "G4ios.hh"
34
35//
36// G4 Gheisha friend class G4GHEVector
37// J.L. Chuma, TRIUMF, 22-Feb-1996
38// last modified: H. Fesefeldt 02-July--1998
39// Fesefeldt, bug fixed in Defs1, 14 August 2000
40
41#include "G4HEVector.hh"
42#include "G4ParticleDefinition.hh"
43
44G4HEVector::G4HEVector(const G4HadProjectile * aParticle)
45 {
46 G4ThreeVector aMom = 1./GeV*aParticle->Get4Momentum().vect();
47 px = aMom.x();
48 py = aMom.y();
49 pz = aMom.z();
50 energy = aParticle->GetTotalEnergy()/GeV;
51 kineticEnergy = aParticle->GetKineticEnergy()/GeV;
52 mass = aParticle->GetDefinition()->GetPDGMass()/GeV;
53 charge = aParticle->GetDefinition()->GetPDGCharge()/eplus;
54 timeOfFlight = 0.0;
55 side = 0;
56 flag = false;
57 code = aParticle->GetDefinition()->GetPDGEncoding();
58 baryon = aParticle->GetDefinition()->GetBaryonNumber();
59 particleName = getParticleName(code, baryon);
60 particleType = aParticle->GetDefinition()->GetParticleType();
61 strangeness = aParticle->GetDefinition()->GetQuarkContent(3);
62 }
63
64G4double G4HEVector::Amax(G4double a, G4double b)
65 {
66 G4double c = a;
67 if(b > a) c = b;
68 return c;
69 }
70
71G4String G4HEVector::getParticleName(G4int aCode, G4int aBaryon)
72 {
73 G4String name;
74 if(aCode == 211) name = "PionPlus";
75 else if(aCode == 111) name = "PionZero";
76 else if(aCode == -211) name = "PionMinus";
77 else if(aCode == 321) name = "KaonPlus";
78 else if(aCode == 311) name = "KaonZero";
79 else if(aCode == -311) name = "AntiKaonZero";
80 else if(aCode == -321) name = "KaonMinus";
81 else if(aCode == 310) name = "KaonZeroShort";
82 else if(aCode == 130) name = "KaonZeroLong";
83 else if(aCode == 2212) name = "Proton";
84 else if(aCode == -2212) name = "AntiProton";
85 else if(aCode == 2112) name = "Neutron";
86 else if(aCode == -2112) name = "AntiNeutron";
87 else if(aCode == 3122) name = "Lambda";
88 else if(aCode == -3122) name = "AntiLambda";
89 else if(aCode == 3222) name = "SigmaPlus";
90 else if(aCode == 3212) name = "SigmaZero";
91 else if(aCode == 3112) name = "SigmaMinus";
92 else if(aCode == -3222) name = "AntiSigmaPlus";
93 else if(aCode == -3212) name = "AntiSigmaZero";
94 else if(aCode == -3112) name = "AntiSigmaMinus";
95 else if(aCode == 3322) name = "XiZero";
96 else if(aCode == 3312) name = "XiMinus";
97 else if(aCode == -3322) name = "AntiXiZero";
98 else if(aCode == -3312) name = "AntiXiMinus";
99 else if(aCode == 3334) name = "OmegaMinus";
100 else if(aCode == -3334) name = "AntiOmegaMinus";
101 else if(aCode == 0)
102 {
103 if(aBaryon==2) name = "Deuteron";
104 else if(aBaryon==3) name = "Triton";
105 else if(aBaryon==4) name = "Alpha";
106 }
107 else if(aCode == 22) name = "Gamma";
108 else
109 {
110 G4cout << "particle " << aCode << " " <<aBaryon<< " not known in this generator!!" << G4endl;
111 }
112 return name;
113 }
114
115
116void
117G4HEVector::setMomentum(const G4ParticleMomentum mom )
118 {
119 px = mom.x();
120 py = mom.y();
121 pz = mom.z();
122 return;
123 }
124
125void
126G4HEVector::setMomentum(const G4ParticleMomentum * mom )
127 {
128 px = mom->x();
129 py = mom->y();
130 pz = mom->z();
131 return;
132 }
133
134void
135G4HEVector::setMomentumAndUpdate( const G4ParticleMomentum mom )
136 {
137 px = mom.x();
138 py = mom.y();
139 pz = mom.z();
140 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
141 kineticEnergy = Amax(0.,energy - mass);
142 return;
143 }
144
145void
146G4HEVector::setMomentumAndUpdate( const G4ParticleMomentum * mom )
147 {
148 px = mom->x();
149 py = mom->y();
150 pz = mom->z();
151 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
152 kineticEnergy = Amax(0.,energy - mass);
153 return;
154 }
155
156const G4ParticleMomentum
157G4HEVector::getMomentum() const
158 {
159 G4ParticleMomentum mom;
160 mom.setX(px);
161 mom.setY(py);
162 mom.setZ(pz);
163 return mom;
164 }
165
166G4double
167G4HEVector::getTotalMomentum()
168 {
169 return std::sqrt(px*px + py*py + pz*pz);
170 }
171
172void
173G4HEVector::setMomentum( G4double x, G4double y, G4double z)
174 {
175 px = x;
176 py = y;
177 pz = z;
178 return;
179 }
180
181void
182G4HEVector::setMomentumAndUpdate( G4double x, G4double y, G4double z )
183 {
184 px = x;
185 py = y;
186 pz = z;
187 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
188 kineticEnergy = Amax(0.,energy-mass);
189 return;
190 }
191
192void
193G4HEVector::setMomentum( G4double x, G4double y )
194 {
195 px = x;
196 py = y;
197 return;
198 }
199
200void
201G4HEVector::setMomentumAndUpdate( G4double x, G4double y )
202 {
203 px = x;
204 py = y;
205 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
206 kineticEnergy = Amax(0.,energy-mass);
207 return;
208 }
209
210void
211G4HEVector::setMomentum( G4double z )
212 {
213 pz = z;
214 return;
215 }
216
217void
218G4HEVector::setMomentumAndUpdate( G4double z )
219 {
220 pz = z;
221 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
222 kineticEnergy = Amax(0.,energy-mass);
223 return;
224 }
225
226void
227G4HEVector::setEnergy( G4double e )
228 {
229 energy = e;
230 return;
231 }
232
233void
234G4HEVector::setEnergyAndUpdate( G4double e )
235 {
236 if (e <= mass)
237 {
238 energy = mass;
239 kineticEnergy = 0.;
240 px = 0.;
241 py = 0.;
242 pz = 0.;
243 }
244 else
245 {
246 energy = e;
247 kineticEnergy = energy - mass;
248 G4double momold = std::sqrt(px*px + py*py + pz*pz);
249 G4double momnew = std::sqrt(energy*energy - mass*mass);
250 if (momold == 0.)
251 {
252 G4double cost = 1.0- 2.0*G4UniformRand();
253 G4double sint = std::sqrt(1. - cost*cost);
254 G4double phi = twopi* G4UniformRand();
255 px = momnew * sint * std::cos(phi);
256 py = momnew * sint * std::sin(phi);
257 pz = momnew * cost;
258 }
259 else
260 {
261 momnew /= momold;
262 px *= momnew;
263 py *= momnew;
264 pz *= momnew;
265 }
266 }
267 return;
268 }
269
270void
271G4HEVector::setKineticEnergy( G4double ekin )
272 {
273 kineticEnergy = ekin;
274 return;
275 }
276
277void
278G4HEVector::setKineticEnergyAndUpdate(G4double ekin)
279 {
280 if (ekin <= 0.)
281 {
282 energy = mass;
283 kineticEnergy = 0.;
284 px = 0.;
285 py = 0.;
286 pz = 0.;
287 }
288 else
289 {
290 energy = ekin + mass;
291 kineticEnergy = ekin;
292 G4double momold = std::sqrt(px*px + py*py + pz*pz);
293 G4double momnew = std::sqrt(energy*energy - mass*mass);
294 if (momold == 0.)
295 {
296 G4double cost = 1.0-2.0*G4UniformRand();
297 G4double sint = std::sqrt(1. - cost*cost);
298 G4double phi = twopi* G4UniformRand();
299 px = momnew * sint * std::cos(phi);
300 py = momnew * sint * std::sin(phi);
301 pz = momnew * cost;
302 }
303 else
304 {
305 momnew /= momold;
306 px *= momnew;
307 py *= momnew;
308 pz *= momnew;
309 }
310 }
311 return;
312 }
313
314G4double
315G4HEVector::getEnergy()
316 {
317 return energy;
318 }
319
320G4double
321G4HEVector::getKineticEnergy()
322 {
323 return kineticEnergy;
324 }
325
326void
327G4HEVector::setMass( G4double m )
328 {
329 mass = m;
330 return;
331 }
332
333void
334G4HEVector::setMassAndUpdate( G4double m )
335 {
336 kineticEnergy = Amax(0., energy - mass);
337 mass = m;
338 energy = kineticEnergy + mass;
339 G4double momnew = std::sqrt(Amax(0., energy*energy - mass*mass));
340 if ( momnew == 0.0)
341 {
342 px = 0.;
343 py = 0.;
344 pz = 0.;
345 }
346 else
347 {
348 G4double momold = std::sqrt(px*px + py*py + pz*pz);
349 if (momold == 0.)
350 {
351 G4double cost = 1.-2.*G4UniformRand();
352 G4double sint = std::sqrt(1.-cost*cost);
353 G4double phi = twopi*G4UniformRand();
354 px = momnew*sint*std::cos(phi);
355 py = momnew*sint*std::sin(phi);
356 pz = momnew*cost;
357 }
358 else
359 {
360 momnew /= momold;
361 px *= momnew ;
362 py *= momnew ;
363 pz *= momnew ;
364 }
365 }
366 return;
367 }
368
369G4double
370G4HEVector::getMass()
371 {
372 return mass;
373 }
374
375void
376G4HEVector::setCharge( G4double c )
377 {
378 charge = c;
379 return;
380 }
381
382G4double
383G4HEVector::getCharge()
384 {
385 return charge;
386 }
387
388void
389G4HEVector::setTOF( G4double t )
390 {
391 timeOfFlight = t;
392 return;
393 }
394
395G4double
396G4HEVector::getTOF()
397 {
398 return timeOfFlight;
399 }
400
401void
402G4HEVector::setSide( G4int s )
403 {
404 side = s;
405 return;
406 }
407
408G4int
409G4HEVector::getSide()
410 {
411 return side;
412 }
413
414
415void
416G4HEVector::setFlag( G4bool f )
417 {
418 flag = f;
419 return;
420 }
421
422G4bool
423G4HEVector::getFlag()
424 {
425 return flag;
426 }
427
428void
429G4HEVector::setCode( G4int c )
430 {
431 code = c;
432 return;
433 }
434
435G4int
436G4HEVector::getCode()
437 {
438 return code;
439 }
440
441G4String
442G4HEVector::getName()
443 {
444 return particleName;
445 }
446
447G4String
448G4HEVector::getType()
449 {
450 return particleType;
451 }
452
453G4int
454G4HEVector::getBaryonNumber()
455 {
456 return baryon;
457 }
458
459G4int
460G4HEVector::getStrangenessNumber()
461 {
462 return strangeness;
463 }
464G4int
465G4HEVector::getQuarkContent(G4int flavor)
466 {
467 if(flavor > 0 && flavor < 8)
468 {
469 G4int check;
470 check = FillQuarkContent();
471 if(check != code)
472 {
473 return 0;
474 }
475 else
476 {
477 return theQuarkContent[flavor-1];
478 }
479 }
480 else
481 {
482 return 0;
483 }
484 }
485
486G4int
487G4HEVector::getAntiQuarkContent(G4int flavor)
488 {
489 if(flavor > 0 && flavor < 8)
490 {
491 G4int check;
492 check = FillQuarkContent();
493 if(check != code)
494 {
495 return 0;
496 }
497 else
498 {
499 return theAntiQuarkContent[flavor-1];
500 }
501 }
502 else
503 {
504 return 0;
505 }
506 }
507
508void
509G4HEVector::setZero()
510 {
511 px = 0.0;
512 py = 0.0;
513 pz = 0.0;
514 energy = 0.0;
515 kineticEnergy = 0.0;
516 mass = 0.0;
517 charge = 0.0;
518 timeOfFlight = 0.0;
519 side = 0;
520 flag = false;
521 code = 0;
522 particleName = "";
523 particleType = "";
524 baryon = 0;
525 strangeness = 0;
526 }
527
528void
529G4HEVector::Add( const G4HEVector & p1, const G4HEVector & p2 )
530 {
531 px = p1.px + p2.px;
532 py = p1.py + p2.py;
533 pz = p1.pz + p2.pz;
534 energy = p1.energy + p2.energy;
535 G4double b = energy*energy - px*px - py*py - pz*pz;
536 if( b < 0 )
537 mass = -1. * std::sqrt( -b );
538 else
539 mass = std::sqrt( b );
540 kineticEnergy = Amax(0.,energy - mass);
541 charge = p1.charge + p2.charge;
542 code = 0;
543 particleName = "";
544 particleType = "";
545 baryon = 0;
546 strangeness = 0;
547 }
548
549void
550G4HEVector::Sub( const G4HEVector & p1, const G4HEVector & p2 )
551 {
552 px = p1.px - p2.px;
553 py = p1.py - p2.py;
554 pz = p1.pz - p2.pz;
555 energy = p1.energy - p2.energy;
556 G4double b = energy*energy - px*px - py*py - pz*pz;
557 if( b < 0 )
558 mass = -1. * std::sqrt( -b );
559 else
560 mass = std::sqrt( b );
561 kineticEnergy = Amax(0.,energy - mass);
562 charge = p1.charge - p2.charge;
563 code = 0;
564 particleName = "";
565 particleType = "";
566 baryon = 0;
567 strangeness = 0;
568 }
569
570void
571G4HEVector::Lor( const G4HEVector & p1, const G4HEVector & p2 )
572 {
573 G4double a;
574 a = ( Dot(p1,p2)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
575 px = p1.px + a*p2.px;
576 py = p1.py + a*p2.py;
577 pz = p1.pz + a*p2.pz;
578 energy = std::sqrt( sqr(p1.mass) + px*px + py*py + pz*pz);
579 mass = p1.mass;
580 kineticEnergy = Amax(0.,energy - mass);
581 timeOfFlight = p1.timeOfFlight;
582 side = p1.side;
583 flag = p1.flag;
584 code = p1.code;
585 particleName = p1.particleName;
586 particleType = p1.particleType;
587 baryon = p1.baryon;
588 strangeness = p1.strangeness;
589 }
590
591G4double
592G4HEVector::CosAng( const G4HEVector & p )
593 {
594 G4double a = std::sqrt( (px*px + py*py + pz*pz)*(p.px*p.px + p.py*p.py + p.pz*p.pz) );
595 if( a != 0.0 )
596 {
597 a = (px*p.px + py*p.py + pz*p.pz)/a;
598 if( std::fabs(a) > 1.0 )
599 {
600 if(a<0.0) a=-1.0;
601 else a=1.0;
602 }
603 }
604 return a;
605 }
606
607G4double
608G4HEVector::Ang(const G4HEVector & p )
609 {
610 G4double a = std::sqrt( (px*px + py*py + pz*pz)*(p.px*p.px + p.py*p.py + p.pz*p.pz) );
611 if( a != 0.0 )
612 {
613 a = (px*p.px + py*p.py + pz*p.pz)/a;
614 if( std::fabs(a) > 1.0 )
615 {
616 if(a<0.0) a=-1.0;
617 else a=1.0;
618 }
619 }
620 return std::acos(a);
621 }
622
623G4double
624G4HEVector::Dot4( const G4HEVector & p1, const G4HEVector & p2)
625 {
626 return ( p1.energy*p2.energy - p1.px*p2.px - p1.py*p2.py - p1.pz*p2.pz );
627 }
628
629G4double
630G4HEVector::Impu( const G4HEVector & p1, const G4HEVector & p2)
631 {
632 return ( - sqr( p1.energy - p2.energy)
633 + sqr( p1.px - p2.px)
634 + sqr( p1.py - p2.py)
635 + sqr( p1.pz - p2.pz) );
636 }
637
638void
639G4HEVector::Add3( const G4HEVector & p1, const G4HEVector & p2)
640 {
641 px = p1.px + p2.px;
642 py = p1.py + p2.py;
643 pz = p1.pz + p2.pz;
644 return;
645 }
646
647void
648G4HEVector::Sub3( const G4HEVector & p1, const G4HEVector & p2)
649 {
650 px = p1.px - p2.px;
651 py = p1.py - p2.py;
652 pz = p1.pz - p2.pz;
653 return;
654 }
655
656void
657G4HEVector::Cross( const G4HEVector & p1, const G4HEVector & p2)
658 {
659 G4double tpx = p1.py * p2.pz - p1.pz * p2.py;
660 G4double tpy = p1.pz * p2.px - p1.px * p2.pz;
661 G4double tpz = p1.px * p2.py - p1.py * p2.px;
662 px=tpx;
663 py=tpy;
664 pz=tpz;
665 return;
666 }
667
668G4double
669G4HEVector::Dot( const G4HEVector & p1, const G4HEVector & p2)
670 {
671 return ( p1.px * p2.px + p1.py * p2.py + p1.pz * p2.pz );
672 }
673
674void
675G4HEVector::Smul( const G4HEVector & p, G4double h)
676 {
677 px = h * p.px;
678 py = h * p.py;
679 pz = h * p.pz;
680 return;
681 }
682
683void
684G4HEVector::SmulAndUpdate( const G4HEVector & p, G4double h)
685 {
686 px = h * p.px;
687 py = h * p.py;
688 pz = h * p.pz;
689 mass = p.mass;
690 energy = std::sqrt(px*px + py*py + pz*pz + mass*mass);
691 kineticEnergy = energy - mass;
692 charge = p.charge;
693 timeOfFlight = p.timeOfFlight;
694 side = p.side;
695 flag = p.flag;
696 code = p.code;
697 particleName = p.particleName;
698 particleType = p.particleType;
699 baryon = p.baryon;
700 strangeness = p.strangeness;
701 return;
702 }
703
704void
705G4HEVector::Norz( const G4HEVector & p )
706 {
707 G4double a = p.px*p.px + p.py*p.py + p.pz*p.pz;
708 if (a > 0.0) a = 1./std::sqrt(a);
709 px = a * p.px;
710 py = a * p.py;
711 pz = a * p.pz;
712 mass = p.mass;
713 energy = std::sqrt(px*px + py*py + pz*pz + mass*mass);
714 kineticEnergy = energy - mass;
715 charge = p.charge;
716 timeOfFlight = p.timeOfFlight;
717 side = p.side;
718 flag = p.flag;
719 code = p.code;
720 particleName = p.particleName;
721 particleType = p.particleType;
722 baryon = p.baryon;
723 strangeness = p.strangeness;
724 return;
725 }
726
727G4double
728G4HEVector::Length()
729 {
730 return std::sqrt(px*px + py*py + pz*pz);
731 }
732
733void
734G4HEVector::Exch( G4HEVector & p1)
735 {
736 G4HEVector mx = *this;
737 *this = p1;
738 p1 = mx;
739 return;
740 }
741
742void
743G4HEVector::Defs1( const G4HEVector & p1, const G4HEVector & p2)
744 {
745 G4double pt2 = p2.px*p2.px + p2.py*p2.py;
746 if (pt2 > 0.0)
747 {
748 G4double ph, qx, qy, qz;
749 G4double a = std::sqrt(p2.px*p2.px + p2.py*p2.py + p2.pz*p2.pz);
750 G4double cost = p2.pz/a;
751 G4double sint = 0.5 * (std::sqrt(std::fabs((1.-cost)*(1.+cost))) + std::sqrt(pt2)/a);
752 if(p2.py < 0.) ph = 1.5*pi;
753 else ph = halfpi;
754 if( p2.px != 0.0)
755 ph = std::atan2(p2.py,p2.px);
756 qx = cost*std::cos(ph)*p1.px - std::sin(ph)*p1.py
757 + sint*std::cos(ph)*p1.pz;
758 qy = cost*std::sin(ph)*p1.px + std::cos(ph)*p1.py
759 + sint*std::sin(ph)*p1.pz;
760 qz = - sint *p1.px
761 + cost *p1.pz;
762 px = qx;
763 py = qy;
764 pz = qz;
765 }
766 else
767 {
768 px = p1.px;
769 py = p1.py;
770 pz = p1.pz;
771 }
772 }
773
774void
775G4HEVector::Defs( const G4HEVector & p1, const G4HEVector & p2,
776 G4HEVector & my, G4HEVector & mz )
777 {
778 my = p1;
779 mz = p2;
780 px = my.py*mz.pz - my.pz*mz.py;
781 py = my.pz*mz.px - my.px*mz.pz;
782 pz = my.px*mz.py - my.py*mz.px;
783 my.px = mz.py*pz - mz.pz*py;
784 my.py = mz.pz*px - mz.px*pz;
785 my.pz = mz.px*py - mz.py*px;
786 G4double pp;
787 pp = std::sqrt(px*px + py*py + pz*pz);
788 if (pp > 0.)
789 {
790 pp = 1./pp;
791 px = px*pp ;
792 py = py*pp ;
793 pz = pz*pp ;
794 }
795 pp = std::sqrt(my.px*my.px + my.py*my.py + my.pz*my.pz);
796 if (pp > 0.)
797 {
798 pp = 1./pp;
799 my.px = my.px*pp ;
800 my.py = my.py*pp ;
801 my.pz = my.pz*pp ;
802 }
803 pp = std::sqrt(mz.px*mz.px + mz.py*mz.py + mz.pz*mz.pz);
804 if (pp > 0.)
805 {
806 pp = 1./pp;
807 mz.px = mz.px*pp ;
808 mz.py = mz.py*pp ;
809 mz.pz = mz.pz*pp ;
810 }
811 return;
812 }
813
814void
815G4HEVector::Trac( const G4HEVector & p1, const G4HEVector & mx,
816 const G4HEVector & my, const G4HEVector & mz)
817 {
818 G4double qx, qy, qz;
819 qx = mx.px*p1.px + mx.py*p1.py + mx.pz*p1.pz;
820 qy = my.px*p1.px + my.py*p1.py + my.pz*p1.pz;
821 qz = mz.px*p1.px + mz.py*p1.py + mz.pz*p1.pz;
822 px = qx ;
823 py = qy ;
824 pz = qz ;
825 return;
826 }
827
828void
829G4HEVector::setDefinition(G4String name)
830 {
831 if(name == "PionPlus")
832 {
833 mass = 0.1395700;
834 charge = 1.;
835 code = 211;
836 particleType = "Meson";
837 particleName = name;
838 baryon = 0;
839 strangeness = 0;
840 }
841 else if(name == "PionZero")
842 {
843 mass = 0.1349764;
844 charge = 0.;
845 code = 111;
846 particleType = "Meson";
847 particleName = name;
848 baryon = 0;
849 strangeness = 0;
850 }
851 else if(name == "PionMinus")
852 {
853 mass = 0.1395700;
854 charge = -1.;
855 code = -211;
856 particleType = "Meson";
857 particleName = name;
858 baryon = 0;
859 strangeness = 0;
860 }
861 else if(name == "KaonPlus")
862 {
863 mass = 0.493677;
864 charge = 1.;
865 code = 321;
866 particleType = "Meson";
867 particleName = name;
868 baryon = 0;
869 strangeness = 1;
870 }
871 else if(name == "KaonZero")
872 {
873 mass = 0.497672;
874 charge = 0.;
875 code = 311;
876 particleType = "Meson";
877 particleName = name;
878 baryon = 0;
879 strangeness = 1;
880 }
881 else if(name == "AntiKaonZero")
882 {
883 mass = 0.497672;
884 charge = 0.;
885 code = -311;
886 particleType = "Meson";
887 particleName = name;
888 baryon = 0;
889 strangeness =-1;
890 }
891 else if(name == "KaonMinus")
892 {
893 mass = 0.493677;
894 charge = -1.;
895 code = -321;
896 particleType = "Meson";
897 particleName = name;
898 baryon = 0;
899 strangeness = -1;
900 }
901 else if(name == "KaonZeroShort")
902 {
903 mass = 0.497672;
904 charge = 0.;
905 code = 310;
906 particleType = "Meson";
907 particleName = name;
908 baryon = 0;
909 strangeness = 0;
910 }
911 else if(name == "KaonZeroLong")
912 {
913 mass = 0.497672;
914 charge = 0.;
915 code = 130;
916 particleType = "Meson";
917 particleName = name;
918 baryon = 0;
919 strangeness = 0;
920 }
921 else if(name == "Proton")
922 {
923 mass = 0.9382723;
924 charge = 1.;
925 code = 2212;
926 particleType = "Baryon";
927 particleName = name;
928 baryon = 1;
929 strangeness = 0;
930 }
931 else if(name == "AntiProton")
932 {
933 mass = 0.9382723;
934 charge = -1.;
935 code = -2212;
936 particleType = "AntiBaryon";
937 particleName = name;
938 baryon = -1;
939 strangeness = 0;
940 }
941 else if(name == "Neutron")
942 {
943 mass = 0.93956563;
944 charge = 0.;
945 code = 2112;
946 particleType = "Baryon";
947 particleName = name;
948 baryon = 1;
949 strangeness = 0;
950 }
951 else if(name == "AntiNeutron")
952 {
953 mass = 0.93956563;
954 charge = 0.;
955 code = -2112;
956 particleType = "AntiBaryon";
957 particleName = name;
958 baryon = -1;
959 strangeness = 0;
960 }
961 else if(name == "Lambda")
962 {
963 mass = 1.115684;
964 charge = 0.;
965 code = 3122;
966 particleType = "Baryon";
967 particleName = name;
968 baryon = 1;
969 strangeness = -1;
970 }
971 else if(name == "AntiLambda")
972 {
973 mass = 1.115684;
974 charge = 0.;
975 code = -3122;
976 particleType = "AntiBaryon";
977 particleName = name;
978 baryon = -1;
979 strangeness = 1;
980 }
981 else if(name == "SigmaPlus")
982 {
983 mass = 1.18937;
984 charge = 1.;
985 code = 3222;
986 particleType = "Baryon";
987 particleName = name;
988 baryon = 1;
989 strangeness = -1;
990 }
991 else if(name == "SigmaZero")
992 {
993 mass = 1.19255;
994 charge = 0.;
995 code = 3212;
996 particleType = "Baryon";
997 particleName = name;
998 baryon = 1;
999 strangeness = -1;
1000 }
1001 else if(name == "SigmaMinus")
1002 {
1003 mass = 1.19744;
1004 charge = -1.;
1005 code = 3112;
1006 particleType = "Baryon";
1007 particleName = name;
1008 baryon = 1;
1009 strangeness = -1;
1010 }
1011 else if(name == "AntiSigmaPlus")
1012 {
1013 mass = 1.18937;
1014 charge = -1.;
1015 code = -3222;
1016 particleType = "AntiBaryon";
1017 particleName = name;
1018 baryon = -1;
1019 strangeness = 1;
1020 }
1021 else if(name == "AntiSigmaZero")
1022 {
1023 mass = 1.19255;
1024 charge = 0.;
1025 code = -3212;
1026 particleType = "AntiBaryon";
1027 particleName = name;
1028 baryon = -1;
1029 strangeness = 1;
1030 }
1031 else if(name == "AntiSigmaMinus")
1032 {
1033 mass = 1.19744;
1034 charge = 1.;
1035 code = -3112;
1036 particleType = "AntiBaryon";
1037 particleName = name;
1038 baryon = -1;
1039 strangeness = 1;
1040 }
1041 else if(name == "XiZero")
1042 {
1043 mass = 1.3149;
1044 charge = 0.;
1045 code = 3322;
1046 particleType = "Baryon";
1047 particleName = name;
1048 baryon = 1;
1049 strangeness = -2;
1050 }
1051 else if(name == "XiMinus")
1052 {
1053 mass = 1.32132;
1054 charge = -1.;
1055 code = 3312;
1056 particleType = "Baryon";
1057 particleName = name;
1058 baryon = 1;
1059 strangeness = -2;
1060 }
1061 else if(name == "AntiXiZero")
1062 {
1063 mass = 1.3149;
1064 charge = 0.;
1065 code = -3322;
1066 particleType = "AntiBaryon";
1067 particleName = name;
1068 baryon = -1;
1069 strangeness = 2;
1070 }
1071 else if(name == "AntiXiMinus")
1072 {
1073 mass = 1.32132;
1074 charge = 1.;
1075 code = -3312;
1076 particleType = "AntiBaryon";
1077 particleName = name;
1078 baryon = -1;
1079 strangeness = 2;
1080 }
1081 else if(name == "OmegaMinus")
1082 {
1083 mass = 1.67245;
1084 charge = -1.;
1085 code = 3334;
1086 particleType = "Baryon";
1087 particleName = name;
1088 baryon = 1;
1089 strangeness = -3;
1090 }
1091 else if(name == "AntiOmegaMinus")
1092 {
1093 mass = 1.67245;
1094 charge = 1.;
1095 code = -3334;
1096 particleType = "AntiBaryon";
1097 particleName = name;
1098 baryon = -1;
1099 strangeness = 3;
1100 }
1101 else if(name == "Deuteron")
1102 {
1103 mass = 1.875613;
1104 charge = 1.;
1105 code = 0;
1106 particleType = "Nucleus";
1107 particleName = name;
1108 baryon = 2;
1109 strangeness = 0;
1110 }
1111 else if(name == "Triton")
1112 {
1113 mass = 2.80925;
1114 charge = 1.;
1115 code = 0;
1116 particleType = "Nucleus";
1117 particleName = name;
1118 baryon = 3;
1119 strangeness = 0;
1120 }
1121 else if(name == "Alpha")
1122 {
1123 mass = 3.727417;
1124 charge = 2.;
1125 code = 0;
1126 particleType = "Nucleus";
1127 particleName = name;
1128 baryon = 4;
1129 strangeness = 0;
1130 }
1131 else if(name == "Gamma")
1132 {
1133 mass = 0.;
1134 charge = 0.;
1135 code = 22;
1136 particleType = "Boson";
1137 particleName = name;
1138 baryon = 0;
1139 strangeness = 0;
1140 }
1141 else
1142 {
1143 G4cout << "particle " << name << " not known in this generator!!" << G4endl;
1144 return;
1145 }
1146 px = 0.;
1147 py = 0.;
1148 pz = 0.;
1149 kineticEnergy = 0.;
1150 energy = mass;
1151 timeOfFlight = 0.;
1152 side = 0;
1153 flag = false;
1154 return;
1155 }
1156
1157G4int G4HEVector::FillQuarkContent()
1158 // calculate quark and anti-quark contents
1159 // return value is PDG encoding for this particle.
1160 // It means error if the return value is differnt from
1161 // this->thePDGEncoding.
1162{
1163 G4int tempPDGcode = code;
1164 G4double eplus = 1.;
1165
1166 for (G4int flavor=0; flavor<NumberOfQuarkFlavor; flavor++){
1167 theQuarkContent[flavor] =0;
1168 theAntiQuarkContent[flavor] =0;
1169 }
1170
1171 G4int temp = std::abs(tempPDGcode);
1172 G4int multiplet = temp/10000;
1173 temp -= G4int(multiplet*10000);
1174 G4int quark1 = temp/1000;
1175 temp -= G4int(quark1*1000);
1176 G4int quark2 = temp/100;
1177 temp -= G4int(quark2*100);
1178 G4int quark3 = temp/10;
1179 temp -= G4int(quark3*10);
1180 G4int spin= (temp-1);
1181
1182 if (particleType =="quark") {
1183 if (tempPDGcode>0){
1184 if (tempPDGcode<=NumberOfQuarkFlavor){
1185 theQuarkContent[tempPDGcode-1] =1;
1186 } else {
1187 // --- thePDGEncoding is wrong
1188 tempPDGcode = 0;
1189 }
1190 } else {
1191 G4int temp = -1*tempPDGcode;
1192 if (temp<=NumberOfQuarkFlavor){
1193 theAntiQuarkContent[temp-1] =1;
1194 } else {
1195 // --- thePDGEncoding is wrong
1196 tempPDGcode = 0;
1197 }
1198 }
1199
1200 } else if (particleType == "Meson") {
1201 // -- exceptions --
1202 if (tempPDGcode == 310) spin = 0; //K0s
1203 if (tempPDGcode == 130) { //K0l
1204 spin = 0;
1205 quark2 = 3;
1206 quark3 = 1;
1207 }
1208
1209 if (quark1 !=0)
1210 {
1211 tempPDGcode = 0;
1212 }
1213 if ((quark2==0)||(quark3==0)){
1214 tempPDGcode = 0;
1215 }
1216 if (quark2<quark3) {
1217 tempPDGcode = 0;
1218 }
1219 // check quark flavor
1220 if (quark2>=NumberOfQuarkFlavor){
1221 tempPDGcode = 0;
1222 }
1223 // check heavier quark type
1224 if (quark2 & 1) {
1225 // down type qurak
1226 if (tempPDGcode >0) {
1227 theQuarkContent[quark3-1] =1;
1228 theAntiQuarkContent[quark2-1] =1;
1229 } else {
1230 theQuarkContent[quark2-1] =1;
1231 theAntiQuarkContent[quark3-1] =1;
1232 }
1233 } else {
1234 // up type quark
1235 if (tempPDGcode >0) {
1236 theQuarkContent[quark2-1] =1;
1237 theAntiQuarkContent[quark3-1] =1;
1238 } else {
1239 theQuarkContent[quark3-1] =1;
1240 theAntiQuarkContent[quark2-1] =1;
1241 }
1242 }
1243 // check charge
1244 G4double totalCharge = 0.0;
1245 for (G4int flavor= 0; flavor<NumberOfQuarkFlavor-1; flavor+=2){
1246 totalCharge += (-1./3.)*eplus*theQuarkContent[flavor];
1247 totalCharge += 1./3.*eplus*theAntiQuarkContent[flavor];
1248 totalCharge += 2./3.*eplus*theQuarkContent[flavor+1];
1249 totalCharge += (-2./3.)*eplus*theAntiQuarkContent[flavor+1];
1250 }
1251 if (std::abs(totalCharge-charge)>0.1*eplus) {
1252 tempPDGcode = 0;
1253 }
1254 } else if (particleType == "baryon"){
1255 // check Meson or not
1256 if ((quark1==0)||(quark2==0)||(quark3==0)){
1257 tempPDGcode = 0;
1258 }
1259 //exceptions
1260 if (std::abs(tempPDGcode) == 3122) {
1261 // Lambda
1262 quark2=2; quark3 = 1; spin = 1;
1263 } else if (std::abs(tempPDGcode) == 4122) {
1264 // Lambda_c
1265 quark2=2; quark3 = 1; spin = 1;
1266 }
1267 // check quark flavor
1268 if ((quark1<quark2)||(quark2<quark3)||(quark1<quark3)) {
1269 tempPDGcode = 0;
1270 }
1271 if (quark1>=NumberOfQuarkFlavor) {
1272 tempPDGcode = 0;
1273 }
1274 if (tempPDGcode >0) {
1275 theQuarkContent[quark1-1] ++;
1276 theQuarkContent[quark2-1] ++;
1277 theQuarkContent[quark3-1] ++;
1278 } else {
1279 theAntiQuarkContent[quark1-1] ++;
1280 theAntiQuarkContent[quark2-1] ++;
1281 theAntiQuarkContent[quark3-1] ++;
1282 }
1283 // check charge
1284 G4double totalCharge = 0.0;
1285 for (G4int flavor= 0; flavor<NumberOfQuarkFlavor-1; flavor+=2){
1286 totalCharge += (-1./3.)*eplus*theQuarkContent[flavor];
1287 totalCharge += 1./3.*eplus*theAntiQuarkContent[flavor];
1288 totalCharge += 2./3.*eplus*theQuarkContent[flavor+1];
1289 totalCharge += (-2./3.)*eplus*theAntiQuarkContent[flavor+1];
1290 }
1291 if (std::abs(totalCharge-charge)>0.1*eplus) {
1292 tempPDGcode = 0;
1293 }
1294 } else {
1295 }
1296 return tempPDGcode;
1297}
1298
1299void
1300G4HEVector::Print( G4int L)
1301 {
1302 G4cout << "HEV: "
1303 << L << " " << px << " " << py << " " << pz << " "
1304 << energy << " " << mass << " " << charge << " "
1305 << timeOfFlight << " " << side << " " << flag << " "
1306 << code << " " << baryon << " " << particleName << G4endl;
1307 /*
1308 printf("HEV: %3d %6f.2 %6f.2 %6f.2 %6f.2 %6f.2 %3f.0 %4f.1 %3d
1309 %3d %6d %6d %s \n", L, px, py, pz, energy, mass, charge,
1310 timeOfFlight, side, flag, code, baryon, particleName);
1311 */
1312 return;
1313 }
1314
Note: See TracBrowser for help on using the repository browser.