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

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

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