CustomPhysics/src/CMSmuDarkBremsstrahlungModel.cc

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