Belle II Software  release-08-01-10
Teegg.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/teegg/Teegg.h>
10 #include <framework/gearbox/Unit.h>
11 #include <framework/logging/Logger.h>
12 #include <framework/dataobjects/EventMetaData.h>
13 #include <framework/dataobjects/EventExtraInfo.h>
14 
15 #include <framework/datastore/StoreObjPtr.h>
16 
17 #include <TDatabasePDG.h>
18 #include <Math/Vector4D.h>
19 #include <TMath.h>
20 #include <TRandom3.h>
21 
22 using namespace std;
23 using namespace Belle2;
24 
25 extern "C" {
26 
27  //same convention as BHWIDE
28  extern struct {
29  double tp1[4];
30  double tq1[4];
31  double tp2[4];
32  double tq2[4];
33  double tphot[4][2];
34  int tnphot;
35  } momset_;
36 
37 // results for cross section and photon multiplicity
38  extern struct {
39  int resngen;
40  int resntrials;
41  double reseff;
42  double rescross;
43  double rescrosserr;
44  double avgq2;
45  } teeggresults_;
46 
47 // extrainfo
48  extern struct {
49  double tt;
50  double tw2;
51  double tweight;
52  double tvp2;
53  } teeggextra_;
54 
55 
57  void teegg_rndmarr_(double* drvec, const int* lengt)
58  {
59  for (int i = 0; i < *lengt; ++i) {
60  double rr = gRandom->Rndm();
61  drvec[i] = rr;
62  }
63  }
64 
65  double teegg_rndm_(int*)
66  {
67  double r = gRandom->Rndm();
68  return r;
69  }
70 
72  void teeggm_(int* mode, double* xpar, int* npar);
73 
75  void teegg_warning_generic_(const double* weight, const double* max)
76  {
77  B2WARNING("TEEGG: Maximum weight " << *max << " to small, increase fmax to at least " << *weight);
78  }
79 
81  void report_() { B2FATAL("Teegg: report_() is not implemented"); }
83  void golife_() { B2FATAL("Teegg: golife_() is not implemented"); }
85  void golint_() { B2FATAL("Teegg: golint_() is not implemented"); }
86 }
87 
88 Teegg::Teegg()
89 {
90  for (int i = 0; i < 100; ++i) {
91  m_npar[i] = 0;
92  m_xpar[i] = 0.0;
93  }
94 
95  setDefaultSettings();
96 }
97 
98 Teegg::~Teegg()
99 {
100 
101 }
102 
103 void Teegg::setDefaultSettings()
104 {
105  m_cmsEnergy = 5.28695 * 2.0 * Unit::GeV;
106 
107  m_pi = 3.1415926535897932384626433832795029;
108  m_conversionFactor = 0.389379660e6;
109  m_alphaQED0 = 1.0 / 137.0359895;
110  m_massElectron = 0.51099906 * Unit::MeV;
111 
112  m_sVACPOL = "HLMNT";
113  m_sRADCOR = "NONE";
114  m_sCONFIG = "EGAMMA";
115  m_sMATRIX = "BKM2";
116  m_sMTRXGG = "EPADC";
117 
118  m_TEVETO = 0.1;
119  m_TEMIN = 20.0;// TMath::ACos(0.750);
120  m_TGMIN = 20.0;// TMath::ACos(0.750);
121  m_TGVETO = 0.05;
122  m_EEMIN = 2.00;
123  m_EGMIN = 2.00;
124  m_PEGMIN = m_pi / 4.00;
125  m_EEVETO = 0.00;
126  m_EGVETO = 0.00;
127  m_PHVETO = m_pi / 4.00;
128  m_CUTOFF = 0.25;
129  m_EPS = 0.01;
130  m_FRAPHI = 0.00;
131  m_EPSPHI = 1.0e-4;
132  m_WGHT1M = 1.001;
133  m_WGHTMX = 1.000;
134 
135  m_RADCOR = -1;
136  m_CONFIG = -1;
137  m_MATRIX = -1;
138  m_MTRXGG = -1;
139  m_VACPOL = -1;
140  m_UNWGHT = 1;
141  m_t = 0.0;
142  m_w2 = 1.0;
143  m_weight = 1.0;
144  m_vp2 = 1.0;
145 }
146 
147 void Teegg::initExtraInfo()
148 {
149  StoreObjPtr<EventExtraInfo> eventextrainfo;
150  eventextrainfo.registerInDataStore();
151 }
152 
153 void Teegg::init()
154 {
155  applySettings();
156 }
157 
158 
159 void Teegg::generateEvent(MCParticleGraph& mcGraph, ROOT::Math::XYZVector vertex, ROOT::Math::LorentzRotation boost)
160 {
161  //Generate event
162  int mode = 1;
163  teeggm_(&mode, m_xpar, m_npar);
164 
165  //Store the final state fermions as real particle into the MCParticleGraph
166  storeParticle(mcGraph, momset_.tq2, 11, vertex, boost, false, false);
167  storeParticle(mcGraph, momset_.tp2, -11, vertex, boost, false, false);
168 
169  //Store the real photon(s) into the MCParticleGraph
170  for (int iPhot = 0; iPhot < momset_.tnphot; ++iPhot) {
171  double photMom[4] = {momset_.tphot[0][iPhot], momset_.tphot[1][iPhot], momset_.tphot[2][iPhot], momset_.tphot[3][iPhot]};
172  storeParticle(mcGraph, photMom, 22, vertex, boost, false, false);
173  }
174 
175  // Get extra information
176  m_t = teeggextra_.tt;
177  m_w2 = teeggextra_.tw2;
178  m_weight = teeggextra_.tweight;
179  m_vp2 = teeggextra_.tvp2;
180 
181  // Fill event weight
182  StoreObjPtr<EventMetaData> eventMetaDataPtr("EventMetaData", DataStore::c_Event);
183  eventMetaDataPtr->setGeneratedWeight(m_weight);
184 
185  // Fill vacuum polarization correction
186  StoreObjPtr<EventExtraInfo> eventExtraInfo;
187  if (not eventExtraInfo.isValid())
188  eventExtraInfo.create();
189  if (eventExtraInfo->hasExtraInfo("GeneratorVP2")) {
190  B2WARNING("EventExtraInfo with given name is already set! I won't set it again!");
191  } else {
192  eventExtraInfo->addExtraInfo("GeneratorVP2", m_vp2);
193  }
194 
195 }
196 
197 void Teegg::term()
198 {
199  int mode = 2;
200  teeggm_(&mode, m_xpar, m_npar);
201 
202  B2RESULT("Cross-section (nb) = " << teeggresults_.rescross / 1.e3 << " +/- " << teeggresults_.rescrosserr / 1.e3 << "");
203  B2RESULT("Events (unweighted) = " << teeggresults_.resngen);
204  B2RESULT("Trials = " << teeggresults_.resntrials);
205  B2RESULT("Efficiency = " << 100.*teeggresults_.reseff << " %");
206  B2RESULT("Average Q2 = " << teeggresults_.avgq2);
207 
208 }
209 
210 //=========================================================================
211 // Protected methods
212 //=========================================================================
213 
214 void Teegg::applySettings()
215 {
216  //--------------------
217  // Integer parameters
218  //--------------------
219  m_npar[0] = m_UNWGHT; //producing strange results...
220 
221  if (m_sRADCOR == "NONE") m_RADCOR = 3; //DEFAULT
222  else if (m_sRADCOR == "SOFT") m_RADCOR = 2;
223  else if (m_sRADCOR == "HARD") m_RADCOR = 1;
224  else m_sRADCOR = 3;
225  m_npar[1] = m_RADCOR;
226 
227  if (m_sCONFIG == "EGAMMA") m_CONFIG = 11; //DEFAULT
228  else if (m_sCONFIG == "ETRON") m_CONFIG = 12;
229  else if (m_sCONFIG == "GAMMA") m_CONFIG = 13;
230  else if (m_sCONFIG == "GAMMAE") m_CONFIG = 14;
231  else m_CONFIG = 11;
232  m_npar[2] = m_CONFIG;
233 
234  if (m_sMATRIX == "BK") m_MATRIX = 21;
235  else if (m_sMATRIX == "BKM2") m_MATRIX = 22; //DEFAULT
236  else if (m_sMATRIX == "TCHAN") m_MATRIX = 23;
237  else if (m_sMATRIX == "EPA") m_MATRIX = 24;
238  else m_MATRIX = 22;
239  m_npar[3] = m_MATRIX;
240 
241  if (m_sMTRXGG == "EPADC") m_MTRXGG = 31; //DEFAULT
242  else if (m_sMTRXGG == "BEEGG") m_MTRXGG = 32;
243  else if (m_sMTRXGG == "MEEGG") m_MTRXGG = 33;
244  else if (m_sMTRXGG == "HEEGG") m_MTRXGG = 34;
245  else m_MTRXGG = 31;
246  m_npar[4] = m_MTRXGG;
247 
248  if (m_sVACPOL == "OFF") m_VACPOL = 41;
249  else if (m_sVACPOL == "HLMNT") m_VACPOL = 42; //DEFAULT
250  else if (m_sVACPOL == "NSK") m_VACPOL = 43;
251  else m_VACPOL = 42;
252  m_npar[5] = m_VACPOL;
253 
254  //--------------------
255  // Double parameters
256  //--------------------
257  double toRad = TMath::DegToRad();
258  m_xpar[0] = m_TEVETO * toRad;
259  m_xpar[1] = m_TEMIN * toRad;
260  m_xpar[2] = m_TGMIN * toRad;
261  m_xpar[3] = m_TGVETO * toRad;
262  m_xpar[4] = m_EEMIN;
263  m_xpar[5] = m_EGMIN;
264  m_xpar[6] = m_PEGMIN * toRad;
265  m_xpar[7] = m_EEVETO;
266  m_xpar[8] = m_EGVETO;
267  m_xpar[9] = m_PHVETO * toRad;
268  m_xpar[10] = m_CUTOFF;
269  m_xpar[11] = m_EPS;
270  m_xpar[12] = m_FRAPHI;
271  m_xpar[13] = m_EPSPHI;
272  m_xpar[14] = m_WGHT1M;
273  m_xpar[15] = m_WGHTMX;
274  m_xpar[30] = m_cmsEnergy;
275 
276  int mode = -1; //use mode to control init/generation/finalize in FORTRAN code
277  teeggm_(&mode, m_xpar, m_npar);
278 }
279 
280 
281 void Teegg::storeParticle(MCParticleGraph& mcGraph, const double* mom, int pdg, ROOT::Math::XYZVector vertex,
282  ROOT::Math::LorentzRotation boost,
283  bool isVirtual, bool isInitial)
284 {
285  // Create particle
286  MCParticleGraph::GraphParticle& part = mcGraph.addParticle();
287  if (isVirtual) {
288  part.addStatus(MCParticle::c_IsVirtual);
289  } else if (isInitial) {
290  part.addStatus(MCParticle::c_Initial);
291  }
292 
293  // All particles from a generator are primary
294  part.addStatus(MCParticle::c_PrimaryParticle);
295 
296  // All particles from TEEGG are stable
297  part.addStatus(MCParticle::c_StableInGenerator);
298 
299  // All gammas from TEEGG are ISR or FSR (impossible to distinguish due to IFI)
300  if (pdg == 22) {
301  part.addStatus(MCParticleGraph::GraphParticle::c_IsISRPhoton);
302  part.addStatus(MCParticleGraph::GraphParticle::c_IsFSRPhoton);
303  }
304 
305  part.setPDG(pdg);
306  part.setFirstDaughter(0);
307  part.setLastDaughter(0);
308  part.setMomentum(ROOT::Math::XYZVector(mom[0], mom[1], mom[2]));
309  part.setMass(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
310  part.setEnergy(mom[3]);
311 
312  //boost
313  ROOT::Math::PxPyPzEVector p4 = part.get4Vector();
314  p4.SetPz(-1.0 * p4.Pz()); //TEEGG uses other direction convention
315  p4 = boost * p4;
316  part.set4Vector(p4);
317 
318  //set vertex
319  if (!isInitial) {
320  ROOT::Math::XYZVector v3 = part.getProductionVertex();
321  v3 = v3 + vertex;
322  part.setProductionVertex(v3);
323  part.setValidVertex(true);
324  }
325 }
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:75
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:37
Belle2::StoreAccessorBase::registerInDataStore
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Definition: StoreAccessorBase.h:52
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
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