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.5 by dasu, Tue Feb 15 06:49:11 2011 UTC vs.
Revision 1.16 by dasu, Wed Mar 2 21:57:09 2011 UTC

# Line 14 | Line 14 | using namespace std;
14   using namespace Pythia8;
15  
16   #include "UltraFastSim.h"
17 < #include "Analysis.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 24 | Line 45 | int main(int argc, char **argv) {
45    int meanPileupEventCount = 0;
46    if(argc < 3 || argc > 5)
47      {
48 <      cerr << "Command syntax: " << argv[0] << " BBH[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
49 <      cerr << "Command syntax: " << argv[0] << " BBA[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
50 <      cerr << "or              " << argv[0] << " ZLL[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51 <      cerr << "or              " << argv[0] << " ZMM[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
52 <      cerr << "or              " << argv[0] << " TLN[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
53 <      cerr << "or              " << argv[0] << " WLN[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
48 >      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
49 >      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
50 >      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51 >      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
52 >      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
53 >      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
54 >      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
55 >      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
56 >      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
57 >      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
58 >      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
59 >      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
60 >      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
61 >      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
62 >      cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
63        exit(1);
64      }
65    else {
66      nEvents = atoi(argv[2]);
67      if(argc >= 4) runNumber = atoi(argv[3]);
68      if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
69 <    if(strncmp(argv[1], "ZLL", 3) == 0)
69 >    if(strncmp(argv[1], "Zee", 3) == 0)
70        {
71          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
72          pythia.readString("PhaseSpace:mHatMin = 60.");    
73          pythia.readString("PhaseSpace:mHatMax = -1");    
74 <        pythia.readString("23:onMode = 0");       // Z decays to leptons only
75 <        pythia.readString("23:onIfAny = 11 13 15");
46 <        pythia.readString("15:onMode = 2");       // tau- decays to electrons and muons - second tau decays in all modes
47 <        pythia.readString("15:onIfAny = 11 13");
74 >        pythia.readString("23:onMode = 0");       // Z decays to electrons only
75 >        pythia.readString("23:onIfAny = 11");
76        }
77 <    else if(strncmp(argv[1], "ZMM", 3) == 0)
77 >    else if(strncmp(argv[1], "Zmm", 3) == 0)
78        {
79          pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
80          pythia.readString("PhaseSpace:mHatMin = 60.");    
# Line 54 | Line 82 | int main(int argc, char **argv) {
82          pythia.readString("23:onMode = 0");       // Z decays to muons only
83          pythia.readString("23:onIfAny = 13");
84        }
85 <    else if(strncmp(argv[1], "WLN", 3) == 0)
85 >    else if(strncmp(argv[1], "Zeh", 3) == 0)
86 >      {
87 >        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
88 >        pythia.readString("PhaseSpace:mHatMin = 60.");    
89 >        pythia.readString("PhaseSpace:mHatMax = -1");    
90 >        pythia.readString("23:onMode = 0");       // Z decays to taus only
91 >        pythia.readString("23:onIfAny = 15");
92 >        pythia.readString("15:onMode = 2");       // tau- decays to electrons - second tau decays in all modes
93 >        pythia.readString("15:onIfAny = 11");
94 >      }
95 >    else if(strncmp(argv[1], "Zmh", 3) == 0)
96 >      {
97 >        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
98 >        pythia.readString("PhaseSpace:mHatMin = 60.");    
99 >        pythia.readString("PhaseSpace:mHatMax = -1");    
100 >        pythia.readString("23:onMode = 0");       // Z decays to taus only
101 >        pythia.readString("23:onIfAny = 15");
102 >        pythia.readString("15:onMode = 2");       // tau- decays to muons - second tau decays in all modes
103 >        pythia.readString("15:onIfAny = 13");
104 >      }
105 >    else if(strncmp(argv[1], "Zhh", 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 taus only
111 >        pythia.readString("23:onIfAny = 15");
112 >      }
113 >    else if(strncmp(argv[1], "Wen", 3) == 0)
114 >      {
115 >        pythia.readString("WeakSingleBoson:ffbar2W = on");
116 >        pythia.readString("24:onMode = 0");       // W decays to electron
117 >        pythia.readString("24:onIfAny = 11");
118 >      }
119 >    else if(strncmp(argv[1], "Wmn", 3) == 0)
120        {
121          pythia.readString("WeakSingleBoson:ffbar2W = on");
122 <        pythia.readString("24:onMode = 0");       // W decays to light leptons only
123 <        pythia.readString("24:onIfAny = 11 13");
122 >        pythia.readString("24:onMode = 0");       // W decays to muon
123 >        pythia.readString("24:onIfAny = 13");
124 >      }
125 >    else if(strncmp(argv[1], "Ten", 3) == 0)
126 >      {
127 >        pythia.readString("Top:all = on");
128 >        pythia.readString("24:onMode = 2");       // 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], "Tmn", 3) == 0)
132        {
133          pythia.readString("Top:all = on");
134 <        pythia.readString("24:onMode = 2");       // W+ decays to light leptons only
135 <        pythia.readString("24:onIfAny = 11 13");
134 >        pythia.readString("24:onMode = 2");       // W+ decays to muon
135 >        pythia.readString("24:onIfAny = 13");
136        }
137 <    else if(strncmp(argv[1], "ZH", 2) == 0)
137 >    else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
138        {
139          pythia.readString("HiggsSM:ffbar2HZ = on");
140          pythia.readString("25:m0 = 120");
# Line 75 | Line 143 | int main(int argc, char **argv) {
143          pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
144          pythia.readString("23:onIfAny = 13");
145        }
146 +    else if(strncmp(argv[1], "ZHeebb", 6) == 0)
147 +      {
148 +        pythia.readString("HiggsSM:ffbar2HZ = on");
149 +        pythia.readString("25:m0 = 120");
150 +        pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
151 +        pythia.readString("25:onIfAny = 5");
152 +        pythia.readString("23:onMode = 0");      // Z decays to e+,e- only (for trigger)
153 +        pythia.readString("23:onIfAny = 11");
154 +      }
155      else if(strncmp(argv[1], "BB", 2) == 0)
156        {
157          pythia.readString("Higgs:useBSM = on");
81 –        if(strncmp(argv[1], "BBH", 3) == 0)
82 –          pythia.readString("HiggsBSM:gg2H2bbbar = on");
83 –        else
84 –          pythia.readString("HiggsBSM:gg2A3bbbar = on");
158          pythia.readString("25:m0 = 129");
159          pythia.readString("35:m0 = 499.7");
160          pythia.readString("36:m0 = 500");
161          pythia.readString("37:m0 = 506");
162          pythia.readString("HiggsHchg:tanBeta = 30");
163 <        pythia.readString("35:onMode = 0");      // H2(H_2) decays only to tau-pairs
164 <        pythia.readString("35:onIfAny = 15");
165 <        pythia.readString("36:onMode = 0");      // A0(H_3) decays only to tau-pairs
166 <        pythia.readString("36:onIfAny = 15");
163 >        if(strncmp(argv[1], "BBH", 3) == 0) {
164 >          pythia.readString("HiggsBSM:gg2H2bbbar = on");
165 >          pythia.readString("35:onMode = 0");      // H2(H_2) decays only to tau-pairs
166 >          pythia.readString("35:onIfAny = 15");
167 >        }
168 >        else {
169 >          pythia.readString("HiggsBSM:gg2A3bbbar = on");
170 >          pythia.readString("36:onMode = 0");      // A0(H_3) decays only to tau-pairs
171 >          pythia.readString("36:onIfAny = 15");
172 >        }
173 >        if(strncmp(argv[1], "BBHmh", 5) == 0 ||
174 >           strncmp(argv[1], "BBAmh", 5) == 0) {
175 >          pythia.readString("15:onMode = 2");       // tau- decays to muons - second tau decays in all modes
176 >          pythia.readString("15:onIfAny = 13");
177 >        }
178        }
179 +     else if(strncmp(argv[1], "ZZmmbb", 6) == 0)
180 +       {
181 +         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
182 +         pythia.readString("23:onMode = 0");
183 +         pythia.readString("23:onIfAny = 5 13");   // Let the Z decay to bQuarks or muons
184 +       }
185      else
186        {
187          cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
# Line 115 | Line 205 | int main(int argc, char **argv) {
205  
206    // Ultra Fast Simulator
207  
208 <  UltraFastSim ufs(rndmPtr);
209 <  Analysis analysis(argv[1], &ufs);
208 >  UltraFastSim ufs;
209 >  char jobName[256];
210 >  sprintf(jobName, "%s-PU%3.3d-%4.4d", argv[1], meanPileupEventCount, runNumber);
211 >  UFSDataStore dataStore(jobName, &ufs);
212  
213    // Begin event loop
214    for (int iEvent = 0; iEvent < nEvents; ) {
# Line 129 | Line 221 | int main(int argc, char **argv) {
221      // Add pileup
222  
223      int pileupEventCount = 0;
224 <    if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
224 >    if(meanPileupEventCount > 0) pileupEventCount = poisson(meanPileupEventCount, rndmPtr);
225      for (int puEvent = 0; puEvent < pileupEventCount; ) {
226        if(!pileupPythia.next()) continue;
227 +      double vx = rndmPtr->gauss()*2.; // Mean beam spread is 2 mm in x-y and 75 mm in z
228 +      double vy = rndmPtr->gauss()*2.;
229 +      double vz = rndmPtr->gauss()*75.;
230        for(int i = 0; i < pileupPythia.event.size(); i++) {
231          Particle& particle = pileupPythia.event[i];
232 +        particle.xProd(vx+particle.xProd());
233 +        particle.yProd(vy+particle.yProd());
234 +        particle.zProd(vz+particle.zProd());
235          if(particle.status() > 0) {
236            if(particle.isVisible()) {
237              if(particle.pT() > 0.5) {
# Line 147 | Line 245 | int main(int argc, char **argv) {
245  
246      // Ultra fast simulation
247  
248 <    if(!ufs.run(pythia.event))
248 >    if(!ufs.run(pythia.event, rndmPtr))
249        {
250          cerr << "Ultra fast simulation failed - aborting" << endl;
251          exit(1);
252        }
253  
254 <    cout << "Number of Gen Chrgd = " << ufs.chargedHadronList().size() << endl;
157 <    cout << "Number of Gen Neutr = " << ufs.neutralHadronList().size() << endl;
158 <    cout << "Number of Gen Phtns = " << ufs.photonList().size() << endl;
159 <    cout << "Number of Gen Elcns = " << ufs.electronList().size() << endl;
160 <    cout << "Number of Gen Muons = " << ufs.muonList().size() << endl;
161 <    cout << "Number of Gen Taus =  " << ufs.genTauList().size() << endl;
162 <    cout << "Number of taus =      " << ufs.tauList().size() << endl;
163 <    cout << "Number of c Quarks =  " << ufs.cQuarkList().size() << endl;
164 <    cout << "Number of b Quarks =  " << ufs.bQuarkList().size() << endl;
165 <    cout << "Number of Jets =      " << ufs.jetList().size() << endl;
166 <    cout << "Number of bJets =     " << ufs.bJetList().size() << endl;
167 <    cout << endl;
168 <
169 <    // Run analysis
254 >    // Store data
255  
256 <    if(!analysis.run())
256 >    if(!dataStore.run())
257        {
258 <        cerr << "Analysis failed - aborting" << endl;
259 <        exit(2);
258 >        cerr << "Failed to store data" << endl;
259 >        exit(3);
260        }
261  
262      // End of event loop. Statistics. Histogram. Done.

Diff Legend

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