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.11 by dasu, Thu Feb 24 07:34:29 2011 UTC vs.
Revision 1.23 by dasu, Sun Sep 9 17:43:30 2012 UTC

# Line 1 | Line 1
1   // This simple program generates various signal processes for
2 < // MSSM bbPhi process and its backgrounds.
2 > // various physics processes.
3   // It includes the mixing of required level of pileup.
4   // It then takes pythia output and runs it through very simplistic
5   // detector simulation using parameterizations, including making
# 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 >  bool lhaFile = false;
49 >  if(argc < 3 || argc > 6)
50      {
51 <      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
52 <      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
53 <      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
54 <      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
55 <      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
56 <      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
57 <      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
58 <      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
59 <      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
60 <      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
61 <      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
62 <      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
63 <      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
64 <      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51 >      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
52 >      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
53 >      cerr << "or              " << argv[0] << " ZHmmXX[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
54 >      cerr << "or              " << argv[0] << " ZHeeXX[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
55 >      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
56 >      cerr << "or              " << argv[0] << " ZZmmnn[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
57 >      cerr << "or              " << argv[0] << " ZZeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
58 >      cerr << "or              " << argv[0] << " ZZeenn[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
59 >      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
60 >      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
61 >      cerr << "or              " << argv[0] << " VBFHXX[..] " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
62 >      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
63 >      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
64 >      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
65 >      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
66 >      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
67 >      cerr << "or              " << argv[0] << " Znn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
68 >      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
69 >      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
70 >      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
71 >      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
72 >      cerr << "or              " << argv[0] << " QCD[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
73 >      cerr << "or              " << argv[0] << " QED[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
74 >      cerr << "or              " << argv[0] << " EWK[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
75 >      cerr << "or              " << argv[0] << " Top[...]   " << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
76 >      cerr << "or              " << argv[0] << " QCDMET[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
77 >      cerr << "or              " << argv[0] << " ZnnMET[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
78 >      cerr << "or              " << argv[0] << " ZZTo4L[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
79 >      cerr << "or              " << argv[0] << " HZZTo4L[...]" << " <nEvents> [RandomSeed] [meanPileupEventcount] [cmEnergy]" << endl;
80        cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
81 +      cerr << "QCD, QED, EWK and Top do not specify the decays -- they are useful for generic BG production" << endl;
82 +      cerr << "QCDMET has MET > 50 GeV to make efficient to study QCD MET BG" << endl;
83 +      cerr << "ZnnMET has MET > 50 GeV to make efficient to study Znn MET BG" << endl;
84        exit(1);
85      }
86    else {
87      nEvents = atoi(argv[2]);
88      if(argc >= 4) runNumber = atoi(argv[3]);
89      if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
90 +    if(argc >= 6) cmEnergy = atof(argv[5]);
91      if(strncmp(argv[1], "Zee", 3) == 0)
92        {
93          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
# Line 65 | Line 104 | int main(int argc, char **argv) {
104          pythia.readString("23:onMode = 0");       // Z decays to muons only
105          pythia.readString("23:onIfAny = 13");
106        }
107 +    else if(strncmp(argv[1], "Znn", 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 neutrinos only
113 +        pythia.readString("23:onIfAny = 12, 14, 16");
114 +        if(strncmp(argv[1], "ZnnMET", 6) == 0) {
115 +          cutOnMET = true;
116 +        }
117 +      }
118      else if(strncmp(argv[1], "Zeh", 3) == 0)
119        {
120          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
# Line 117 | Line 167 | int main(int argc, char **argv) {
167          pythia.readString("24:onMode = 2");       // W+ decays to muon
168          pythia.readString("24:onIfAny = 13");
169        }
170 +    else if(strncmp(argv[1], "HZGamma", 7) == 0)
171 +      {
172 +        pythia.readString("HiggsSM:all = on");
173 +        pythia.readString("25:m0 = 125");
174 +        pythia.readString("25:onMode = 0");      // Higgs decays to Zgamma only
175 +        pythia.readString("25:onIfMatch = 22 23");
176 +        pythia.readString("23:onMode = 0");      // Z decays to light leptons
177 +        pythia.readString("23:onIfAny = 11 13");
178 +      }
179      else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
180        {
181          pythia.readString("HiggsSM:ffbar2HZ = on");
182 <        pythia.readString("25:m0 = 120");
182 >        pythia.readString("25:m0 = 125");
183          pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
184          pythia.readString("25:onIfAny = 5");
185          pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
186          pythia.readString("23:onIfAny = 13");
187        }
188 +     else if(strncmp(argv[1], "ZHmmXX", 6) == 0)
189 +       {
190 +        pythia.readString("HiggsSM:ffbar2HZ = on");
191 +        pythia.readString("25:m0 = 125");
192 +        pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
193 +        pythia.readString("25:isVisible = false");
194 +        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
195 +        pythia.readString("23:onIfAny = 13");
196 +       }
197      else if(strncmp(argv[1], "ZHeebb", 6) == 0)
198        {
199          pythia.readString("HiggsSM:ffbar2HZ = on");
200 <        pythia.readString("25:m0 = 120");
200 >        pythia.readString("25:m0 = 125");
201          pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
202          pythia.readString("25:onIfAny = 5");
203          pythia.readString("23:onMode = 0");      // Z decays to e+,e- only (for trigger)
204          pythia.readString("23:onIfAny = 11");
205        }
206 +     else if(strncmp(argv[1], "ZHeeXX", 6) == 0)
207 +       {
208 +        pythia.readString("HiggsSM:ffbar2HZ = on");
209 +        pythia.readString("25:m0 = 125");
210 +        pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
211 +        pythia.readString("25:isVisible = false");
212 +        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
213 +        pythia.readString("23:onIfAny = 11");
214 +       }
215 +     else if(strncmp(argv[1], "H+tb", 4) == 0)
216 +      {
217 +        pythia.readString("Higgs:useBSM = on");        // Select these consistently for (MA,tanBeta)
218 +        pythia.readString("25:m0 = 129");              // Light higgs
219 +        pythia.readString("35:m0 = 499.7");            // Heavy CP-Even higgs
220 +        pythia.readString("36:m0 = 500");              // Heavy CP-Odd higgs
221 +        pythia.readString("37:m0 = 506");              // Heavy charged higgs
222 +        pythia.readString("HiggsHchg:tanBeta = 30");   // tanBeta
223 +        pythia.readString("HiggsBSM:allH+- = on");     // Turn on charged higgs produciton
224 +        pythia.readString("37:onMode = 0");            // Turn on all decays
225 +        pythia.readString("37:onIfAny = 6");           // Turn on decays to top (and b)
226 +      }
227      else if(strncmp(argv[1], "BB", 2) == 0)
228        {
229          pythia.readString("Higgs:useBSM = on");
# Line 165 | Line 254 | int main(int argc, char **argv) {
254           pythia.readString("23:onMode = 0");
255           pythia.readString("23:onIfAny = 5 13");   // Let the Z decay to bQuarks or muons
256         }
257 <    else
257 >     else if(strncmp(argv[1], "ZZTo4L", 6) == 0)
258 >       {
259 >         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
260 >         pythia.readString("23:onMode = 0");
261 >         pythia.readString("23:onIfAny = 11 13 15");   // Let the Z decay to charged lepton pairs
262 >       }
263 >     else if(strncmp(argv[1], "HZZTo4L", 7) == 0)
264        {
265 <        cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
266 <        exit(1);
265 >        pythia.readString("HiggsSM:all = on");
266 >        pythia.readString("25:m0 = 125");
267 >        pythia.readString("25:onMode = 0");      // Higgs decays to ZZ only
268 >        pythia.readString("25:onIfMatch = 23 23");
269 >        pythia.readString("23:onMode = 0");      // Z decays to light leptons
270 >        pythia.readString("23:onIfAny = 11 13");
271        }
272 +     else if(strncmp(argv[1], "ZZmmnn", 6) == 0)
273 +       {
274 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
275 +         pythia.readString("23:onMode = 0");
276 +         pythia.readString("23:onIfAny = 13 12 14 16");   // Let the Z decay to muons or neutrinos
277 +       }
278 +     else if(strncmp(argv[1], "ZZeebb", 6) == 0)
279 +       {
280 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
281 +         pythia.readString("23:onMode = 0");
282 +         pythia.readString("23:onIfAny = 5 11");   // Let the Z decay to bQuarks or electrons
283 +       }
284 +     else if(strncmp(argv[1], "ZZeenn", 6) == 0)
285 +       {
286 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
287 +         pythia.readString("23:onMode = 0");
288 +         pythia.readString("23:onIfAny = 11 12 14 16");   // Let the Z decay to electrons or neutrinos
289 +       }
290 +     else if(strncmp(argv[1], "QCD", 3) == 0)
291 +       {
292 +         pythia.readString("HardQCD:all = on");    
293 +         pythia.readString("PhaseSpace:pTHatMin = 50.");
294 +         if(strncmp(argv[1], "QCDMET", 6) == 0) {
295 +           cutOnMET = true;
296 +         }
297 +       }
298 +     else if(strcmp(argv[1], "QED") == 0)
299 +       {
300 +         pythia.readString("PromptPhoton:all = on");
301 +         pythia.readString("PhaseSpace:mHatMin = 60.");    
302 +         pythia.readString("PhaseSpace:mHatMax = -1");
303 +       }
304 +     else if(strncmp(argv[1], "EWK", 3) == 0)
305 +       {
306 +         pythia.readString("WeakSingleBoson:all = on");
307 +         pythia.readString("PhaseSpace:mHatMin = 60.");    
308 +         pythia.readString("PhaseSpace:mHatMax = -1");    
309 +         pythia.readString("WeakDoubleBoson:all = on");
310 +       }
311 +     else if(strncmp(argv[1], "Top", 3) == 0)
312 +       {
313 +         pythia.readString("Top:all = on");
314 +       }
315 +     else if(strncmp(argv[1], "VBFHXX", 6) == 0)
316 +       {
317 +         pythia.readString("HiggsSM:ff2Hff(t:ZZ) = on");
318 +         pythia.readString("HiggsSM:ff2Hff(t:WW) = on");
319 +         pythia.readString("25:m0 = 125");
320 +         pythia.readString("25:mayDecay = false");      // Higgs does not decay -- becomes invisible
321 +         pythia.readString("25:isVisible = false");
322 +       }
323 +     else if(fopen(argv[1], "r") != NULL)
324 +       {
325 +         lhaFile = true;
326 +       }
327 +     else
328 +       {
329 +         cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
330 +         exit(1);
331 +       }
332    }
333  
175  string name(argv[1]);
176  name += ".root";
177  TFile *outFile = TFile::Open(name.c_str(), "recreate");
178
334    // Initialize pythia
335  
336 <  pythia.init( 2212, 2212, 14000.);
336 >  if(lhaFile) {
337 >    pythia.init(argv[1]);
338 >  }
339 >  else {
340 >    pythia.init(2212, 2212, cmEnergy);
341 >  }
342    pythia.particleData.listAll();
343    Rndm * rndmPtr = &pythia.rndm;
344    rndmPtr->init(runNumber);
# Line 187 | Line 347 | int main(int argc, char **argv) {
347  
348    Pythia pileupPythia;
349    pileupPythia.readString("SoftQCD:minBias = on");
350 <  pileupPythia.init( 2212, 2212, 14000.);
350 >  pileupPythia.init( 2212, 2212, cmEnergy);
351    pileupPythia.rndm.init(runNumber);
352  
353    // Ultra Fast Simulator
354  
355 <  UltraFastSim ufs(rndmPtr);
356 <  bbHAnalysis bbH(outFile, &pythia, &ufs, true);
357 <  ZHAnalysis ZH; // outFile, &pythia, &ufs, true);
358 <  ZH.BookTree();
355 >  UltraFastSim ufs;
356 >  char jobName[256];
357 >  sprintf(jobName, "%s-PU%3.3d-%4.4d", argv[1], meanPileupEventCount, runNumber);
358 >  UFSDataStore dataStore(jobName, &ufs);
359  
360    // Begin event loop
361    for (int iEvent = 0; iEvent < nEvents; ) {
# Line 208 | Line 368 | int main(int argc, char **argv) {
368      // Add pileup
369  
370      int pileupEventCount = 0;
371 <    if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
371 >    if(meanPileupEventCount > 0) pileupEventCount = poisson(meanPileupEventCount, rndmPtr);
372      for (int puEvent = 0; puEvent < pileupEventCount; ) {
373        if(!pileupPythia.next()) continue;
374        double vx = rndmPtr->gauss()*2.; // Mean beam spread is 2 mm in x-y and 75 mm in z
# Line 232 | Line 392 | int main(int argc, char **argv) {
392  
393      // Ultra fast simulation
394  
395 <    if(!ufs.run(pythia.event))
395 >    if(!ufs.run(pythia.event, rndmPtr))
396        {
397          cerr << "Ultra fast simulation failed - aborting" << endl;
398          exit(1);
399        }
400  
401 <    // Run bbHAnalysis
401 >    // If necessary, check MET and decide to store or not
402  
403 <    if(!bbH.run())
404 <      {
405 <        cerr << "bbH Analysis failed - aborting" << endl;
406 <        exit(2);
403 >    bool doStore = true;
404 >    if(cutOnMET) {
405 >      double MET = ufs.getMET().Pt();
406 >      if(MET < 50.) {
407 >        doStore = false;
408        }
409 +    }
410 +    // Store data
411  
412 <    // Run ZHAnalysis
250 <
251 <    if(!ZH.Run(&ufs))
412 >    if(doStore && !dataStore.run())
413        {
414 <        cerr << "ZH Analysis failed - aborting" << endl;
414 >        cerr << "Failed to store data" << endl;
415          exit(3);
416        }
417  
# Line 266 | Line 427 | int main(int argc, char **argv) {
427    pythia.statistics();
428    pileupPythia.statistics();
429  
269  bbH.end();
270
271  outFile->cd();
272  outFile->Write();
273  outFile->Close();
274
430    return 0;
431  
432   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines