ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dasu/UltraFastSim/generateData.cc
Revision: 1.5
Committed: Tue Feb 15 06:49:11 2011 UTC (14 years, 3 months ago) by dasu
Content type: text/plain
Branch: MAIN
Changes since 1.4: +18 -5 lines
Log Message: 
Added Mike's analysis -- which doesn't work yet :(

File Contents

# User Rev Content
1 dasu 1.1 // This simple program generates various signal processes for
2     // MSSM bbPhi process and its backgrounds.
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
6     // isolated light leptons, photons, composite tau, jet and b-tagged
7     // jet objects.
8     // It saves the output in root format for later analysis.
9    
10     using namespace std;
11     #include <math.h>
12    
13     #include "Pythia.h"
14     using namespace Pythia8;
15    
16     #include "UltraFastSim.h"
17 dasu 1.5 #include "Analysis.h"
18 dasu 1.1
19     int main(int argc, char **argv) {
20     // Generator. Process selection. LHC initialization. Histogram.
21     Pythia pythia;
22     int nEvents = 0;
23     int runNumber = 0;
24     int meanPileupEventCount = 0;
25     if(argc < 3 || argc > 5)
26     {
27 dasu 1.3 cerr << "Command syntax: " << argv[0] << " BBH[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
28 dasu 1.1 cerr << "Command syntax: " << argv[0] << " BBA[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
29     cerr << "or " << argv[0] << " ZLL[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
30 dasu 1.3 cerr << "or " << argv[0] << " ZMM[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
31 dasu 1.1 cerr << "or " << argv[0] << " TLN[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
32     cerr << "or " << argv[0] << " WLN[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
33     exit(1);
34     }
35     else {
36     nEvents = atoi(argv[2]);
37     if(argc >= 4) runNumber = atoi(argv[3]);
38     if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
39     if(strncmp(argv[1], "ZLL", 3) == 0)
40     {
41     pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
42     pythia.readString("PhaseSpace:mHatMin = 60.");
43     pythia.readString("PhaseSpace:mHatMax = -1");
44     pythia.readString("23:onMode = 0"); // Z decays to leptons only
45     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");
48     }
49 dasu 1.3 else if(strncmp(argv[1], "ZMM", 3) == 0)
50     {
51     pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
52     pythia.readString("PhaseSpace:mHatMin = 60.");
53     pythia.readString("PhaseSpace:mHatMax = -1");
54     pythia.readString("23:onMode = 0"); // Z decays to muons only
55     pythia.readString("23:onIfAny = 13");
56     }
57 dasu 1.1 else if(strncmp(argv[1], "WLN", 3) == 0)
58     {
59     pythia.readString("WeakSingleBoson:ffbar2W = on");
60     pythia.readString("24:onMode = 0"); // W decays to light leptons only
61     pythia.readString("24:onIfAny = 11 13");
62     }
63     else if(strncmp(argv[1], "TLN", 3) == 0)
64     {
65     pythia.readString("Top:all = on");
66     pythia.readString("24:onMode = 2"); // W+ decays to light leptons only
67     pythia.readString("24:onIfAny = 11 13");
68     }
69 dasu 1.3 else if(strncmp(argv[1], "ZH", 2) == 0)
70     {
71     pythia.readString("HiggsSM:ffbar2HZ = on");
72     pythia.readString("25:m0 = 120");
73     pythia.readString("25:onMode = 0"); // Higgs decays to bbBar only
74     pythia.readString("25:onIfAny = 5");
75     pythia.readString("23:onMode = 0"); // Z decays to mu+,mu- only (for trigger)
76     pythia.readString("23:onIfAny = 13");
77     }
78     else if(strncmp(argv[1], "BB", 2) == 0)
79 dasu 1.1 {
80     pythia.readString("Higgs:useBSM = on");
81 dasu 1.3 if(strncmp(argv[1], "BBH", 3) == 0)
82     pythia.readString("HiggsBSM:gg2H2bbbar = on");
83     else
84     pythia.readString("HiggsBSM:gg2A3bbbar = on");
85     pythia.readString("25:m0 = 129");
86     pythia.readString("35:m0 = 499.7");
87     pythia.readString("36:m0 = 500");
88     pythia.readString("37:m0 = 506");
89     pythia.readString("HiggsHchg:tanBeta = 30");
90     pythia.readString("35:onMode = 0"); // H2(H_2) decays only to tau-pairs
91     pythia.readString("35:onIfAny = 15");
92 dasu 1.1 pythia.readString("36:onMode = 0"); // A0(H_3) decays only to tau-pairs
93     pythia.readString("36:onIfAny = 15");
94     }
95     else
96     {
97     cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
98     exit(1);
99     }
100     }
101    
102     // Initialize pythia
103    
104     pythia.init( 2212, 2212, 14000.);
105     pythia.particleData.listAll();
106     Rndm * rndmPtr = &pythia.rndm;
107     rndmPtr->init(runNumber);
108    
109     // Pileup pythia
110    
111     Pythia pileupPythia;
112     pileupPythia.readString("SoftQCD:minBias = on");
113     pileupPythia.init( 2212, 2212, 14000.);
114     pileupPythia.rndm.init(runNumber);
115    
116     // Ultra Fast Simulator
117    
118     UltraFastSim ufs(rndmPtr);
119 dasu 1.5 Analysis analysis(argv[1], &ufs);
120 dasu 1.1
121     // Begin event loop
122     for (int iEvent = 0; iEvent < nEvents; ) {
123    
124     // Generate event. Skip if error. List first one.
125     if (!pythia.next()) continue;
126    
127 dasu 1.4 if(iEvent < 10) pythia.event.list();
128    
129 dasu 1.1 // Add pileup
130    
131     int pileupEventCount = 0;
132     if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
133     for (int puEvent = 0; puEvent < pileupEventCount; ) {
134     if(!pileupPythia.next()) continue;
135     for(int i = 0; i < pileupPythia.event.size(); i++) {
136     Particle& particle = pileupPythia.event[i];
137     if(particle.status() > 0) {
138     if(particle.isVisible()) {
139     if(particle.pT() > 0.5) {
140     pythia.event.append(particle);
141     }
142     }
143     }
144     }
145     puEvent++;
146     }
147    
148     // Ultra fast simulation
149    
150     if(!ufs.run(pythia.event))
151     {
152     cerr << "Ultra fast simulation failed - aborting" << endl;
153     exit(1);
154     }
155    
156 dasu 1.5 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 dasu 1.2 cout << "Number of Gen Muons = " << ufs.muonList().size() << endl;
161     cout << "Number of Gen Taus = " << ufs.genTauList().size() << endl;
162 dasu 1.5 cout << "Number of taus = " << ufs.tauList().size() << endl;
163     cout << "Number of c Quarks = " << ufs.cQuarkList().size() << endl;
164 dasu 1.2 cout << "Number of b Quarks = " << ufs.bQuarkList().size() << endl;
165 dasu 1.5 cout << "Number of Jets = " << ufs.jetList().size() << endl;
166 dasu 1.2 cout << "Number of bJets = " << ufs.bJetList().size() << endl;
167     cout << endl;
168 dasu 1.5
169     // Run analysis
170    
171     if(!analysis.run())
172     {
173     cerr << "Analysis failed - aborting" << endl;
174     exit(2);
175     }
176    
177 dasu 1.1 // End of event loop. Statistics. Histogram. Done.
178    
179     if(!(iEvent % 100)) cout << "Processed event " << iEvent << endl;
180     iEvent++;
181    
182     }
183    
184     // Save output
185    
186     pythia.statistics();
187     pileupPythia.statistics();
188    
189     return 0;
190    
191     }