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