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.14 by dasu, Fri Feb 25 17:10:39 2011 UTC vs.
Revision 1.18 by dasu, Wed Nov 2 21:18:08 2011 UTC

# Line 16 | Line 16 | using namespace Pythia8;
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.
42    Pythia pythia;
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] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
50 <      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51 <      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
52 <      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
53 <      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
54 <      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
55 <      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
56 <      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
57 <      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
58 <      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
59 <      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
60 <      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
61 <      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
62 <      cerr << "or              " << argv[0] << " Tmn[...]   " << " <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(argc >= 6) cmEnergy = atof(argv[5]);
81      if(strncmp(argv[1], "Zee", 3) == 0)
82        {
83          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
# Line 122 | 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");
# Line 131 | Line 173 | int main(int argc, char **argv) {
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");
# Line 161 | Line 212 | int main(int argc, char **argv) {
212           pythia.readString("23:onMode = 0");
213           pythia.readString("23:onIfAny = 5 13");   // Let the Z decay to bQuarks or muons
214         }
215 <    else
216 <      {
217 <        cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
218 <        exit(1);
219 <      }
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 185 | Line 276 | int main(int argc, char **argv) {
276    // Ultra Fast Simulator
277  
278    UltraFastSim ufs;
279 <  UFSDataStore dataStore(argv[1], &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 198 | Line 291 | int main(int argc, char **argv) {
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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines