11 #include <generators/phokhara/Phokhara.h>
12 #include <framework/logging/Logger.h>
14 #include <TDatabasePDG.h>
15 #include <TLorentzVector.h>
34 double phokhara_rndm_(
int*)
36 double r = gRandom->Rndm();
45 void phokhara_rndmarray_(
double* drvec,
int* lengt)
47 for (
int i = 0; i < *lengt; ++i) {
48 drvec[i] = gRandom->Rndm();
53 void phokhara_(
int* mode,
double* xpar,
int* npar);
57 void phokhara_setparamfile_(
const char* paramfilename,
size_t* length);
60 void phokhara_error_trials_(
double* trials)
62 B2ERROR(
"PHOKHARA: Event rejected! Number of maximal trials (" << *trials <<
" << exceeded");
66 void phokhara_warning_weight_(
int* photmult,
double* weight,
double* max)
68 B2WARNING(
"PHOKHARA: Event weight (" << *weight <<
") exceeds limit (" << *max <<
"), photon multiplicity: " << *photmult);
72 void phokhara_warning_negweight_(
int* photmult,
double* weight)
74 B2WARNING(
"PHOKHARA: Event weight (" << *weight <<
") is negative, photon multiplicity: " << *photmult);
81 for (
int i = 0; i < 100; ++i) {
94 void Phokhara::setDefaultSettings()
124 m_replaceMuonsByVirtualPhoton =
false;
134 m_pionstructure = -1;
138 m_ScatteringAngleRangePhoton = make_pair(0.0, 180.0);
139 m_ScatteringAngleRangeFinalStates = make_pair(0.0, 180.0);
140 m_MinInvMassHadronsGamma = -1;
141 m_MinInvMassHadrons = 0.;
142 m_ForceMinInvMassHadronsCut =
false;
143 m_MaxInvMassHadrons = 0.;
144 m_MinEnergyGamma = 0.;
147 m_conversionFactor = 0.;
157 void Phokhara::init(
const std::string& paramFile)
159 B2INFO(
"Phokhara::init, using paramater file: " << paramFile);
161 if (paramFile.empty()) B2FATAL(
"Phokhara: The specified param file is empty!");
162 size_t fileLength = paramFile.size();
163 phokhara_setparamfile_(paramFile.c_str(), &fileLength);
173 void Phokhara::generateEvent(
MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
178 if (m_ForceMinInvMassHadronsCut) {
180 phokhara_(&mode, m_xpar, m_npar);
181 double partMom[4] = {0, 0, 0, 0};
182 for (
int iPart = 0; iPart < momset_.bnhad; ++iPart) {
183 for (
int j = 0; j < 4; ++j)
184 partMom[j] += momset_.bp2[j][iPart];
186 double m2 = partMom[0] * partMom[0] - partMom[1] * partMom[1]
187 - partMom[2] * partMom[2] - partMom[3] * partMom[3];
188 if (m2 > m_MinInvMassHadrons)
192 phokhara_(&mode, m_xpar, m_npar);
195 double eMom[4] = {momset_.bp1[1], momset_.bp1[2], momset_.bp1[3], momset_.bp1[0]};
196 double pMom[4] = {momset_.bq1[1], momset_.bq1[2], momset_.bq1[3], momset_.bq1[0]};
198 storeParticle(mcGraph, eMom, 11, vertex, boost,
false,
true);
199 storeParticle(mcGraph, pMom, -11, vertex, boost,
false,
true);
202 for (
int iPhot = 0; iPhot < momset_.bnphot; ++iPhot) {
203 double photMom[4] = {momset_.bphot[1][iPhot], momset_.bphot[2][iPhot], momset_.bphot[3][iPhot], momset_.bphot[0][iPhot]};
204 storeParticle(mcGraph, photMom, 22, vertex, boost,
false,
false);
208 if ((m_finalState == 0) && m_replaceMuonsByVirtualPhoton) {
209 if (momset_.bnhad != 2)
210 B2FATAL(
"Number of particles generated by PHOKHARA does not match the "
211 "requested final state (mu+ mu-).");
212 double partMom[4] = {momset_.bp2[1][0] + momset_.bp2[1][1],
213 momset_.bp2[2][0] + momset_.bp2[2][1],
214 momset_.bp2[3][0] + momset_.bp2[3][1],
215 momset_.bp2[0][0] + momset_.bp2[0][1]
217 storeParticle(mcGraph, partMom, 10022, vertex, boost,
false,
false);
219 for (
int iPart = 0; iPart < momset_.bnhad; ++iPart) {
220 double partMom[4] = {momset_.bp2[1][iPart], momset_.bp2[2][iPart], momset_.bp2[3][iPart], momset_.bp2[0][iPart]};
222 storeParticle(mcGraph, partMom, momset_.bp2[4][iPart], vertex, boost,
false,
false);
227 if (m_finalState == 9) {
228 int id = 2 + momset_.bnphot;
238 daughter1 = &mcGraph[
id + 4];
240 daughter2 = &mcGraph[
id + 5];
243 lambdabar->
removeStatus(MCParticle::c_StableInGenerator);
250 void Phokhara::term()
254 phokhara_(&mode, m_xpar, m_npar);
277 void Phokhara::applySettings()
283 m_npar[1] = m_nMaxTrials;
284 m_npar[2] = m_nSearchMax;
285 m_npar[20] = m_finalState;
289 m_npar[33] = m_NLOIFI;
290 m_npar[34] = m_alpha;
291 m_npar[35] = m_pionff;
292 m_npar[36] = m_pionstructure;
293 m_npar[37] = m_kaonff;
294 m_npar[38] = m_narres;
295 m_npar[39] = m_protonff;
300 m_xpar[0] = m_cmsEnergy;
301 m_xpar[11] = m_epsilon;
303 m_xpar[15] = m_MinInvMassHadronsGamma;
304 m_xpar[16] = m_MinInvMassHadrons;
305 m_xpar[17] = m_MaxInvMassHadrons;
306 m_xpar[18] = m_MinEnergyGamma;
308 m_xpar[20] = m_ScatteringAngleRangePhoton.first;
309 m_xpar[21] = m_ScatteringAngleRangePhoton.second;
310 m_xpar[22] = m_ScatteringAngleRangeFinalStates.first;
311 m_xpar[23] = m_ScatteringAngleRangeFinalStates.second;
314 phokhara_(&mode, m_xpar, m_npar);
318 void Phokhara::storeParticle(
MCParticleGraph& mcGraph,
const double* mom,
int pdg, TVector3 vertex, TLorentzRotation boost,
319 bool isVirtual,
bool isInitial)
326 part.addStatus(MCParticle::c_PrimaryParticle);
329 part.addStatus(MCParticle::c_StableInGenerator);
332 part.addStatus(MCParticle::c_IsVirtual);
333 }
else if (isInitial) {
334 part.addStatus(MCParticle::c_Initial);
339 part.addStatus(MCParticle::c_IsISRPhoton);
340 part.addStatus(MCParticle::c_IsFSRPhoton);
344 part.setFirstDaughter(0);
345 part.setLastDaughter(0);
346 part.setMomentum(TVector3(mom[0], mom[1], mom[2]));
348 if ((m_finalState == 0) && m_replaceMuonsByVirtualPhoton && (pdg == 10022))
349 part.setMass(sqrt(mom[3] * mom[3] - mom[0] * mom[0] - mom[1] * mom[1] -
352 part.setMass(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
353 part.setEnergy(mom[3]);
356 TLorentzVector p4 = part.get4Vector();
362 TVector3 v3 = part.getProductionVertex();
364 part.setProductionVertex(v3);
365 part.setValidVertex(
true);