G4HEVector.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //
00030 
00031 //
00032 // G4 Gheisha friend class G4GHEVector
00033 // J.L. Chuma, TRIUMF, 22-Feb-1996
00034 // last modified: H. Fesefeldt 02-July--1998
00035 // Fesefeldt, bug fixed in Defs1, 14 August 2000
00036 
00037 #include "G4HEVector.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4ParticleDefinition.hh"
00043 
00044 G4HEVector::G4HEVector(const G4HadProjectile * aParticle)
00045   {
00046      G4ThreeVector aMom = 1./GeV*aParticle->Get4Momentum().vect();
00047      px               = aMom.x();
00048      py               = aMom.y();
00049      pz               = aMom.z();
00050      energy           = aParticle->GetTotalEnergy()/GeV;
00051      kineticEnergy    = aParticle->GetKineticEnergy()/GeV;
00052      mass             = aParticle->GetDefinition()->GetPDGMass()/GeV;
00053      charge           = aParticle->GetDefinition()->GetPDGCharge()/eplus;
00054      timeOfFlight     = 0.0;
00055      side             = 0;
00056      flag             = false;
00057      code             = aParticle->GetDefinition()->GetPDGEncoding();
00058      baryon           = aParticle->GetDefinition()->GetBaryonNumber();
00059      particleName     = getParticleName(code, baryon);
00060      particleType     = aParticle->GetDefinition()->GetParticleType();
00061      strangeness       = aParticle->GetDefinition()->GetQuarkContent(3);
00062   }
00063   
00064 G4double G4HEVector::Amax(G4double a, G4double b)
00065   {
00066      G4double c = a;
00067      if(b > a) c = b;
00068      return c;
00069   } 
00070 
00071 G4String G4HEVector::getParticleName(G4int aCode, G4int aBaryon)
00072    {
00073         G4String name;
00074         if(aCode == 211) name = "PionPlus";
00075         else if(aCode == 111) name = "PionZero";
00076         else if(aCode == -211) name = "PionMinus";
00077         else if(aCode == 321) name = "KaonPlus";
00078         else if(aCode == 311) name = "KaonZero";
00079         else if(aCode == -311) name = "AntiKaonZero";
00080         else if(aCode == -321) name = "KaonMinus";
00081         else if(aCode == 310) name = "KaonZeroShort";
00082         else if(aCode == 130) name = "KaonZeroLong";
00083         else if(aCode == 2212) name = "Proton";
00084         else if(aCode == -2212) name = "AntiProton";
00085         else if(aCode == 2112) name = "Neutron";
00086         else if(aCode == -2112) name = "AntiNeutron";
00087         else if(aCode == 3122) name = "Lambda";
00088         else if(aCode == -3122) name = "AntiLambda";
00089         else if(aCode == 3222) name = "SigmaPlus";
00090         else if(aCode == 3212) name = "SigmaZero";
00091         else if(aCode == 3112) name = "SigmaMinus";
00092         else if(aCode == -3222) name = "AntiSigmaPlus";
00093         else if(aCode == -3212) name = "AntiSigmaZero";
00094         else if(aCode == -3112) name = "AntiSigmaMinus";
00095         else if(aCode == 3322) name = "XiZero";
00096         else if(aCode == 3312) name = "XiMinus";
00097         else if(aCode == -3322) name = "AntiXiZero";
00098         else if(aCode == -3312) name = "AntiXiMinus";
00099         else if(aCode == 3334) name = "OmegaMinus";
00100         else if(aCode == -3334) name = "AntiOmegaMinus";
00101         else if(aCode == 0)
00102         { 
00103           if(aBaryon==2) name = "Deuteron";
00104           else if(aBaryon==3) name = "Triton";
00105           else if(aBaryon==4) name = "Alpha";
00106         }
00107         else if(aCode == 22) name = "Gamma";
00108         else
00109           {
00110                G4cout << "particle " << aCode << "  "  <<aBaryon<< " not known in this generator!!" << G4endl;
00111           }
00112         return name;
00113    } 
00114 
00115 
00116 void 
00117 G4HEVector::setMomentum(const G4ParticleMomentum mom ) 
00118    {
00119       px  = mom.x();
00120       py  = mom.y();
00121       pz  = mom.z(); 
00122       return; 
00123    }
00124 
00125 void 
00126 G4HEVector::setMomentum(const G4ParticleMomentum * mom ) 
00127    {
00128       px  = mom->x();
00129       py  = mom->y();
00130       pz  = mom->z(); 
00131       return; 
00132    }
00133 
00134 void 
00135 G4HEVector::setMomentumAndUpdate( const G4ParticleMomentum mom )
00136    {
00137      px = mom.x();
00138      py = mom.y();
00139      pz = mom.z();
00140      energy        = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00141      kineticEnergy = Amax(0.,energy - mass);
00142      return;
00143    }
00144 
00145 void 
00146 G4HEVector::setMomentumAndUpdate( const G4ParticleMomentum * mom )
00147    {
00148      px = mom->x();
00149      py = mom->y();
00150      pz = mom->z();
00151      energy        = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00152      kineticEnergy = Amax(0.,energy - mass);
00153      return;
00154    }
00155 
00156 const G4ParticleMomentum 
00157 G4HEVector::getMomentum() const 
00158    { 
00159      G4ParticleMomentum mom;
00160      mom.setX(px);
00161      mom.setY(py);
00162      mom.setZ(pz);
00163      return mom; 
00164    }
00165 
00166 G4double G4HEVector::getTotalMomentum() const
00167 {
00168   return std::sqrt(px*px + py*py + pz*pz);
00169 }
00170 
00171 void
00172 G4HEVector::setMomentum(G4double x, G4double y, G4double z)
00173 { 
00174   px = x;
00175   py = y;
00176   pz = z;
00177   return;
00178 } 
00179 
00180 void 
00181 G4HEVector::setMomentumAndUpdate(G4double x, G4double y, G4double z)
00182 {
00183   px = x;
00184   py = y;
00185   pz = z;
00186   energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00187   kineticEnergy = Amax(0.,energy-mass);
00188   return;
00189 }
00190 
00191 void 
00192 G4HEVector::setMomentum(G4double x, G4double y)
00193 {
00194   px = x;
00195   py = y;
00196   return;
00197 }
00198 
00199 void 
00200 G4HEVector::setMomentumAndUpdate( G4double x, G4double y )
00201    {
00202      px = x;
00203      py = y;
00204      energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00205      kineticEnergy = Amax(0.,energy-mass);
00206      return;
00207    }
00208 
00209 void 
00210 G4HEVector::setMomentum( G4double z )
00211    {
00212      pz = z;
00213      return;
00214    }
00215 
00216 void 
00217 G4HEVector::setMomentumAndUpdate( G4double z )
00218    {
00219      pz = z;
00220      energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00221      kineticEnergy = Amax(0.,energy-mass);
00222      return;
00223    }
00224 
00225 void 
00226 G4HEVector::setEnergy( G4double e ) 
00227    { 
00228      energy = e; 
00229      return; 
00230    }
00231 
00232 void 
00233 G4HEVector::setEnergyAndUpdate( G4double e )
00234    {
00235      if (e <= mass)
00236        { 
00237          energy        = mass;
00238          kineticEnergy = 0.;
00239          px            = 0.;
00240          py            = 0.;
00241          pz            = 0.;
00242        }
00243      else
00244        {
00245          energy = e;
00246          kineticEnergy   = energy - mass;
00247          G4double momold = std::sqrt(px*px + py*py + pz*pz);
00248          G4double momnew = std::sqrt(energy*energy - mass*mass);
00249          if (momold == 0.)
00250            {
00251              G4double cost = 1.0- 2.0*G4UniformRand();
00252              G4double sint = std::sqrt(1. - cost*cost);
00253              G4double phi  = twopi* G4UniformRand();
00254              px            = momnew * sint * std::cos(phi);
00255              py            = momnew * sint * std::sin(phi);
00256              pz            = momnew * cost;
00257            }
00258          else
00259            {
00260              momnew /= momold;
00261              px     *= momnew;
00262              py     *= momnew;
00263              pz     *= momnew;
00264            }
00265        }    
00266      return;
00267    }
00268 
00269 void 
00270 G4HEVector::setKineticEnergy( G4double ekin ) 
00271    { 
00272      kineticEnergy = ekin;
00273      return; 
00274    }
00275  
00276 void 
00277 G4HEVector::setKineticEnergyAndUpdate(G4double ekin) 
00278    {
00279      if (ekin <= 0.)
00280        { 
00281          energy        = mass;
00282          kineticEnergy = 0.;
00283          px            = 0.;
00284          py            = 0.;
00285          pz            = 0.;
00286        }
00287      else
00288        {
00289          energy = ekin + mass;
00290          kineticEnergy   = ekin;
00291          G4double momold = std::sqrt(px*px + py*py + pz*pz);
00292          G4double momnew = std::sqrt(energy*energy - mass*mass);
00293          if (momold == 0.)
00294            {
00295              G4double cost = 1.0-2.0*G4UniformRand();
00296              G4double sint = std::sqrt(1. - cost*cost);
00297              G4double phi  = twopi* G4UniformRand();
00298              px            = momnew * sint * std::cos(phi);
00299              py            = momnew * sint * std::sin(phi);
00300              pz            = momnew * cost;
00301            }
00302          else
00303            {
00304              momnew /= momold;
00305              px     *= momnew;
00306              py     *= momnew;
00307              pz     *= momnew;
00308            }
00309        }    
00310      return;
00311    }
00312 
00313 G4double G4HEVector::getEnergy() const 
00314 {
00315   return energy;
00316 }
00317 
00318 G4double G4HEVector::getKineticEnergy() const 
00319 {
00320   return kineticEnergy;
00321 }
00322 
00323 
00324 void G4HEVector::setMass(G4double amass) 
00325 { 
00326   mass = amass; 
00327   return; 
00328 }
00329 
00330 
00331 void G4HEVector::setMassAndUpdate(G4double amass)
00332 {
00333   kineticEnergy = Amax(0., energy - mass);
00334   mass = amass;
00335   energy = kineticEnergy + mass;
00336   G4double momnew = std::sqrt(Amax(0., energy*energy - mass*mass));
00337   if (momnew == 0.0) {
00338     px = 0.;
00339     py = 0.;
00340     pz = 0.;
00341   } else {
00342     G4double momold = std::sqrt(px*px + py*py + pz*pz);
00343     if (momold == 0.) { 
00344       G4double cost = 1.-2.*G4UniformRand();
00345       G4double sint = std::sqrt(1.-cost*cost);
00346       G4double phi = twopi*G4UniformRand();
00347       px = momnew*sint*std::cos(phi);
00348       py = momnew*sint*std::sin(phi);
00349       pz = momnew*cost;
00350     } else {
00351       momnew /= momold;
00352       px *= momnew ;
00353       py *= momnew ;
00354       pz *= momnew ;
00355     }
00356   }     
00357   return;
00358 }
00359 
00360 
00361 G4double G4HEVector::getMass() const 
00362 { 
00363   return mass; 
00364 }
00365 
00366 void 
00367 G4HEVector::setCharge( G4double c ) 
00368    { 
00369      charge = c; 
00370      return; 
00371    }
00372 
00373 G4double G4HEVector::getCharge() const 
00374 {
00375   return charge; 
00376 }
00377 
00378 void 
00379 G4HEVector::setTOF( G4double t ) 
00380    { 
00381      timeOfFlight = t; 
00382      return; 
00383    }
00384 
00385 G4double 
00386 G4HEVector::getTOF() 
00387    { 
00388      return timeOfFlight; 
00389    }
00390 
00391 
00392 void G4HEVector::setSide(G4int aside) 
00393 { 
00394   side = aside; 
00395   return; 
00396 }
00397 
00398 
00399 G4int 
00400 G4HEVector::getSide() 
00401    { 
00402      return side; 
00403    }
00404 
00405 
00406 void 
00407 G4HEVector::setFlag( G4bool f ) 
00408    { 
00409      flag = f; 
00410      return; 
00411    }
00412 
00413 G4bool 
00414 G4HEVector::getFlag() 
00415    { 
00416      return flag; 
00417    }
00418 
00419 void 
00420 G4HEVector::setCode( G4int c ) 
00421    { 
00422      code = c; 
00423      return; 
00424    }
00425 
00426 G4int G4HEVector::getCode() const 
00427 { 
00428   return code; 
00429 } 
00430 
00431 G4String G4HEVector::getName() const
00432 {
00433   return particleName;
00434 }
00435 
00436 G4String G4HEVector::getType() const
00437 {
00438   return particleType;
00439 }
00440  
00441 G4int G4HEVector::getBaryonNumber() const
00442 {
00443   return baryon;
00444 }
00445 
00446 G4int G4HEVector::getStrangenessNumber() const
00447 {
00448   return strangeness;
00449 }
00450 
00451 G4int 
00452 G4HEVector::getQuarkContent(G4int flavor)
00453    {
00454      if(flavor > 0 && flavor < 8)
00455      { 
00456        G4int check;
00457        check = FillQuarkContent();
00458        if(check  != code) 
00459        { 
00460          return 0;
00461        }
00462        else
00463        {
00464          return theQuarkContent[flavor-1];
00465        }
00466      }
00467      else
00468      {
00469        return 0;
00470      }
00471    } 
00472 
00473 G4int 
00474 G4HEVector::getAntiQuarkContent(G4int flavor)
00475    {
00476      if(flavor > 0 && flavor < 8)
00477      { 
00478        G4int check;
00479        check = FillQuarkContent();
00480        if(check  != code) 
00481        { 
00482          return 0;
00483        }
00484        else
00485        {
00486          return theAntiQuarkContent[flavor-1];
00487        }
00488      }
00489      else
00490      {
00491        return 0;
00492      }
00493    } 
00494 
00495 void 
00496 G4HEVector::setZero()
00497    {
00498      px            = 0.0;
00499      py            = 0.0;
00500      pz            = 0.0;
00501      energy        = 0.0;
00502      kineticEnergy = 0.0;
00503      mass          = 0.0;
00504      charge        = 0.0;
00505      timeOfFlight  = 0.0;
00506      side          = 0;
00507      flag          = false;
00508      code          = 0;
00509      particleName  = "";
00510      particleType  = "";
00511      baryon        = 0;
00512      strangeness   = 0;
00513    }
00514 
00515 void 
00516 G4HEVector::Add( const G4HEVector & p1, const G4HEVector & p2 )
00517    {
00518      px  = p1.px + p2.px;
00519      py  = p1.py + p2.py;
00520      pz  = p1.pz + p2.pz;
00521      energy = p1.energy + p2.energy;
00522      G4double b = energy*energy - px*px - py*py - pz*pz;
00523      if( b < 0 )
00524        mass = -1. * std::sqrt( -b );
00525      else
00526        mass = std::sqrt( b );
00527      kineticEnergy = Amax(0.,energy - mass);
00528      charge        = p1.charge + p2.charge;
00529      code          = 0;
00530      particleName  = "";
00531      particleType  = "";
00532      baryon        = 0;
00533      strangeness   = 0;
00534    }
00535 
00536 void 
00537 G4HEVector::Sub( const G4HEVector & p1, const G4HEVector & p2 )
00538    {
00539      px  = p1.px - p2.px;
00540      py  = p1.py - p2.py;
00541      pz  = p1.pz - p2.pz;
00542      energy = p1.energy - p2.energy;
00543      G4double b = energy*energy - px*px - py*py - pz*pz;
00544      if( b < 0 )
00545        mass = -1. * std::sqrt( -b );
00546      else
00547        mass = std::sqrt( b );
00548      kineticEnergy = Amax(0.,energy - mass);
00549      charge        = p1.charge - p2.charge;
00550      code          = 0;
00551      particleName  = "";
00552      particleType  = "";
00553      baryon        = 0;
00554      strangeness   = 0;
00555    }
00556 
00557 void 
00558 G4HEVector::Lor( const G4HEVector & p1, const G4HEVector & p2 )
00559    {
00560      G4double a;
00561      a  = ( Dot(p1,p2)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
00562      px = p1.px + a*p2.px;
00563      py = p1.py + a*p2.py;
00564      pz = p1.pz + a*p2.pz; 
00565      energy = std::sqrt( sqr(p1.mass) + px*px + py*py + pz*pz);
00566      mass = p1.mass;
00567      kineticEnergy = Amax(0.,energy - mass);
00568      timeOfFlight  = p1.timeOfFlight;
00569      side          = p1.side;
00570      flag          = p1.flag;
00571      code          = p1.code;
00572      particleName  = p1.particleName;
00573      particleType  = p1.particleType; 
00574      baryon        = p1.baryon;
00575      strangeness   = p1.strangeness;
00576    }
00577 
00578 G4double G4HEVector::CosAng(const G4HEVector& p) const
00579 {
00580   G4double a = std::sqrt( (px*px + py*py + pz*pz)*(p.px*p.px + p.py*p.py + p.pz*p.pz) );
00581   if (a != 0.0) {
00582     a = (px*p.px + py*p.py + pz*p.pz)/a;
00583     if (std::fabs(a) > 1.0) {
00584       if (a < 0.0) a = -1.0;
00585       else a = 1.0;
00586     }   
00587   }
00588   return a;
00589 }
00590 
00591 G4double 
00592 G4HEVector::Ang(const G4HEVector & p )
00593    {
00594      G4double a = std::sqrt( (px*px + py*py + pz*pz)*(p.px*p.px + p.py*p.py + p.pz*p.pz) );
00595      if( a != 0.0 ) 
00596        {
00597          a = (px*p.px + py*p.py + pz*p.pz)/a;
00598          if( std::fabs(a) > 1.0 ) 
00599          {
00600            if(a<0.0) a=-1.0;
00601            else a=1.0;
00602          } 
00603        }
00604      return std::acos(a);
00605    }      
00606 
00607 G4double 
00608 G4HEVector::Dot4( const G4HEVector & p1, const G4HEVector & p2)
00609    {
00610      return ( p1.energy*p2.energy - p1.px*p2.px - p1.py*p2.py - p1.pz*p2.pz );
00611    } 
00612   
00613 G4double 
00614 G4HEVector::Impu( const G4HEVector & p1, const G4HEVector & p2)
00615    {
00616      return ( - sqr( p1.energy  - p2.energy)  
00617               + sqr( p1.px      - p2.px)
00618               + sqr( p1.py      - p2.py) 
00619               + sqr( p1.pz      - p2.pz) );
00620    }
00621 
00622 void 
00623 G4HEVector::Add3( const G4HEVector & p1, const G4HEVector & p2)
00624    {
00625      px =  p1.px + p2.px;
00626      py =  p1.py + p2.py;
00627      pz =  p1.pz + p2.pz;    
00628      return;
00629    }
00630 
00631 void 
00632 G4HEVector::Sub3( const G4HEVector & p1, const G4HEVector & p2)
00633    {
00634      px =  p1.px - p2.px;
00635      py =  p1.py - p2.py;
00636      pz =  p1.pz - p2.pz;
00637      return;
00638    }
00639 
00640 void 
00641 G4HEVector::Cross( const G4HEVector & p1, const G4HEVector & p2)
00642    {
00643      G4double tpx = p1.py * p2.pz - p1.pz * p2.py;
00644      G4double tpy = p1.pz * p2.px - p1.px * p2.pz;
00645      G4double tpz = p1.px * p2.py - p1.py * p2.px;
00646      px=tpx;
00647      py=tpy;
00648      pz=tpz;
00649      return;
00650    }  
00651  
00652 G4double 
00653 G4HEVector::Dot( const G4HEVector & p1, const G4HEVector & p2)
00654    {
00655      return ( p1.px * p2.px + p1.py * p2.py + p1.pz * p2.pz );
00656    }
00657 
00658 void 
00659 G4HEVector::Smul( const G4HEVector & p, G4double h)
00660    {
00661      px =  h * p.px;
00662      py =  h * p.py;
00663      pz =  h * p.pz;
00664      return;
00665    }   
00666 
00667 void 
00668 G4HEVector::SmulAndUpdate( const G4HEVector & p, G4double h)
00669    {
00670      px = h * p.px;
00671      py = h * p.py;
00672      pz = h * p.pz;
00673      mass          = p.mass;
00674      energy        = std::sqrt(px*px + py*py + pz*pz + mass*mass);
00675      kineticEnergy = energy - mass;
00676      charge        = p.charge;
00677      timeOfFlight  = p.timeOfFlight;
00678      side          = p.side;
00679      flag          = p.flag;
00680      code          = p.code;
00681      particleName  = p.particleName;
00682      particleType  = p.particleType;
00683      baryon        = p.baryon;
00684      strangeness   = p.strangeness;
00685      return;
00686    }
00687  
00688 void 
00689 G4HEVector::Norz( const G4HEVector & p )
00690    {
00691      G4double a =   p.px*p.px + p.py*p.py + p.pz*p.pz;
00692      if (a > 0.0) a = 1./std::sqrt(a);
00693      px = a * p.px;
00694      py = a * p.py;
00695      pz = a * p.pz;
00696      mass          = p.mass;
00697      energy        = std::sqrt(px*px + py*py + pz*pz + mass*mass);
00698      kineticEnergy = energy - mass;
00699      charge        = p.charge;
00700      timeOfFlight  = p.timeOfFlight;
00701      side          = p.side;
00702      flag          = p.flag;
00703      code          = p.code; 
00704      particleName  = p.particleName;
00705      particleType  = p.particleType;
00706      baryon        = p.baryon;
00707      strangeness   = p.strangeness;
00708      return;
00709    }
00710 
00711 G4double G4HEVector::Length() const
00712 {
00713   return std::sqrt(px*px + py*py + pz*pz);
00714 }
00715 
00716 void 
00717 G4HEVector::Exch( G4HEVector & p1)
00718    {
00719      G4HEVector mx = *this;
00720      *this = p1;
00721      p1 = mx;
00722      return; 
00723    }      
00724 
00725 void 
00726 G4HEVector::Defs1( const G4HEVector & p1, const G4HEVector & p2)
00727    {
00728      G4double pt2 = p2.px*p2.px + p2.py*p2.py;
00729      if (pt2 > 0.0)
00730         {
00731           G4double ph, qx, qy, qz;
00732           G4double a     = std::sqrt(p2.px*p2.px + p2.py*p2.py + p2.pz*p2.pz);
00733           G4double cost  = p2.pz/a; 
00734           G4double sint  = 0.5 * (std::sqrt(std::fabs((1.-cost)*(1.+cost))) + std::sqrt(pt2)/a);
00735           if(p2.py < 0.) ph = 1.5*pi;
00736           else ph = halfpi;
00737           if( p2.px != 0.0) 
00738              ph = std::atan2(p2.py,p2.px);             
00739           qx =   cost*std::cos(ph)*p1.px - std::sin(ph)*p1.py
00740                + sint*std::cos(ph)*p1.pz;
00741           qy =   cost*std::sin(ph)*p1.px + std::cos(ph)*p1.py
00742                + sint*std::sin(ph)*p1.pz;
00743           qz = - sint        *p1.px 
00744                + cost        *p1.pz;
00745           px = qx; 
00746           py = qy;
00747           pz = qz;     
00748         }
00749      else
00750         {
00751           px = p1.px;
00752           py = p1.py;
00753           pz = p1.pz;               
00754         } 
00755    }
00756  
00757 void 
00758 G4HEVector::Defs( const G4HEVector & p1, const G4HEVector & p2,
00759                   G4HEVector & my,       G4HEVector & mz )
00760    {
00761      my = p1;
00762      mz = p2;
00763      px = my.py*mz.pz - my.pz*mz.py;
00764      py = my.pz*mz.px - my.px*mz.pz;
00765      pz = my.px*mz.py - my.py*mz.px;
00766      my.px = mz.py*pz - mz.pz*py;
00767      my.py = mz.pz*px - mz.px*pz;
00768      my.pz = mz.px*py - mz.py*px;
00769      G4double pp;
00770      pp = std::sqrt(px*px + py*py + pz*pz);
00771      if (pp > 0.)
00772         {
00773           pp = 1./pp; 
00774           px = px*pp ;
00775           py = py*pp ;
00776           pz = pz*pp ;
00777         }
00778      pp = std::sqrt(my.px*my.px + my.py*my.py + my.pz*my.pz);
00779      if (pp > 0.)
00780         {
00781           pp = 1./pp;
00782           my.px = my.px*pp ;
00783           my.py = my.py*pp ;
00784           my.pz = my.pz*pp ;
00785         }
00786      pp = std::sqrt(mz.px*mz.px + mz.py*mz.py + mz.pz*mz.pz);
00787      if (pp > 0.)
00788         {
00789           pp = 1./pp;
00790           mz.px = mz.px*pp ;
00791           mz.py = mz.py*pp ;
00792           mz.pz = mz.pz*pp ;
00793         }
00794      return; 
00795    }  
00796              
00797 void 
00798 G4HEVector::Trac( const G4HEVector & p1, const G4HEVector & mx,
00799               const G4HEVector & my, const G4HEVector & mz)
00800    {
00801      G4double qx, qy, qz;
00802      qx =   mx.px*p1.px + mx.py*p1.py + mx.pz*p1.pz;
00803      qy =   my.px*p1.px + my.py*p1.py + my.pz*p1.pz;
00804      qz =   mz.px*p1.px + mz.py*p1.py + mz.pz*p1.pz;
00805      px = qx ;
00806      py = qy ;
00807      pz = qz ;
00808      return;
00809    }                 
00810 
00811 void 
00812 G4HEVector::setDefinition(G4String name)
00813    {
00814         if(name == "PionPlus")
00815           { 
00816             mass = 0.1395700;
00817             charge = 1.;
00818             code = 211;
00819             particleType = "Meson";
00820             particleName = name;
00821             baryon       = 0;
00822             strangeness  = 0;
00823           }
00824         else if(name == "PionZero")
00825           {
00826              mass = 0.1349764;
00827              charge = 0.;
00828              code = 111;
00829              particleType = "Meson";
00830              particleName = name;
00831              baryon       = 0;
00832              strangeness  = 0; 
00833           }
00834         else if(name == "PionMinus")
00835           {
00836               mass = 0.1395700;
00837               charge = -1.;
00838               code = -211;
00839               particleType = "Meson";
00840               particleName = name;
00841               baryon       = 0;
00842               strangeness  = 0;
00843           }        
00844         else if(name == "KaonPlus")
00845           {
00846               mass = 0.493677;
00847               charge = 1.;
00848               code = 321;
00849               particleType = "Meson";
00850               particleName = name;
00851               baryon       = 0;
00852               strangeness  = 1;
00853           }
00854         else if(name == "KaonZero")
00855           {
00856               mass = 0.497672;
00857               charge = 0.;
00858               code = 311;
00859               particleType = "Meson";
00860               particleName = name;
00861               baryon       = 0;
00862               strangeness  = 1;
00863           }
00864         else if(name == "AntiKaonZero")
00865           {
00866               mass = 0.497672;
00867               charge = 0.;
00868               code = -311;
00869               particleType = "Meson";
00870               particleName = name;
00871               baryon       = 0;
00872               strangeness  =-1;
00873           }
00874         else if(name == "KaonMinus")
00875           {
00876               mass = 0.493677;
00877               charge = -1.;
00878               code = -321;
00879               particleType = "Meson";
00880               particleName = name;
00881               baryon       = 0;
00882               strangeness  = -1;          
00883           }
00884         else if(name == "KaonZeroShort")
00885           {
00886               mass = 0.497672;
00887               charge = 0.;
00888               code = 310;
00889               particleType = "Meson";
00890               particleName = name;
00891               baryon       = 0;
00892               strangeness  = 0;
00893           }
00894         else if(name == "KaonZeroLong")
00895           {
00896               mass = 0.497672;
00897               charge = 0.;
00898               code = 130;
00899               particleType = "Meson";
00900               particleName = name;
00901               baryon       = 0;
00902               strangeness  = 0;
00903           }
00904         else if(name == "Proton")
00905           {
00906               mass = 0.9382723;
00907               charge = 1.;
00908               code = 2212;
00909               particleType = "Baryon";
00910               particleName = name;
00911               baryon       = 1;
00912               strangeness  = 0;
00913           }
00914         else if(name == "AntiProton")
00915           {
00916               mass = 0.9382723;
00917               charge = -1.;
00918               code = -2212;
00919               particleType = "AntiBaryon";
00920               particleName = name;
00921               baryon       = -1;
00922               strangeness  = 0;
00923           }
00924         else if(name == "Neutron")
00925           {
00926               mass = 0.93956563;
00927               charge = 0.;
00928               code = 2112;
00929               particleType = "Baryon";
00930               particleName = name;
00931               baryon       = 1;
00932               strangeness  = 0;
00933           }
00934         else if(name == "AntiNeutron")
00935           {
00936               mass = 0.93956563;
00937               charge = 0.;
00938               code = -2112;
00939               particleType = "AntiBaryon";
00940               particleName = name;
00941               baryon = -1;
00942               strangeness  = 0;
00943           }
00944         else if(name == "Lambda")
00945           {
00946               mass = 1.115684; 
00947               charge = 0.;
00948               code = 3122;
00949               particleType = "Baryon";
00950               particleName = name;
00951               baryon = 1;
00952               strangeness  = -1;
00953           }
00954         else if(name == "AntiLambda")
00955           {
00956               mass = 1.115684;
00957               charge = 0.;
00958               code = -3122;
00959               particleType = "AntiBaryon";
00960               particleName = name;
00961               baryon = -1;
00962               strangeness  = 1;
00963           }
00964         else if(name == "SigmaPlus")
00965           {
00966               mass = 1.18937;
00967               charge = 1.;
00968               code = 3222;
00969               particleType = "Baryon";
00970               particleName = name;
00971               baryon       = 1;
00972               strangeness  = -1;
00973           }
00974         else if(name == "SigmaZero") 
00975           {
00976               mass = 1.19255;
00977               charge = 0.;
00978               code = 3212;
00979               particleType = "Baryon";
00980               particleName = name;
00981               baryon       = 1;
00982               strangeness  = -1;
00983           }
00984         else if(name == "SigmaMinus")
00985           {
00986               mass = 1.19744;
00987               charge = -1.;
00988               code = 3112;
00989               particleType = "Baryon";
00990               particleName = name;
00991               baryon       = 1;
00992               strangeness  = -1;
00993           }
00994         else if(name == "AntiSigmaPlus") 
00995           {
00996               mass = 1.18937;
00997               charge = -1.;
00998               code = -3222;
00999               particleType = "AntiBaryon";
01000               particleName = name;
01001               baryon = -1;
01002               strangeness  = 1;
01003           }
01004         else if(name == "AntiSigmaZero")
01005           {
01006               mass = 1.19255;
01007               charge = 0.;
01008               code = -3212;
01009               particleType = "AntiBaryon";
01010               particleName = name;
01011               baryon       = -1;
01012               strangeness  = 1;
01013           }        
01014         else if(name == "AntiSigmaMinus")
01015           {
01016               mass = 1.19744;
01017               charge = 1.;
01018               code = -3112;
01019               particleType = "AntiBaryon";
01020               particleName = name;
01021               baryon       = -1;
01022               strangeness  = 1;
01023           }
01024         else if(name == "XiZero")
01025           {
01026               mass = 1.3149;
01027               charge = 0.;
01028               code = 3322;
01029               particleType = "Baryon";
01030               particleName = name;
01031               baryon       = 1;
01032               strangeness  = -2;
01033           }        
01034         else if(name == "XiMinus")
01035           {
01036               mass = 1.32132;
01037               charge = -1.;
01038               code = 3312;
01039               particleType = "Baryon";
01040               particleName = name;
01041               baryon       = 1;
01042               strangeness  = -2;
01043           }
01044         else if(name == "AntiXiZero")
01045           {
01046                mass = 1.3149;
01047                charge = 0.;
01048                code = -3322;
01049                particleType = "AntiBaryon";
01050                particleName = name;
01051                baryon = -1;
01052                strangeness  = 2;
01053           }        
01054         else if(name == "AntiXiMinus")
01055           {
01056                mass = 1.32132;
01057                charge = 1.;
01058                code = -3312;
01059                particleType = "AntiBaryon";
01060                particleName = name;
01061                baryon       = -1;
01062                strangeness  = 2;
01063           }
01064         else if(name == "OmegaMinus")
01065           {
01066                mass = 1.67245;
01067                charge = -1.;
01068                code = 3334;
01069                particleType = "Baryon";
01070                particleName = name;
01071                baryon       = 1;
01072                strangeness  = -3;
01073           }        
01074         else if(name == "AntiOmegaMinus")
01075           {
01076                mass = 1.67245;
01077                charge = 1.;
01078                code = -3334;
01079                particleType = "AntiBaryon";
01080                particleName = name;
01081                baryon       = -1;
01082                strangeness  = 3;
01083           }
01084         else if(name == "Deuteron")
01085           {
01086                mass = 1.875613;
01087                charge = 1.;
01088                code = 0;
01089                particleType = "Nucleus";
01090                particleName = name;
01091                baryon       = 2;
01092                strangeness  = 0;
01093           }        
01094         else if(name == "Triton")
01095           {
01096                mass = 2.80925;
01097                charge = 1.;
01098                code = 0;
01099                particleType = "Nucleus";
01100                particleName = name;
01101                baryon       = 3;
01102                strangeness  = 0;
01103           }
01104         else if(name == "Alpha")
01105           {
01106                mass = 3.727417;
01107                charge = 2.;
01108                code = 0;
01109                particleType = "Nucleus";
01110                particleName = name;
01111                baryon       = 4;
01112                strangeness  = 0;
01113           }
01114         else if(name == "Gamma")
01115           {
01116                mass = 0.;
01117                charge = 0.;
01118                code = 22;
01119                particleType = "Boson";
01120                particleName = name;
01121                baryon = 0;
01122                strangeness  = 0;
01123           }
01124         else
01125           {
01126                G4cout << "particle " << name << " not known in this generator!!" << G4endl;
01127                return;
01128           }
01129         px = 0.;
01130         py = 0.;
01131         pz = 0.;
01132         kineticEnergy = 0.;
01133         energy = mass;
01134         timeOfFlight = 0.;
01135         side = 0;
01136         flag = false;
01137         return;
01138    } 
01139 
01140 G4int G4HEVector::FillQuarkContent()
01141 {
01142   // Calculate quark and anti-quark content
01143   // Return value is PDG encoding for this particle
01144   // Error if return value is differnt from this->thePDGEncoding
01145 
01146   G4int tempPDGcode = code;
01147   G4double echarge = 1.;
01148 
01149   for (G4int flavor=0; flavor<NumberOfQuarkFlavor; flavor++){
01150     theQuarkContent[flavor] =0;
01151     theAntiQuarkContent[flavor] =0;
01152   }
01153 
01154   G4int temp = std::abs(tempPDGcode);
01155   G4int multiplet = temp/10000;
01156   temp -= G4int(multiplet*10000);
01157   G4int quark1 = temp/1000;
01158   temp -= G4int(quark1*1000);
01159   G4int quark2 = temp/100;
01160   temp -= G4int(quark2*100);
01161   G4int quark3 = temp/10;
01162   temp -= G4int(quark3*10);
01163 
01164   // G4int spin= (temp-1);  DHW 19 May 2011: variable set but not used
01165 
01166   if (particleType =="quark") {
01167     if (tempPDGcode>0){
01168       if (tempPDGcode<=NumberOfQuarkFlavor){
01169         theQuarkContent[tempPDGcode-1] =1;
01170       } else {
01171         //  --- thePDGEncoding is wrong 
01172         tempPDGcode = 0;
01173       }
01174     } else {
01175       G4int temp0 = -1*tempPDGcode; 
01176       if (temp0 <= NumberOfQuarkFlavor) {
01177         theAntiQuarkContent[temp0-1] =1;
01178       } else {
01179         //  --- thePDGEncoding is wrong 
01180         tempPDGcode = 0;
01181       }
01182     }
01183 
01184   } else if (particleType == "Meson") {
01185     //   -- exceptions --
01186     // if (tempPDGcode == 310) spin = 0;        //K0s
01187     // DHW 19 May 2011: variable set but not used
01188 
01189     if (tempPDGcode == 130) {     //K0l
01190       // spin = 0;   DHW 19 May 2011: variable set but not used        
01191       quark2 = 3;
01192       quark3 = 1;
01193     }
01194 
01195     if (quark1 !=0) 
01196      { 
01197        tempPDGcode = 0; 
01198      } 
01199     if ((quark2==0)||(quark3==0)){ 
01200       tempPDGcode = 0;
01201      }
01202     if (quark2<quark3) { 
01203       tempPDGcode = 0;
01204     }
01205     // check quark flavor
01206     if (quark2>=NumberOfQuarkFlavor){
01207       tempPDGcode = 0;
01208     }
01209     // check heavier quark type
01210     if (quark2 & 1) {
01211       // down type qurak
01212       if (tempPDGcode >0) {
01213         theQuarkContent[quark3-1] =1;
01214         theAntiQuarkContent[quark2-1] =1;
01215       } else {
01216         theQuarkContent[quark2-1] =1;
01217         theAntiQuarkContent[quark3-1] =1;
01218       }
01219     } else {
01220       // up type quark
01221       if (tempPDGcode >0) {
01222         theQuarkContent[quark2-1] =1;
01223         theAntiQuarkContent[quark3-1] =1;
01224       } else {
01225         theQuarkContent[quark3-1] =1;
01226         theAntiQuarkContent[quark2-1] =1;
01227       }
01228     }
01229     // check charge
01230     G4double totalCharge = 0.0;
01231     for (G4int flavor= 0; flavor<NumberOfQuarkFlavor-1; flavor+=2){
01232       totalCharge += (-1./3.)*echarge*theQuarkContent[flavor];
01233       totalCharge += 1./3.*echarge*theAntiQuarkContent[flavor];
01234       totalCharge += 2./3.*echarge*theQuarkContent[flavor+1];
01235       totalCharge += (-2./3.)*echarge*theAntiQuarkContent[flavor+1];
01236     }
01237     if (std::abs(totalCharge-charge)>0.1*echarge) { 
01238       tempPDGcode = 0;
01239     }
01240   } else if (particleType == "baryon"){
01241     // check Meson or not
01242     if ((quark1==0)||(quark2==0)||(quark3==0)){ 
01243       tempPDGcode = 0;
01244     }
01245     //exceptions
01246     if (std::abs(tempPDGcode) == 3122) { 
01247       // Lambda
01248       quark2 = 2;
01249       quark3 = 1;
01250       // spin = 1;  DHW 19 May 2011: variable set but not used
01251     } else if (std::abs(tempPDGcode) == 4122) { 
01252       // Lambda_c
01253       quark2 = 2;
01254       quark3 = 1;
01255       // spin = 1;  DHW 19 May 2011: variable set but not used
01256     } 
01257     // check quark flavor
01258     if ((quark1<quark2)||(quark2<quark3)||(quark1<quark3)) { 
01259       tempPDGcode = 0;
01260     }
01261     if (quark1>=NumberOfQuarkFlavor) {
01262       tempPDGcode = 0;
01263     }
01264     if (tempPDGcode >0) {
01265       theQuarkContent[quark1-1] ++;
01266       theQuarkContent[quark2-1] ++;
01267       theQuarkContent[quark3-1] ++;
01268     } else {
01269       theAntiQuarkContent[quark1-1] ++;
01270       theAntiQuarkContent[quark2-1] ++;
01271       theAntiQuarkContent[quark3-1] ++;
01272     }
01273     // check charge 
01274     G4double totalCharge = 0.0;
01275     for (G4int flavor= 0; flavor<NumberOfQuarkFlavor-1; flavor+=2){
01276       totalCharge += (-1./3.)*echarge*theQuarkContent[flavor];
01277       totalCharge += 1./3.*echarge*theAntiQuarkContent[flavor];
01278       totalCharge += 2./3.*echarge*theQuarkContent[flavor+1];
01279       totalCharge += (-2./3.)*echarge*theAntiQuarkContent[flavor+1];
01280     }
01281     if (std::abs(totalCharge - charge) > 0.1*echarge) tempPDGcode = 0;
01282 
01283   } else {
01284   }
01285   return tempPDGcode;
01286 }
01287 
01288 
01289 void G4HEVector::Print(G4int L) const
01290 {
01291   G4cout << "HEV: " 
01292          << L << " " << px << " " <<  py << " " <<  pz << " "
01293          << energy << " " << mass << " " << charge << " " 
01294          << timeOfFlight << " " << side << " " << flag << " " 
01295          << code << " " << baryon << " " << particleName << G4endl;
01296   return;                         
01297 }

Generated on Mon May 27 17:48:31 2013 for Geant4 by  doxygen 1.4.7