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