Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-04 04:35:19

0001 
0002 #include "G4PhysicalConstants.hh"
0003 #include "G4SystemOfUnits.hh"
0004 #include "G4ParticleTable.hh"
0005 #include "G4ParticleDefinition.hh"
0006 #include "Randomize.hh"
0007 #include "G4NucleiProperties.hh"
0008 #include <math.h>
0009 #include "TMath.h"
0010 
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 
0013 #include "SimG4Core/CustomPhysics/interface/CMSSQNeutronAnnih.h"
0014 #include "SimG4Core/CustomPhysics/interface/CMSSQ.h"
0015 
0016 CMSSQNeutronAnnih::CMSSQNeutronAnnih(double mass) : G4HadronicInteraction("SexaQuark-neutron annihilation") {
0017   SetMinEnergy(0.0 * GeV);
0018   SetMaxEnergy(100. * TeV);
0019 
0020   theSQ = CMSSQ::SQ(mass);
0021   theK0S = G4KaonZeroShort::KaonZeroShort();
0022   theAntiL = G4AntiLambda::AntiLambda();
0023   theProton = G4Proton::
0024       Proton();  //proton only used when the particle which the sexaquark hits is a deutereon and the neutron dissapears, so what stays behind is a proton
0025 }
0026 
0027 CMSSQNeutronAnnih::~CMSSQNeutronAnnih() {}
0028 
0029 //9Be momentum distribution from Jan Ryckebusch
0030 G4double CMSSQNeutronAnnih::momDistr(G4double x_in) {
0031   const int n_entries = 50;
0032 
0033   G4double CDF_k[n_entries] = {0,   0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,   1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
0034                                1.7, 1.8, 1.9, 2,   2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3,   3.1, 3.2, 3.3,
0035                                3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.,  4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9};
0036 
0037   G4double x[n_entries] = {0,
0038                            0.0038033182,
0039                            0.0187291764,
0040                            0.0510409777,
0041                            0.1048223609,
0042                            0.1807862863,
0043                            0.2756514534,
0044                            0.3825832103,
0045                            0.4926859745,
0046                            0.5970673837,
0047                            0.6887542272,
0048                            0.7637748784,
0049                            0.8212490273,
0050                            0.8627259608,
0051                            0.8911605331,
0052                            0.9099115186,
0053                            0.9220525854,
0054                            0.9300190818,
0055                            0.9355376091,
0056                            0.9397242185,
0057                            0.9432387722,
0058                            0.946438928,
0059                            0.9495023924,
0060                            0.9525032995,
0061                            0.9554669848,
0062                            0.9583936672,
0063                            0.9612770117,
0064                            0.9641067202,
0065                            0.9668727859,
0066                            0.9695676121,
0067                            0.9721815799,
0068                            0.9747092981,
0069                            0.9771426396,
0070                            0.9794740235,
0071                            0.9816956807,
0072                            0.9838003583,
0073                            0.9857816165,
0074                            0.9876331761,
0075                            0.9893513365,
0076                            0.9909333198,
0077                            0.992378513,
0078                            0.9936885054,
0079                            0.9948665964,
0080                            0.9959179448,
0081                            0.9968491104,
0082                            0.9976680755,
0083                            0.9983832508,
0084                            0.9990041784,
0085                            0.9995400073,
0086                            1};
0087 
0088   //now interpolate the above points for x_in
0089   G4double result = 9999;
0090   for (int i = 0; i < n_entries; i++) {
0091     if (x[i] > x_in) {
0092       result = (CDF_k[i] - CDF_k[i - 1]) / (x[i] - x[i - 1]) * (x_in - x[i - 1]) + CDF_k[i - 1];
0093       break;
0094     }
0095   }
0096   //ROOT::Math::Interpolator inter(n_entries, ROOT::Math::Interpolation::kAKIMA);
0097   //inter.SetData(n_entries,x,CDF_k);
0098   //result = inter.Eval(x_in);
0099 
0100   return result;
0101   //return 1;
0102 }
0103 
0104 G4HadFinalState* CMSSQNeutronAnnih::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus) {
0105   theParticleChange.Clear();
0106   const G4HadProjectile* aParticle = &aTrack;
0107   G4double ekin = aParticle->GetKineticEnergy();
0108 
0109   G4int A = targetNucleus.GetA_asInt();
0110   G4int Z = targetNucleus.GetZ_asInt();
0111 
0112   G4double m_K0S = G4KaonZeroShort::KaonZeroShort()->GetPDGMass();
0113   G4double m_L = G4AntiLambda::AntiLambda()->GetPDGMass();
0114 
0115   //G4double plab = aParticle->GetTotalMomentum();
0116 
0117   //    edm::LogVerbatim("CMSSWNeutronAnnih") << "CMSSQNeutronAnnih: Incident particle p (GeV), total Energy (GeV), particle name, eta ="
0118   //       << plab/GeV << "  "
0119   //       << aParticle->GetTotalEnergy()/GeV << "  "
0120   //       << aParticle->GetDefinition()->GetParticleName() << " "
0121   //       << aParticle->Get4Momentum();
0122 
0123   // Scattered particle referred to axis of incident particle
0124   //const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
0125 
0126   //G4int projPDG = theParticle->GetPDGEncoding();
0127   //    edm::LogVerbatim("CMSSWNeutronAnnih") << "CMSSQNeutronAnnih: for " << theParticle->GetParticleName()
0128   //           << " PDGcode= " << projPDG << " on nucleus Z= " << Z
0129   //           << " A= " << A << " N= " << N;
0130 
0131   G4LorentzVector lv1 = aParticle->Get4Momentum();
0132   edm::LogVerbatim("CMSSWNeutronAnnih") << "The neutron Fermi momentum (mag, x, y, z) "
0133                                         << targetNucleus.GetFermiMomentum().mag() / MeV << " "
0134                                         << targetNucleus.GetFermiMomentum().x() / MeV << " "
0135                                         << targetNucleus.GetFermiMomentum().y() / MeV << " "
0136                                         << targetNucleus.GetFermiMomentum().z() / MeV;
0137 
0138   //calculate fermi momentum
0139 
0140   G4double k_neutron = momDistr(G4UniformRand());
0141   G4double momentum_neutron = 0.1973 * GeV * k_neutron;
0142 
0143   G4double theta_neutron = TMath::ACos(2 * G4UniformRand() - 1);
0144   G4double phi_neutron = 2. * TMath::Pi() * G4UniformRand();
0145 
0146   G4double p_neutron_x = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Cos(phi_neutron);
0147   G4double p_neutron_y = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Sin(phi_neutron);
0148   G4double p_neutron_z = momentum_neutron * TMath::Cos(theta_neutron);
0149 
0150   //G4LorentzVector lv0(targetNucleus.GetFermiMomentum(), sqrt( pow(G4Neutron::Neutron()->GetPDGMass(),2) + targetNucleus.GetFermiMomentum().mag2()  ) );
0151   G4LorentzVector lv0(p_neutron_x,
0152                       p_neutron_y,
0153                       p_neutron_z,
0154                       sqrt(pow(G4Neutron::Neutron()->GetPDGMass(), 2) + momentum_neutron * momentum_neutron));
0155 
0156   //const G4Nucleus* aNucleus = &targetNucleus;
0157   G4double BENeutronInNucleus = 0;
0158   if (A != 0)
0159     BENeutronInNucleus = G4NucleiProperties::GetBindingEnergy(A, Z) / (A);
0160 
0161   edm::LogVerbatim("CMSSWNeutronAnnih") << "BE of nucleon in the nucleus (GeV): " << BENeutronInNucleus / GeV;
0162 
0163   G4LorentzVector lvBE(0, 0, 0, BENeutronInNucleus / GeV);
0164   G4LorentzVector lv = lv0 + lv1 - lvBE;
0165 
0166   // kinematiacally impossible ?
0167   G4double etot = lv0.e() + lv1.e() - lvBE.e();
0168   if (etot < theK0S->GetPDGMass() + theAntiL->GetPDGMass()) {
0169     theParticleChange.SetEnergyChange(ekin);
0170     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
0171     return &theParticleChange;
0172   }
0173 
0174   float newIonMass = targetNucleus.AtomicMass(A - 1, Z) * 931.5 * MeV;
0175   ;
0176   G4LorentzVector nlvIon(0, 0, 0, newIonMass);
0177 
0178   G4double theta_KS0_star = TMath::ACos(2 * G4UniformRand() - 1);
0179   G4double phi_KS0_star = 2. * TMath::Pi() * G4UniformRand();
0180 
0181   G4double p_K0S_star_x = TMath::Sin(theta_KS0_star) * TMath::Cos(phi_KS0_star);
0182   G4double p_K0S_star_y = TMath::Sin(theta_KS0_star) * TMath::Sin(phi_KS0_star);
0183   G4double p_K0S_star_z = TMath::Cos(theta_KS0_star);
0184 
0185   G4ThreeVector p(p_K0S_star_x, p_K0S_star_y, p_K0S_star_z);
0186   double m0 = lv.m();
0187   double m0_2 = m0 * m0;
0188   double m1_2 = m_K0S * m_K0S;
0189   double m2_2 = m_L * m_L;
0190 
0191   p *= 0.5 / m0 * sqrt(m0_2 * m0_2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m0_2 * m1_2 - 2 * m0_2 * m2_2 - 2 * m1_2 * m2_2);
0192   double p2 = p.mag2();
0193 
0194   G4LorentzVector nlvK0S(p, sqrt(p2 + m1_2));
0195   G4LorentzVector nlvAntiL(-p, sqrt(p2 + m2_2));
0196 
0197   // Boost out of the rest frame.
0198   nlvK0S.boost(lv.boostVector());
0199   nlvAntiL.boost(lv.boostVector());
0200 
0201   // now move to implement the interaction
0202   theParticleChange.SetStatusChange(stopAndKill);
0203   //theParticleChange.SetEnergyChange(ekin); // was 0.0
0204 
0205   G4DynamicParticle* aSec1 = new G4DynamicParticle(theK0S, nlvK0S);
0206   theParticleChange.AddSecondary(aSec1);
0207   G4DynamicParticle* aSec2 = new G4DynamicParticle(theAntiL, nlvAntiL);
0208   theParticleChange.AddSecondary(aSec2);
0209 
0210   const G4ParticleDefinition* theRemainingNucleusDef = theProton;
0211   if (A != 1)
0212     theRemainingNucleusDef = G4IonTable::GetIonTable()->GetIon(Z, A - 1);
0213   G4DynamicParticle* aSec3 = new G4DynamicParticle(theRemainingNucleusDef, nlvIon);
0214   theParticleChange.AddSecondary(aSec3);
0215 
0216   // return as is; we don't care about what happens to the nucleus
0217   return &theParticleChange;
0218 }