ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dasu/UltraFastSim/generateData.cc
(Generate patch)

Comparing UserCode/dasu/UltraFastSim/generateData.cc (file contents):
Revision 1.3 by dasu, Thu Feb 10 14:15:45 2011 UTC vs.
Revision 1.18 by dasu, Wed Nov 2 21:18:08 2011 UTC

# Line 14 | Line 14 | using namespace std;
14   using namespace Pythia8;
15  
16   #include "UltraFastSim.h"
17 + #include "UFSDataStore.h"
18 +
19 + int poisson(double mean, Rndm *rndmPtr) {
20 +  static double oldMean = -1;
21 +  static double g;
22 +  if(mean != oldMean) {
23 +    oldMean = mean;
24 +    if(mean == 0) {
25 +      g = 0;
26 +    }
27 +    else {
28 +      g = exp(-mean);
29 +    }
30 +  }    
31 +  double em = -1;
32 +  double t = 1;
33 +  do {
34 +    em++;
35 +    t *= rndmPtr->flat();
36 +  } while(t > g);
37 +  return em;
38 + }
39  
40   int main(int argc, char **argv) {
41    // Generator. Process selection. LHC initialization. Histogram.
# Line 21 | Line 43 | int main(int argc, char **argv) {
43    int nEvents = 0;
44    int runNumber = 0;
45    int meanPileupEventCount = 0;
46 <  if(argc < 3 || argc > 5)
46 >  float cmEnergy = 14000.;
47 >  if(argc < 3 || argc > 6)
48      {
49 <      cerr << "Command syntax: " << argv[0] << " BBH[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
50 <      cerr << "Command syntax: " << argv[0] << " BBA[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51 <      cerr << "or              " << argv[0] << " ZLL[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
52 <      cerr << "or              " << argv[0] << " ZMM[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
53 <      cerr << "or              " << argv[0] << " TLN[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
54 <      cerr << "or              " << argv[0] << " WLN[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
49 >      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
50 >      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
51 >      cerr << "or              " << argv[0] << " ZHmmXX[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
52 >      cerr << "or              " << argv[0] << " ZHeeXX[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
53 >      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
54 >      cerr << "or              " << argv[0] << " ZZmmnn[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
55 >      cerr << "or              " << argv[0] << " ZZeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
56 >      cerr << "or              " << argv[0] << " ZZeenn[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
57 >      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
58 >      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
59 >      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
60 >      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
61 >      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
62 >      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
63 >      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
64 >      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
65 >      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
66 >      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
67 >      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
68 >      cerr << "or              " << argv[0] << " QCD[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
69 >      cerr << "or              " << argv[0] << " QED[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
70 >      cerr << "or              " << argv[0] << " EWK[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
71 >      cerr << "or              " << argv[0] << " Top[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
72 >      cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
73 >      cerr << "QCD, QED, EWK and Top do not specify the decays -- they are useful for generic BG production" << endl;
74        exit(1);
75      }
76    else {
77      nEvents = atoi(argv[2]);
78      if(argc >= 4) runNumber = atoi(argv[3]);
79      if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
80 <    if(strncmp(argv[1], "ZLL", 3) == 0)
80 >    if(argc >= 6) cmEnergy = atof(argv[5]);
81 >    if(strncmp(argv[1], "Zee", 3) == 0)
82        {
83          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
84          pythia.readString("PhaseSpace:mHatMin = 60.");    
85          pythia.readString("PhaseSpace:mHatMax = -1");    
86 <        pythia.readString("23:onMode = 0");       // Z decays to leptons only
87 <        pythia.readString("23:onIfAny = 11 13 15");
45 <        pythia.readString("15:onMode = 2");       // tau- decays to electrons and muons - second tau decays in all modes
46 <        pythia.readString("15:onIfAny = 11 13");
86 >        pythia.readString("23:onMode = 0");       // Z decays to electrons only
87 >        pythia.readString("23:onIfAny = 11");
88        }
89 <    else if(strncmp(argv[1], "ZMM", 3) == 0)
89 >    else if(strncmp(argv[1], "Zmm", 3) == 0)
90        {
91          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
92          pythia.readString("PhaseSpace:mHatMin = 60.");    
# Line 53 | Line 94 | int main(int argc, char **argv) {
94          pythia.readString("23:onMode = 0");       // Z decays to muons only
95          pythia.readString("23:onIfAny = 13");
96        }
97 <    else if(strncmp(argv[1], "WLN", 3) == 0)
97 >    else if(strncmp(argv[1], "Zeh", 3) == 0)
98 >      {
99 >        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
100 >        pythia.readString("PhaseSpace:mHatMin = 60.");    
101 >        pythia.readString("PhaseSpace:mHatMax = -1");    
102 >        pythia.readString("23:onMode = 0");       // Z decays to taus only
103 >        pythia.readString("23:onIfAny = 15");
104 >        pythia.readString("15:onMode = 2");       // tau- decays to electrons - second tau decays in all modes
105 >        pythia.readString("15:onIfAny = 11");
106 >      }
107 >    else if(strncmp(argv[1], "Zmh", 3) == 0)
108 >      {
109 >        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
110 >        pythia.readString("PhaseSpace:mHatMin = 60.");    
111 >        pythia.readString("PhaseSpace:mHatMax = -1");    
112 >        pythia.readString("23:onMode = 0");       // Z decays to taus only
113 >        pythia.readString("23:onIfAny = 15");
114 >        pythia.readString("15:onMode = 2");       // tau- decays to muons - second tau decays in all modes
115 >        pythia.readString("15:onIfAny = 13");
116 >      }
117 >    else if(strncmp(argv[1], "Zhh", 3) == 0)
118 >      {
119 >        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
120 >        pythia.readString("PhaseSpace:mHatMin = 60.");    
121 >        pythia.readString("PhaseSpace:mHatMax = -1");    
122 >        pythia.readString("23:onMode = 0");       // Z decays to taus only
123 >        pythia.readString("23:onIfAny = 15");
124 >      }
125 >    else if(strncmp(argv[1], "Wen", 3) == 0)
126        {
127          pythia.readString("WeakSingleBoson:ffbar2W = on");
128 <        pythia.readString("24:onMode = 0");       // W decays to light leptons only
129 <        pythia.readString("24:onIfAny = 11 13");
128 >        pythia.readString("24:onMode = 0");       // W decays to electron
129 >        pythia.readString("24:onIfAny = 11");
130        }
131 <    else if(strncmp(argv[1], "TLN", 3) == 0)
131 >    else if(strncmp(argv[1], "Wmn", 3) == 0)
132 >      {
133 >        pythia.readString("WeakSingleBoson:ffbar2W = on");
134 >        pythia.readString("24:onMode = 0");       // W decays to muon
135 >        pythia.readString("24:onIfAny = 13");
136 >      }
137 >    else if(strncmp(argv[1], "Ten", 3) == 0)
138 >      {
139 >        pythia.readString("Top:all = on");
140 >        pythia.readString("24:onMode = 2");       // W+ decays to electron
141 >        pythia.readString("24:onIfAny = 11");
142 >      }
143 >    else if(strncmp(argv[1], "Tmn", 3) == 0)
144        {
145          pythia.readString("Top:all = on");
146 <        pythia.readString("24:onMode = 2");       // W+ decays to light leptons only
147 <        pythia.readString("24:onIfAny = 11 13");
146 >        pythia.readString("24:onMode = 2");       // W+ decays to muon
147 >        pythia.readString("24:onIfAny = 13");
148        }
149 <    else if(strncmp(argv[1], "ZH", 2) == 0)
149 >    else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
150        {
151          pythia.readString("HiggsSM:ffbar2HZ = on");
152          pythia.readString("25:m0 = 120");
# Line 74 | Line 155 | int main(int argc, char **argv) {
155          pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
156          pythia.readString("23:onIfAny = 13");
157        }
158 +     else if(strncmp(argv[1], "ZHmmXX", 6) == 0)
159 +       {
160 +        pythia.readString("HiggsSM:ffbar2HZ = on");
161 +        pythia.readString("25:m0 = 120");
162 +        pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
163 +        pythia.readString("25:isVisible = false");
164 +        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
165 +        pythia.readString("23:onIfAny = 13");
166 +       }
167 +    else if(strncmp(argv[1], "ZHeebb", 6) == 0)
168 +      {
169 +        pythia.readString("HiggsSM:ffbar2HZ = on");
170 +        pythia.readString("25:m0 = 120");
171 +        pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
172 +        pythia.readString("25:onIfAny = 5");
173 +        pythia.readString("23:onMode = 0");      // Z decays to e+,e- only (for trigger)
174 +        pythia.readString("23:onIfAny = 11");
175 +      }
176 +     else if(strncmp(argv[1], "ZHeeXX", 6) == 0)
177 +       {
178 +        pythia.readString("HiggsSM:ffbar2HZ = on");
179 +        pythia.readString("25:m0 = 120");
180 +        pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
181 +        pythia.readString("25:isVisible = false");
182 +        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
183 +        pythia.readString("23:onIfAny = 11");
184 +       }
185      else if(strncmp(argv[1], "BB", 2) == 0)
186        {
187          pythia.readString("Higgs:useBSM = on");
80        if(strncmp(argv[1], "BBH", 3) == 0)
81          pythia.readString("HiggsBSM:gg2H2bbbar = on");
82        else
83          pythia.readString("HiggsBSM:gg2A3bbbar = on");
188          pythia.readString("25:m0 = 129");
189          pythia.readString("35:m0 = 499.7");
190          pythia.readString("36:m0 = 500");
191          pythia.readString("37:m0 = 506");
192          pythia.readString("HiggsHchg:tanBeta = 30");
193 <        pythia.readString("35:onMode = 0");      // H2(H_2) decays only to tau-pairs
194 <        pythia.readString("35:onIfAny = 15");
195 <        pythia.readString("36:onMode = 0");      // A0(H_3) decays only to tau-pairs
196 <        pythia.readString("36:onIfAny = 15");
197 <      }
198 <    else
199 <      {
200 <        cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
201 <        exit(1);
193 >        if(strncmp(argv[1], "BBH", 3) == 0) {
194 >          pythia.readString("HiggsBSM:gg2H2bbbar = on");
195 >          pythia.readString("35:onMode = 0");      // H2(H_2) decays only to tau-pairs
196 >          pythia.readString("35:onIfAny = 15");
197 >        }
198 >        else {
199 >          pythia.readString("HiggsBSM:gg2A3bbbar = on");
200 >          pythia.readString("36:onMode = 0");      // A0(H_3) decays only to tau-pairs
201 >          pythia.readString("36:onIfAny = 15");
202 >        }
203 >        if(strncmp(argv[1], "BBHmh", 5) == 0 ||
204 >           strncmp(argv[1], "BBAmh", 5) == 0) {
205 >          pythia.readString("15:onMode = 2");       // tau- decays to muons - second tau decays in all modes
206 >          pythia.readString("15:onIfAny = 13");
207 >        }
208        }
209 +     else if(strncmp(argv[1], "ZZmmbb", 6) == 0)
210 +       {
211 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
212 +         pythia.readString("23:onMode = 0");
213 +         pythia.readString("23:onIfAny = 5 13");   // Let the Z decay to bQuarks or muons
214 +       }
215 +     else if(strncmp(argv[1], "ZZmmnn", 6) == 0)
216 +       {
217 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
218 +         pythia.readString("23:onMode = 0");
219 +         pythia.readString("23:onIfAny = 13 12 14 16");   // Let the Z decay to muons or neutrinos
220 +       }
221 +     else if(strncmp(argv[1], "ZZeebb", 6) == 0)
222 +       {
223 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
224 +         pythia.readString("23:onMode = 0");
225 +         pythia.readString("23:onIfAny = 5 11");   // Let the Z decay to bQuarks or electrons
226 +       }
227 +     else if(strncmp(argv[1], "ZZeenn", 6) == 0)
228 +       {
229 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
230 +         pythia.readString("23:onMode = 0");
231 +         pythia.readString("23:onIfAny = 11 12 14 16");   // Let the Z decay to electrons or neutrinos
232 +       }
233 +     else if(strncmp(argv[1], "QCD", 3) == 0)
234 +       {
235 +         pythia.readString("HardQCD:all = on");    
236 +         pythia.readString("PhaseSpace:pTHatMin = 50.");  
237 +       }
238 +     else if(strcmp(argv[1], "QED") == 0)
239 +       {
240 +         pythia.readString("PromptPhoton:all = on");
241 +         pythia.readString("PhaseSpace:mHatMin = 60.");    
242 +         pythia.readString("PhaseSpace:mHatMax = -1");
243 +       }
244 +     else if(strncmp(argv[1], "EWK", 3) == 0)
245 +       {
246 +         pythia.readString("WeakSingleBoson:all = on");
247 +         pythia.readString("PhaseSpace:mHatMin = 60.");    
248 +         pythia.readString("PhaseSpace:mHatMax = -1");    
249 +         pythia.readString("WeakDoubleBoson:all = on");
250 +       }
251 +     else if(strncmp(argv[1], "Top", 3) == 0)
252 +       {
253 +         pythia.readString("Top:all = on");
254 +       }
255 +     else
256 +       {
257 +         cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
258 +         exit(1);
259 +       }
260    }
261  
262    // Initialize pythia
263  
264 <  pythia.init( 2212, 2212, 14000.);
264 >  pythia.init(2212, 2212, cmEnergy);
265    pythia.particleData.listAll();
266    Rndm * rndmPtr = &pythia.rndm;
267    rndmPtr->init(runNumber);
# Line 114 | Line 275 | int main(int argc, char **argv) {
275  
276    // Ultra Fast Simulator
277  
278 <  UltraFastSim ufs(rndmPtr);
278 >  UltraFastSim ufs;
279 >  char jobName[256];
280 >  sprintf(jobName, "%s-PU%3.3d-%4.4d", argv[1], meanPileupEventCount, runNumber);
281 >  UFSDataStore dataStore(jobName, &ufs);
282  
283    // Begin event loop
284    for (int iEvent = 0; iEvent < nEvents; ) {
# Line 122 | Line 286 | int main(int argc, char **argv) {
286      // Generate event. Skip if error. List first one.
287      if (!pythia.next()) continue;
288  
289 +    if(iEvent < 10) pythia.event.list();
290 +
291      // Add pileup
292  
293      int pileupEventCount = 0;
294 <    if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
294 >    if(meanPileupEventCount > 0) pileupEventCount = poisson(meanPileupEventCount, rndmPtr);
295      for (int puEvent = 0; puEvent < pileupEventCount; ) {
296        if(!pileupPythia.next()) continue;
297 +      double vx = rndmPtr->gauss()*2.; // Mean beam spread is 2 mm in x-y and 75 mm in z
298 +      double vy = rndmPtr->gauss()*2.;
299 +      double vz = rndmPtr->gauss()*75.;
300        for(int i = 0; i < pileupPythia.event.size(); i++) {
301          Particle& particle = pileupPythia.event[i];
302 +        particle.xProd(vx+particle.xProd());
303 +        particle.yProd(vy+particle.yProd());
304 +        particle.zProd(vz+particle.zProd());
305          if(particle.status() > 0) {
306            if(particle.isVisible()) {
307              if(particle.pT() > 0.5) {
# Line 143 | Line 315 | int main(int argc, char **argv) {
315  
316      // Ultra fast simulation
317  
318 <    if(!ufs.run(pythia.event))
318 >    if(!ufs.run(pythia.event, rndmPtr))
319        {
320          cerr << "Ultra fast simulation failed - aborting" << endl;
321          exit(1);
322        }
323  
324 <    cout << "Number of Particles = " << ufs.particleList().size() << endl;
325 <    cout << "Number of Gen Elecs = " << ufs.electronList().size() << endl;
326 <    cout << "Number of Gen Muons = " << ufs.muonList().size() << endl;
327 <    cout << "Number of Gen Taus =  " << ufs.genTauList().size() << endl;
328 <    cout << "Number of b Quarks =  " << ufs.bQuarkList().size() << endl;
329 <    cout << "Number of Jets =      " << ufs.sortedJetList().size() << endl;
330 <    cout << "Number of bJets =     " << ufs.bJetList().size() << endl;
331 <    cout << "Number of taus =      " << ufs.tauList().size() << endl;
160 <    cout << endl;
161 <  
324 >    // Store data
325 >
326 >    if(!dataStore.run())
327 >      {
328 >        cerr << "Failed to store data" << endl;
329 >        exit(3);
330 >      }
331 >
332      // End of event loop. Statistics. Histogram. Done.
333  
334      if(!(iEvent % 100)) cout << "Processed event " << iEvent << endl;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines