DIS-Pythia6/eventgen/generate_ep_dis.C

0001 /*======Electron-Proton Collision Mark V========
0002  *Project: Electron Ion Collider Simulations
0003  *Author: Thomas Krahulik (Stony Brook University)
0004  *Date Updated: February 25, 2017
0005  *Purpose: To simulate electron proton collisions
0006  *To run: root -l generate_ep_dis.C
0007  ===============================================*/
0008
0009
0010 #ifndef __CINT__
0011 #include "TPythia6.h"
0012 #include "TTree.h"
0013 #include "TClonesArray.h"
0014 #include "TH1.h"
0015 #include "TStyle.h"
0016 #include "TCanvas.h"
0017 using namespace std;
0018 #endif
0019
0020 void loadLibraries()
0021 {
0022 #ifdef __CINT__
0023   // Load the Event Generator abstraction library, Pythia 6
0024   // library, and the Pythia 6 interface library.
0025   gSystem->Load("libEG");
0026   gSystem->Load("libPythia6");
0027   gSystem->Load("libEGPythia6");
0028 #endif
0029 }
0030
0031 /**
0032  * eELEC = Electron Energy
0033  * ePROT = Proton Energy
0034  * nEvents = Number of Events
0035  */
0036 int generate_ep_dis(   double eELEC = 10.,
0037                double ePROT = 250.,
0038                int nEvents = 1000 )
0039 {
0040   loadLibraries();
0041
0042   /* build output file name */
0043   TString filename("pythia6.ep.dis.e");
0044   filename += eELEC;
0045   filename += "xp";
0046   filename += ePROT;
0047   filename += ".";
0048   filename += nEvents;
0049   filename += ".root";
0050
0051   /* create output file */
0052   TFile *f = new TFile(filename, "recreate");
0053
0054   TPythia6* pythia = new TPythia6;
0055
0056   //Set Process
0057   pythia->SetMSEL(2);
0058
0059   //Switches
0060   pythia->SetMSTP(14, 30);
0061   pythia->SetMSTP(15, 0);
0062   pythia->SetMSTP(16, 1);
0063   pythia->SetMSTP(17, 4);
0064   pythia->SetMSTP(18, 3);
0065   pythia->SetMSTP(19, 1);
0066   pythia->SetMSTP(20, 0);
0067   pythia->SetMSTP(32, 8);
0068   pythia->SetMSTP(38, 4);
0069   // pythia->SetMSTP(51, 10050);
0070   // pythia->SetMSTP(52, 2);
0071   pythia->SetMSTP(53, 3);
0072   pythia->SetMSTP(54, 1);
0073   pythia->SetMSTP(55, 5);
0074   pythia->SetMSTP(56, 1);
0075   pythia->SetMSTP(57, 1);
0076   pythia->SetMSTP(58, 5);
0077   pythia->SetMSTP(59, 1);
0078   pythia->SetMSTP(60, 7);
0079   pythia->SetMSTP(61, 2);
0080   pythia->SetMSTP(71, 1);
0081   pythia->SetMSTP(81, 0);
0082   pythia->SetMSTP(82, 0);
0083   pythia->SetMSTP(91, 1);
0084   pythia->SetMSTP(92, 3);
0085   pythia->SetMSTP(93, 1);
0086   pythia->SetMSTP(101, 3);
0087   pythia->SetMSTP(102, 1);
0088   pythia->SetMSTP(111, 1);
0089   pythia->SetMSTP(121, 0);
0090   pythia->SetMSTP(127, 1);
0091
0092   //Parameters
0093   pythia->SetPARP(13, 1);
0094   pythia->SetPARP(18, 0.40);
0095   pythia->SetPARP(81, 1.9);
0096   pythia->SetPARP(89, 1800);
0097   pythia->SetPARP(90, 0.16);
0098   pythia->SetPARP(91, 0.40);
0099   pythia->SetPARP(93, 5.);
0100   pythia->SetPARP(99, 0.40);
0101   pythia->SetPARP(100 ,5);
0102   pythia->SetPARP(102, 0.28);
0103   pythia->SetPARP(103, 1.0);
0104   pythia->SetPARP(104, 0.8);
0105   pythia->SetPARP(111, 2.);
0106   pythia->SetPARP(161, 3.00);
0107   pythia->SetPARP(162, 24.6);
0108   pythia->SetPARP(163, 18.8);
0109   pythia->SetPARP(164, 11.5);
0110   pythia->SetPARP(165, 0.47679);
0111   pythia->SetPARP(166, 0.67597);
0112
0113   //Switches for JetSet
0114   pythia->SetMSTJ(1, 1);
0115   pythia->SetMSTJ(12, 1);
0116   pythia->SetMSTJ(45, 5);
0117
0118   pythia->SetMSTU(16, 2);
0119   pythia->SetMSTU(112, 5);
0120   pythia->SetMSTU(113, 5);
0121   pythia->SetMSTU(114, 5);
0122
0123   //Parameters for JetSet
0124   pythia->SetPARJ(1, 0.100);
0125   pythia->SetPARJ(2, 0.300);
0126   pythia->SetPARJ(11, 0.5);
0127   pythia->SetPARJ(12, 0.6);
0128   pythia->SetPARJ(21, 0.40);
0129   pythia->SetPARJ(32, 1.0);
0130   pythia->SetPARJ(33, 0.80);
0131   pythia->SetPARJ(41, 0.30);
0132   pythia->SetPARJ(42, 0.58);
0133   pythia->SetPARJ(45, 0.5);
0134
0135   //Kinematic Variables
0136   pythia->SetCKIN(1, 1.);
0137   pythia->SetCKIN(2, -1.);
0138   pythia->SetCKIN(3, 0.);
0139   pythia->SetCKIN(4, -1.);
0140   pythia->SetCKIN(5, 1.00);
0141   pythia->SetCKIN(6, 1.00);
0142   pythia->SetCKIN(7, -10.);
0143   pythia->SetCKIN(8, 10.);
0144   pythia->SetCKIN(9, -40.);
0145   pythia->SetCKIN(10, 40.);
0146   pythia->SetCKIN(11, -40.);
0147   pythia->SetCKIN(12, 40.);
0148   pythia->SetCKIN(13, -40.);
0149   pythia->SetCKIN(14, 40.);
0150   pythia->SetCKIN(15, -40.);
0151   pythia->SetCKIN(16, 40.);
0152   pythia->SetCKIN(17, -1.);
0153   pythia->SetCKIN(18, 1.);
0154   pythia->SetCKIN(19, -1.);
0155   pythia->SetCKIN(20, 1.);
0156   pythia->SetCKIN(21, 0.);
0157   pythia->SetCKIN(22, 1.);
0158   pythia->SetCKIN(23, 0.);
0159   pythia->SetCKIN(24, 1.);
0160   pythia->SetCKIN(25, -1.);
0161   pythia->SetCKIN(26, 1.);
0162   pythia->SetCKIN(27, -1.);
0163   pythia->SetCKIN(28, 1.);
0164   pythia->SetCKIN(31, 2.);
0165   pythia->SetCKIN(32, -1.);
0166   pythia->SetCKIN(35, 0.);
0167   pythia->SetCKIN(36, -1.);
0168   pythia->SetCKIN(37, 0.);
0169   pythia->SetCKIN(38, -1.);
0170   pythia->SetCKIN(39, 4.);
0171   pythia->SetCKIN(40, -1.);
0172   pythia->SetCKIN(65, 1.e-09);
0173   pythia->SetCKIN(66, -1.);
0174   pythia->SetCKIN(67, 0.);
0175   pythia->SetCKIN(68, -1.);
0176   pythia->SetCKIN(77, 2.0);
0177   pythia->SetCKIN(78, -1.);
0178
0179   //Initialize
0180   pythia->SetP(1, 1, 0.); // x Momentum of Electron
0181   pythia->SetP(1, 2, 0.); // y Momentum of Electron
0182   pythia->SetP(1, 3, -eELEC); // z Momentum of Electron
0183   pythia->SetP(2, 1, 0.); // x Momentum of Proton
0184   pythia->SetP(2, 2, 0.); // y Momentum of Proton
0185   pythia->SetP(2, 3, ePROT); // z Momentum of Proton
0186   pythia->Initialize("3MOM", "e-", "p", 0.);
0187
0188   TTree* tree = new TTree("tree", "Pythia 6 tree");
0189
0190   Int_t nPart;
0191   Double_t Q2, y, x;
0192   TClonesArray* p = (TClonesArray*)pythia->GetListOfParticles();
0193   tree->Branch("p", &p);
0194   TBranch *b_nPart = tree->Branch("nPart", &nPart, "nPart/I");
0195   TBranch *b_Q2 = tree->Branch("Q2", &Q2, "Q2/D");
0196   TBranch *b_y = tree->Branch("y", &y, "y/D");
0197   TBranch *b_x = tree->Branch("x", &x, "x/D");
0198   // TBranch *b_W2 = tree->Branch("W2", &W2, "W2/D");
0199   // TBranch *b_nu = tree->Branch("nu", &nu, "nu/D");
0200
0201   // Now we generate some events
0202   for (Int_t i = 0; i < nEvents; i++) {
0203     // Show how far we got every 100'th event.
0204     if (i % 1000 == 0) cout << "Event # " << i << endl;
0205     // Make one event and fill the tree
0206     pythia->GenerateEvent();
0207
0208     nPart = pythia->GetN();
0209     Q2 = pythia->GetPARI(22);
0210     x = pythia->GetPARI(32);
0211     y = Q2 / (4.*x*eELEC*ePROT);
0212
0213     tree->Fill();
0214     // nParticles->Fill();
0215   }
0216
0217   // Show tree structure
0218   tree->Print();
0219
0220   f->Write();
0221
0222   return 0;
0223 }