StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main46.cc
1 // main46.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 //
6 // Author: S.Chekanov (ANL).
7 // An example of how to write ProMC files.
8 //
9 // The ProMC library is described at http://atlaswww.hep.anl.gov/asc/promc/
10 // A makefile can be found in the ProMC package, in examples/pythia.
11 
12 #include <map>
13 #include "Pythia.h"
14 // ProMC file
15 #include "ProMCHeader.pb.h"
16 #include "ProMC.pb.h"
17 #include "ProMCBook.h"
18 
19 using namespace Pythia8;
20 
21 //--------------------------------------------------------------------------
22 
23 string getEnvVar( std::string const & key ) {
24  char * val = getenv( key.c_str() );
25  return (val == NULL) ? std::string("") : std::string(val);
26 }
27 
28 //--------------------------------------------------------------------------
29 
30 void readPDG( ProMCHeader * header ) {
31 
32  string temp_string;
33  istringstream curstring;
34  string PdgTableFilename = getEnvVar("PROMC") + "/data/particle.tbl";
35  if (PdgTableFilename.size() < 2) {
36  cout << "** ERROR: PROMC variable not set. Did you run source.sh"
37  << " **" << endl;
38  exit(1);
39  }
40 
41  ifstream file_to_read(PdgTableFilename.c_str());
42  if (!file_to_read.good()) {
43  cout << "** ERROR: PDG Table (" << PdgTableFilename
44  << ") not found! exit. **" << endl;
45  exit(1);
46  return;
47  }
48 
49  // First three lines of the file are useless.
50  getline(file_to_read,temp_string);
51  getline(file_to_read,temp_string);
52  getline(file_to_read,temp_string);
53 
54  while (getline(file_to_read,temp_string)) {
55  // Needed when using several times istringstream::str(string).
56  curstring.clear();
57  curstring.str(temp_string);
58  long int ID; std::string name; int charge; float mass; float width;
59  float lifetime;
60  // ID name chg mass total width lifetime
61  // 1 d -1 0.33000 0.00000 0.00000E+00
62  // in the table, the charge is in units of e+/3
63  // the total width is in GeV
64  // the lifetime is ctau in mm
65  curstring >> ID >> name >> charge >> mass >> width >> lifetime;
66  ProMCHeader_ParticleData* pp= header->add_particledata();
67  pp->set_id(ID);
68  pp->set_mass(mass);
69  pp->set_name(name);
70  pp->set_width(width);
71  pp->set_lifetime(lifetime);
72  cout << ID << " " << name << " " << mass << endl;
73  }
74 
75 }
76 
77 //--------------------------------------------------------------------------
78 
79 int main() {
80 
81  int Ntot = 1000; // Total number of events.
82 
83  // Generator. Process selection. LHC initialization.
84  Pythia pythia;
85  pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
86  pythia.readString("PhaseSpace:mHatMin = 80.");
87  pythia.readString("PhaseSpace:mHatMax = 120.");
88  pythia.readString("Random:setSeed = on");
89  pythia.readString("Random:seed = 0");
90  pythia.init( 2212, 2212, 14000.);
91 
92  // **************** book ProMC file **********************
93  ProMCBook* epbook = new ProMCBook("Pythia8.promc","w");
94  epbook->setDescription(Ntot,"PYTHIA8");
95 
96  // Info on incoming beams and CM energy.
97  ProMCHeader header;
98  header.set_id1( pythia.info.idA() );
99  header.set_id2( pythia.info.idB() );
100  header.set_ecm( pythia.info.eCM() );
101  header.set_s( pythia.info.s() );
102 
103  // Use the range 0.01 MeV to 20 TeV using varints (integers).
104  // With particle in GeV, we multiply it by kEV, to get 0.01 MeV = 1 unit.
105  const double kEV = 1000*100;
106  // With particle in mm, we multiply it by kL to get 0.01 mm = 1 unit.
107  const double kL = 100;
108 
109  // Set units.
110  header.set_momentumunit( (int)kEV );
111  header.set_lengthunit( (int)kL );
112 
113  // Store a map with PDG information (stored in the header).
114  readPDG( &header );
115  epbook->setHeader(header); // write header
116 
117  // Begin event loop. Generate event. Skip if error.
118  for (int n = 0; n < Ntot; n++) {
119  if (!pythia.next()) {
120  // If failure because reached end of file then exit event loop.
121  if (pythia.info.atEndOfFile()) {
122  cout << " Aborted since reached end of Les Houches Event File\n";
123  break;
124  }
125  continue;
126  }
127 
128  //************ ProMC file ***************//
129  ProMCEvent promc;
130 
131  // Fill event information.
132  ProMCEvent_Event *eve = promc.mutable_event();
133  eve->set_number(n);
134  eve->set_process_id( pythia.info.code() ); // process ID
135  eve->set_scale( pythia.info.pTHat( )); // relevant for 2 -> 2 only
136  eve->set_alpha_qed( pythia.info.alphaEM() );
137  eve->set_alpha_qcd( pythia.info.alphaS() );
138  eve->set_scale_pdf( pythia.info.QFac() );
139  eve->set_x1( pythia.info.x1pdf() );
140  eve->set_x2( pythia.info.x2pdf() );
141  eve->set_id1( pythia.info.id1pdf() );
142  eve->set_id2( pythia.info.id2pdf() );
143  eve->set_pdf1( pythia.info.pdf1() );
144  eve->set_pdf2( pythia.info.pdf2() );
145  eve->set_weight( pythia.info.weight() );
146 
147  // Fill truth particle information, looping over all particles in event.
148  ProMCEvent_Particles *pa= promc.mutable_particles();
149  for (int i = 0; i < pythia.event.size(); i++) {
150 
151  // Fill information particle by particle.
152  pa->add_id( i );
153  pa->add_pdg_id( pythia.event[i].id() );
154  // Particle status in HEPMC style.
155  pa->add_status( pythia.event.statusHepMC(i) );
156  pa->add_mother1( pythia.event[i].mother1() );
157  pa->add_mother2( pythia.event[i].mother2() );
158  pa->add_daughter1( pythia.event[i].daughter1() );
159  pa->add_daughter2( pythia.event[i].daughter2() );
160  // Only store three-momentum and mass, so need to calculate energy.
161  pa->add_px( (int)(pythia.event[i].px()*kEV) );
162  pa->add_py( (int)(pythia.event[i].py()*kEV) );
163  pa->add_pz( (int)(pythia.event[i].pz()*kEV) );
164  pa->add_mass( (int)(pythia.event[i].m()*kEV) );
165  // Store production vertex; will often be the origin.
166  pa->add_x( int(pythia.event[i].xProd()*kL) );
167  pa->add_y( int(pythia.event[i].yProd()*kL) );
168  pa->add_z( int(pythia.event[i].zProd()*kL) );
169  pa->add_t( int(pythia.event[i].tProd()*kL) );
170 
171  } // end loop over particles in the event
172  epbook->write(promc); // write event
173 
174  } // end loop over events
175 
176  // Print statistics.
177  pythia.stat();
178  double sigmapb = pythia.info.sigmaGen() * 1.0E9;
179  cout << "== Cross section for this run = " << sigmapb << " pb" << endl;
180  cout << "== Events for this run = " << Ntot << endl;
181  double lumi = (Ntot/sigmapb)/1000;
182  cout << "== Luminosity for this run = " << lumi << " fb-1" << endl;
183  cout << "\n\n";
184 
185  // Save post-generation statistics for ProMC.
186  ProMCStat stat;
187  stat.set_cross_section_accumulated( sigmapb ); // in pb
188  stat.set_cross_section_error_accumulated( pythia.info.sigmaErr() * 1e9 );
189  stat.set_luminosity_accumulated( Ntot/sigmapb );
190  stat.set_ntried( pythia.info.nTried() );
191  stat.set_nselected( pythia.info.nSelected() );
192  stat.set_naccepted( pythia.info.nAccepted() );
193  epbook->setStatistics( stat );
194 
195  // Close the ProMC file.
196  epbook->close(); // close
197 
198  return 0;
199 }
Generated by   doxygen 1.8.5