Belle II Software  release-08-01-10
Phokhara.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <generators/phokhara/Phokhara.h>
10 #include <framework/logging/Logger.h>
11 
12 #include <TDatabasePDG.h>
13 #include <Math/Vector4D.h>
14 #include <TRandom3.h>
15 
16 using namespace std;
17 using namespace Belle2;
18 
19 extern "C" {
20 
22  extern struct {
23  double bp1[4];
24  double bq1[4];
25  double bp2[6][10];
26  double bphot[4][2];
27  double wgt;
28  int bnphot;
29  int bnhad;
30  } belle2_phokhara_particles;
31 
33  double phokhara_rndm()
34  {
35  double r = gRandom->Rndm();
36  return r;
37  }
38 
44  void phokhara_rndmarray(double* drvec, const int* lengt)
45  {
46  for (int i = 0; i < *lengt; ++i) {
47  drvec[i] = gRandom->Rndm();
48  }
49  }
50 
52  void phokhara(int* mode, double* xpar, int* npar);
53 
55 // void phokhara_setinputfile_(const char* inputfilename, size_t* length);
56  void phokhara_set_parameter_file(const char* paramfilename);
57 
59  void phokhara_error_trials_(const double* trials)
60  {
61  B2ERROR("PHOKHARA: Event rejected! Number of maximal trials (" << *trials << " << exceeded");
62  }
63 
65  void phokhara_warning_weight_(const int* photmult, const double* weight, const double* max)
66  {
67  B2WARNING("PHOKHARA: Event weight (" << *weight << ") exceeds limit (" << *max << "), photon multiplicity: " << *photmult);
68  }
69 
71  void phokhara_warning_negweight_(const int* photmult, const double* weight)
72  {
73  B2WARNING("PHOKHARA: Event weight (" << *weight << ") is negative, photon multiplicity: " << *photmult);
74  }
75 
76 }
77 
78 Phokhara::Phokhara()
79 {
80  for (int i = 0; i < 100; ++i) {
81  m_npar[i] = 0;
82  m_xpar[i] = 0.0;
83  }
84 
85  setDefaultSettings();
86 }
87 
88 Phokhara::~Phokhara()
89 {
90 
91 }
92 
93 void Phokhara::setDefaultSettings()
94 {
95 // std::cout << "Phokhara::setDefaultSettings()" << std::endl;
96 
97 // these are the settings to reproduce the phokhara9.1 standalone result
98 // m_cmsEnergy = 10.580 * Unit::GeV;
99 // m_applyBoost = true;
100 // m_finalState = 1;
101 // m_nMaxTrials = 10000;
102 // m_nSearchMax = 50000;
103 // m_epsilon = 1.e-4;
104 // m_LO = 1;
105 // m_NLO = 0;
106 // m_QED = 0;
107 // m_IFSNLO = 0;
108 // m_alpha = 1;
109 // m_pionff = 0;
110 // m_pionstructure = 0;
111 // m_kaonff = 0;
112 // m_narres = 0;
113 // m_protonff = 1;
114 // m_ScatteringAngleRangePhoton = make_pair(0.0, 180.0); //in [deg]
115 // m_ScatteringAngleRangeFinalStates = make_pair(0.0, 180.0); //in [deg]
116 // m_MinInvMassHadronsGamma = 0.;
117 // m_MinInvMassHadrons = 0.2;
118 // m_MaxInvMassHadrons = 0.5;
119 // m_MinEnergyGamma = 4.0;
120 
121  m_cmsEnergy = 0.0;
122  m_finalState = -1;
123  m_replaceMuonsByVirtualPhoton = false;
124  m_nMaxTrials = -1;
125  m_nSearchMax = -1;
126  m_epsilon = 1.e-4;
127  m_weighted = 0;
128  m_LO = -1;
129  m_NLO = -1;
130  m_fullNLO = -1;
131  m_QED = -1;
132  m_IFSNLO = -1;
133  m_alpha = -1;
134  m_pionff = -1;
135  m_pionstructure = -1;
136  m_kaonff = -1;
137  m_narres = -1;
138  m_protonff = -1;
139  m_ScatteringAngleRangePhoton = make_pair(0.0, 180.0); //in [deg]
140  m_ScatteringAngleRangeFinalStates = make_pair(0.0, 180.0); //in [deg]
141  m_MinInvMassHadronsGamma = -1;
142  m_MinInvMassHadrons = 0.;
143  m_ForceMinInvMassHadronsCut = false;
144  m_MaxInvMassHadrons = 0.;
145  m_MinEnergyGamma = 0.;
146  m_chi_sw = 0;
147  m_be_r = 0;
148  m_beamres = 0.;
149 
150  m_pi = 0.;
151  m_conversionFactor = 0.;
152  m_alphaQED0 = 0.;
153  m_massElectron = 0.;
154  m_massMuon = 0.;
155  m_massW = 0.;
156  m_massZ = 0.;
157  m_widthZ = 0.;
158 
159 }
160 
161 void Phokhara::init(const std::string& paramFile)
162 {
163  B2INFO("Phokhara::init, using paramater file: " << paramFile);
164 
165  if (paramFile.empty()) B2FATAL("Phokhara: The specified param file is empty!");
166  phokhara_set_parameter_file(paramFile.c_str());
167 
168 // if (inputFile.empty()) B2FATAL("Phokhara: The specified input file is empty!")
169 // fileLength = inputFile.size();
170 // phokhara_setinputfile_(inputFile.c_str(), &fileLength);
171 
172  applySettings();
173 }
174 
175 
176 double Phokhara::generateEvent(MCParticleGraph& mcGraph, ROOT::Math::XYZVector vertex, ROOT::Math::LorentzRotation boost)
177 {
178 
179  //Generate event
180  int mode = 1;
181  if (m_ForceMinInvMassHadronsCut) {
182  while (1) {
183  phokhara(&mode, m_xpar, m_npar);
184  double partMom[4] = {0, 0, 0, 0};
185  for (int iPart = 0; iPart < belle2_phokhara_particles.bnhad; ++iPart) {
186  for (int j = 0; j < 4; ++j)
187  partMom[j] += belle2_phokhara_particles.bp2[j][iPart];
188  }
189  double m2 = partMom[0] * partMom[0] - partMom[1] * partMom[1]
190  - partMom[2] * partMom[2] - partMom[3] * partMom[3];
191  if (m2 > m_MinInvMassHadrons)
192  break;
193  }
194  } else
195  phokhara(&mode, m_xpar, m_npar);
196 
197  //Store the initial particles as virtual particles into the MCParticleGraph
198  double eMom[4] = {belle2_phokhara_particles.bp1[1], belle2_phokhara_particles.bp1[2], belle2_phokhara_particles.bp1[3], belle2_phokhara_particles.bp1[0]};
199  double pMom[4] = {belle2_phokhara_particles.bq1[1], belle2_phokhara_particles.bq1[2], belle2_phokhara_particles.bq1[3], belle2_phokhara_particles.bq1[0]};
200 
201  storeParticle(mcGraph, eMom, 11, vertex, boost, false, true);
202  storeParticle(mcGraph, pMom, -11, vertex, boost, false, true);
203 
204  //Store the real photons into the MCParticleGraph
205  for (int iPhot = 0; iPhot < belle2_phokhara_particles.bnphot; ++iPhot) {
206  double photMom[4] = {belle2_phokhara_particles.bphot[1][iPhot], belle2_phokhara_particles.bphot[2][iPhot], belle2_phokhara_particles.bphot[3][iPhot], belle2_phokhara_particles.bphot[0][iPhot]};
207  storeParticle(mcGraph, photMom, 22, vertex, boost, false, false);
208  }
209 
210  //Store the other final state particles into the MCParticleGraph
211  if ((m_finalState == 0) && m_replaceMuonsByVirtualPhoton) {
212  if (belle2_phokhara_particles.bnhad != 2)
213  B2FATAL("Number of particles generated by PHOKHARA does not match the "
214  "requested final state (mu+ mu-).");
215  double partMom[4] = {belle2_phokhara_particles.bp2[1][0] + belle2_phokhara_particles.bp2[1][1],
216  belle2_phokhara_particles.bp2[2][0] + belle2_phokhara_particles.bp2[2][1],
217  belle2_phokhara_particles.bp2[3][0] + belle2_phokhara_particles.bp2[3][1],
218  belle2_phokhara_particles.bp2[0][0] + belle2_phokhara_particles.bp2[0][1]
219  };
220  storeParticle(mcGraph, partMom, 10022, vertex, boost, false, false);
221  } else {
222  for (int iPart = 0; iPart < belle2_phokhara_particles.bnhad; ++iPart) {
223  double partMom[4] = {belle2_phokhara_particles.bp2[1][iPart], belle2_phokhara_particles.bp2[2][iPart], belle2_phokhara_particles.bp2[3][iPart], belle2_phokhara_particles.bp2[0][iPart]};
224 
225  storeParticle(mcGraph, partMom, belle2_phokhara_particles.bp2[4][iPart], vertex, boost, false, false);
226  }
227  }
228 
229  //some PHOKHARA final states contain unstable particle
230  if (m_finalState == 9) { //Lambda, Lambdabar
231  int id = 2 + belle2_phokhara_particles.bnphot;
232 
233  //get lambdabar -> p+ anti-p and lambda -> pi p
234  MCParticleGraph::GraphParticle* lambdabar = &mcGraph[id];
235  MCParticleGraph::GraphParticle* daughter1 = &mcGraph[id + 2];
236  daughter1->comesFrom(*lambdabar);
237  MCParticleGraph::GraphParticle* daughter2 = &mcGraph[id + 3];
238  daughter2->comesFrom(*lambdabar);
239 
240  MCParticleGraph::GraphParticle* lambda = &mcGraph[id + 1];
241  daughter1 = &mcGraph[id + 4];
242  daughter1->comesFrom(*lambda);
243  daughter2 = &mcGraph[id + 5];
244  daughter2->comesFrom(*lambda);
245 
246  lambdabar->removeStatus(MCParticle::c_StableInGenerator);
247  lambda->removeStatus(MCParticle::c_StableInGenerator);
248  }
249 
250  if (m_weighted)
251  return belle2_phokhara_particles.wgt;
252  return 1.0;
253 
254 }
255 
256 
257 void Phokhara::term()
258 {
259 
260  int mode = 2;
261  phokhara(&mode, m_xpar, m_npar);
262 
263 // B2INFO("> Crosssection ")
264 // B2INFO(">> xsec (weighted) = (" << bresults_.rescross << " +/- " << bresults_.rescrosserr << ") nb")
265 // for (int i = 0; i < 3; i++) {
266 // B2INFO(">>> " << i << " photon(s), xsec (weighted) = (" << bresults_.rescrossphot[i] << " +/- " << bresults_.rescrossphoterr[i] << ") nb, fraction = " << bresults_.rescrossphotfrac[i]*100.<< "%")
267 // }
268 //
269 // B2INFO(">>> xsec (unweighted) = (" << bhitnmiss_.hnmcross << " +/- " << bhitnmiss_.hnmcrosserr << ") nb")
270 // for (int i = 0; i < 3; i++) {
271 // B2INFO(">>> " << i << " photon(s), xsec (unweighted) = (" << bhitnmiss_.hnmcrossphot[i] << " +/- " << bhitnmiss_.hnmcrossphoterr[i] << ") nb, fraction.= " << bhitnmiss_.hnmcrossphotfrac[i]*100.<< "%")
272 // }
273 //
274 // B2INFO("> hit/miss efficiency")
275 // for (int i = 0; i < 3; i++) {
276 // B2INFO(">>> " << i << " photon(s),eff.= " << bhitnmiss_.hnmeff[i]*100.<< "%")
277 // }
278 }
279 
280 //=========================================================================
281 // Protected methods
282 //=========================================================================
283 
284 void Phokhara::applySettings()
285 {
286 
287  //--------------------
288  // Integer parameters
289  //--------------------
290  m_npar[1] = m_nMaxTrials;
291  m_npar[2] = m_nSearchMax;
292  m_npar[3] = m_weighted;
293  m_npar[20] = m_finalState;
294  m_npar[30] = m_LO;
295  m_npar[31] = m_NLO;
296  m_npar[32] = m_QED;
297  m_npar[33] = m_IFSNLO;
298  m_npar[34] = m_alpha;
299  m_npar[35] = m_pionff;
300  m_npar[36] = m_pionstructure;
301  m_npar[37] = m_kaonff;
302  m_npar[38] = m_narres;
303  m_npar[39] = m_protonff;
304  m_npar[40] = m_fullNLO;
305  m_npar[50] = m_chi_sw;
306  m_npar[51] = m_be_r;
307  m_npar[60] = m_weighted;
308 
309  //--------------------
310  // Double parameters
311  //--------------------
312  m_xpar[0] = m_cmsEnergy;
313  m_xpar[11] = m_epsilon;
314 
315  m_xpar[15] = m_MinInvMassHadronsGamma;
316  m_xpar[16] = m_MinInvMassHadrons;
317  m_xpar[17] = m_MaxInvMassHadrons;
318  m_xpar[18] = m_MinEnergyGamma;
319 
320  m_xpar[20] = m_ScatteringAngleRangePhoton.first;
321  m_xpar[21] = m_ScatteringAngleRangePhoton.second;
322  m_xpar[22] = m_ScatteringAngleRangeFinalStates.first;
323  m_xpar[23] = m_ScatteringAngleRangeFinalStates.second;
324 
325  m_xpar[30] = m_beamres;
326 
327  int mode = -1; //use mode to control init/generation/finalize in FORTRAN code
328  phokhara(&mode, m_xpar, m_npar);
329 }
330 
331 
332 void Phokhara::storeParticle(MCParticleGraph& mcGraph, const double* mom, int pdg, ROOT::Math::XYZVector vertex,
333  ROOT::Math::LorentzRotation boost, bool isVirtual, bool isInitial)
334 {
335 
336  // Create particle
337  MCParticleGraph::GraphParticle& part = mcGraph.addParticle();
338 
339  //all particles are primary!
340  part.addStatus(MCParticle::c_PrimaryParticle);
341 
342  //all particles are stable (if not, will be changed later)!
343  part.addStatus(MCParticle::c_StableInGenerator);
344 
345  if (isVirtual) {
346  part.addStatus(MCParticle::c_IsVirtual);
347  } else if (isInitial) {
348  part.addStatus(MCParticle::c_Initial);
349  }
350 
351  //set photon flags
352  if (pdg == 22) {
353  part.addStatus(MCParticle::c_IsISRPhoton);
354  part.addStatus(MCParticle::c_IsFSRPhoton);
355  }
356 
357  part.setPDG(pdg);
358  part.setFirstDaughter(0);
359  part.setLastDaughter(0);
360  part.setMomentum(ROOT::Math::XYZVector(mom[0], mom[1], mom[2]));
361  // part.get4Vector() uses mass, need to set invariant mass for virtual photons
362  if ((m_finalState == 0) && m_replaceMuonsByVirtualPhoton && (pdg == 10022))
363  part.setMass(sqrt(mom[3] * mom[3] - mom[0] * mom[0] - mom[1] * mom[1] -
364  mom[2] * mom[2]));
365  else
366  part.setMass(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
367  part.setEnergy(mom[3]);
368 
369  //boost
370  ROOT::Math::PxPyPzEVector p4 = part.get4Vector();
371  p4 = boost * p4;
372  part.set4Vector(p4);
373 
374  //set vertex
375  if (!isInitial) {
376  ROOT::Math::XYZVector v3 = part.getProductionVertex();
377  v3 = v3 + vertex;
378  part.setProductionVertex(v3);
379  part.setValidVertex(true);
380  }
381 }
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:75
Belle2::MCParticleGraph::GraphParticle::comesFrom
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
Definition: MCParticleGraph.h:105
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:37
Belle2::MCParticle::removeStatus
void removeStatus(unsigned short int bitmask)
Remove bitmask from current status.
Definition: MCParticle.h:360
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:294
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17