Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
generate_ep_dis.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file generate_ep_dis.C
1 /*======Electron-Proton Collision Mark V========
2  *Project: Electron Ion Collider Simulations
3  *Author: Thomas Krahulik (Stony Brook University)
4  *Date Updated: February 25, 2017
5  *Purpose: To simulate electron proton collisions
6  *To run: root -l generate_ep_dis.C
7  ===============================================*/
8 
9 
10 #ifndef __CINT__
11 #include "TPythia6.h"
12 #include "TTree.h"
13 #include "TClonesArray.h"
14 #include "TH1.h"
15 #include "TStyle.h"
16 #include "TCanvas.h"
17 using namespace std;
18 #endif
19 
20 void loadLibraries()
21 {
22 #ifdef __CINT__
23  // Load the Event Generator abstraction library, Pythia 6
24  // library, and the Pythia 6 interface library.
25  gSystem->Load("libEG");
26  gSystem->Load("libPythia6");
27  gSystem->Load("libEGPythia6");
28 #endif
29 }
30 
36 int generate_ep_dis( double eELEC = 10.,
37  double ePROT = 250.,
38  int nEvents = 1000 )
39 {
40  loadLibraries();
41 
42  /* build output file name */
43  TString filename("pythia6.ep.dis.e");
44  filename += eELEC;
45  filename += "xp";
46  filename += ePROT;
47  filename += ".";
48  filename += nEvents;
49  filename += ".root";
50 
51  /* create output file */
52  TFile *f = new TFile(filename, "recreate");
53 
54  TPythia6* pythia = new TPythia6;
55 
56  //Set Process
57  pythia->SetMSEL(2);
58 
59  //Switches
60  pythia->SetMSTP(14, 30);
61  pythia->SetMSTP(15, 0);
62  pythia->SetMSTP(16, 1);
63  pythia->SetMSTP(17, 4);
64  pythia->SetMSTP(18, 3);
65  pythia->SetMSTP(19, 1);
66  pythia->SetMSTP(20, 0);
67  pythia->SetMSTP(32, 8);
68  pythia->SetMSTP(38, 4);
69  // pythia->SetMSTP(51, 10050);
70  // pythia->SetMSTP(52, 2);
71  pythia->SetMSTP(53, 3);
72  pythia->SetMSTP(54, 1);
73  pythia->SetMSTP(55, 5);
74  pythia->SetMSTP(56, 1);
75  pythia->SetMSTP(57, 1);
76  pythia->SetMSTP(58, 5);
77  pythia->SetMSTP(59, 1);
78  pythia->SetMSTP(60, 7);
79  pythia->SetMSTP(61, 2);
80  pythia->SetMSTP(71, 1);
81  pythia->SetMSTP(81, 0);
82  pythia->SetMSTP(82, 0);
83  pythia->SetMSTP(91, 1);
84  pythia->SetMSTP(92, 3);
85  pythia->SetMSTP(93, 1);
86  pythia->SetMSTP(101, 3);
87  pythia->SetMSTP(102, 1);
88  pythia->SetMSTP(111, 1);
89  pythia->SetMSTP(121, 0);
90  pythia->SetMSTP(127, 1);
91 
92  //Parameters
93  pythia->SetPARP(13, 1);
94  pythia->SetPARP(18, 0.40);
95  pythia->SetPARP(81, 1.9);
96  pythia->SetPARP(89, 1800);
97  pythia->SetPARP(90, 0.16);
98  pythia->SetPARP(91, 0.40);
99  pythia->SetPARP(93, 5.);
100  pythia->SetPARP(99, 0.40);
101  pythia->SetPARP(100 ,5);
102  pythia->SetPARP(102, 0.28);
103  pythia->SetPARP(103, 1.0);
104  pythia->SetPARP(104, 0.8);
105  pythia->SetPARP(111, 2.);
106  pythia->SetPARP(161, 3.00);
107  pythia->SetPARP(162, 24.6);
108  pythia->SetPARP(163, 18.8);
109  pythia->SetPARP(164, 11.5);
110  pythia->SetPARP(165, 0.47679);
111  pythia->SetPARP(166, 0.67597);
112 
113  //Switches for JetSet
114  pythia->SetMSTJ(1, 1);
115  pythia->SetMSTJ(12, 1);
116  pythia->SetMSTJ(45, 5);
117 
118  pythia->SetMSTU(16, 2);
119  pythia->SetMSTU(112, 5);
120  pythia->SetMSTU(113, 5);
121  pythia->SetMSTU(114, 5);
122 
123  //Parameters for JetSet
124  pythia->SetPARJ(1, 0.100);
125  pythia->SetPARJ(2, 0.300);
126  pythia->SetPARJ(11, 0.5);
127  pythia->SetPARJ(12, 0.6);
128  pythia->SetPARJ(21, 0.40);
129  pythia->SetPARJ(32, 1.0);
130  pythia->SetPARJ(33, 0.80);
131  pythia->SetPARJ(41, 0.30);
132  pythia->SetPARJ(42, 0.58);
133  pythia->SetPARJ(45, 0.5);
134 
135  //Kinematic Variables
136  pythia->SetCKIN(1, 1.);
137  pythia->SetCKIN(2, -1.);
138  pythia->SetCKIN(3, 0.);
139  pythia->SetCKIN(4, -1.);
140  pythia->SetCKIN(5, 1.00);
141  pythia->SetCKIN(6, 1.00);
142  pythia->SetCKIN(7, -10.);
143  pythia->SetCKIN(8, 10.);
144  pythia->SetCKIN(9, -40.);
145  pythia->SetCKIN(10, 40.);
146  pythia->SetCKIN(11, -40.);
147  pythia->SetCKIN(12, 40.);
148  pythia->SetCKIN(13, -40.);
149  pythia->SetCKIN(14, 40.);
150  pythia->SetCKIN(15, -40.);
151  pythia->SetCKIN(16, 40.);
152  pythia->SetCKIN(17, -1.);
153  pythia->SetCKIN(18, 1.);
154  pythia->SetCKIN(19, -1.);
155  pythia->SetCKIN(20, 1.);
156  pythia->SetCKIN(21, 0.);
157  pythia->SetCKIN(22, 1.);
158  pythia->SetCKIN(23, 0.);
159  pythia->SetCKIN(24, 1.);
160  pythia->SetCKIN(25, -1.);
161  pythia->SetCKIN(26, 1.);
162  pythia->SetCKIN(27, -1.);
163  pythia->SetCKIN(28, 1.);
164  pythia->SetCKIN(31, 2.);
165  pythia->SetCKIN(32, -1.);
166  pythia->SetCKIN(35, 0.);
167  pythia->SetCKIN(36, -1.);
168  pythia->SetCKIN(37, 0.);
169  pythia->SetCKIN(38, -1.);
170  pythia->SetCKIN(39, 4.);
171  pythia->SetCKIN(40, -1.);
172  pythia->SetCKIN(65, 1.e-09);
173  pythia->SetCKIN(66, -1.);
174  pythia->SetCKIN(67, 0.);
175  pythia->SetCKIN(68, -1.);
176  pythia->SetCKIN(77, 2.0);
177  pythia->SetCKIN(78, -1.);
178 
179  //Initialize
180  pythia->SetP(1, 1, 0.); // x Momentum of Electron
181  pythia->SetP(1, 2, 0.); // y Momentum of Electron
182  pythia->SetP(1, 3, -eELEC); // z Momentum of Electron
183  pythia->SetP(2, 1, 0.); // x Momentum of Proton
184  pythia->SetP(2, 2, 0.); // y Momentum of Proton
185  pythia->SetP(2, 3, ePROT); // z Momentum of Proton
186  pythia->Initialize("3MOM", "e-", "p", 0.);
187 
188  TTree* tree = new TTree("tree", "Pythia 6 tree");
189 
190  Int_t nPart;
191  Double_t Q2, y, x;
192  TClonesArray* p = (TClonesArray*)pythia->GetListOfParticles();
193  tree->Branch("p", &p);
194  TBranch *b_nPart = tree->Branch("nPart", &nPart, "nPart/I");
195  TBranch *b_Q2 = tree->Branch("Q2", &Q2, "Q2/D");
196  TBranch *b_y = tree->Branch("y", &y, "y/D");
197  TBranch *b_x = tree->Branch("x", &x, "x/D");
198  // TBranch *b_W2 = tree->Branch("W2", &W2, "W2/D");
199  // TBranch *b_nu = tree->Branch("nu", &nu, "nu/D");
200 
201  // Now we generate some events
202  for (Int_t i = 0; i < nEvents; i++) {
203  // Show how far we got every 100'th event.
204  if (i % 1000 == 0) cout << "Event # " << i << endl;
205  // Make one event and fill the tree
206  pythia->GenerateEvent();
207 
208  nPart = pythia->GetN();
209  Q2 = pythia->GetPARI(22);
210  x = pythia->GetPARI(32);
211  y = Q2 / (4.*x*eELEC*ePROT);
212 
213  tree->Fill();
214  // nParticles->Fill();
215  }
216 
217  // Show tree structure
218  tree->Print();
219 
220  f->Write();
221 
222  return 0;
223 }