EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
runpythia.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file runpythia.cpp
1 // Example for generating EventPythia objects directly from PYTHIA
2 // output, without intermediate text output.
3 
4 #include<iostream>
5 #include<string>
6 #include<TString.h>
7 #include<TPythia6.h>
8 #include<TRandom3.h>
9 #include<TFile.h>
10 #include<TStopwatch.h>
11 
12 #include <eicsmear/erhic/Pythia6.h>
13 #include <eicsmear/erhic/Pythia6EventBuilder.h>
14 
15 using namespace std;
16 
17 
18 void runpythia(TString outFile,
19  int nEvents,
20  double pElectron,
21  double pProton,
22  double minQ2 = 1.,
23  double maxQ2 = -1.,
24  int messageInterval = 1000){
25  // bool doFiltering = false,
26  // double minz = 0.,
27  // double maxcorrelation = 2.)
28 
29  cout << "Input arguments:" << endl;
30  cout << "Output file: " << outFile << endl;
31  cout << "Number of events: " << nEvents << endl;
32  cout << "Electron momentum: " << pElectron << endl;
33  cout << "Proton momentum: " << pProton << endl;
34  cout << "Minimum Q2: " << minQ2 << endl;
35  cout << "Maximum Q2: " << maxQ2 << endl;
36  cout << "Message interval: " << messageInterval << endl;
37  // cout << "Filtering: " << doFiltering << endl;
38  // cout << "Minimum z: " << minz << endl;
39  // cout << "Max pair correlation: " << maxcorrelation << endl;
40 
41  // The behaviour of PYTHIA 6 is configured via the ROOT interface
42  // class TPythia6.
43  TPythia6* pythia = TPythia6::Instance();
44 
45  // Configure for charm production
46  // I copied these parameters from Elke's charm log files
47  pythia->SetMSEL(4);
48  pythia->SetMSTP(14, 30);
49  pythia->SetMSTP(15, 0);
50  pythia->SetMSTP(16, 1);
51  pythia->SetMSTP(17, 4); // MSTP 17=6 is the R-rho measured as by hermes, =4 Default
52  pythia->SetMSTP(18, 3);
53  pythia->SetMSTP(19, 1); // Hermes MSTP-19=1 different Q2 suppression, default = 4
54  pythia->SetMSTP(20, 0); // Hermes MSTP(20, 0 , default MSTP(20, 3
55  pythia->SetMSTP(32, 8);
56  pythia->SetMSTP(38, 4);
57  pythia->SetMSTP(51, 10800); // PDFLIB/LHAPDF 10800 = CT10
58  pythia->SetMSTP(52, 2); // MSTP(52) = 2 indicates to use PDFLIB
59  pythia->SetMSTP(53, 3);
60  pythia->SetMSTP(54, 1);
61  pythia->SetMSTP(55, 5);
62  pythia->SetMSTP(56, 1);
63  pythia->SetMSTP(57, 1);
64  pythia->SetMSTP(58, 5);
65  pythia->SetMSTP(59, 1);
66  pythia->SetMSTP(60, 7);
67  pythia->SetMSTP(61, 2);
68  pythia->SetMSTP(71, 1);
69  pythia->SetMSTP(81, 0);
70  pythia->SetMSTP(82, 1);
71  pythia->SetMSTP(91, 1);
72  pythia->SetMSTP(92, 3); // hermes MSTP(92, 4
73  pythia->SetMSTP(93, 1);
74  pythia->SetMSTP(101, 3);
75  pythia->SetMSTP(102, 1);
76  pythia->SetMSTP(111, 1);
77  pythia->SetMSTP(121, 0);
78  pythia->SetPARP(13, 1);
79  pythia->SetPARP(18, 0.4); // hermes pythia->SetPARP(18, 0.17
80  pythia->SetPARP(81, 1.9);
81  pythia->SetPARP(89, 1800);
82  pythia->SetPARP(90, 0.16);
83  pythia->SetPARP(91, 0.40);
84  pythia->SetPARP(93, 5.);
85  pythia->SetPARP(99, 0.40);
86  pythia->SetPARP(100, 5);
87  pythia->SetPARP(102, 0.28);
88  pythia->SetPARP(103, 1.0);
89  pythia->SetPARP(104, 0.8);
90  pythia->SetPARP(111, 2.);
91  pythia->SetPARP(161, 3.00);
92  pythia->SetPARP(162, 24.6);
93  pythia->SetPARP(163, 18.8);
94  pythia->SetPARP(164, 11.5);
95  pythia->SetPARP(165, 0.47679);
96  pythia->SetPARP(166, 0.67597); // PARP165/166 are linked to MSTP17 as R_rho of HERMES is used
97  pythia->SetPARJ(1, 0.100);
98  pythia->SetPARJ(2, 0.300);
99  pythia->SetPARJ(11, 0.5);
100  pythia->SetPARJ(12, 0.6);
101  pythia->SetPARJ(21, 0.40);
102  pythia->SetPARJ(32, 1.0);
103  pythia->SetPARJ(33, 0.80);
104  pythia->SetPARJ(41, 0.30);
105  pythia->SetPARJ(42, 0.58);
106  pythia->SetPARJ(45, 0.5);
107  pythia->SetMSTJ(1, 1);
108  pythia->SetMSTJ(12, 1);
109  pythia->SetMSTJ(45, 5);
110  pythia->SetMSTU(112, 5);
111  pythia->SetMSTU(113, 5);
112  pythia->SetMSTU(114, 5);
113  pythia->SetCKIN(1, 1.);
114  pythia->SetCKIN(2, -1.);
115  pythia->SetCKIN(3, 0.);
116  pythia->SetCKIN(4, -1.);
117  pythia->SetCKIN(5, 1.00);
118  pythia->SetCKIN(6, 1.00);
119  pythia->SetCKIN(7, -10.);
120  pythia->SetCKIN(8, 10.);
121  pythia->SetCKIN(9, -40.);
122  pythia->SetCKIN(10, 40.);
123  pythia->SetCKIN(11, -40.);
124  pythia->SetCKIN(12, 40.);
125  pythia->SetCKIN(13, -40.);
126  pythia->SetCKIN(14, 40.);
127  pythia->SetCKIN(15, -40.);
128  pythia->SetCKIN(16, 40.);
129  pythia->SetCKIN(17, -1.);
130  pythia->SetCKIN(18, 1.);
131  pythia->SetCKIN(19, -1.);
132  pythia->SetCKIN(20, 1.);
133  pythia->SetCKIN(21, 0.);
134  pythia->SetCKIN(22, 1.);
135  pythia->SetCKIN(23, 0.);
136  pythia->SetCKIN(24, 1.);
137  pythia->SetCKIN(25, -1.);
138  pythia->SetCKIN(26, 1.);
139  pythia->SetCKIN(27, -1.);
140  pythia->SetCKIN(28, 1.);
141  pythia->SetCKIN(31, 2.);
142  pythia->SetCKIN(32, -1.);
143  pythia->SetCKIN(35, 0.);
144  pythia->SetCKIN(36, -1);
145  pythia->SetCKIN(37, 0.);
146  pythia->SetCKIN(38, -1.);
147  pythia->SetCKIN(39, 4.);
148  pythia->SetCKIN(40, -1.);
149  pythia->SetCKIN(65, minQ2); // Min for Q^2
150  pythia->SetCKIN(66, maxQ2); // Max for Q^2
151  pythia->SetCKIN(67, 0.);
152  pythia->SetCKIN(68, -1.);
153  pythia->SetCKIN(77, 2.0);
154  pythia->SetCKIN(78, -1.);
155 
156  // Set beam momenta.
157  // The fourth argument is a dummy value when calling with frame "3MOM".
158  // Set the electron 3-momenta via SetP(1, ...) (beam 1)
159  // and the proton 3-momenta via SetP(2, ...) (beam 2)
160  pythia->SetP(1, 1, 0.); // px
161  pythia->SetP(1, 2, 0.); // py
162  pythia->SetP(1, 3, -pElectron); // pz
163  pythia->SetP(2, 1, 0.); // px
164  pythia->SetP(2, 2, 0.); // py
165  pythia->SetP(2, 3, pProton); // pz
166 
167  // Set the random number generator seed.
168  // Seeds can be in the range [0, 900,000,000].
169  // TRandom::Integer(N) returns number in the range [0, N).
170  // Set the TRandom seed to zero so the seed is set to a random value.
171  TRandom3 random(0);
172  pythia->SetMRPY(1, random.Integer(900000001));
173  std::cout << "PYTHIA random number seed MRPY(1) = " << pythia->GetMRPY(1)
174  << std::endl;
175 
176  // Intialise PYTHIA for e-p collisions
177  char* target = (char*) "p+";
178  char* beam = (char*) "gamma/e-";
179  char* frame = (char*) "3MOM";
180  float WIN = pElectron;
181  pythia->Pyinit(frame, beam, target, WIN);
182 
183  // Open the output
184  TFile* file = new TFile(outFile, "RECREATE");
185  erhic::VirtualEventFactory* factory = new erhic::Pythia6EventBuilder();
186  erhic::Pythia6 pythia6(file, factory, nEvents, "EICTree", "event", messageInterval);
187 
188  // The event filter code uses the particle filter
189  // classes for the dihadron code.
190  // The file is provided but untested.
191  // gROOT->LoadMacro("filter.cpp+g"); // Defines class DiHadronEventFilter
192  // DiHadronEventFilter filter(minz, maxcorrelation);
193  // if(doFiltering) pythia6.SetFilter(&filter);
194 
195  // Run PYTHIA and write output
196  TStopwatch timer;
197  pythia6.Run();
198 
199  std::cout << "Completed in " << timer.RealTime() << " seconds" << std::endl;
200  file->Close();
201 
202  return;
203 }