Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main27.cc
1 // main27.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 // Kaluza-Klein gamma*/Z resonances in TeV-sized extra dimensions.
7
8 #include <assert.h>
9 #include <time.h>
10 #include <sstream>
11
12 #include "Pythia.h"
13 using namespace Pythia8; 
14
15 int main() {
16
17   // Generator.
18   Pythia pythia;
19   ParticleData& pdt = pythia.particleData;
20
21   // Pick new random number seed for each run, based on clock.
22   pythia.readString("Random:setSeed = on");
23   pythia.readString("Random:seed = 0");
24  
25   // Process selection.
26   // ANY COMBINATION OF THE PROCESSES FLAGS BELOW IS ALLOWED
27   // HERE WE SWITCH ON ONLY THE MU+MU- FINAL STATE.
28   // TO SWITCH ALL POSSIBLE FINAL STATES ON, UNCOMMENT ALL
29   // THE RELEVANT LINES BELOW:
30   //pythia.readString("ExtraDimensionsTEV:ffbar2e+e- = on");
31   pythia.readString("ExtraDimensionsTEV:ffbar2mu+mu- = on");
32   //pythia.readString("ExtraDimensionsTEV:ffbar2tau+tau- = on");
33   //pythia.readString("ExtraDimensionsTEV:ffbar2uubar = on");
34   //pythia.readString("ExtraDimensionsTEV:ffbar2ddbar = on");
35   //pythia.readString("ExtraDimensionsTEV:ffbar2ccbar = on");
36   //pythia.readString("ExtraDimensionsTEV:ffbar2ssbar = on");
37   //pythia.readString("ExtraDimensionsTEV:ffbar2bbbar = on");
38   //pythia.readString("ExtraDimensionsTEV:ffbar2ttbar = on");
39   //pythia.readString("ExtraDimensionsTEV:ffbar2nuenuebar = on");
40   //pythia.readString("ExtraDimensionsTEV:ffbar2numunumubar = on");
41   //pythia.readString("ExtraDimensionsTEV:ffbar2nutaunutaubar = on");
42
43   // Pick KK mass.  
44   double newMass = 4000.; // GeV
45   cout << "|-------------------------" << endl;
46   cout << "| KK mass is: " << newMass << endl;
47   cout << "|-------------------------" << endl;
48   stringstream strm;
49   string sNewMass, sNewWidth, sNewLowBound, sNewHighBound;
50
51   // Manually set the mass and therefore the width 
52   // and the phase space for the sampling
53   strm.clear();
54   strm << newMass;
55   strm >> sNewMass;
56   strm.clear();
57   strm << newMass / pdt.m0(5000023) * pdt.mWidth(5000023);
58   strm >> sNewWidth;
59   strm.clear();
60   strm << newMass/4.;
61   strm >> sNewLowBound;
62   strm.clear();
63   strm << newMass*2.;
64   strm >> sNewHighBound;
65
66   // Feed in KK state information and other generation specifics.  
67   pythia.readString("5000023:m0 = " + sNewMass);
68   pythia.readString("5000023:mWidth = " + sNewWidth);
69   pythia.readString("5000023:mMin = " + sNewLowBound);
70   pythia.readString("5000023:mMax = " + sNewHighBound);
71   //////////////////////////////////////////////////////////////////////////
72   pythia.readString("5000023:isResonance = false"); // THIS IS MANDATORY  //
73   //////////////////////////////////////////////////////////////////////////
74   // 0=(gm+Z), 1=(gm), 2=(Z), 3=(gm+Z+gmKK+ZKK), 4=(m+Z+gmKK), 5=(m+Z+ZKK)
75   pythia.readString("ExtraDimensionsTEV:gmZmode = 3"); 
76   // min=0, max=100, default=10
77   pythia.readString("ExtraDimensionsTEV:nMax = 100"); 
78   pythia.readString("ExtraDimensionsTEV:mStar = " + sNewMass);
79   pythia.readString("PhaseSpace:mHatMin = " + sNewLowBound);
80   pythia.readString("PhaseSpace:mHatMax = " + sNewHighBound);
81
82   // Initialize for LHC.
83   pythia.readString("Beams:eCM = 14000.");    
84   pythia.init();
85
86   // Histograms.
87   Hist mHatHisto("dN/dmHat", 50, newMass/4., newMass*2.);
88   Hist pTmuHisto("(dN/dpT)_mu^-", 50, 1., 2501.);
89
90   vector<int> moms;
91   
92   // Measure the cpu runtime.
93   clock_t start, stop;
94   double t = 0.0;
95   // Depending on operating system, either of lines below gives warning.
96   //assert((start = clock()) != -1); // Start timer; clock_t signed.
97   //assert((start = clock()) != -1u); // Start timer; clock_t unsigned.
98   // Simpler option, not using assert.
99   start = clock();
100   
101   // Begin event loop. Generate event. Skip if error. List first one.
102   for (int iEvent = 0 ; iEvent < 500 ; ++iEvent) {
103     if (!pythia.next()) continue;
104
105     // Begin event analysis.
106     bool isZ = false;
107     bool ismu = false;
108     int iZ = 0;
109     int imu = 0;
110     for (int i = 0 ; i < pythia.event.size() ; ++i) {
111       // find the most recent Z
112       if (pythia.event[i].id() == 5000023) {
113         iZ = i;
114         isZ = true;
115       }
116       // find the final muon who's first mother is the Z
117       if (pythia.event[i].id() == 13 && pythia.event[i].isFinal()) {
118         moms.clear();
119         moms = pythia.event.motherList(i);
120         for (int m = 0 ; m < int(moms.size()) ; m++) {
121           if( pythia.event[ moms[m] ].id() == 5000023 ) {
122             imu = i;
123             ismu = true;
124             break;
125           } // end if 5000023
126         } // end for moms.size()
127       } // end if final muon
128     } // end for event.size()
129     if(isZ && ismu) {
130       mHatHisto.fill( pythia.event[iZ].m() );
131       pTmuHisto.fill( pythia.event[imu].pT() );
132     }
133     if(iEvent%10 == 0) cout << "Event: " << iEvent << endl << flush;
134   } // end for iEvent<500
135
136   // Done. Print results.
137   stop = clock(); // Stop timer
138   t = (double) (stop-start)/CLOCKS_PER_SEC;
139
140   pythia.stat();
141   cout << mHatHisto;
142   cout << pTmuHisto;
143
144   cout << "\n" << "|----------------------------------------|" << endl;
145   cout << "| CPU Runtime = " << t << " sec" << endl;
146   cout << "|----------------------------------------|" << "\n" << endl;
147   
148   return 0;
149 }