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.10 by dasu, Fri Feb 18 05:06:02 2011 UTC vs.
Revision 1.22 by dasu, Wed Aug 8 00:12:30 2012 UTC

# Line 13 | Line 13 | using namespace std;
13   #include "Pythia.h"
14   using namespace Pythia8;
15  
16 – #include "TROOT.h"
17 – #include "TFile.h"
18 –
16   #include "UltraFastSim.h"
17 < #include "bbHAnalysis.h"
18 < #include "ZHAnalysis.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 26 | 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 >  bool cutOnMET = false;
48 >  if(argc < 3 || argc > 6)
49      {
50 <      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51 <      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
52 <      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
53 <      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
54 <      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
55 <      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
56 <      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
57 <      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
58 <      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
59 <      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
60 <      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
61 <      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
62 <      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
63 <      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
50 >      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
51 >      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
52 >      cerr << "or              " << argv[0] << " ZHmmXX[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
53 >      cerr << "or              " << argv[0] << " ZHeeXX[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
54 >      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
55 >      cerr << "or              " << argv[0] << " ZZmmnn[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
56 >      cerr << "or              " << argv[0] << " ZZeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
57 >      cerr << "or              " << argv[0] << " ZZeenn[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
58 >      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
59 >      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
60 >      cerr << "or              " << argv[0] << " VBFHXX[..] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
61 >      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
62 >      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
63 >      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
64 >      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
65 >      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
66 >      cerr << "or              " << argv[0] << " Znn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
67 >      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
68 >      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
69 >      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
70 >      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
71 >      cerr << "or              " << argv[0] << " QCD[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
72 >      cerr << "or              " << argv[0] << " QED[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
73 >      cerr << "or              " << argv[0] << " EWK[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
74 >      cerr << "or              " << argv[0] << " Top[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
75 >      cerr << "or              " << argv[0] << " QCDMET[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
76 >      cerr << "or              " << argv[0] << " ZnnMET[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
77 >      cerr << "or              " << argv[0] << " ZZTo4L[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
78        cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
79 +      cerr << "QCD, QED, EWK and Top do not specify the decays -- they are useful for generic BG production" << endl;
80 +      cerr << "QCDMET has MET > 50 GeV to make efficient to study QCD MET BG" << endl;
81 +      cerr << "ZnnMET has MET > 50 GeV to make efficient to study Znn MET BG" << endl;
82        exit(1);
83      }
84    else {
85      nEvents = atoi(argv[2]);
86      if(argc >= 4) runNumber = atoi(argv[3]);
87      if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
88 +    if(argc >= 6) cmEnergy = atof(argv[5]);
89      if(strncmp(argv[1], "Zee", 3) == 0)
90        {
91          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
# Line 65 | Line 102 | int main(int argc, char **argv) {
102          pythia.readString("23:onMode = 0");       // Z decays to muons only
103          pythia.readString("23:onIfAny = 13");
104        }
105 +    else if(strncmp(argv[1], "Znn", 3) == 0)
106 +      {
107 +        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
108 +        pythia.readString("PhaseSpace:mHatMin = 60.");    
109 +        pythia.readString("PhaseSpace:mHatMax = -1");    
110 +        pythia.readString("23:onMode = 0");       // Z decays to neutrinos only
111 +        pythia.readString("23:onIfAny = 12, 14, 16");
112 +        if(strncmp(argv[1], "ZnnMET", 6) == 0) {
113 +          cutOnMET = true;
114 +        }
115 +      }
116      else if(strncmp(argv[1], "Zeh", 3) == 0)
117        {
118          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
# Line 117 | Line 165 | int main(int argc, char **argv) {
165          pythia.readString("24:onMode = 2");       // W+ decays to muon
166          pythia.readString("24:onIfAny = 13");
167        }
168 +    else if(strncmp(argv[1], "HZGamma", 7) == 0)
169 +      {
170 +        pythia.readString("HiggsSM:all = on");
171 +        pythia.readString("25:m0 = 125");
172 +        pythia.readString("25:onMode = 0");      // Higgs decays to Zgamma only
173 +        pythia.readString("25:onIfMatch = 22 23");
174 +        pythia.readString("23:onMode = 0");      // Z decays to light leptons
175 +        pythia.readString("23:onIfAny = 11 13");
176 +      }
177      else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
178        {
179          pythia.readString("HiggsSM:ffbar2HZ = on");
180 <        pythia.readString("25:m0 = 120");
180 >        pythia.readString("25:m0 = 125");
181          pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
182          pythia.readString("25:onIfAny = 5");
183          pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
184          pythia.readString("23:onIfAny = 13");
185        }
186 +     else if(strncmp(argv[1], "ZHmmXX", 6) == 0)
187 +       {
188 +        pythia.readString("HiggsSM:ffbar2HZ = on");
189 +        pythia.readString("25:m0 = 125");
190 +        pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
191 +        pythia.readString("25:isVisible = false");
192 +        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
193 +        pythia.readString("23:onIfAny = 13");
194 +       }
195      else if(strncmp(argv[1], "ZHeebb", 6) == 0)
196        {
197          pythia.readString("HiggsSM:ffbar2HZ = on");
198 <        pythia.readString("25:m0 = 120");
198 >        pythia.readString("25:m0 = 125");
199          pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
200          pythia.readString("25:onIfAny = 5");
201          pythia.readString("23:onMode = 0");      // Z decays to e+,e- only (for trigger)
202          pythia.readString("23:onIfAny = 11");
203        }
204 +     else if(strncmp(argv[1], "ZHeeXX", 6) == 0)
205 +       {
206 +        pythia.readString("HiggsSM:ffbar2HZ = on");
207 +        pythia.readString("25:m0 = 125");
208 +        pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
209 +        pythia.readString("25:isVisible = false");
210 +        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
211 +        pythia.readString("23:onIfAny = 11");
212 +       }
213      else if(strncmp(argv[1], "BB", 2) == 0)
214        {
215          pythia.readString("Higgs:useBSM = on");
# Line 165 | Line 240 | int main(int argc, char **argv) {
240           pythia.readString("23:onMode = 0");
241           pythia.readString("23:onIfAny = 5 13");   // Let the Z decay to bQuarks or muons
242         }
243 <    else
244 <      {
245 <        cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
246 <        exit(1);
247 <      }
243 >     else if(strncmp(argv[1], "ZZTo4L", 6) == 0)
244 >       {
245 >         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
246 >         pythia.readString("23:onMode = 0");
247 >         pythia.readString("23:onIfAny = 11 13 15");   // Let the Z decay to charged lepton pairs
248 >       }
249 >     else if(strncmp(argv[1], "ZZmmnn", 6) == 0)
250 >       {
251 >         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
252 >         pythia.readString("23:onMode = 0");
253 >         pythia.readString("23:onIfAny = 13 12 14 16");   // Let the Z decay to muons or neutrinos
254 >       }
255 >     else if(strncmp(argv[1], "ZZeebb", 6) == 0)
256 >       {
257 >         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
258 >         pythia.readString("23:onMode = 0");
259 >         pythia.readString("23:onIfAny = 5 11");   // Let the Z decay to bQuarks or electrons
260 >       }
261 >     else if(strncmp(argv[1], "ZZeenn", 6) == 0)
262 >       {
263 >         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
264 >         pythia.readString("23:onMode = 0");
265 >         pythia.readString("23:onIfAny = 11 12 14 16");   // Let the Z decay to electrons or neutrinos
266 >       }
267 >     else if(strncmp(argv[1], "QCD", 3) == 0)
268 >       {
269 >         pythia.readString("HardQCD:all = on");    
270 >         pythia.readString("PhaseSpace:pTHatMin = 50.");
271 >         if(strncmp(argv[1], "QCDMET", 6) == 0) {
272 >           cutOnMET = true;
273 >         }
274 >       }
275 >     else if(strcmp(argv[1], "QED") == 0)
276 >       {
277 >         pythia.readString("PromptPhoton:all = on");
278 >         pythia.readString("PhaseSpace:mHatMin = 60.");    
279 >         pythia.readString("PhaseSpace:mHatMax = -1");
280 >       }
281 >     else if(strncmp(argv[1], "EWK", 3) == 0)
282 >       {
283 >         pythia.readString("WeakSingleBoson:all = on");
284 >         pythia.readString("PhaseSpace:mHatMin = 60.");    
285 >         pythia.readString("PhaseSpace:mHatMax = -1");    
286 >         pythia.readString("WeakDoubleBoson:all = on");
287 >       }
288 >     else if(strncmp(argv[1], "Top", 3) == 0)
289 >       {
290 >         pythia.readString("Top:all = on");
291 >       }
292 >     else if(strncmp(argv[1], "VBFHXX", 6) == 0)
293 >       {
294 >         pythia.readString("HiggsSM:ff2Hff(t:ZZ) = on");
295 >         pythia.readString("HiggsSM:ff2Hff(t:WW) = on");
296 >         pythia.readString("25:m0 = 125");
297 >         pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
298 >         pythia.readString("25:isVisible = false");
299 >       }
300 >     else
301 >       {
302 >         cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
303 >         exit(1);
304 >       }
305    }
306  
175 –  string name(argv[1]);
176 –  name += ".root";
177 –  TFile *outFile = TFile::Open(name.c_str(), "recreate");
178 –
307    // Initialize pythia
308  
309 <  pythia.init( 2212, 2212, 14000.);
309 >  pythia.init(2212, 2212, cmEnergy);
310    pythia.particleData.listAll();
311    Rndm * rndmPtr = &pythia.rndm;
312    rndmPtr->init(runNumber);
# Line 192 | Line 320 | int main(int argc, char **argv) {
320  
321    // Ultra Fast Simulator
322  
323 <  UltraFastSim ufs(rndmPtr);
324 <  bbHAnalysis bbH(outFile, &pythia, &ufs, true);
325 <  ZHAnalysis ZH(outFile, &pythia, &ufs, true);
323 >  UltraFastSim ufs;
324 >  char jobName[256];
325 >  sprintf(jobName, "%s-PU%3.3d-%4.4d", argv[1], meanPileupEventCount, runNumber);
326 >  UFSDataStore dataStore(jobName, &ufs);
327  
328    // Begin event loop
329    for (int iEvent = 0; iEvent < nEvents; ) {
# Line 207 | Line 336 | int main(int argc, char **argv) {
336      // Add pileup
337  
338      int pileupEventCount = 0;
339 <    if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
339 >    if(meanPileupEventCount > 0) pileupEventCount = poisson(meanPileupEventCount, rndmPtr);
340      for (int puEvent = 0; puEvent < pileupEventCount; ) {
341        if(!pileupPythia.next()) continue;
342 +      double vx = rndmPtr->gauss()*2.; // Mean beam spread is 2 mm in x-y and 75 mm in z
343 +      double vy = rndmPtr->gauss()*2.;
344 +      double vz = rndmPtr->gauss()*75.;
345        for(int i = 0; i < pileupPythia.event.size(); i++) {
346          Particle& particle = pileupPythia.event[i];
347 +        particle.xProd(vx+particle.xProd());
348 +        particle.yProd(vy+particle.yProd());
349 +        particle.zProd(vz+particle.zProd());
350          if(particle.status() > 0) {
351            if(particle.isVisible()) {
352              if(particle.pT() > 0.5) {
# Line 225 | Line 360 | int main(int argc, char **argv) {
360  
361      // Ultra fast simulation
362  
363 <    if(!ufs.run(pythia.event))
363 >    if(!ufs.run(pythia.event, rndmPtr))
364        {
365          cerr << "Ultra fast simulation failed - aborting" << endl;
366          exit(1);
367        }
368  
369 <    // Run bbHAnalysis
369 >    // If necessary, check MET and decide to store or not
370  
371 <    if(!bbH.run())
372 <      {
373 <        cerr << "bbH Analysis failed - aborting" << endl;
374 <        exit(2);
371 >    bool doStore = true;
372 >    if(cutOnMET) {
373 >      double MET = ufs.getMET().Pt();
374 >      if(MET < 50.) {
375 >        doStore = false;
376        }
377 +    }
378 +    // Store data
379  
380 <    // Run ZHAnalysis
243 <
244 <    if(!ZH.run())
380 >    if(doStore && !dataStore.run())
381        {
382 <        cerr << "ZH Analysis failed - aborting" << endl;
382 >        cerr << "Failed to store data" << endl;
383          exit(3);
384        }
385  
# Line 259 | Line 395 | int main(int argc, char **argv) {
395    pythia.statistics();
396    pileupPythia.statistics();
397  
262 –  bbH.end();
263 –  ZH.end();
264 –
265 –  outFile->cd();
266 –  outFile->Write();
267 –  outFile->Close();
268 –
398    return 0;
399  
400   }

Diff Legend

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