CustomPhysics/src/CMSSQNeutronAnnih.cc

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   if (x_in <= 0.0)
0090     return 0.;
0091   if (x_in >= 1.0)
0092     return CDF_k[n_entries - 1];
0093   for (int i = 1; i < n_entries; i++) {
0094     if (x[i] >= x_in) {
0095       return (CDF_k[i] - CDF_k[i - 1]) / (x[i] - x[i - 1]) * (x_in - x[i - 1]) + CDF_k[i - 1];
0096     }
0097   }
0098   return 0.0;
0099 }
0100
0101 G4HadFinalState* CMSSQNeutronAnnih::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus) {
0102   theParticleChange.Clear();
0103   const G4HadProjectile* aParticle = &aTrack;
0104   G4double ekin = aParticle->GetKineticEnergy();
0105
0106   G4int A = targetNucleus.GetA_asInt();
0107   G4int Z = targetNucleus.GetZ_asInt();
0108
0109   G4double m_K0S = G4KaonZeroShort::KaonZeroShort()->GetPDGMass();
0110   G4double m_L = G4AntiLambda::AntiLambda()->GetPDGMass();
0111
0112   //G4double plab = aParticle->GetTotalMomentum();
0113
0114   //    edm::LogVerbatim("CMSSWNeutronAnnih") << "CMSSQNeutronAnnih: Incident particle p (GeV), total Energy (GeV), particle name, eta ="
0115   //       << plab/GeV << "  "
0116   //       << aParticle->GetTotalEnergy()/GeV << "  "
0117   //       << aParticle->GetDefinition()->GetParticleName() << " "
0118   //       << aParticle->Get4Momentum();
0119
0120   // Scattered particle referred to axis of incident particle
0121   //const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
0122
0123   //G4int projPDG = theParticle->GetPDGEncoding();
0124   //    edm::LogVerbatim("CMSSWNeutronAnnih") << "CMSSQNeutronAnnih: for " << theParticle->GetParticleName()
0125   //           << " PDGcode= " << projPDG << " on nucleus Z= " << Z
0126   //           << " A= " << A << " N= " << N;
0127
0128   G4LorentzVector lv1 = aParticle->Get4Momentum();
0129   edm::LogVerbatim("CMSSWNeutronAnnih") << "The neutron Fermi momentum (mag, x, y, z) "
0130                                         << targetNucleus.GetFermiMomentum().mag() / MeV << " "
0131                                         << targetNucleus.GetFermiMomentum().x() / MeV << " "
0132                                         << targetNucleus.GetFermiMomentum().y() / MeV << " "
0133                                         << targetNucleus.GetFermiMomentum().z() / MeV;
0134
0135   //calculate fermi momentum
0136
0137   G4double k_neutron = momDistr(G4UniformRand());
0138   G4double momentum_neutron = 0.1973 * GeV * k_neutron;
0139
0140   G4double theta_neutron = TMath::ACos(2 * G4UniformRand() - 1);
0141   G4double phi_neutron = 2. * TMath::Pi() * G4UniformRand();
0142
0143   G4double p_neutron_x = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Cos(phi_neutron);
0144   G4double p_neutron_y = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Sin(phi_neutron);
0145   G4double p_neutron_z = momentum_neutron * TMath::Cos(theta_neutron);
0146
0147   //G4LorentzVector lv0(targetNucleus.GetFermiMomentum(), sqrt( pow(G4Neutron::Neutron()->GetPDGMass(),2) + targetNucleus.GetFermiMomentum().mag2()  ) );
0148   G4LorentzVector lv0(p_neutron_x,
0149                       p_neutron_y,
0150                       p_neutron_z,
0151                       sqrt(pow(G4Neutron::Neutron()->GetPDGMass(), 2) + momentum_neutron * momentum_neutron));
0152
0153   //const G4Nucleus* aNucleus = &targetNucleus;
0154   G4double BENeutronInNucleus = 0;
0155   if (A != 0)
0156     BENeutronInNucleus = G4NucleiProperties::GetBindingEnergy(A, Z) / (A);
0157
0158   edm::LogVerbatim("CMSSWNeutronAnnih") << "BE of nucleon in the nucleus (GeV): " << BENeutronInNucleus / GeV;
0159
0160   G4LorentzVector lvBE(0, 0, 0, BENeutronInNucleus / GeV);
0161   G4LorentzVector lv = lv0 + lv1 - lvBE;
0162
0163   // kinematiacally impossible ?
0164   G4double etot = lv0.e() + lv1.e() - lvBE.e();
0165   if (etot < theK0S->GetPDGMass() + theAntiL->GetPDGMass()) {
0166     theParticleChange.SetEnergyChange(ekin);
0167     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
0168     return &theParticleChange;
0169   }
0170
0171   float newIonMass = targetNucleus.AtomicMass(A - 1, Z) * 931.5 * MeV;
0172   ;
0173   G4LorentzVector nlvIon(0, 0, 0, newIonMass);
0174
0175   G4double theta_KS0_star = TMath::ACos(2 * G4UniformRand() - 1);
0176   G4double phi_KS0_star = 2. * TMath::Pi() * G4UniformRand();
0177
0178   G4double p_K0S_star_x = TMath::Sin(theta_KS0_star) * TMath::Cos(phi_KS0_star);
0179   G4double p_K0S_star_y = TMath::Sin(theta_KS0_star) * TMath::Sin(phi_KS0_star);
0180   G4double p_K0S_star_z = TMath::Cos(theta_KS0_star);
0181
0182   G4ThreeVector p(p_K0S_star_x, p_K0S_star_y, p_K0S_star_z);
0183   double m0 = lv.m();
0184   double m0_2 = m0 * m0;
0185   double m1_2 = m_K0S * m_K0S;
0186   double m2_2 = m_L * m_L;
0187
0188   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);
0189   double p2 = p.mag2();
0190
0191   G4LorentzVector nlvK0S(p, sqrt(p2 + m1_2));
0192   G4LorentzVector nlvAntiL(-p, sqrt(p2 + m2_2));
0193
0194   // Boost out of the rest frame.
0195   nlvK0S.boost(lv.boostVector());
0196   nlvAntiL.boost(lv.boostVector());
0197
0198   // now move to implement the interaction
0199   theParticleChange.SetStatusChange(stopAndKill);
0200   //theParticleChange.SetEnergyChange(ekin); // was 0.0
0201
0202   G4DynamicParticle* aSec1 = new G4DynamicParticle(theK0S, nlvK0S);
0203   theParticleChange.AddSecondary(aSec1);
0204   G4DynamicParticle* aSec2 = new G4DynamicParticle(theAntiL, nlvAntiL);
0205   theParticleChange.AddSecondary(aSec2);
0206
0207   const G4ParticleDefinition* theRemainingNucleusDef = theProton;
0208   if (A != 1)
0209     theRemainingNucleusDef = G4IonTable::GetIonTable()->GetIon(Z, A - 1);
0210   G4DynamicParticle* aSec3 = new G4DynamicParticle(theRemainingNucleusDef, nlvIon);
0211   theParticleChange.AddSecondary(aSec3);
0212
0213   // return as is; we don't care about what happens to the nucleus
0214   return &theParticleChange;
0215 }