dasu/UltraFastSim/generateData.cc

// This simple program generates various signal processes for 
// MSSM bbPhi process and its backgrounds.
// It includes the mixing of required level of pileup.
// It then takes pythia output and runs it through very simplistic
// detector simulation using parameterizations, including making
// isolated light leptons, photons, composite tau, jet and b-tagged 
// jet objects.
// It saves the output in root format for later analysis.

using namespace std;
#include <math.h>

#include "Pythia.h"
using namespace Pythia8; 

#include "TROOT.h"
#include "TFile.h"

#include "UltraFastSim.h"
#include "bbHAnalysis.h"
#include "ZHAnalysis.h"

int main(int argc, char **argv) {
  // Generator. Process selection. LHC initialization. Histogram.
  Pythia pythia;
  int nEvents = 0;
  int runNumber = 0;
  int meanPileupEventCount = 0;
  if(argc < 3 || argc > 5)
    {
      cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Zee[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Zmm[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Zeh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Zmh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Zhh[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Wen[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Wmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Ten[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "or              " << argv[0] << " Tmn[...]   " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
      cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
      exit(1);
    }
  else {
    nEvents = atoi(argv[2]);
    if(argc >= 4) runNumber = atoi(argv[3]);
    if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
    if(strncmp(argv[1], "Zee", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("PhaseSpace:mHatMin = 60.");    
        pythia.readString("PhaseSpace:mHatMax = -1");    
        pythia.readString("23:onMode = 0");       // Z decays to electrons only
        pythia.readString("23:onIfAny = 11");
      }
    else if(strncmp(argv[1], "Zmm", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("PhaseSpace:mHatMin = 60.");    
        pythia.readString("PhaseSpace:mHatMax = -1");    
        pythia.readString("23:onMode = 0");       // Z decays to muons only
        pythia.readString("23:onIfAny = 13");
      }
    else if(strncmp(argv[1], "Zeh", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("PhaseSpace:mHatMin = 60.");    
        pythia.readString("PhaseSpace:mHatMax = -1");    
        pythia.readString("23:onMode = 0");       // Z decays to taus only
        pythia.readString("23:onIfAny = 15");
        pythia.readString("15:onMode = 2");       // tau- decays to electrons - second tau decays in all modes
        pythia.readString("15:onIfAny = 11");
      }
    else if(strncmp(argv[1], "Zmh", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("PhaseSpace:mHatMin = 60.");    
        pythia.readString("PhaseSpace:mHatMax = -1");    
        pythia.readString("23:onMode = 0");       // Z decays to taus only
        pythia.readString("23:onIfAny = 15");
        pythia.readString("15:onMode = 2");       // tau- decays to muons - second tau decays in all modes
        pythia.readString("15:onIfAny = 13");
      }
    else if(strncmp(argv[1], "Zhh", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("PhaseSpace:mHatMin = 60.");    
        pythia.readString("PhaseSpace:mHatMax = -1");    
        pythia.readString("23:onMode = 0");       // Z decays to taus only
        pythia.readString("23:onIfAny = 15");
      }
    else if(strncmp(argv[1], "Wen", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2W = on");
        pythia.readString("24:onMode = 0");       // W decays to electron
        pythia.readString("24:onIfAny = 11");
      }
    else if(strncmp(argv[1], "Wmn", 3) == 0)
      {
        pythia.readString("WeakSingleBoson:ffbar2W = on");
        pythia.readString("24:onMode = 0");       // W decays to muon
        pythia.readString("24:onIfAny = 13");
      }
    else if(strncmp(argv[1], "Ten", 3) == 0)
      {
        pythia.readString("Top:all = on");
        pythia.readString("24:onMode = 2");       // W+ decays to electron
        pythia.readString("24:onIfAny = 11");
      }
    else if(strncmp(argv[1], "Tmn", 3) == 0)
      {
        pythia.readString("Top:all = on");
        pythia.readString("24:onMode = 2");       // W+ decays to muon
        pythia.readString("24:onIfAny = 13");
      }
    else if(strncmp(argv[1], "ZHmmbb", 6) == 0) 
      {
        pythia.readString("HiggsSM:ffbar2HZ = on");
        pythia.readString("25:m0 = 120");
        pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
        pythia.readString("25:onIfAny = 5");
        pythia.readString("23:onMode = 0");      // Z decays to mu+,mu- only (for trigger)
        pythia.readString("23:onIfAny = 13");
      }
    else if(strncmp(argv[1], "ZHeebb", 6) == 0) 
      {
        pythia.readString("HiggsSM:ffbar2HZ = on");
        pythia.readString("25:m0 = 120");
        pythia.readString("25:onMode = 0");      // Higgs decays to bbBar only
        pythia.readString("25:onIfAny = 5");
        pythia.readString("23:onMode = 0");      // Z decays to e+,e- only (for trigger)
        pythia.readString("23:onIfAny = 11");
      }
    else if(strncmp(argv[1], "BB", 2) == 0)
      {
        pythia.readString("Higgs:useBSM = on");
        pythia.readString("25:m0 = 129");
        pythia.readString("35:m0 = 499.7");
        pythia.readString("36:m0 = 500");
        pythia.readString("37:m0 = 506");
        pythia.readString("HiggsHchg:tanBeta = 30");
        if(strncmp(argv[1], "BBH", 3) == 0) {
          pythia.readString("HiggsBSM:gg2H2bbbar = on");
          pythia.readString("35:onMode = 0");      // H2(H_2) decays only to tau-pairs
          pythia.readString("35:onIfAny = 15");
        }
        else {
          pythia.readString("HiggsBSM:gg2A3bbbar = on");
          pythia.readString("36:onMode = 0");      // A0(H_3) decays only to tau-pairs
          pythia.readString("36:onIfAny = 15");
        }
        if(strncmp(argv[1], "BBHmh", 5) == 0 ||
           strncmp(argv[1], "BBAmh", 5) == 0) {
          pythia.readString("15:onMode = 2");       // tau- decays to muons - second tau decays in all modes
          pythia.readString("15:onIfAny = 13");
        }
      }
     else if(strncmp(argv[1], "ZZmmbb", 6) == 0)
       {
         pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
         pythia.readString("23:onMode = 0");
         pythia.readString("23:onIfAny = 5 13");   // Let the Z decay to bQuarks or muons 
       }
    else
      {
        cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
        exit(1);
      }
  }

  string name(argv[1]);
  name += ".root";
  TFile *outFile = TFile::Open(name.c_str(), "recreate");

  // Initialize pythia

  pythia.init( 2212, 2212, 14000.);
  pythia.particleData.listAll();
  Rndm * rndmPtr = &pythia.rndm;
  rndmPtr->init(runNumber);

  // Pileup pythia

  Pythia pileupPythia;
  pileupPythia.readString("SoftQCD:minBias = on");
  pileupPythia.init( 2212, 2212, 14000.);
  pileupPythia.rndm.init(runNumber);

  // Ultra Fast Simulator

  UltraFastSim ufs(rndmPtr);
  bbHAnalysis bbH(outFile, &pythia, &ufs, true);
  ZHAnalysis ZH; // outFile, &pythia, &ufs, true);
  ZH.BookTree();

  // Begin event loop
  for (int iEvent = 0; iEvent < nEvents; ) {

    // Generate event. Skip if error. List first one.
    if (!pythia.next()) continue;

    if(iEvent < 10) pythia.event.list();

    // Add pileup

    int pileupEventCount = 0;
    if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
    for (int puEvent = 0; puEvent < pileupEventCount; ) {
      if(!pileupPythia.next()) continue;
      double vx = rndmPtr->gauss()*2.; // Mean beam spread is 2 mm in x-y and 75 mm in z 
      double vy = rndmPtr->gauss()*2.;
      double vz = rndmPtr->gauss()*75.;
      for(int i = 0; i < pileupPythia.event.size(); i++) {
        Particle& particle = pileupPythia.event[i];
        particle.xProd(vx+particle.xProd());
        particle.yProd(vy+particle.yProd());
        particle.zProd(vz+particle.zProd());
        if(particle.status() > 0) {
          if(particle.isVisible()) {
            if(particle.pT() > 0.5) {
              pythia.event.append(particle);
            }
          }
        }
      }
      puEvent++;
    }

    // Ultra fast simulation

    if(!ufs.run(pythia.event))
      {
        cerr << "Ultra fast simulation failed - aborting" << endl;
        exit(1);
      }

    // Run bbHAnalysis

    if(!bbH.run())
      {
        cerr << "bbH Analysis failed - aborting" << endl;
        exit(2);
      }

    // Run ZHAnalysis

    if(!ZH.Run(&ufs))
      {
        cerr << "ZH Analysis failed - aborting" << endl;
        exit(3);
      }

    // End of event loop. Statistics. Histogram. Done.

    if(!(iEvent % 100)) cout << "Processed event " << iEvent << endl;
    iEvent++;
    
  }

  // Save output

  pythia.statistics();
  pileupPythia.statistics();

  bbH.end();

  outFile->cd();
  outFile->Write();
  outFile->Close();

  return 0;

}
Revision:	1.11
Committed:	Thu Feb 24 07:34:29 2011 UTC (14 years, 2 months ago) by dasu
Content type:	text/plain
Branch:	MAIN
Changes since 1.10:	+9 -3 lines
Log Message:	Added Majid's changes
#	Content
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 "TROOT.h"
17	#include "TFile.h"
18
19	#include "UltraFastSim.h"
20	#include "bbHAnalysis.h"
21	#include "ZHAnalysis.h"
22
23	int main(int argc, char **argv) {
24	// Generator. Process selection. LHC initialization. Histogram.
25	Pythia pythia;
26	int nEvents = 0;
27	int runNumber = 0;
28	int meanPileupEventCount = 0;
29	if(argc < 3 \|\| argc > 5)
30	{
31	cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
32	cerr << "or " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
33	cerr << "or " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
34	cerr << "or " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
35	cerr << "or " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
36	cerr << "or " << argv[0] << " Zee[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
37	cerr << "or " << argv[0] << " Zmm[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
38	cerr << "or " << argv[0] << " Zeh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
39	cerr << "or " << argv[0] << " Zmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
40	cerr << "or " << argv[0] << " Zhh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
41	cerr << "or " << argv[0] << " Wen[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
42	cerr << "or " << argv[0] << " Wmn[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
43	cerr << "or " << argv[0] << " Ten[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
44	cerr << "or " << argv[0] << " Tmn[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
45	cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
46	exit(1);
47	}
48	else {
49	nEvents = atoi(argv[2]);
50	if(argc >= 4) runNumber = atoi(argv[3]);
51	if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
52	if(strncmp(argv[1], "Zee", 3) == 0)
53	{
54	pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
55	pythia.readString("PhaseSpace:mHatMin = 60.");
56	pythia.readString("PhaseSpace:mHatMax = -1");
57	pythia.readString("23:onMode = 0"); // Z decays to electrons only
58	pythia.readString("23:onIfAny = 11");
59	}
60	else if(strncmp(argv[1], "Zmm", 3) == 0)
61	{
62	pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
63	pythia.readString("PhaseSpace:mHatMin = 60.");
64	pythia.readString("PhaseSpace:mHatMax = -1");
65	pythia.readString("23:onMode = 0"); // Z decays to muons only
66	pythia.readString("23:onIfAny = 13");
67	}
68	else if(strncmp(argv[1], "Zeh", 3) == 0)
69	{
70	pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
71	pythia.readString("PhaseSpace:mHatMin = 60.");
72	pythia.readString("PhaseSpace:mHatMax = -1");
73	pythia.readString("23:onMode = 0"); // Z decays to taus only
74	pythia.readString("23:onIfAny = 15");
75	pythia.readString("15:onMode = 2"); // tau- decays to electrons - second tau decays in all modes
76	pythia.readString("15:onIfAny = 11");
77	}
78	else if(strncmp(argv[1], "Zmh", 3) == 0)
79	{
80	pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
81	pythia.readString("PhaseSpace:mHatMin = 60.");
82	pythia.readString("PhaseSpace:mHatMax = -1");
83	pythia.readString("23:onMode = 0"); // Z decays to taus only
84	pythia.readString("23:onIfAny = 15");
85	pythia.readString("15:onMode = 2"); // tau- decays to muons - second tau decays in all modes
86	pythia.readString("15:onIfAny = 13");
87	}
88	else if(strncmp(argv[1], "Zhh", 3) == 0)
89	{
90	pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
91	pythia.readString("PhaseSpace:mHatMin = 60.");
92	pythia.readString("PhaseSpace:mHatMax = -1");
93	pythia.readString("23:onMode = 0"); // Z decays to taus only
94	pythia.readString("23:onIfAny = 15");
95	}
96	else if(strncmp(argv[1], "Wen", 3) == 0)
97	{
98	pythia.readString("WeakSingleBoson:ffbar2W = on");
99	pythia.readString("24:onMode = 0"); // W decays to electron
100	pythia.readString("24:onIfAny = 11");
101	}
102	else if(strncmp(argv[1], "Wmn", 3) == 0)
103	{
104	pythia.readString("WeakSingleBoson:ffbar2W = on");
105	pythia.readString("24:onMode = 0"); // W decays to muon
106	pythia.readString("24:onIfAny = 13");
107	}
108	else if(strncmp(argv[1], "Ten", 3) == 0)
109	{
110	pythia.readString("Top:all = on");
111	pythia.readString("24:onMode = 2"); // W+ decays to electron
112	pythia.readString("24:onIfAny = 11");
113	}
114	else if(strncmp(argv[1], "Tmn", 3) == 0)
115	{
116	pythia.readString("Top:all = on");
117	pythia.readString("24:onMode = 2"); // W+ decays to muon
118	pythia.readString("24:onIfAny = 13");
119	}
120	else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
121	{
122	pythia.readString("HiggsSM:ffbar2HZ = on");
123	pythia.readString("25:m0 = 120");
124	pythia.readString("25:onMode = 0"); // Higgs decays to bbBar only
125	pythia.readString("25:onIfAny = 5");
126	pythia.readString("23:onMode = 0"); // Z decays to mu+,mu- only (for trigger)
127	pythia.readString("23:onIfAny = 13");
128	}
129	else if(strncmp(argv[1], "ZHeebb", 6) == 0)
130	{
131	pythia.readString("HiggsSM:ffbar2HZ = on");
132	pythia.readString("25:m0 = 120");
133	pythia.readString("25:onMode = 0"); // Higgs decays to bbBar only
134	pythia.readString("25:onIfAny = 5");
135	pythia.readString("23:onMode = 0"); // Z decays to e+,e- only (for trigger)
136	pythia.readString("23:onIfAny = 11");
137	}
138	else if(strncmp(argv[1], "BB", 2) == 0)
139	{
140	pythia.readString("Higgs:useBSM = on");
141	pythia.readString("25:m0 = 129");
142	pythia.readString("35:m0 = 499.7");
143	pythia.readString("36:m0 = 500");
144	pythia.readString("37:m0 = 506");
145	pythia.readString("HiggsHchg:tanBeta = 30");
146	if(strncmp(argv[1], "BBH", 3) == 0) {
147	pythia.readString("HiggsBSM:gg2H2bbbar = on");
148	pythia.readString("35:onMode = 0"); // H2(H_2) decays only to tau-pairs
149	pythia.readString("35:onIfAny = 15");
150	}
151	else {
152	pythia.readString("HiggsBSM:gg2A3bbbar = on");
153	pythia.readString("36:onMode = 0"); // A0(H_3) decays only to tau-pairs
154	pythia.readString("36:onIfAny = 15");
155	}
156	if(strncmp(argv[1], "BBHmh", 5) == 0 \|\|
157	strncmp(argv[1], "BBAmh", 5) == 0) {
158	pythia.readString("15:onMode = 2"); // tau- decays to muons - second tau decays in all modes
159	pythia.readString("15:onIfAny = 13");
160	}
161	}
162	else if(strncmp(argv[1], "ZZmmbb", 6) == 0)
163	{
164	pythia.readString("WeakDoubleBoson:ffbar2gmZgmZ = on");
165	pythia.readString("23:onMode = 0");
166	pythia.readString("23:onIfAny = 5 13"); // Let the Z decay to bQuarks or muons
167	}
168	else
169	{
170	cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
171	exit(1);
172	}
173	}
174
175	string name(argv[1]);
176	name += ".root";
177	TFile *outFile = TFile::Open(name.c_str(), "recreate");
178
179	// Initialize pythia
180
181	pythia.init( 2212, 2212, 14000.);
182	pythia.particleData.listAll();
183	Rndm * rndmPtr = &pythia.rndm;
184	rndmPtr->init(runNumber);
185
186	// Pileup pythia
187
188	Pythia pileupPythia;
189	pileupPythia.readString("SoftQCD:minBias = on");
190	pileupPythia.init( 2212, 2212, 14000.);
191	pileupPythia.rndm.init(runNumber);
192
193	// Ultra Fast Simulator
194
195	UltraFastSim ufs(rndmPtr);
196	bbHAnalysis bbH(outFile, &pythia, &ufs, true);
197	ZHAnalysis ZH; // outFile, &pythia, &ufs, true);
198	ZH.BookTree();
199
200	// Begin event loop
201	for (int iEvent = 0; iEvent < nEvents; ) {
202
203	// Generate event. Skip if error. List first one.
204	if (!pythia.next()) continue;
205
206	if(iEvent < 10) pythia.event.list();
207
208	// Add pileup
209
210	int pileupEventCount = 0;
211	if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
212	for (int puEvent = 0; puEvent < pileupEventCount; ) {
213	if(!pileupPythia.next()) continue;
214	double vx = rndmPtr->gauss()*2.; // Mean beam spread is 2 mm in x-y and 75 mm in z
215	double vy = rndmPtr->gauss()*2.;
216	double vz = rndmPtr->gauss()*75.;
217	for(int i = 0; i < pileupPythia.event.size(); i++) {
218	Particle& particle = pileupPythia.event[i];
219	particle.xProd(vx+particle.xProd());
220	particle.yProd(vy+particle.yProd());
221	particle.zProd(vz+particle.zProd());
222	if(particle.status() > 0) {
223	if(particle.isVisible()) {
224	if(particle.pT() > 0.5) {
225	pythia.event.append(particle);
226	}
227	}
228	}
229	}
230	puEvent++;
231	}
232
233	// Ultra fast simulation
234
235	if(!ufs.run(pythia.event))
236	{
237	cerr << "Ultra fast simulation failed - aborting" << endl;
238	exit(1);
239	}
240
241	// Run bbHAnalysis
242
243	if(!bbH.run())
244	{
245	cerr << "bbH Analysis failed - aborting" << endl;
246	exit(2);
247	}
248
249	// Run ZHAnalysis
250
251	if(!ZH.Run(&ufs))
252	{
253	cerr << "ZH Analysis failed - aborting" << endl;
254	exit(3);
255	}
256
257	// End of event loop. Statistics. Histogram. Done.
258
259	if(!(iEvent % 100)) cout << "Processed event " << iEvent << endl;
260	iEvent++;
261
262	}
263
264	// Save output
265
266	pythia.statistics();
267	pileupPythia.statistics();
268
269	bbH.end();
270
271	outFile->cd();
272	outFile->Write();
273	outFile->Close();
274
275	return 0;
276
277	}