Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main91.cc
1 // main91.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 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 // This is a simple test program illustrating how to link in Pythia 6.4 
7 // for the generation of hard processes. That part is now considered
8 // obsolete, but for some debug work this example has been collected
9 // from various code pieces that were separately available up until 
10 // Pythia 8.125. In addition to modifying the below code to fit your
11 // needs, you also have to modify examples/Makefile to link properly
12 // for main91, including the Pythia 6.4xx library to be used (where
13 // xx is the current subversion number). 
14
15 // All hard PYTHIA 6.4 processes should be available for full generation
16 // in PYTHIA 8, at least to the extent that they are defined for p p,
17 // pbar p or e+ e-collisions. Soft processes, i.e. elastic and diffractive 
18 // scattering, as well as minimum-bias events, require a different 
19 // kinematics machinery, and can only be generated with the internal 
20 // PYTHIA 8 processes.
21
22 // PYTHIA 6.4 does its own rejection of events internally, according to
23 // the strategy option 3. However, the current runtime interface does not 
24 // take cross-section information from PYTHIA 6.4. This means that both
25 // the initial maxima and the final cross sections printed by the PYTHIA 8
26 // routines are irrelevant in this case. Instead you have to study the
27 // standard PYTHIA 6.4 initialization printout and call on pystat(...)
28 // at the end of the run to obtain this information. It also means that 
29 // you cannot mix with internally generated PYTHIA 8 events.
30
31 //==========================================================================
32
33 #include "Pythia.h"
34 #include "LHAFortran.h"
35
36 using namespace Pythia8; 
37
38 //==========================================================================
39
40 // Declare the Fortran subroutines that may be used.
41 // This code section is generic.
42
43 #ifdef _WIN32
44   #define pygive_ PYGIVE
45   #define pyinit_ PYINIT
46   #define pyupin_ PYUPIN
47   #define pyupev_ PYUPEV
48   #define pylist_ PYLIST
49   #define pystat_ PYSTAT
50 #endif
51
52 extern "C" {
53 #ifdef _WIN32
54   extern void pyinit_(const char*, int, const char*, int, const char*, 
55     int, double&);
56 #else
57   extern void pyinit_(const char*, const char*, const char*, double&, 
58     int, int, int);
59 #endif
60 }
61
62 extern "C" {
63   extern void pygive_(const char*, int);
64   extern void pyupin_();
65   extern void pyupev_();
66   extern void pylist_(int&);
67   extern void pystat_(int&);
68 }
69
70 //==========================================================================
71
72 // Interfaces to the above routines, to make the C++ calls similar to Fortran.
73 // This code section is generic.
74
75 class Pythia6Interface {
76
77 public:
78
79   // Give in a command to change a setting.
80   static void pygive(const string cmnd) { 
81     const char* cstring = cmnd.c_str(); int len = cmnd.length(); 
82     pygive_(cstring, len);
83   }
84
85   // Initialize the generation for the given beam confiuration.
86   static void pyinit(const string frame, const string beam, 
87     const string target, double wIn) { 
88     const char* cframe = frame.c_str(); int lenframe = frame.length();
89     const char* cbeam = beam.c_str(); int lenbeam = beam.length();
90     const char* ctarget = target.c_str(); int lentarget = target.length();
91 #ifdef _WIN32
92     pyinit_(cframe, lenframe, cbeam, lenbeam, ctarget, lentarget, wIn);
93 #else
94     pyinit_(cframe, cbeam, ctarget, wIn, lenframe, lenbeam, lentarget); 
95 #endif
96   }
97   
98   // Fill the initialization information in the HEPRUP commonblock.
99   static void pyupin() {pyupin_();}
100
101   // Generate the next hard process and 
102   // fill the event information in the HEPEUP commonblock
103   static void pyupev() {pyupev_();}
104
105   // List the event at the process level.
106   static void pylist(int mode) {pylist_(mode);}
107
108   // Print statistics on the event generation process.
109   static void pystat(int mode) {pystat_(mode);}
110
111 };
112
113 //==========================================================================
114
115 // Implement initialization fillHepRup method for Pythia6 example.
116 // This code section is specific to the kind of precesses you generate.
117
118 // Of all parameters that could be set with pygive, only those that 
119 // influence the generation of the hard processes have any impact, 
120 // since this is the only part of the Fortran code that is used. 
121
122 bool LHAupFortran::fillHepRup() { 
123
124   // Set process to generate.
125   // Example 1: QCD production; must set pTmin.  
126   //Pythia6Interface::pygive("msel = 1"); 
127   //Pythia6Interface::pygive("ckin(3) = 20.");
128   // Example 2: t-tbar production.  
129   //Pythia6Interface::pygive("msel = 6"); 
130   // Example 3: Z0 production; must set mMin.
131   Pythia6Interface::pygive("msel = 11"); 
132   Pythia6Interface::pygive("ckin(1) = 50."); 
133
134   // Speed up initialization: multiparton interactions only in C++ code.
135   Pythia6Interface::pygive("mstp(81)=0");
136     
137   // Initialize for 14 TeV pp collider.
138   Pythia6Interface::pyinit("cms","p","p",14000.);   
139
140   // Fill initialization information in HEPRUP.
141   Pythia6Interface::pyupin();
142
143   // Done.
144   return true;
145
146 }
147
148 //==========================================================================
149
150 // Implement event generation fillHepEup method for Pythia 6 example.
151 // This code section is generic.
152
153 bool LHAupFortran::fillHepEup() { 
154
155   // Generate and fill the next Pythia6 event in HEPEUP.
156   Pythia6Interface::pyupev();
157
158   // Done.
159   return true;
160
161 }
162
163 //==========================================================================
164
165 // The main program.
166 // This code section is specific to the physics study you want to do.
167
168 int main() {
169
170   // Generator. Shorthand for the event and for settings.
171   Pythia pythia;
172   Event& event = pythia.event;
173   Settings& settings = pythia.settings;
174
175   // Set Pythia8 generation aspects. Examples only. 
176   pythia.readString("BeamRemnants:primordialKThard = 2.");    
177   pythia.readString("MultipartonInteractions:bProfile = 3");    
178   pythia.readString("Next:numberShowInfo = 0"); 
179   pythia.readString("Next:numberShowProcess = 0"); 
180   pythia.readString("Next:numberShowEvent = 0"); 
181
182   // Initialize to access Pythia6 generator by Les Houches interface.
183   pythia.readString("Beams:frameType = 5"); 
184   LHAupFortran pythia6;
185   pythia.setLHAupPtr( &pythia6);
186   pythia.init();    
187
188   // Set some generation values.
189   int nEvent = 100;
190   int nList  = 1;
191   int nAbort = 10;
192
193   // List changed settings data.
194   settings.listChanged();
195
196   // Histograms.
197   double eCM = 14000.;
198   double epTol = 1e-7 * eCM;
199   Hist epCons("deviation from energy-momentum conservation",100,0.,epTol);
200   Hist nFinal("final particle multiplicity",100,-0.5,1599.5);
201   Hist nChg("final charged multiplicity",100,-0.5,799.5);
202
203   // Begin event loop.
204   int iAbort = 0; 
205   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
206
207     // Generate events. Quit if too many failures.
208     if (!pythia.next()) {
209       if (++iAbort < nAbort) continue;
210       cout << " Event generation aborted prematurely, owing to error!\n"; 
211       break;
212     }
213  
214     // List first few events, both hard process and complete events.
215     if (iEvent < nList) { 
216       pythia.info.list();
217       // This call to Pythia6 is superfluous, but shows it can be done.
218       Pythia6Interface::pylist(1);
219       pythia.process.list();
220       event.list();
221     }
222
223     // Reset quantities to be summed over event.
224     int nfin = 0;
225     int nch = 0;
226     Vec4 pSum = - (event[1].p() + event[2].p());
227
228     // Loop over final particles in the event. 
229     for (int i = 0; i < event.size(); ++i) 
230     if (event[i].isFinal()) {
231       ++nfin;
232       if (event[i].isCharged()) ++nch;
233       pSum += event[i].p();
234     }
235
236     // Fill summed quantities. 
237     double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
238       + abs(pSum.pz());
239     epCons.fill(epDev);
240     nFinal.fill(nfin);
241     nChg.fill(nch);
242
243   // End of event loop.
244   }
245
246   // Final statistics. Must do call to Pythia6 explicitly.
247   pythia.stat();
248   Pythia6Interface::pystat(1);  
249
250   // Histogram output.
251   cout << epCons << nFinal<< nChg; 
252
253   // Done.
254   return 0;
255 }