]>
Commit | Line | Data |
---|---|---|
c6b60c38 | 1 | // main28.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2013 Peter Skands, 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 | // Example of of R-hadron production. | |
7 | // Several of the possibilities shown here, like displaced vertices, | |
8 | // are extras that need not be used for the basic setup. | |
9 | ||
10 | #include "Pythia.h" | |
11 | ||
12 | using namespace Pythia8; | |
13 | ||
14 | int main() { | |
15 | ||
16 | // Key settings to be used in the main program. | |
17 | // nGluino = 0, 1, 2 give stop pair, single gluino or gluino pair. | |
18 | int nGluino = 2; | |
19 | int nEvent = 200; | |
20 | int nAbort = 3; | |
21 | int nList = 0; | |
22 | double eCM = 7000.; | |
23 | ||
24 | // Generator. Shorthand for the event. | |
25 | Pythia pythia; | |
26 | Event& event = pythia.event; | |
27 | ||
28 | // Set up beams: p p is default so only need set energy. | |
29 | pythia.settings.parm("Beams:eCM", eCM); | |
30 | ||
31 | // Squark pair: use stop-antistop as example. | |
32 | if (nGluino == 0) { | |
33 | pythia.readString("SUSY:gg2squarkantisquark = on"); | |
34 | pythia.readString("SUSY:idA = 1000006"); | |
35 | pythia.readString("SUSY:idB = 1000006"); | |
36 | // Squark-gluino pair: also supersymmetric u has been made long-lived. | |
37 | // Stop does not work since then one would need inoming top PDF. | |
38 | // Nevertheless R-hadrons are numbered/named as if containing a stop. | |
39 | } else if (nGluino == 1) { | |
40 | pythia.readString("SUSY:qg2squarkgluino = on"); | |
41 | pythia.readString("SUSY:idA = 1000002"); | |
42 | pythia.readString("RHadrons:idStop = 1000002"); | |
43 | pythia.readString("SUSY:idB = 1000021"); | |
44 | // Gluino pair. | |
45 | } else { | |
46 | pythia.readString("SUSY:gg2gluinogluino = on"); | |
47 | } | |
48 | ||
49 | // Use hacked sps1a file, with stop (+su) and gluino made long-lived. | |
50 | // This is based on the width being less than 0.2 GeV by default. | |
51 | pythia.readString("SLHA:file = sps1aNarrowStopGluino.spc"); | |
52 | // Further hacked file, to test R-parity violating gluino decay. | |
53 | //pythia.readString("SLHA:file = sps1aNarrowStopGluinoRPV.spc"); | |
54 | ||
55 | // Allow R-hadron formation. | |
56 | pythia.readString("Rhadrons:allow = on"); | |
57 | ||
58 | // If you want to do the decay separately later, | |
59 | // you need to switch off automatic decays. | |
60 | pythia.readString("RHadrons:allowDecay = off"); | |
61 | ||
62 | // Fraction of gluinoballs. | |
63 | pythia.readString("RHadrons:probGluinoball = 0.1"); | |
64 | ||
65 | // Switch off key components. | |
66 | //pythia.readString("PartonLevel:MPI = off"); | |
67 | //pythia.readString("PartonLevel:ISR = off"); | |
68 | //pythia.readString("PartonLevel:FSR = off"); | |
69 | //pythia.readString("HadronLevel:Hadronize = off"); | |
70 | ||
71 | // Allow the R-hadrons to have secondary vertices: set c*tau in mm. | |
72 | // Note that width and lifetime can be set independently. | |
73 | // (Nonzero small widths are needed e.g. to select branching ratios.) | |
74 | pythia.readString("1000002:tau0 = 200."); | |
75 | pythia.readString("1000006:tau0 = 250."); | |
76 | pythia.readString("1000021:tau0 = 300."); | |
77 | ||
78 | // Checks. Optionally relax E-p-conservation. | |
79 | pythia.readString("Check:nErrList = 2"); | |
80 | //pythia.readString("Check:epTolErr = 2e-3"); | |
81 | ||
82 | // Possibility to switch off particle data and event listings. | |
83 | // Also to shop location of displaced vertices. | |
84 | pythia.readString("Init:showChangedSettings = on"); | |
85 | pythia.readString("Init:showChangedParticleData = off"); | |
86 | pythia.readString("Next:numberShowInfo = 1"); | |
87 | pythia.readString("Next:numberShowProcess = 1"); | |
88 | pythia.readString("Next:numberShowEvent = 0"); | |
89 | pythia.readString("Next:showScaleAndVertex = on"); | |
90 | ||
91 | // Initialize. | |
92 | pythia.init(); | |
93 | ||
94 | // Histograms. | |
95 | Hist nChargedH("charged multiplicity", 100, -0.5, 799.5); | |
96 | Hist dndyChargedH("dn/dy charged", 100, -10., 10.); | |
97 | Hist dndyRH("dn/dy R-hadrons", 100, -5., 5.); | |
98 | Hist pTRH("pT R-hadrons", 100, 0., 1000.); | |
99 | Hist xRH("p_RHadron / p_sparticle", 100, 0.9, 1.1); | |
100 | Hist mDiff("m(Rhadron) - m(sparticle)", 100, 0., 5.); | |
101 | Hist decVtx("R-hadron decay vertex (mm from origin)", 100, 0., 1000.); | |
102 | ||
103 | // R-hadron flavour composition. | |
104 | map<int, int> flavours; | |
105 | ||
106 | // Begin event loop. | |
107 | int iAbort = 0; | |
108 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { | |
109 | ||
110 | // Generate events. Quit if failure. | |
111 | if (!pythia.next()) { | |
112 | if (++iAbort < nAbort) continue; | |
113 | cout << " Event generation aborted prematurely, owing to error!\n"; | |
114 | break; | |
115 | } | |
116 | ||
117 | // Loop over final charged particles in the event. | |
118 | // The R-hadrons may not yet have decayed here. | |
119 | int nCharged = 0; | |
120 | Vec4 pSum; | |
121 | for (int i = 0; i < event.size(); ++i) { | |
122 | if (event[i].isFinal()) { | |
123 | pSum += event[i].p(); | |
124 | if (event[i].isCharged()) { | |
125 | ++nCharged; | |
126 | dndyChargedH.fill( event[i].y() ); | |
127 | } | |
128 | } | |
129 | } | |
130 | nChargedH.fill( nCharged ); | |
131 | ||
132 | // Loop over final R-hadrons in the event: kinematic distribution | |
133 | for (int i = 0; i < event.size(); ++i) { | |
134 | int idAbs = event[i].idAbs(); | |
135 | if (idAbs > 1000100 && idAbs < 2000000 && idAbs != 1009002) { | |
136 | ++flavours[ event[i].id() ]; | |
137 | dndyRH.fill( event[i].y() ); | |
138 | pTRH.fill( event[i].pT() ); | |
139 | // Trace back to mother; compare momenta and masses. | |
140 | int iMother = i; | |
141 | while( event[iMother].statusAbs() > 100) | |
142 | iMother = event[iMother].mother1(); | |
143 | double xFrac = event[i].pAbs() / event[iMother].pAbs(); | |
144 | xRH.fill( xFrac); | |
145 | double mShift = event[i].m() - event[iMother].m(); | |
146 | mDiff.fill( mShift ); | |
147 | // Separation of R-hadron decay vertex from origin. | |
148 | // Don't be fooled by pAbs(); it gives the three-vector length | |
149 | // of any Vec4, also one representing spatial coordinates. | |
150 | double dist = event[i].vDec().pAbs(); | |
151 | decVtx.fill( dist); | |
152 | ||
153 | // This is a place where you could allow a R-hadron shift of | |
154 | // identity, momentum and decay vertex to allow for detector effects. | |
155 | // Identity not illustrated here; requires a change of mass as well. | |
156 | // Toy model: assume an exponential energy loss, < > = 1 GeV, | |
157 | // but at most half of kinetic energy. Unchanged direction. | |
158 | // Note that event will no longer conserve energy and momentum. | |
159 | double eLossAvg = 1.; | |
160 | double eLoss = 0.; | |
161 | do { eLoss = eLossAvg * pythia.rndm.exp(); } | |
162 | while (eLoss > 0.5 * (event[i].e() - event[i].m())); | |
163 | double eNew = event[i].e() - eLoss; | |
164 | Vec4 pNew = event[i].p() * sqrt( pow2(eNew) - pow2(event[i].m()) ) | |
165 | / event[i].pAbs(); | |
166 | pNew.e( eNew); | |
167 | event[i].p( pNew); | |
168 | // The decay vertex will be calculated based on the production vertex, | |
169 | // the proper lifetime tau and the NEW four-momentum, rather than | |
170 | // e.g. some average momentum, if you do not set it by hand. | |
171 | // This commented-out piece illustrates brute-force setting, | |
172 | // but you should provide real numbers from some tracking program. | |
173 | // With tau = 0 the decay is right at the chosen point. | |
174 | //event[i].tau( 0.); | |
175 | //event[i].vProd( 132., 155., 233., 177.); | |
176 | ||
177 | // End of loop over final R-hadrons. | |
178 | } | |
179 | } | |
180 | ||
181 | // If you have set R-hadrons stable above, | |
182 | // you can still force them to decay at this stage. | |
183 | pythia.forceRHadronDecays(); | |
184 | if (iEvent < nList) pythia.event.list(true); | |
185 | ||
186 | // End of event loop. | |
187 | } | |
188 | ||
189 | // Final statistics, flavour composition and histogram output. | |
190 | pythia.stat(); | |
191 | cout << "\n Composition of produced R-hadrons \n code " | |
192 | << "name times " << endl; | |
193 | for (map<int, int>::iterator flavNow = flavours.begin(); | |
194 | flavNow != flavours.end(); ++flavNow) cout << setw(8) | |
195 | << flavNow->first << setw(16) << pythia.particleData.name(flavNow->first) | |
196 | << setw(8) << flavNow->second << endl; | |
197 | cout << nChargedH << dndyChargedH << dndyRH << pTRH << xRH << mDiff | |
198 | << decVtx; | |
199 | ||
200 | // Done. | |
201 | return 0; | |
202 | } |