CustomPhysics/src/G4muDarkBremsstrahlungModel.cc

0001 #include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h"
0002 #include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h"
0003 #include "SimG4Core/CustomPhysics/interface/G4APrime.h"
0004
0005 //Geant 4
0006 #include "G4ProcessTable.hh"
0007 #include "G4ProcTblElement.hh"
0008 #include "G4MuonMinus.hh"
0009 #include "G4MuonPlus.hh"
0010 #include "G4ProductionCutsTable.hh"
0011 #include <CLHEP/Units/SystemOfUnits.h>
0012 //Root
0013 #include "TFile.h"
0014 #include "TTree.h"
0015 //gsl
0016 #include <gsl/gsl_math.h>
0017 #include <gsl/gsl_integration.h>
0018
0019 using namespace std;
0020
0021 G4muDarkBremsstrahlungModel::G4muDarkBremsstrahlungModel(const G4String& scalefile,
0022                                                          const G4double biasFactor,
0023                                                          const G4ParticleDefinition* p,
0024                                                          const G4String& nam)
0025     : G4VEmModel(nam),
0026       mgfile(scalefile),
0027       cxBias(biasFactor),
0028       particle(nullptr),
0029       isMuon(true),
0030       probsup(1.0),
0031       isInitialised(false),
0032       method("forward_only"),
0033       mg_loaded(false) {
0034   if (p) {
0035     SetParticle(p);
0036   }  //Verify that the particle is a muon.
0037   theAPrime = G4APrime::APrime();
0038   MA = G4APrime::APrime()->GetPDGMass() / CLHEP::GeV;                        //Get the A' mass.
0039   muonMass = G4MuonMinus::MuonMinusDefinition()->GetPDGMass() / CLHEP::GeV;  //Get the muon mass
0040   highKinEnergy = HighEnergyLimit();
0041   lowKinEnergy = LowEnergyLimit();
0042   fParticleChange = nullptr;
0043 }
0044
0045 G4muDarkBremsstrahlungModel::~G4muDarkBremsstrahlungModel() {
0046   size_t n = partialSumSigma.size();
0047   if (n > 0) {
0048     for (size_t i = 0; i < n; i++) {
0049       delete partialSumSigma[i];
0050     }
0051   }
0052 }
0053
0054 void G4muDarkBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p) {
0055   particle = p;
0056
0057   if ((p == G4MuonPlus::MuonPlus()) || (p == G4MuonMinus::MuonMinus())) {
0058     isMuon = true;
0059   } else {
0060     isMuon = false;
0061   }
0062 }
0063
0064 void G4muDarkBremsstrahlungModel::Initialise(const G4ParticleDefinition* p, const G4DataVector& cuts) {
0065   if (p) {
0066     SetParticle(p);
0067   }
0068
0069   highKinEnergy = HighEnergyLimit();
0070   lowKinEnergy = LowEnergyLimit();
0071   const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
0072
0073   if (theCoupleTable) {
0074     G4int numOfCouples = theCoupleTable->GetTableSize();
0075     G4int nn = partialSumSigma.size();
0076     G4int nc = cuts.size();
0077     if (nn > 0) {
0078       for (G4int ii = 0; ii < nn; ii++) {
0079         G4DataVector* a = partialSumSigma[ii];
0080         if (a)
0081           delete a;
0082       }
0083       partialSumSigma.clear();
0084     }
0085     if (numOfCouples > 0) {
0086       for (G4int i = 0; i < numOfCouples; i++) {
0087         G4double cute = DBL_MAX;
0088         if (i < nc)
0089           cute = cuts[i];
0090         const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
0091         const G4Material* material = couple->GetMaterial();
0092         G4DataVector* dv = ComputePartialSumSigma(material, 0.5 * highKinEnergy, std::min(cute, 0.25 * highKinEnergy));
0093         partialSumSigma.push_back(dv);
0094       }
0095     }
0096   }
0097
0098   if (isInitialised)
0099     return;
0100   fParticleChange = GetParticleChangeForLoss();
0101   isInitialised = true;
0102 }
0103
0104 void G4muDarkBremsstrahlungModel::LoadMG()
0105 //Parses a Madgraph root file to extract the kinetic energy fraction and pt of the outgoing electron in each event. Loads the two numbers from every event into a map of vectors of pairs (mgdata). Map is keyed by energy, vector pairs are energy fraction + pt. Also creates an list of energies and placeholders (energies), so that different energies can be looped separately.
0106 {
0107   TFile* f = TFile::Open(mgfile.c_str());
0108   TTree* tree = (TTree*)f->Get("Events");
0109   TLorentzVector* mvec = new TLorentzVector();
0110   TLorentzVector* avec = new TLorentzVector();
0111   TLorentzVector* nvec = new TLorentzVector();
0112   tree->SetBranchAddress("IncidentParticle", &mvec);
0113   tree->SetBranchAddress("APrime", &avec);
0114   tree->SetBranchAddress("Nucleus", &nvec);
0115   int entries = tree->GetEntries();
0116   int start = int(G4UniformRand() * entries);
0117   for (int i = start; i < int(start + entries / 10.); i++) {
0118     if (i < entries) {
0119       tree->GetEntry(i);
0120     } else {
0121       tree->GetEntry(i - entries);
0122     }
0123     frame evnt;
0124     evnt.fEl = (TLorentzVector*)mvec->Clone();
0125     evnt.cm = (TLorentzVector*)avec->Clone();
0126     *evnt.cm = *evnt.cm + *evnt.fEl;
0127     TLorentzVector* ebeam = (TLorentzVector*)nvec->Clone();
0128     *ebeam = *ebeam + *evnt.cm;
0129     evnt.E = round(ebeam->Z() * 10.) / 10.;
0130     if (mgdata.count(evnt.E) == 0) {
0131       mgdata[evnt.E];
0132     }
0133     mgdata[evnt.E].push_back(evnt);
0134   }
0135   f->Close();
0136 }
0137
0138 void G4muDarkBremsstrahlungModel::MakePlaceholders() {
0139   //Need to do this to set up random sampling of mg distributions
0140   for (const auto& iter : mgdata) {
0141     energies.push_back(std::make_pair(iter.first, iter.second.size()));
0142   }
0143
0144   for (uint64_t i = 0; i < energies.size(); i++) {
0145     energies[i].second = int(G4UniformRand() * mgdata[energies[i].first].size());
0146   }
0147 }
0148
0149 frame G4muDarkBremsstrahlungModel::GetMadgraphData(double E0)
0150 //Gets the energy fraction and Pt from the imported LHE data. E0 should be in GeV, returns the total energy and Pt in GeV. Scales from the closest imported beam energy above the given value (scales down to avoid biasing issues).
0151 {
0152   frame cmdata;
0153   uint64_t i = 0;
0154   bool pass = false;
0155   //G4double Mel = 5.1E-04;
0156   //Cycle through imported beam energies until the closest one above is found, or the max is reached.
0157   while (!pass) {
0158     i++;
0159     if ((E0 <= energies[i].first) || (i >= energies.size())) {
0160       pass = true;
0161     }
0162   }
0163
0164   if (i == energies.size()) {
0165     i--;
0166   }  //Decrement if the loop hit the maximum, to prevent a segfault.
0167   //energies is a vector of pairs. energies[i].first holds the imported beam energy,
0168   //energies[i].second holds the place as we loop through different energy sampling files.
0169   //Need to loop around if we hit the end, when the size of mgdata is smaller than our index
0170   //on energies[i].second.
0171   if (energies[i].second >= double(mgdata[energies[i].first].size())) {
0172     energies[i].second = 0;
0173   }
0174
0175   //Get the lorentz vectors from the index given by the placeholder.
0176   cmdata = mgdata[energies[i].first].at(energies[i].second);
0177
0178   //Increment the placeholder.
0179   energies[i].second = energies[i].second + 1;
0180
0181   return cmdata;
0182 }
0183
0184 G4double G4muDarkBremsstrahlungModel::DsigmaDx(double x, void* pp) {
0185   ParamsForChi* params = (ParamsForChi*)pp;
0186
0187   G4double beta = sqrt(1 - (params->MMA) * (params->MMA) / (params->EE0) / (params->EE0));
0188   G4double num = 1. - x + x * x / 3.;
0189   G4double denom = (params->MMA) * (params->MMA) * (1. - x) / x + (params->MMu) * (params->MMu) * x;
0190   G4double DsDx = beta * num / denom;
0191
0192   return DsDx;
0193 }
0194
0195 G4double G4muDarkBremsstrahlungModel::chi(double t, void* pp) {
0196   ParamsForChi* params = (ParamsForChi*)pp;
0197   /* Reminder II:
0198  * params->AA;
0199  * params->ZZ;
0200  * params->MMA;
0201  * params->MMu;
0202  * params->EE0;
0203  */
0204   G4double Mel = 5.1E-04;
0205   G4double MUp = 2.79;
0206   G4double Mpr = 0.938;
0207   G4double d = 0.164 / pow((params->AA), 2. / 3.);
0208   G4double ap = 773.0 / Mel / pow((params->ZZ), 2. / 3.);
0209   G4double a = 111.0 / Mel / pow((params->ZZ), 1. / 3.);
0210   G4double G2el = (params->ZZ) * (params->ZZ) * a * a * a * a * t * t / (1.0 + a * a * t) / (1.0 + a * a * t) /
0211                   (1.0 + t / d) / (1.0 + t / d);
0212   G4double G2in = (params->ZZ) * ap * ap * ap * ap * t * t / (1.0 + ap * ap * t) / (1.0 + ap * ap * t) /
0213                   (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) /
0214                   (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) *
0215                   (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr) * (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr);
0216   G4double G2 = G2el + G2in;
0217   G4double ttmin = (params->MMA) * (params->MMA) * (params->MMA) * (params->MMA) / 4.0 / (params->EE0) / (params->EE0);
0218   //   G4double ttmin = lowerLimit(x,theta,p);
0219   G4double Under = G2 * (t - ttmin) / t / t;
0220
0221   return Under;
0222 }
0223
0224 G4double G4muDarkBremsstrahlungModel::ComputeCrossSectionPerAtom(
0225     const G4ParticleDefinition*, G4double E0, G4double Z, G4double A, G4double cut, G4double)
0226 // Calculates the cross section per atom in GEANT4 internal units. Uses WW approximation to find the total cross section, performing numerical integrals over x and theta.
0227 {
0228   G4double cross = 0.0;
0229   if (E0 < CLHEP::keV || E0 < cut) {
0230     return cross;
0231   }
0232
0233   E0 = E0 / CLHEP::GeV;  //Change energy to GeV.
0234   if (E0 < 2. * MA)
0235     return 0.;
0236
0237   //begin: chi-formfactor calculation
0238   gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
0239
0240   G4double result, error;
0241   G4double tmin = MA * MA * MA * MA / (4. * E0 * E0);
0242   G4double tmax = MA * MA;
0243
0244   gsl_function F;
0245   ParamsForChi alpha = {1.0, 1.0, 1.0, 1.0, 1.0};
0246   F.function = &G4muDarkBremsstrahlungModel::chi;
0247   F.params = &alpha;
0248   alpha.AA = A;
0249   alpha.ZZ = Z;
0250   alpha.MMA = MA;
0251   alpha.MMu = muonMass;
0252   alpha.EE0 = E0;
0253
0254   //Integrate over chi.
0255   gsl_integration_qags(&F, tmin, tmax, 0, 1e-7, 1000, w, &result, &error);
0256
0257   G4double ChiRes = result;
0258   gsl_integration_workspace_free(w);
0259
0260   //Integrate over x. Can use log approximation instead, which falls off at high A' mass.
0261   gsl_integration_workspace* dxspace = gsl_integration_workspace_alloc(1000);
0262   gsl_function G;
0263   G.function = &DsigmaDx;
0264   G.params = &alpha;
0265   G4double xmin = 0;
0266   G4double xmax = muonMass > MA ? 1. - muonMass / E0 : 1. - MA / E0;
0267   G4double res, err;
0268
0269   gsl_integration_qags(&G, xmin, xmax, 0, 1e-7, 1000, dxspace, &res, &err);
0270
0271   G4double DsDx = res;
0272   gsl_integration_workspace_free(dxspace);
0273
0274   G4double GeVtoPb = 3.894E08;
0275   G4double alphaEW = 1.0 / 137.0;
0276   G4double epsilBench = 1;
0277
0278   cross = GeVtoPb * 4. * alphaEW * alphaEW * alphaEW * epsilBench * epsilBench * ChiRes * DsDx * CLHEP::picobarn;
0279   if (cross < 0.) {
0280     cross = 0.;
0281   }
0282
0283   cross = cross * cxBias;
0284   return cross;
0285 }
0286
0287 G4DataVector* G4muDarkBremsstrahlungModel::ComputePartialSumSigma(const G4Material* material,
0288                                                                   G4double kineticEnergy,
0289                                                                   G4double cut)
0290
0291 // Build the table of cross section per element.
0292 //The table is built for MATERIALS.
0293 // This table is used by DoIt to select randomly an element in the material.
0294 {
0295   G4int nElements = material->GetNumberOfElements();
0296   const G4ElementVector* theElementVector = material->GetElementVector();
0297   const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
0298   G4DataVector* dv = new G4DataVector();
0299   G4double cross = 0.0;
0300
0301   for (G4int i = 0; i < nElements; i++) {
0302     cross += theAtomNumDensityVector[i] *
0303              ComputeCrossSectionPerAtom(
0304                  particle, kineticEnergy, (*theElementVector)[i]->GetZ(), (*theElementVector)[i]->GetN(), cut);
0305     dv->push_back(cross);
0306   }
0307   return dv;
0308 }
0309
0310 void G4muDarkBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
0311                                                     const G4MaterialCutsCouple* couple,
0312                                                     const G4DynamicParticle* dp,
0313                                                     G4double tmin,
0314                                                     G4double maxEnergy)
0315 //Simulates the emission of a dark photon + electron. Gets an energy fraction and Pt from madgraph files. Scales the energy so that the fraction of kinectic energy is constant, keeps the Pt constant. If the Pt is larger than the new energy, that event is skipped, and a new one is taken from the file.
0316 {
0317   //Deactivate the process after one dark brem. Needs to be reactivated in the end of event action. If this is in the stepping action instead, more than one brem can occur within each step.
0318   G4bool state = false;
0319   G4String pname = "muDBrem";
0320   G4ProcessTable* ptable = G4ProcessTable::GetProcessTable();
0321   ptable->SetProcessActivation(pname, state);
0322   if (mg_loaded == false) {
0323     LoadMG();
0324     MakePlaceholders();  //Setup the placeholder offsets for getting data.
0325     mg_loaded = true;
0326   }
0327   G4double E0 = dp->GetTotalEnergy();
0328   G4double tmax = min(maxEnergy, E0);
0329   if (tmin >= tmax) {
0330     return;
0331   }  // limits of the energy sampling
0332
0333   E0 = E0 / CLHEP::GeV;  //Convert the energy to GeV, the units used in the sampling files.
0334
0335   frame data = GetMadgraphData(E0);
0336   double EAcc, Pt, P, PhiAcc;
0337   if (method == "forward_only") {
0338     EAcc = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
0339     EAcc = muonMass + EAcc;
0340     Pt = data.fEl->Pt();
0341     P = sqrt(EAcc * EAcc - muonMass * muonMass);
0342     PhiAcc = data.fEl->Phi();
0343     int i = 0;
0344     while (Pt * Pt + muonMass * muonMass > EAcc * EAcc)  //Skip events until the Pt is less than the energy.
0345     {
0346       i++;
0347       data = GetMadgraphData(E0);
0348       EAcc = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
0349       EAcc = muonMass + EAcc;
0350       Pt = data.fEl->Pt();
0351       P = sqrt(EAcc * EAcc - muonMass * muonMass);
0352       PhiAcc = data.fEl->Phi();
0353
0354       if (i > 10000) {
0355         break;
0356       }
0357     }
0358   } else if (method == "cm_scaling") {
0359     TLorentzVector* el = new TLorentzVector(data.fEl->X(), data.fEl->Y(), data.fEl->Z(), data.fEl->E());
0360     double ediff = data.E - E0;
0361     TLorentzVector* newcm = new TLorentzVector(data.cm->X(), data.cm->Y(), data.cm->Z() - ediff, data.cm->E() - ediff);
0362     el->Boost(-1. * data.cm->BoostVector());
0363     el->Boost(newcm->BoostVector());
0364     double newE = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
0365     el->SetE(newE);
0366     EAcc = el->E();
0367     Pt = el->Pt();
0368     P = el->P();
0369     PhiAcc = el->Phi();
0370   } else {
0371     EAcc = E0;
0372     P = dp->GetTotalMomentum();
0373     Pt = sqrt(dp->Get4Momentum().px() * dp->Get4Momentum().px() + dp->Get4Momentum().py() * dp->Get4Momentum().py());
0374     PhiAcc = dp->Get4Momentum().phi();
0375   }
0376
0377   EAcc = EAcc * CLHEP::GeV;  //Change the energy back to MeV, the internal GEANT unit.
0378
0379   G4double muon_mass_mev = muonMass * CLHEP::GeV;
0380   G4double momentum = sqrt(EAcc * EAcc - muon_mass_mev * muon_mass_mev);  //Muon momentum in MeV.
0381   G4ThreeVector newDirection;
0382   double ThetaAcc = std::asin(Pt / P);
0383   newDirection.set(std::sin(ThetaAcc) * std::cos(PhiAcc), std::sin(ThetaAcc) * std::sin(PhiAcc), std::cos(ThetaAcc));
0384   newDirection.rotateUz(dp->GetMomentumDirection());
0385   newDirection.setMag(momentum);
0386   // create g4dynamicparticle object for the dark photon.
0387   G4ThreeVector direction = (dp->GetMomentum() - newDirection);
0388   G4DynamicParticle* dphoton = new G4DynamicParticle(theAPrime, direction);
0389   vdp->push_back(dphoton);
0390
0391   // energy of primary
0392   G4double finalKE = EAcc - muon_mass_mev;
0393
0394   // stop tracking and create new secondary instead of primary
0395   bool makeSecondary = true;
0396   if (makeSecondary) {
0397     fParticleChange->ProposeTrackStatus(fStopAndKill);
0398     fParticleChange->SetProposedKineticEnergy(0.0);
0399     G4DynamicParticle* mu =
0400         new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle), newDirection.unit(), finalKE);
0401     vdp->push_back(mu);
0402     // continue tracking
0403   } else {
0404     fParticleChange->SetProposedMomentumDirection(newDirection.unit());
0405     fParticleChange->SetProposedKineticEnergy(finalKE);
0406   }
0407 }
0408
0409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0410
0411 const G4Element* G4muDarkBremsstrahlungModel::SelectRandomAtom(const G4MaterialCutsCouple* couple) {
0412   // select randomly 1 element within the material
0413
0414   const G4Material* material = couple->GetMaterial();
0415   G4int nElements = material->GetNumberOfElements();
0416   const G4ElementVector* theElementVector = material->GetElementVector();
0417
0418   const G4Element* elm = nullptr;
0419
0420   if (1 < nElements) {
0421     --nElements;
0422     G4DataVector* dv = partialSumSigma[couple->GetIndex()];
0423     G4double rval = G4UniformRand() * ((*dv)[nElements]);
0424
0425     elm = (*theElementVector)[nElements];
0426     for (G4int i = 0; i < nElements; ++i) {
0427       if (rval <= (*dv)[i]) {
0428         elm = (*theElementVector)[i];
0429         break;
0430       }
0431     }
0432   } else {
0433     elm = (*theElementVector)[0];
0434   }
0435
0436   SetCurrentElement(elm);
0437   return elm;
0438 }
0439
0440 void G4muDarkBremsstrahlungModel::SetMethod(std::string method_in) {
0441   method = method_in;
0442   return;
0443 }