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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
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-03-ref-09 $
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.