Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-04 01:54:17

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 "G4SystemOfUnits.hh"
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   double samplingE = energies[0].first;
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     samplingE = energies[i].first;
0161     if ((E0 <= samplingE) || (i >= energies.size())) {
0162       pass = true;
0163     }
0164   }
0165 
0166   if (i == energies.size()) {
0167     i--;
0168   }  //Decrement if the loop hit the maximum, to prevent a segfault.
0169   //energies is a vector of pairs. energies[i].first holds the imported beam energy,
0170   //energies[i].second holds the place as we loop through different energy sampling files.
0171   //Need to loop around if we hit the end, when the size of mgdata is smaller than our index
0172   //on energies[i].second.
0173   if (energies[i].second >= double(mgdata[energies[i].first].size())) {
0174     energies[i].second = 0;
0175   }
0176 
0177   //Get the lorentz vectors from the index given by the placeholder.
0178   cmdata = mgdata[energies[i].first].at(energies[i].second);
0179 
0180   //Increment the placeholder.
0181   energies[i].second = energies[i].second + 1;
0182 
0183   return cmdata;
0184 }
0185 
0186 G4double G4muDarkBremsstrahlungModel::DsigmaDx(double x, void* pp) {
0187   ParamsForChi* params = (ParamsForChi*)pp;
0188 
0189   G4double beta = sqrt(1 - (params->MMA) * (params->MMA) / (params->EE0) / (params->EE0));
0190   G4double num = 1. - x + x * x / 3.;
0191   G4double denom = (params->MMA) * (params->MMA) * (1. - x) / x + (params->MMu) * (params->MMu) * x;
0192   G4double DsDx = beta * num / denom;
0193 
0194   return DsDx;
0195 }
0196 
0197 G4double G4muDarkBremsstrahlungModel::chi(double t, void* pp) {
0198   ParamsForChi* params = (ParamsForChi*)pp;
0199   /* Reminder II:
0200  * params->AA;
0201  * params->ZZ;
0202  * params->MMA;
0203  * params->MMu;
0204  * params->EE0;
0205  */
0206   G4double Mel = 5.1E-04;
0207   G4double MUp = 2.79;
0208   G4double Mpr = 0.938;
0209   G4double d = 0.164 / pow((params->AA), 2. / 3.);
0210   G4double ap = 773.0 / Mel / pow((params->ZZ), 2. / 3.);
0211   G4double a = 111.0 / Mel / pow((params->ZZ), 1. / 3.);
0212   G4double G2el = (params->ZZ) * (params->ZZ) * a * a * a * a * t * t / (1.0 + a * a * t) / (1.0 + a * a * t) /
0213                   (1.0 + t / d) / (1.0 + t / d);
0214   G4double G2in = (params->ZZ) * ap * ap * ap * ap * t * t / (1.0 + ap * ap * t) / (1.0 + ap * ap * t) /
0215                   (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) /
0216                   (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) *
0217                   (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr) * (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr);
0218   G4double G2 = G2el + G2in;
0219   G4double ttmin = (params->MMA) * (params->MMA) * (params->MMA) * (params->MMA) / 4.0 / (params->EE0) / (params->EE0);
0220   //   G4double ttmin = lowerLimit(x,theta,p);
0221   G4double Under = G2 * (t - ttmin) / t / t;
0222 
0223   return Under;
0224 }
0225 
0226 G4double G4muDarkBremsstrahlungModel::ComputeCrossSectionPerAtom(
0227     const G4ParticleDefinition*, G4double E0, G4double Z, G4double A, G4double cut, G4double)
0228 // 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.
0229 {
0230   G4double cross = 0.0;
0231   if (E0 < keV || E0 < cut) {
0232     return cross;
0233   }
0234 
0235   E0 = E0 / CLHEP::GeV;  //Change energy to GeV.
0236   if (E0 < 2. * MA)
0237     return 0.;
0238 
0239   //begin: chi-formfactor calculation
0240   gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
0241 
0242   G4double result, error;
0243   G4double tmin = MA * MA * MA * MA / (4. * E0 * E0);
0244   G4double tmax = MA * MA;
0245 
0246   gsl_function F;
0247   ParamsForChi alpha = {1.0, 1.0, 1.0, 1.0, 1.0};
0248   F.function = &G4muDarkBremsstrahlungModel::chi;
0249   F.params = &alpha;
0250   alpha.AA = A;
0251   alpha.ZZ = Z;
0252   alpha.MMA = MA;
0253   alpha.MMu = muonMass;
0254   alpha.EE0 = E0;
0255 
0256   //Integrate over chi.
0257   gsl_integration_qags(&F, tmin, tmax, 0, 1e-7, 1000, w, &result, &error);
0258 
0259   G4double ChiRes = result;
0260   gsl_integration_workspace_free(w);
0261 
0262   //Integrate over x. Can use log approximation instead, which falls off at high A' mass.
0263   gsl_integration_workspace* dxspace = gsl_integration_workspace_alloc(1000);
0264   gsl_function G;
0265   G.function = &DsigmaDx;
0266   G.params = &alpha;
0267   G4double xmin = 0;
0268   G4double xmax = 1;
0269   if ((muonMass / E0) > (MA / E0))
0270     xmax = 1 - muonMass / E0;
0271   else
0272     xmax = 1 - MA / E0;
0273   G4double res, err;
0274 
0275   gsl_integration_qags(&G, xmin, xmax, 0, 1e-7, 1000, dxspace, &res, &err);
0276 
0277   G4double DsDx = res;
0278   gsl_integration_workspace_free(dxspace);
0279 
0280   G4double GeVtoPb = 3.894E08;
0281   G4double alphaEW = 1.0 / 137.0;
0282   G4double epsilBench = 1;
0283 
0284   cross = GeVtoPb * 4. * alphaEW * alphaEW * alphaEW * epsilBench * epsilBench * ChiRes * DsDx * CLHEP::picobarn;
0285   if (cross < 0.) {
0286     cross = 0.;
0287   }
0288   E0 = E0 * CLHEP::GeV;
0289 
0290   cross = cross * cxBias;
0291   return cross;
0292 }
0293 
0294 G4DataVector* G4muDarkBremsstrahlungModel::ComputePartialSumSigma(const G4Material* material,
0295                                                                   G4double kineticEnergy,
0296                                                                   G4double cut)
0297 
0298 // Build the table of cross section per element.
0299 //The table is built for MATERIALS.
0300 // This table is used by DoIt to select randomly an element in the material.
0301 {
0302   G4int nElements = material->GetNumberOfElements();
0303   const G4ElementVector* theElementVector = material->GetElementVector();
0304   const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
0305   G4DataVector* dv = new G4DataVector();
0306   G4double cross = 0.0;
0307 
0308   for (G4int i = 0; i < nElements; i++) {
0309     cross += theAtomNumDensityVector[i] *
0310              ComputeCrossSectionPerAtom(
0311                  particle, kineticEnergy, (*theElementVector)[i]->GetZ(), (*theElementVector)[i]->GetN(), cut);
0312     dv->push_back(cross);
0313   }
0314   return dv;
0315 }
0316 
0317 void G4muDarkBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
0318                                                     const G4MaterialCutsCouple* couple,
0319                                                     const G4DynamicParticle* dp,
0320                                                     G4double tmin,
0321                                                     G4double maxEnergy)
0322 //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.
0323 {
0324   //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.
0325   G4bool state = false;
0326   G4String pname = "muDBrem";
0327   G4ProcessTable* ptable = G4ProcessTable::GetProcessTable();
0328   ptable->SetProcessActivation(pname, state);
0329   if (mg_loaded == false) {
0330     LoadMG();
0331     MakePlaceholders();  //Setup the placeholder offsets for getting data.
0332     mg_loaded = true;
0333   }
0334   G4double E0 = dp->GetTotalEnergy();
0335   G4double tmax = min(maxEnergy, E0);
0336   if (tmin >= tmax) {
0337     return;
0338   }  // limits of the energy sampling
0339 
0340   E0 = E0 / CLHEP::GeV;  //Convert the energy to GeV, the units used in the sampling files.
0341 
0342   frame data = GetMadgraphData(E0);
0343   double EAcc, Pt, P, PhiAcc;
0344   if (method == "forward_only") {
0345     EAcc = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
0346     EAcc = muonMass + EAcc;
0347     Pt = data.fEl->Pt();
0348     P = sqrt(EAcc * EAcc - muonMass * muonMass);
0349     PhiAcc = data.fEl->Phi();
0350     int i = 0;
0351     while (Pt * Pt + muonMass * muonMass > EAcc * EAcc)  //Skip events until the Pt is less than the energy.
0352     {
0353       i++;
0354       data = GetMadgraphData(E0);
0355       EAcc = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
0356       EAcc = muonMass + EAcc;
0357       Pt = data.fEl->Pt();
0358       P = sqrt(EAcc * EAcc - muonMass * muonMass);
0359       PhiAcc = data.fEl->Phi();
0360 
0361       if (i > 10000) {
0362         break;
0363       }
0364     }
0365   } else if (method == "cm_scaling") {
0366     TLorentzVector* el = new TLorentzVector(data.fEl->X(), data.fEl->Y(), data.fEl->Z(), data.fEl->E());
0367     double ediff = data.E - E0;
0368     TLorentzVector* newcm = new TLorentzVector(data.cm->X(), data.cm->Y(), data.cm->Z() - ediff, data.cm->E() - ediff);
0369     el->Boost(-1. * data.cm->BoostVector());
0370     el->Boost(newcm->BoostVector());
0371     double newE = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
0372     el->SetE(newE);
0373     EAcc = el->E();
0374     Pt = el->Pt();
0375     P = el->P();
0376     PhiAcc = el->Phi();
0377   } else {
0378     EAcc = E0;
0379     P = dp->GetTotalMomentum();
0380     Pt = sqrt(dp->Get4Momentum().px() * dp->Get4Momentum().px() + dp->Get4Momentum().py() * dp->Get4Momentum().py());
0381     PhiAcc = dp->Get4Momentum().phi();
0382   }
0383 
0384   EAcc = EAcc * CLHEP::GeV;  //Change the energy back to MeV, the internal GEANT unit.
0385 
0386   G4double muon_mass_mev = muonMass * CLHEP::GeV;
0387   G4double momentum = sqrt(EAcc * EAcc - muon_mass_mev * muon_mass_mev);  //Muon momentum in MeV.
0388   G4ThreeVector newDirection;
0389   double ThetaAcc = std::asin(Pt / P);
0390   newDirection.set(std::sin(ThetaAcc) * std::cos(PhiAcc), std::sin(ThetaAcc) * std::sin(PhiAcc), std::cos(ThetaAcc));
0391   newDirection.rotateUz(dp->GetMomentumDirection());
0392   newDirection.setMag(momentum);
0393   // create g4dynamicparticle object for the dark photon.
0394   G4ThreeVector direction = (dp->GetMomentum() - newDirection);
0395   G4DynamicParticle* dphoton = new G4DynamicParticle(theAPrime, direction);
0396   vdp->push_back(dphoton);
0397 
0398   // energy of primary
0399   G4double finalKE = EAcc - muon_mass_mev;
0400 
0401   // stop tracking and create new secondary instead of primary
0402   bool makeSecondary = true;
0403   if (makeSecondary) {
0404     fParticleChange->ProposeTrackStatus(fStopAndKill);
0405     fParticleChange->SetProposedKineticEnergy(0.0);
0406     G4DynamicParticle* mu =
0407         new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle), newDirection.unit(), finalKE);
0408     vdp->push_back(mu);
0409     // continue tracking
0410   } else {
0411     fParticleChange->SetProposedMomentumDirection(newDirection.unit());
0412     fParticleChange->SetProposedKineticEnergy(finalKE);
0413   }
0414 }
0415 
0416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0417 
0418 const G4Element* G4muDarkBremsstrahlungModel::SelectRandomAtom(const G4MaterialCutsCouple* couple) {
0419   // select randomly 1 element within the material
0420 
0421   const G4Material* material = couple->GetMaterial();
0422   G4int nElements = material->GetNumberOfElements();
0423   const G4ElementVector* theElementVector = material->GetElementVector();
0424 
0425   const G4Element* elm = nullptr;
0426 
0427   if (1 < nElements) {
0428     --nElements;
0429     G4DataVector* dv = partialSumSigma[couple->GetIndex()];
0430     G4double rval = G4UniformRand() * ((*dv)[nElements]);
0431 
0432     elm = (*theElementVector)[nElements];
0433     for (G4int i = 0; i < nElements; ++i) {
0434       if (rval <= (*dv)[i]) {
0435         elm = (*theElementVector)[i];
0436         break;
0437       }
0438     }
0439   } else {
0440     elm = (*theElementVector)[0];
0441   }
0442 
0443   SetCurrentElement(elm);
0444   return elm;
0445 }
0446 
0447 void G4muDarkBremsstrahlungModel::SetMethod(std::string method_in) {
0448   method = method_in;
0449   return;
0450 }