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

Last change on this file was 1347, checked in by garnier, 14 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.