Belle II Software  release-05-01-25
BHWide.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/bhwide/BHWide.h>
12 #include <framework/gearbox/Unit.h>
13 
14 #include <TDatabasePDG.h>
15 #include <TLorentzVector.h>
16 #include <TRandom.h>
17 
18 using namespace std;
19 using namespace Belle2;
20 
21 //Declaration of the BHWide FORTRAN methods and common blocks.
22 extern "C" {
23 
24  extern struct {
25  double p1[4];
26  double q1[4];
27  double p2[4];
28  double q2[4];
29  double phit[4][100];
30  int nphot;
31  } momset_;
32 
33  void bhwide_(int* mode, double* xpar, int* npar);
34  void glimit_(int* number);
35 
41  void varran_(double* drvec, int* lengt)
42  {
43  for (int i = 0; i < *lengt; ++i) {
44  do {
45  //BHWide does not want 1 for some reason, at least they are rejected in
46  //the original VARRAN subroutine so we do it here as well
47  drvec[i] = gRandom->Rndm();
48  } while (drvec[i] >= 1.0);
49  }
50  }
55  void ranmar_(double* rvec)
56  {
57  *rvec = gRandom->Rndm();
58  }
59 }
60 
61 
62 BHWide::BHWide()
63 {
64  for (int i = 0; i < 100; ++i) {
65  m_npar[i] = 0;
66  m_xpar[i] = 0.0;
67  }
68 
69  m_crossSection = 0.;
70  m_crossSectionError = 0.;
71 
72  setDefaultSettings();
73 }
74 
75 
76 BHWide::~BHWide()
77 {
78 
79 }
80 
81 
82 void BHWide::setDefaultSettings()
83 {
84 
85  m_zContribution = true;
86  m_channel = CH_BOTH;
87  m_weighted = false;
88  m_randomGenerator = RG_RANMAR;
89  m_weakCorrections = true;
90  m_ewCorrectionLib = EC_ALIBABA;
91  m_hardBremsModel = HM_CALKUL;
92  m_photonVacPol = PP_BURKHARDT;
93 
94  m_cmsEnergy = 10.58 * Unit::GeV;
95  m_ScatteringAngleRangePositron = make_pair(17.0, 150.0); //in [deg]
96  m_ScatteringAngleRangeElectron = make_pair(17.0, 150.0); //in [deg]
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;
103  m_widthZ = 2.4952;
104  m_sinW2 = 0.22225;
105  m_massTop = 174.3 * Unit::GeV;
106  m_massHiggs = 115.0 * Unit::GeV;
107 }
108 
109 
110 void BHWide::init()
111 {
112  int number = 50000;
113  glimit_(&number);
114  applySettings();
115 }
116 
117 
118 void BHWide::generateEvent(MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
119 {
120  //Generate event
121  int mode = 0;
122  bhwide_(&mode, m_xpar, m_npar);
123 
124  //Store the initial particles as virtual particles into the MCParticleGraph
125  storeParticle(mcGraph, momset_.q1, 11, vertex, boost, true);
126  storeParticle(mcGraph, momset_.p1, -11, vertex, boost, true);
127 
128  //Store the final state positron and electron into the MCParticleGraph
129  storeParticle(mcGraph, momset_.q2, 11, vertex, boost);
130  storeParticle(mcGraph, momset_.p2, -11, vertex, boost);
131 
132  //Store the real photons into the MCParticleGraph
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);
136  }
137 }
138 
139 
140 void BHWide::term()
141 {
142  int mode = 2;
143  bhwide_(&mode, m_xpar, m_npar);
144 
145  //Get the cross section
146  m_crossSection = m_xpar[9];
147  m_crossSectionError = m_xpar[10];
148 }
149 
150 //=========================================================================
151 // Protected methods
152 //=========================================================================
153 
154 void BHWide::applySettings()
155 {
156  //--------------------
157  // Integer parameters
158  //--------------------
159  int keyZof;
160  if (m_zContribution) keyZof = 0;
161  else keyZof = 1;
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;
164 
165  //--------------------
166  // Double parameters
167  //--------------------
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;
183 
184  int mode = -1;
185  bhwide_(&mode, m_xpar, m_npar);
186 }
187 
188 
189 void BHWide::storeParticle(MCParticleGraph& mcGraph, const double* mom, int pdg, TVector3 vertex, TLorentzRotation boost,
190  bool isVirtual, bool isInitial)
191 {
192  // //Create particle
193  //MCParticleGraph::GraphParticle& part = mcGraph.addParticle();
194  //if (!isVirtual) {
195  // part.setStatus(MCParticle::c_PrimaryParticle);
196  //} else {
197  // part.setStatus(MCParticle::c_IsVirtual);
198  //}
199  //Create particle
200  // RG 6/25/14 Add new flag for ISR "c_Initial"
202  if (isVirtual) {
203  part.setStatus(MCParticle::c_IsVirtual);
204  } else if (isInitial) {
205  part.setStatus(MCParticle::c_Initial);
206  }
207 
208  // every particle from a generator is primary
209  part.addStatus(MCParticle::c_PrimaryParticle);
210 
211  // all particles produced by BHWIDE are stable
212  part.addStatus(MCParticle::c_StableInGenerator);
213 
214  // all photons from BHWIDE are ISR or FSR
215  if (pdg == 22 && !isVirtual) {
216  part.addStatus(MCParticle::c_IsISRPhoton);
217  part.addStatus(MCParticle::c_IsFSRPhoton);
218  }
219 
220  part.setPDG(pdg);
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]);
226 
227  //boost
228  TLorentzVector p4 = part.get4Vector();
229  p4.SetPz(-1.0 * p4.Pz()); //BHWIDE uses other direction convention
230  p4 = boost * p4;
231  part.set4Vector(p4);
232 
233  //set vertex
234  if (!isInitial) {
235  TVector3 v3 = part.getProductionVertex();
236  v3 = v3 + vertex;
237  part.setProductionVertex(v3);
238  part.setValidVertex(true);
239  }
240 }
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::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86