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