G4LEpp.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  // G4 Low energy model: n-n or p-p scattering
00028  // F.W. Jones, L.G. Greeniaus, H.P. Wellisch
00029 
00030 // FWJ 27-AUG-2010: extended Coulomb-suppressed data to 5 GeV
00031 
00032 #include "G4LEpp.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036 #include "G4ios.hh"
00037 
00038 // Initialization of static data arrays:
00039 #include "G4LEppData.hh"
00040 
00041 G4LEpp::G4LEpp():G4HadronicInteraction("G4LEpp")
00042 {
00043   //    theParticleChange.SetNumberOfSecondaries(1);
00044   //    SetMinEnergy(10.*MeV);
00045   //    SetMaxEnergy(1200.*MeV);
00046 
00047   SetCoulombEffects(0);
00048 
00049   SetMinEnergy(0.);
00050   SetMaxEnergy(5.*GeV);
00051 }
00052 
00053 G4LEpp::~G4LEpp()
00054 {
00055   //    theParticleChange.Clear();
00056 }
00057 
00058 
00059 void
00060 G4LEpp::SetCoulombEffects(G4int State)
00061 {
00062   if (State) {
00063     for(G4int i=0; i<NANGLE; i++)
00064     {
00065       sig[i] = SigCoul[i];
00066     }
00067     elab = ElabCoul;
00068     SetMaxEnergy(1.2*GeV);
00069   }
00070   else {
00071     for(G4int i=0; i<NANGLE; i++)
00072     {
00073       sig[i] = Sig[i];
00074     }
00075     elab = Elab;
00076     SetMaxEnergy(5.*GeV);
00077   }
00078 }
00079 
00080 
00081 G4HadFinalState*
00082 G4LEpp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
00083 {
00084     theParticleChange.Clear();
00085     const G4HadProjectile* aParticle = &aTrack;
00086 
00087     G4double P = aParticle->GetTotalMomentum();
00088     G4double Px = aParticle->Get4Momentum().x();
00089     G4double Py = aParticle->Get4Momentum().y();
00090     G4double Pz = aParticle->Get4Momentum().z();
00091     G4double ek = aParticle->GetKineticEnergy();
00092     G4ThreeVector theInitial = aParticle->Get4Momentum().vect();
00093 
00094     if (verboseLevel > 1) {
00095       G4double E = aParticle->GetTotalEnergy();
00096       G4double E0 = aParticle->GetDefinition()->GetPDGMass();
00097       G4double Q = aParticle->GetDefinition()->GetPDGCharge();
00098       G4int A = targetNucleus.GetA_asInt();
00099       G4int Z = targetNucleus.GetZ_asInt();
00100       G4cout << "G4LEpp:ApplyYourself: incident particle: "
00101              << aParticle->GetDefinition()->GetParticleName() << G4endl;
00102       G4cout << "P = " << P/GeV << " GeV/c"
00103              << ", Px = " << Px/GeV << " GeV/c"
00104              << ", Py = " << Py/GeV << " GeV/c"
00105              << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
00106       G4cout << "E = " << E/GeV << " GeV"
00107              << ", kinetic energy = " << ek/GeV << " GeV"
00108              << ", mass = " << E0/GeV << " GeV"
00109              << ", charge = " << Q << G4endl;
00110       G4cout << "G4LEpp:ApplyYourself: material:" << G4endl;
00111       G4cout << "A = " << A
00112              << ", Z = " << Z
00113              << ", atomic mass " 
00114              <<  G4Proton::Proton()->GetPDGMass()/GeV << "GeV" 
00115              << G4endl;
00116       //
00117       // GHEISHA ADD operation to get total energy, mass, charge
00118       //
00119       E += proton_mass_c2;
00120       G4double E02 = E*E - P*P;
00121       E0 = std::sqrt(std::fabs(E02));
00122       if (E02 < 0)E0 *= -1;
00123       Q += Z;
00124       G4cout << "G4LEpp:ApplyYourself: total:" << G4endl;
00125       G4cout << "E = " << E/GeV << " GeV"
00126              << ", mass = " << E0/GeV << " GeV"
00127              << ", charge = " << Q << G4endl;
00128     }
00129 
00130     // Find energy bin
00131 
00132     G4int je1 = 0;
00133     G4int je2 = NENERGY - 1;
00134     ek = ek/GeV;
00135     do {
00136       G4int midBin = (je1 + je2)/2;
00137       if (ek < elab[midBin])
00138         je2 = midBin;
00139       else
00140         je1 = midBin;
00141     } while (je2 - je1 > 1); 
00142     G4double delab = elab[je2] - elab[je1];
00143 
00144     // Sample the angle
00145 
00146     G4float sample = G4UniformRand();
00147     G4int ke1 = 0;
00148     G4int ke2 = NANGLE - 1;
00149     G4double dsig = sig[je2][0] - sig[je1][0];
00150     G4double rc = dsig/delab;
00151     G4double b = sig[je1][0] - rc*elab[je1];
00152     G4double sigint1 = rc*ek + b;
00153     G4double sigint2 = 0.;
00154 
00155     if (verboseLevel > 1) G4cout << "sample=" << sample << G4endl
00156                                  << ke1 << " " << ke2 << " " 
00157                                  << sigint1 << " " << sigint2 << G4endl;
00158 
00159     do {
00160       G4int midBin = (ke1 + ke2)/2;
00161       dsig = sig[je2][midBin] - sig[je1][midBin];
00162       rc = dsig/delab;
00163       b = sig[je1][midBin] - rc*elab[je1];
00164       G4double sigint = rc*ek + b;
00165       if (sample < sigint) {
00166         ke2 = midBin;
00167         sigint2 = sigint;
00168       }
00169       else {
00170         ke1 = midBin;
00171         sigint1 = sigint;
00172       }
00173       if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " " 
00174                                   << sigint1 << " " << sigint2 << G4endl;
00175     } while (ke2 - ke1 > 1); 
00176 
00177     dsig = sigint2 - sigint1;
00178     rc = 1./dsig;
00179     b = ke1 - rc*sigint1;
00180     G4double kint = rc*sample + b;
00181     G4double theta = (0.5 + kint)*pi/180.;
00182     if (theta < 0.) { theta = 0.; }
00183 
00184     if (verboseLevel > 1) {
00185       G4cout << "   energy bin " << je1 << " energy=" << elab[je1] << G4endl;
00186       G4cout << "   angle bin " << kint << " angle=" << theta/degree << G4endl;
00187     }
00188 
00189     // Get the target particle
00190     G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
00191 
00192     G4double E1 = aParticle->GetTotalEnergy();
00193     G4double M1 = aParticle->GetDefinition()->GetPDGMass();
00194     G4double E2 = targetParticle->GetTotalEnergy();
00195     G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
00196     G4double totalEnergy = E1 + E2;
00197     G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
00198 
00199     // Transform into centre of mass system
00200 
00201     G4double px = (M2/pseudoMass)*Px;
00202     G4double py = (M2/pseudoMass)*Py;
00203     G4double pz = (M2/pseudoMass)*Pz;
00204     G4double p = std::sqrt(px*px + py*py + pz*pz);
00205 
00206     if (verboseLevel > 1) {
00207       G4cout << "  E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
00208       G4cout << "  E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
00209       G4cout << "  particle  1 momentum in CM " << px/GeV << " " << py/GeV << " "
00210              << pz/GeV << " " << p/GeV << G4endl;
00211     }
00212 
00213     // First scatter w.r.t. Z axis
00214     G4double phi = G4UniformRand()*twopi;
00215     G4double pxnew = p*std::sin(theta)*std::cos(phi);
00216     G4double pynew = p*std::sin(theta)*std::sin(phi);
00217     G4double pznew = p*std::cos(theta);
00218 
00219     // Rotate according to the direction of the incident particle
00220     if (px*px + py*py > 0) {
00221       G4double cost, sint, ph, cosp, sinp;
00222       cost = pz/p;
00223       sint = (std::sqrt(std::fabs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2;
00224       py < 0 ? ph = 3*halfpi : ph = halfpi;
00225       if (std::fabs(px) > 0.000001*GeV) ph = std::atan2(py,px);
00226       cosp = std::cos(ph);
00227       sinp = std::sin(ph);
00228       px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
00229       py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
00230       pz = (-sint*pxnew                  + cost*pznew);
00231     }
00232     else {
00233       px = pxnew;
00234       py = pynew;
00235       pz = pznew;
00236     }
00237 
00238     if (verboseLevel > 1) {
00239       G4cout << "  AFTER SCATTER..." << G4endl;
00240       G4cout << "  particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
00241            << pz/GeV << " " << p/GeV << G4endl;
00242     }
00243 
00244     // Transform to lab system
00245 
00246     G4double E1pM2 = E1 + M2;
00247     G4double betaCM  = P/E1pM2;
00248     G4double betaCMx = Px/E1pM2;
00249     G4double betaCMy = Py/E1pM2;
00250     G4double betaCMz = Pz/E1pM2;
00251     G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
00252 
00253     if (verboseLevel > 1) {
00254       G4cout << "  betaCM " << betaCMx << " " << betaCMy << " "
00255              << betaCMz << " " << betaCM << G4endl;
00256       G4cout << "  gammaCM " << gammaCM << G4endl;
00257     }
00258 
00259     // Now following GLOREN...
00260 
00261     G4double BETA[5], PA[5], PB[5];
00262     BETA[1] = -betaCMx;
00263     BETA[2] = -betaCMy;
00264     BETA[3] = -betaCMz;
00265     BETA[4] = gammaCM;
00266 
00267     //The incident particle...
00268 
00269     PA[1] = px;
00270     PA[2] = py;
00271     PA[3] = pz;
00272     PA[4] = std::sqrt(M1*M1 + p*p);
00273 
00274     G4double BETPA  = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
00275     G4double BPGAM  = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
00276 
00277     PB[1] = PA[1] + BPGAM  * BETA[1];
00278     PB[2] = PA[2] + BPGAM  * BETA[2];
00279     PB[3] = PA[3] + BPGAM  * BETA[3];
00280     PB[4] = (PA[4] - BETPA) * BETA[4];
00281 
00282     G4DynamicParticle* newP = new G4DynamicParticle;
00283     newP->SetDefinition(const_cast<G4ParticleDefinition *>(aParticle->GetDefinition()) );
00284     newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
00285 
00286     //The target particle...
00287 
00288     PA[1] = -px;
00289     PA[2] = -py;
00290     PA[3] = -pz;
00291     PA[4] = std::sqrt(M2*M2 + p*p);
00292 
00293     BETPA  = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
00294     BPGAM  = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
00295 
00296     PB[1] = PA[1] + BPGAM  * BETA[1];
00297     PB[2] = PA[2] + BPGAM  * BETA[2];
00298     PB[3] = PA[3] + BPGAM  * BETA[3];
00299     PB[4] = (PA[4] - BETPA) * BETA[4];
00300 
00301     targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
00302 
00303     if (verboseLevel > 1) {
00304       G4cout << "  particle 1 momentum in LAB " 
00305            << newP->GetMomentum()/GeV 
00306            << " " << newP->GetTotalMomentum()/GeV << G4endl;
00307       G4cout << "  particle 2 momentum in LAB " 
00308            << targetParticle->GetMomentum()/GeV 
00309            << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
00310       G4cout << "  TOTAL momentum in LAB " 
00311            << (newP->GetMomentum()+targetParticle->GetMomentum())/GeV 
00312            << " " 
00313            << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
00314            << G4endl;
00315     }
00316 
00317     theParticleChange.SetMomentumChange( newP->GetMomentumDirection());
00318     theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
00319     delete newP;
00320 
00321     // Recoil particle
00322     theParticleChange.AddSecondary(targetParticle);    
00323     return &theParticleChange;
00324 }
00325 
00326  // end of file

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