#include <G4muDarkBremsstrahlungModel.h>

Inheritance diagram for G4muDarkBremsstrahlungModel:

Public Member Functions
G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double tkin, G4double Z, G4double, G4double cut, G4double maxE=DBL_MAX) override

G4DataVector *	ComputePartialSumSigma (const G4Material *material, G4double tkin, G4double cut)

	G4muDarkBremsstrahlungModel (const G4String &scalefile, const G4double biasFactor, const G4ParticleDefinition *p=nullptr, const G4String &nam="eDBrem")

	G4muDarkBremsstrahlungModel (const G4muDarkBremsstrahlungModel &)=delete

frame	GetMadgraphData (double E0)

void	Initialise (const G4ParticleDefinition *, const G4DataVector &) override

void	LoadMG ()

void	MakePlaceholders ()

G4muDarkBremsstrahlungModel &	operator= (const G4muDarkBremsstrahlungModel &right)=delete

void	SampleSecondaries (std::vector< G4DynamicParticle > , const G4MaterialCutsCouple , const G4DynamicParticle , G4double tmin, G4double maxEnergy) override

void	SetMethod (std::string)

	~G4muDarkBremsstrahlungModel () override

Protected Member Functions
const G4Element *	SelectRandomAtom (const G4MaterialCutsCouple *couple)

Protected Attributes
const G4double	cxBias

G4ParticleChangeForLoss *	fParticleChange

G4bool	isMuon

G4double	MA

const G4String &	mgfile

G4double	muonMass

const G4ParticleDefinition *	particle

G4ParticleDefinition *	theAPrime

Private Member Functions
void	SetParticle (const G4ParticleDefinition *p)

Static Private Member Functions
static G4double	chi (double t, void *pp)

static G4double	DsigmaDx (double x, void *pp)

Private Attributes
std::vector< std::pair< double, int > >	energies

G4double	highKinEnergy

G4bool	isInitialised

G4double	lowKinEnergy

std::string	method

G4bool	mg_loaded

std::map< double, std::vector< frame > >	mgdata

std::vector< G4DataVector * >	partialSumSigma

G4double	probsup

Detailed Description

Definition at line 36 of file G4muDarkBremsstrahlungModel.h.

Constructor & Destructor Documentation

◆ G4muDarkBremsstrahlungModel() [1/2]

G4muDarkBremsstrahlungModel::G4muDarkBremsstrahlungModel	(	const G4String &	scalefile,
		const G4double	biasFactor,
		const G4ParticleDefinition *	p = `nullptr`,
		const G4String &	nam = `"eDBrem"`
	)

Definition at line 21 of file G4muDarkBremsstrahlungModel.cc.

References G4APrime::APrime(), fParticleChange, highKinEnergy, lowKinEnergy, MA, muonMass, AlCaHLTBitMon_ParallelJobs::p, SetParticle(), and theAPrime.

     : G4VEmModel(nam),
       mgfile(scalefile),
       cxBias(biasFactor),
       particle(nullptr),
       isMuon(true),
       probsup(1.0),
       isInitialised(false),
       method("forward_only"),
       mg_loaded(false) {
   if (p) {
     SetParticle(p);
   }  //Verify that the particle is a muon.
   theAPrime = G4APrime::APrime();
   MA = G4APrime::APrime()->GetPDGMass() / CLHEP::GeV;                        //Get the A' mass.
   muonMass = G4MuonMinus::MuonMinusDefinition()->GetPDGMass() / CLHEP::GeV;  //Get the muon mass
   highKinEnergy = HighEnergyLimit();
   lowKinEnergy = LowEnergyLimit();
   fParticleChange = nullptr;
 }

◆ ~G4muDarkBremsstrahlungModel()

G4muDarkBremsstrahlungModel::~G4muDarkBremsstrahlungModel ( )

override

Definition at line 45 of file G4muDarkBremsstrahlungModel.cc.

References mps_fire::i, dqmiodumpmetadata::n, and partialSumSigma.

                                                           {
   size_t n = partialSumSigma.size();
   if (n > 0) {
     for (size_t i = 0; i < n; i++) {
       delete partialSumSigma[i];
     }
   }
 }

◆ G4muDarkBremsstrahlungModel() [2/2]

G4muDarkBremsstrahlungModel::G4muDarkBremsstrahlungModel ( const G4muDarkBremsstrahlungModel & )

delete

Member Function Documentation

◆ chi()

G4double G4muDarkBremsstrahlungModel::chi	(	double	t,
		void *	pp
	)

staticprivate

Definition at line 197 of file G4muDarkBremsstrahlungModel.cc.

References a, ztail::d, submitPVValidationJobs::params, funct::pow(), createTree::pp, and submitPVValidationJobs::t.

Referenced by ComputeCrossSectionPerAtom().

                                                             {
   ParamsForChi* params = (ParamsForChi*)pp;
   /* Reminder II:
  * params->AA;
  * params->ZZ;
  * params->MMA;
  * params->MMu;
  * params->EE0;
  */
   G4double Mel = 5.1E-04;
   G4double MUp = 2.79;
   G4double Mpr = 0.938;
   G4double d = 0.164 / pow((params->AA), 2. / 3.);
   G4double ap = 773.0 / Mel / pow((params->ZZ), 2. / 3.);
   G4double a = 111.0 / Mel / pow((params->ZZ), 1. / 3.);
   G4double G2el = (params->ZZ) * (params->ZZ) * a * a * a * a * t * t / (1.0 + a * a * t) / (1.0 + a * a * t) /
                   (1.0 + t / d) / (1.0 + t / d);
   G4double G2in = (params->ZZ) * ap * ap * ap * ap * t * t / (1.0 + ap * ap * t) / (1.0 + ap * ap * t) /
                   (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) /
                   (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) *
                   (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr) * (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr);
   G4double G2 = G2el + G2in;
   G4double ttmin = (params->MMA) * (params->MMA) * (params->MMA) * (params->MMA) / 4.0 / (params->EE0) / (params->EE0);
   //   G4double ttmin = lowerLimit(x,theta,p);
   G4double Under = G2 * (t - ttmin) / t / t;
 
   return Under;
 }

◆ ComputeCrossSectionPerAtom()

G4double G4muDarkBremsstrahlungModel::ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	,
		G4double	tkin,
		G4double	Z,
		G4double	A,
		G4double	cut,
		G4double	maxE = `DBL_MAX`
	)

override

Definition at line 226 of file G4muDarkBremsstrahlungModel.cc.

References A, alpha, chi(), cross(), PA_MinBiasSkim_cff::cut, cxBias, DsigmaDx(), MillePedeFileConverter_cfg::e, submitPVResolutionJobs::err, relativeConstraints::error, F(), cmssw_cycle_finder::G, MA, muonMass, mps_fire::result, tmax, muonTiming_cfi::tmin, w(), TrackerOfflineValidation_Dqm_cff::xmax, TrackerOfflineValidation_Dqm_cff::xmin, and BeamSpotPI::Z.

Referenced by ComputePartialSumSigma().

 {
   G4double cross = 0.0;
   if (E0 < keV || E0 < cut) {
     return cross;
   }
 
   E0 = E0 / CLHEP::GeV;  //Change energy to GeV.
   if (E0 < 2. * MA)
     return 0.;
 
   //begin: chi-formfactor calculation
   gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
 
   G4double result, error;
   G4double tmin = MA * MA * MA * MA / (4. * E0 * E0);
   G4double tmax = MA * MA;
 
   gsl_function F;
   ParamsForChi alpha = {1.0, 1.0, 1.0, 1.0, 1.0};
   F.function = &G4muDarkBremsstrahlungModel::chi;
   F.params = &alpha;
   alpha.AA = A;
   alpha.ZZ = Z;
   alpha.MMA = MA;
   alpha.MMu = muonMass;
   alpha.EE0 = E0;
 
   //Integrate over chi.
   gsl_integration_qags(&F, tmin, tmax, 0, 1e-7, 1000, w, &result, &error);
 
   G4double ChiRes = result;
   gsl_integration_workspace_free(w);
 
   //Integrate over x. Can use log approximation instead, which falls off at high A' mass.
   gsl_integration_workspace* dxspace = gsl_integration_workspace_alloc(1000);
   gsl_function G;
   G.function = &DsigmaDx;
   G.params = &alpha;
   G4double xmin = 0;
   G4double xmax = 1;
   if ((muonMass / E0) > (MA / E0))
     xmax = 1 - muonMass / E0;
   else
     xmax = 1 - MA / E0;
   G4double res, err;
 
   gsl_integration_qags(&G, xmin, xmax, 0, 1e-7, 1000, dxspace, &res, &err);
 
   G4double DsDx = res;
   gsl_integration_workspace_free(dxspace);
 
   G4double GeVtoPb = 3.894E08;
   G4double alphaEW = 1.0 / 137.0;
   G4double epsilBench = 1;
 
   cross = GeVtoPb * 4. * alphaEW * alphaEW * alphaEW * epsilBench * epsilBench * ChiRes * DsDx * CLHEP::picobarn;
   if (cross < 0.) {
     cross = 0.;
   }
   E0 = E0 * CLHEP::GeV;
 
   cross = cross * cxBias;
   return cross;
 }

◆ ComputePartialSumSigma()

G4DataVector * G4muDarkBremsstrahlungModel::ComputePartialSumSigma	(	const G4Material *	material,
		G4double	tkin,
		G4double	cut
	)

Definition at line 294 of file G4muDarkBremsstrahlungModel.cc.

References ComputeCrossSectionPerAtom(), cross(), PA_MinBiasSkim_cff::cut, mps_fire::i, and particle.

Referenced by Initialise().

 {
   G4int nElements = material->GetNumberOfElements();
   const G4ElementVector* theElementVector = material->GetElementVector();
   const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
   G4DataVector* dv = new G4DataVector();
   G4double cross = 0.0;
 
   for (G4int i = 0; i < nElements; i++) {
     cross += theAtomNumDensityVector[i] *
              ComputeCrossSectionPerAtom(
                  particle, kineticEnergy, (*theElementVector)[i]->GetZ(), (*theElementVector)[i]->GetN(), cut);
     dv->push_back(cross);
   }
   return dv;
 }

◆ DsigmaDx()

G4double G4muDarkBremsstrahlungModel::DsigmaDx	(	double	x,
		void *	pp
	)

staticprivate

Definition at line 186 of file G4muDarkBremsstrahlungModel.cc.

References HLT_2022v12_cff::beta, makePileupJSON::denom, EgammaValidation_cff::num, submitPVValidationJobs::params, createTree::pp, mathSSE::sqrt(), and x.

Referenced by ComputeCrossSectionPerAtom().

                                                                  {
   ParamsForChi* params = (ParamsForChi*)pp;
 
   G4double beta = sqrt(1 - (params->MMA) * (params->MMA) / (params->EE0) / (params->EE0));
   G4double num = 1. - x + x * x / 3.;
   G4double denom = (params->MMA) * (params->MMA) * (1. - x) / x + (params->MMu) * (params->MMu) * x;
   G4double DsDx = beta * num / denom;
 
   return DsDx;
 }

◆ GetMadgraphData()

frame G4muDarkBremsstrahlungModel::GetMadgraphData ( double E0 )

Definition at line 149 of file G4muDarkBremsstrahlungModel.cc.

References energies, first, mps_fire::i, mgdata, edm::second(), and findQualityFiles::size.

Referenced by SampleSecondaries().

 {
   double samplingE = energies[0].first;
   frame cmdata;
   uint64_t i = 0;
   bool pass = false;
   //G4double Mel = 5.1E-04;
   //Cycle through imported beam energies until the closest one above is found, or the max is reached.
   while (!pass) {
     i++;
     samplingE = energies[i].first;
     if ((E0 <= samplingE) || (i >= energies.size())) {
       pass = true;
     }
   }
 
   if (i == energies.size()) {
     i--;
   }  //Decrement if the loop hit the maximum, to prevent a segfault.
   //energies is a vector of pairs. energies[i].first holds the imported beam energy,
   //energies[i].second holds the place as we loop through different energy sampling files.
   //Need to loop around if we hit the end, when the size of mgdata is smaller than our index
   //on energies[i].second.
   if (energies[i].second >= double(mgdata[energies[i].first].size())) {
     energies[i].second = 0;
   }
 
   //Get the lorentz vectors from the index given by the placeholder.
   cmdata = mgdata[energies[i].first].at(energies[i].second);
 
   //Increment the placeholder.
   energies[i].second = energies[i].second + 1;
 
   return cmdata;
 }

◆ Initialise()

void G4muDarkBremsstrahlungModel::Initialise	(	const G4ParticleDefinition *	p,
		const G4DataVector &	cuts
	)

override

Definition at line 64 of file G4muDarkBremsstrahlungModel.cc.

References a, ComputePartialSumSigma(), cuts, fParticleChange, highKinEnergy, mps_fire::i, cuy::ii, isInitialised, lowKinEnergy, SiStripPI::min, groupFilesInBlocks::nn, AlCaHLTBitMon_ParallelJobs::p, partialSumSigma, and SetParticle().

                                                                                                     {
   if (p) {
     SetParticle(p);
   }
 
   highKinEnergy = HighEnergyLimit();
   lowKinEnergy = LowEnergyLimit();
   const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
 
   if (theCoupleTable) {
     G4int numOfCouples = theCoupleTable->GetTableSize();
     G4int nn = partialSumSigma.size();
     G4int nc = cuts.size();
     if (nn > 0) {
       for (G4int ii = 0; ii < nn; ii++) {
         G4DataVector* a = partialSumSigma[ii];
         if (a)
           delete a;
       }
       partialSumSigma.clear();
     }
     if (numOfCouples > 0) {
       for (G4int i = 0; i < numOfCouples; i++) {
         G4double cute = DBL_MAX;
         if (i < nc)
           cute = cuts[i];
         const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
         const G4Material* material = couple->GetMaterial();
         G4DataVector* dv = ComputePartialSumSigma(material, 0.5 * highKinEnergy, std::min(cute, 0.25 * highKinEnergy));
         partialSumSigma.push_back(dv);
       }
     }
   }
 
   if (isInitialised)
     return;
   fParticleChange = GetParticleChangeForLoss();
   isInitialised = true;
 }

◆ LoadMG()

void G4muDarkBremsstrahlungModel::LoadMG ( )

Definition at line 104 of file G4muDarkBremsstrahlungModel.cc.

References frame::cm, frame::E, f, frame::fEl, mps_fire::i, createfilelist::int, mgdata, and mgfile.

Referenced by SampleSecondaries().

 {
   TFile* f = TFile::Open(mgfile.c_str());
   TTree* tree = (TTree*)f->Get("Events");
   TLorentzVector* mvec = new TLorentzVector();
   TLorentzVector* avec = new TLorentzVector();
   TLorentzVector* nvec = new TLorentzVector();
   tree->SetBranchAddress("IncidentParticle", &mvec);
   tree->SetBranchAddress("APrime", &avec);
   tree->SetBranchAddress("Nucleus", &nvec);
   int entries = tree->GetEntries();
   int start = int(G4UniformRand() * entries);
   for (int i = start; i < int(start + entries / 10.); i++) {
     if (i < entries) {
       tree->GetEntry(i);
     } else {
       tree->GetEntry(i - entries);
     }
     frame evnt;
     evnt.fEl = (TLorentzVector*)mvec->Clone();
     evnt.cm = (TLorentzVector*)avec->Clone();
     *evnt.cm = *evnt.cm + *evnt.fEl;
     TLorentzVector* ebeam = (TLorentzVector*)nvec->Clone();
     *ebeam = *ebeam + *evnt.cm;
     evnt.E = round(ebeam->Z() * 10.) / 10.;
     if (mgdata.count(evnt.E) == 0) {
       mgdata[evnt.E];
     }
     mgdata[evnt.E].push_back(evnt);
   }
   f->Close();
 }

◆ MakePlaceholders()

void G4muDarkBremsstrahlungModel::MakePlaceholders ( )

Definition at line 138 of file G4muDarkBremsstrahlungModel.cc.

References energies, first, mps_fire::i, createfilelist::int, mgdata, and findQualityFiles::size.

Referenced by SampleSecondaries().

                                                    {
   //Need to do this to set up random sampling of mg distributions
   for (const auto& iter : mgdata) {
     energies.push_back(std::make_pair(iter.first, iter.second.size()));
   }
 
   for (uint64_t i = 0; i < energies.size(); i++) {
     energies[i].second = int(G4UniformRand() * mgdata[energies[i].first].size());
   }
 }

◆ operator=()

G4muDarkBremsstrahlungModel& G4muDarkBremsstrahlungModel::operator= ( const G4muDarkBremsstrahlungModel & right )

delete

◆ SampleSecondaries()

void G4muDarkBremsstrahlungModel::SampleSecondaries	(	std::vector< G4DynamicParticle >	vdp,
		const G4MaterialCutsCouple *	couple,
		const G4DynamicParticle *	dp,
		G4double	tmin,
		G4double	maxEnergy
	)

override

Definition at line 317 of file G4muDarkBremsstrahlungModel.cc.

References funct::cos(), data, Calorimetry_cff::dp, fParticleChange, GetMadgraphData(), mps_fire::i, LoadMG(), MA, MakePlaceholders(), particleFlowClusterECALTimeSelected_cfi::maxEnergy, method, mg_loaded, SiStripPI::min, amptDefaultParameters_cff::mu, muonMass, particle, unpackData-CaloStage2::pname, funct::sin(), mathSSE::sqrt(), theAPrime, tmax, and muonTiming_cfi::tmin.

 {
   //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.
   G4bool state = false;
   G4String pname = "muDBrem";
   G4ProcessTable* ptable = G4ProcessTable::GetProcessTable();
   ptable->SetProcessActivation(pname, state);
   if (mg_loaded == false) {
     LoadMG();
     MakePlaceholders();  //Setup the placeholder offsets for getting data.
     mg_loaded = true;
   }
   G4double E0 = dp->GetTotalEnergy();
   G4double tmax = min(maxEnergy, E0);
   if (tmin >= tmax) {
     return;
   }  // limits of the energy sampling
 
   E0 = E0 / CLHEP::GeV;  //Convert the energy to GeV, the units used in the sampling files.
 
   frame data = GetMadgraphData(E0);
   double EAcc, Pt, P, PhiAcc;
   if (method == "forward_only") {
     EAcc = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
     EAcc = muonMass + EAcc;
     Pt = data.fEl->Pt();
     P = sqrt(EAcc * EAcc - muonMass * muonMass);
     PhiAcc = data.fEl->Phi();
     int i = 0;
     while (Pt * Pt + muonMass * muonMass > EAcc * EAcc)  //Skip events until the Pt is less than the energy.
     {
       i++;
       data = GetMadgraphData(E0);
       EAcc = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
       EAcc = muonMass + EAcc;
       Pt = data.fEl->Pt();
       P = sqrt(EAcc * EAcc - muonMass * muonMass);
       PhiAcc = data.fEl->Phi();
 
       if (i > 10000) {
         break;
       }
     }
   } else if (method == "cm_scaling") {
     TLorentzVector* el = new TLorentzVector(data.fEl->X(), data.fEl->Y(), data.fEl->Z(), data.fEl->E());
     double ediff = data.E - E0;
     TLorentzVector* newcm = new TLorentzVector(data.cm->X(), data.cm->Y(), data.cm->Z() - ediff, data.cm->E() - ediff);
     el->Boost(-1. * data.cm->BoostVector());
     el->Boost(newcm->BoostVector());
     double newE = (data.fEl->E() - muonMass) / (data.E - muonMass - MA) * (E0 - muonMass - MA);
     el->SetE(newE);
     EAcc = el->E();
     Pt = el->Pt();
     P = el->P();
     PhiAcc = el->Phi();
   } else {
     EAcc = E0;
     P = dp->GetTotalMomentum();
     Pt = sqrt(dp->Get4Momentum().px() * dp->Get4Momentum().px() + dp->Get4Momentum().py() * dp->Get4Momentum().py());
     PhiAcc = dp->Get4Momentum().phi();
   }
 
   EAcc = EAcc * CLHEP::GeV;  //Change the energy back to MeV, the internal GEANT unit.
 
   G4double muon_mass_mev = muonMass * CLHEP::GeV;
   G4double momentum = sqrt(EAcc * EAcc - muon_mass_mev * muon_mass_mev);  //Muon momentum in MeV.
   G4ThreeVector newDirection;
   double ThetaAcc = std::asin(Pt / P);
   newDirection.set(std::sin(ThetaAcc) * std::cos(PhiAcc), std::sin(ThetaAcc) * std::sin(PhiAcc), std::cos(ThetaAcc));
   newDirection.rotateUz(dp->GetMomentumDirection());
   newDirection.setMag(momentum);
   // create g4dynamicparticle object for the dark photon.
   G4ThreeVector direction = (dp->GetMomentum() - newDirection);
   G4DynamicParticle* dphoton = new G4DynamicParticle(theAPrime, direction);
   vdp->push_back(dphoton);
 
   // energy of primary
   G4double finalKE = EAcc - muon_mass_mev;
 
   // stop tracking and create new secondary instead of primary
   bool makeSecondary = true;
   if (makeSecondary) {
     fParticleChange->ProposeTrackStatus(fStopAndKill);
     fParticleChange->SetProposedKineticEnergy(0.0);
     G4DynamicParticle* mu =
         new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle), newDirection.unit(), finalKE);
     vdp->push_back(mu);
     // continue tracking
   } else {
     fParticleChange->SetProposedMomentumDirection(newDirection.unit());
     fParticleChange->SetProposedKineticEnergy(finalKE);
   }
 }

◆ SelectRandomAtom()

const G4Element * G4muDarkBremsstrahlungModel::SelectRandomAtom ( const G4MaterialCutsCouple * couple )

protected

Definition at line 418 of file G4muDarkBremsstrahlungModel.cc.

References mps_fire::i, and partialSumSigma.

                                                                                                  {
   // select randomly 1 element within the material
 
   const G4Material* material = couple->GetMaterial();
   G4int nElements = material->GetNumberOfElements();
   const G4ElementVector* theElementVector = material->GetElementVector();
 
   const G4Element* elm = nullptr;
 
   if (1 < nElements) {
     --nElements;
     G4DataVector* dv = partialSumSigma[couple->GetIndex()];
     G4double rval = G4UniformRand() * ((*dv)[nElements]);
 
     elm = (*theElementVector)[nElements];
     for (G4int i = 0; i < nElements; ++i) {
       if (rval <= (*dv)[i]) {
         elm = (*theElementVector)[i];
         break;
       }
     }
   } else {
     elm = (*theElementVector)[0];
   }
 
   SetCurrentElement(elm);
   return elm;
 }

◆ SetMethod()

void G4muDarkBremsstrahlungModel::SetMethod ( std::string method_in )

Definition at line 447 of file G4muDarkBremsstrahlungModel.cc.

References method.

                                                                {
   method = method_in;
   return;
 }

◆ SetParticle()

void G4muDarkBremsstrahlungModel::SetParticle ( const G4ParticleDefinition * p )

private

Definition at line 54 of file G4muDarkBremsstrahlungModel.cc.

References isMuon, AlCaHLTBitMon_ParallelJobs::p, and particle.

Referenced by G4muDarkBremsstrahlungModel(), and Initialise().

                                                                            {
   particle = p;
 
   if ((p == G4MuonPlus::MuonPlus()) || (p == G4MuonMinus::MuonMinus())) {
     isMuon = true;
   } else {
     isMuon = false;
   }
 }