11 #include <generators/bhwide/BHWide.h>
12 #include <framework/gearbox/Unit.h>
14 #include <TDatabasePDG.h>
15 #include <TLorentzVector.h>
33 void bhwide_(
int* mode,
double* xpar,
int* npar);
34 void glimit_(
int* number);
41 void varran_(
double* drvec,
int* lengt)
43 for (
int i = 0; i < *lengt; ++i) {
47 drvec[i] = gRandom->Rndm();
48 }
while (drvec[i] >= 1.0);
55 void ranmar_(
double* rvec)
57 *rvec = gRandom->Rndm();
64 for (
int i = 0; i < 100; ++i) {
70 m_crossSectionError = 0.;
82 void BHWide::setDefaultSettings()
85 m_zContribution =
true;
88 m_randomGenerator = RG_RANMAR;
89 m_weakCorrections =
true;
90 m_ewCorrectionLib = EC_ALIBABA;
91 m_hardBremsModel = HM_CALKUL;
92 m_photonVacPol = PP_BURKHARDT;
94 m_cmsEnergy = 10.58 * Unit::GeV;
95 m_ScatteringAngleRangePositron = make_pair(17.0, 150.0);
96 m_ScatteringAngleRangeElectron = make_pair(17.0, 150.0);
97 m_minEnergyFinalStatePos = 0.2 * Unit::GeV;
98 m_minEnergyFinalStateElc = 0.2 * Unit::GeV;
99 m_maxAcollinearity = 10.0;
100 m_infCutCMSEnergy = 1e-5;
101 m_maxRejectionWeight = 3.0;
102 m_massZ = 91.1882 * Unit::GeV;
105 m_massTop = 174.3 * Unit::GeV;
106 m_massHiggs = 115.0 * Unit::GeV;
118 void BHWide::generateEvent(
MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
122 bhwide_(&mode, m_xpar, m_npar);
125 storeParticle(mcGraph, momset_.q1, 11, vertex, boost,
true);
126 storeParticle(mcGraph, momset_.p1, -11, vertex, boost,
true);
129 storeParticle(mcGraph, momset_.q2, 11, vertex, boost);
130 storeParticle(mcGraph, momset_.p2, -11, vertex, boost);
133 for (
int iPhot = 0; iPhot < momset_.nphot; ++iPhot) {
134 double photMom[4] = {momset_.phit[0][iPhot], momset_.phit[1][iPhot], momset_.phit[2][iPhot], momset_.phit[3][iPhot]};
135 storeParticle(mcGraph, photMom, 22, vertex, boost);
143 bhwide_(&mode, m_xpar, m_npar);
146 m_crossSection = m_xpar[9];
147 m_crossSectionError = m_xpar[10];
154 void BHWide::applySettings()
160 if (m_zContribution) keyZof = 0;
162 m_npar[0] = (1000 * keyZof) + (100 * m_channel) + (10 * m_weighted) + m_randomGenerator;
163 m_npar[1] = (1000 * m_weakCorrections) + (100 * m_ewCorrectionLib) + (10 * m_hardBremsModel) + m_photonVacPol;
168 m_xpar[0] = m_cmsEnergy;
169 m_xpar[1] = m_ScatteringAngleRangePositron.first;
170 m_xpar[2] = m_ScatteringAngleRangePositron.second;
171 m_xpar[3] = m_ScatteringAngleRangeElectron.first;
172 m_xpar[4] = m_ScatteringAngleRangeElectron.second;
173 m_xpar[5] = m_minEnergyFinalStatePos;
174 m_xpar[6] = m_minEnergyFinalStateElc;
175 m_xpar[7] = m_maxAcollinearity;
176 m_xpar[8] = m_infCutCMSEnergy;
177 m_xpar[9] = m_maxRejectionWeight;
178 m_xpar[10] = m_massZ;
179 m_xpar[11] = m_widthZ;
180 m_xpar[12] = m_sinW2;
181 m_xpar[13] = m_massTop;
182 m_xpar[14] = m_massHiggs;
185 bhwide_(&mode, m_xpar, m_npar);
189 void BHWide::storeParticle(
MCParticleGraph& mcGraph,
const double* mom,
int pdg, TVector3 vertex, TLorentzRotation boost,
190 bool isVirtual,
bool isInitial)
203 part.setStatus(MCParticle::c_IsVirtual);
204 }
else if (isInitial) {
205 part.setStatus(MCParticle::c_Initial);
209 part.addStatus(MCParticle::c_PrimaryParticle);
212 part.addStatus(MCParticle::c_StableInGenerator);
215 if (pdg == 22 && !isVirtual) {
216 part.addStatus(MCParticle::c_IsISRPhoton);
217 part.addStatus(MCParticle::c_IsFSRPhoton);
221 part.setFirstDaughter(0);
222 part.setLastDaughter(0);
223 part.setMomentum(TVector3(mom[0], mom[1], mom[2]));
224 part.setMass(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
225 part.setEnergy(mom[3]);
228 TLorentzVector p4 = part.get4Vector();
229 p4.SetPz(-1.0 * p4.Pz());
235 TVector3 v3 = part.getProductionVertex();
237 part.setProductionVertex(v3);
238 part.setValidVertex(
true);