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 "UltraFastSim.h"
#include "Analysis.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] << " 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
      {
        cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
        exit(1);
      }
  }

  // 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);
  Analysis analysis(argv[1], &pythia, &ufs, true);

  // 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;
      for(int i = 0; i < pileupPythia.event.size(); i++) {
        Particle& particle = pileupPythia.event[i];
        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 analysis

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

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

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

  // Save output

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

  return 0;

}
Revision:	1.7
Committed:	Wed Feb 16 00:33:49 2011 UTC (14 years, 3 months ago) by dasu
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+1 -14 lines
Log Message:	Made Luk's code compile and run in Analysis.*
#	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.6	cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
28			cerr << "or " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
29			cerr << "or " << argv[0] << " BBHmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
30			cerr << "or " << argv[0] << " BBAmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
31			cerr << "or " << argv[0] << " Zee[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
32			cerr << "or " << argv[0] << " Zmm[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
33			cerr << "or " << argv[0] << " Zeh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
34			cerr << "or " << argv[0] << " Zmh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
35			cerr << "or " << argv[0] << " Zhh[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
36			cerr << "or " << argv[0] << " Wen[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
37			cerr << "or " << argv[0] << " Wmn[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
38			cerr << "or " << argv[0] << " Ten[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
39			cerr << "or " << argv[0] << " Tmn[...] " << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
40			cerr << "In above e-electron, m-muon, h-hadronic tau. n-neutrino" << endl;
41	dasu	1.1	exit(1);
42			}
43			else {
44			nEvents = atoi(argv[2]);
45			if(argc >= 4) runNumber = atoi(argv[3]);
46			if(argc >= 5) meanPileupEventCount = atoi(argv[4]);
47	dasu	1.6	if(strncmp(argv[1], "Zee", 3) == 0)
48	dasu	1.1	{
49			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
50			pythia.readString("PhaseSpace:mHatMin = 60.");
51			pythia.readString("PhaseSpace:mHatMax = -1");
52	dasu	1.6	pythia.readString("23:onMode = 0"); // Z decays to electrons only
53			pythia.readString("23:onIfAny = 11");
54	dasu	1.1	}
55	dasu	1.6	else if(strncmp(argv[1], "Zmm", 3) == 0)
56	dasu	1.3	{
57			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
58			pythia.readString("PhaseSpace:mHatMin = 60.");
59			pythia.readString("PhaseSpace:mHatMax = -1");
60			pythia.readString("23:onMode = 0"); // Z decays to muons only
61			pythia.readString("23:onIfAny = 13");
62			}
63	dasu	1.6	else if(strncmp(argv[1], "Zeh", 3) == 0)
64			{
65			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
66			pythia.readString("PhaseSpace:mHatMin = 60.");
67			pythia.readString("PhaseSpace:mHatMax = -1");
68			pythia.readString("23:onMode = 0"); // Z decays to taus only
69			pythia.readString("23:onIfAny = 15");
70			pythia.readString("15:onMode = 2"); // tau- decays to electrons - second tau decays in all modes
71			pythia.readString("15:onIfAny = 11");
72			}
73			else if(strncmp(argv[1], "Zmh", 3) == 0)
74			{
75			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
76			pythia.readString("PhaseSpace:mHatMin = 60.");
77			pythia.readString("PhaseSpace:mHatMax = -1");
78			pythia.readString("23:onMode = 0"); // Z decays to taus only
79			pythia.readString("23:onIfAny = 15");
80			pythia.readString("15:onMode = 2"); // tau- decays to muons - second tau decays in all modes
81			pythia.readString("15:onIfAny = 13");
82			}
83			else if(strncmp(argv[1], "Zhh", 3) == 0)
84			{
85			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
86			pythia.readString("PhaseSpace:mHatMin = 60.");
87			pythia.readString("PhaseSpace:mHatMax = -1");
88			pythia.readString("23:onMode = 0"); // Z decays to taus only
89			pythia.readString("23:onIfAny = 15");
90			}
91			else if(strncmp(argv[1], "Wen", 3) == 0)
92			{
93			pythia.readString("WeakSingleBoson:ffbar2W = on");
94			pythia.readString("24:onMode = 0"); // W decays to electron
95			pythia.readString("24:onIfAny = 11");
96			}
97			else if(strncmp(argv[1], "Wmn", 3) == 0)
98	dasu	1.1	{
99			pythia.readString("WeakSingleBoson:ffbar2W = on");
100	dasu	1.6	pythia.readString("24:onMode = 0"); // W decays to muon
101			pythia.readString("24:onIfAny = 13");
102	dasu	1.1	}
103	dasu	1.6	else if(strncmp(argv[1], "Ten", 3) == 0)
104	dasu	1.1	{
105			pythia.readString("Top:all = on");
106	dasu	1.6	pythia.readString("24:onMode = 2"); // W+ decays to electron
107			pythia.readString("24:onIfAny = 11");
108	dasu	1.1	}
109	dasu	1.6	else if(strncmp(argv[1], "Tmn", 3) == 0)
110			{
111			pythia.readString("Top:all = on");
112			pythia.readString("24:onMode = 2"); // W+ decays to muon
113			pythia.readString("24:onIfAny = 13");
114			}
115			else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
116	dasu	1.3	{
117			pythia.readString("HiggsSM:ffbar2HZ = on");
118			pythia.readString("25:m0 = 120");
119			pythia.readString("25:onMode = 0"); // Higgs decays to bbBar only
120			pythia.readString("25:onIfAny = 5");
121			pythia.readString("23:onMode = 0"); // Z decays to mu+,mu- only (for trigger)
122			pythia.readString("23:onIfAny = 13");
123			}
124	dasu	1.6	else if(strncmp(argv[1], "ZHeebb", 6) == 0)
125			{
126			pythia.readString("HiggsSM:ffbar2HZ = on");
127			pythia.readString("25:m0 = 120");
128			pythia.readString("25:onMode = 0"); // Higgs decays to bbBar only
129			pythia.readString("25:onIfAny = 5");
130			pythia.readString("23:onMode = 0"); // Z decays to e+,e- only (for trigger)
131			pythia.readString("23:onIfAny = 11");
132			}
133	dasu	1.3	else if(strncmp(argv[1], "BB", 2) == 0)
134	dasu	1.1	{
135			pythia.readString("Higgs:useBSM = on");
136	dasu	1.3	pythia.readString("25:m0 = 129");
137			pythia.readString("35:m0 = 499.7");
138			pythia.readString("36:m0 = 500");
139			pythia.readString("37:m0 = 506");
140			pythia.readString("HiggsHchg:tanBeta = 30");
141	dasu	1.6	if(strncmp(argv[1], "BBH", 3) == 0) {
142			pythia.readString("HiggsBSM:gg2H2bbbar = on");
143			pythia.readString("35:onMode = 0"); // H2(H_2) decays only to tau-pairs
144			pythia.readString("35:onIfAny = 15");
145			}
146			else {
147			pythia.readString("HiggsBSM:gg2A3bbbar = on");
148			pythia.readString("36:onMode = 0"); // A0(H_3) decays only to tau-pairs
149			pythia.readString("36:onIfAny = 15");
150			}
151			if(strncmp(argv[1], "BBHmh", 5) == 0 \|\|
152			strncmp(argv[1], "BBAmh", 5) == 0) {
153			pythia.readString("15:onMode = 2"); // tau- decays to muons - second tau decays in all modes
154			pythia.readString("15:onIfAny = 13");
155			}
156	dasu	1.1	}
157			else
158			{
159			cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
160			exit(1);
161			}
162			}
163
164			// Initialize pythia
165
166			pythia.init( 2212, 2212, 14000.);
167			pythia.particleData.listAll();
168			Rndm * rndmPtr = &pythia.rndm;
169			rndmPtr->init(runNumber);
170
171			// Pileup pythia
172
173			Pythia pileupPythia;
174			pileupPythia.readString("SoftQCD:minBias = on");
175			pileupPythia.init( 2212, 2212, 14000.);
176			pileupPythia.rndm.init(runNumber);
177
178			// Ultra Fast Simulator
179
180			UltraFastSim ufs(rndmPtr);
181	dasu	1.7	Analysis analysis(argv[1], &pythia, &ufs, true);
182	dasu	1.1
183			// Begin event loop
184			for (int iEvent = 0; iEvent < nEvents; ) {
185
186			// Generate event. Skip if error. List first one.
187			if (!pythia.next()) continue;
188
189	dasu	1.4	if(iEvent < 10) pythia.event.list();
190
191	dasu	1.1	// Add pileup
192
193			int pileupEventCount = 0;
194			if(meanPileupEventCount > 0) pileupEventCount = meanPileupEventCount; // * rmdmPtr->pois(); (not available yet)
195			for (int puEvent = 0; puEvent < pileupEventCount; ) {
196			if(!pileupPythia.next()) continue;
197			for(int i = 0; i < pileupPythia.event.size(); i++) {
198			Particle& particle = pileupPythia.event[i];
199			if(particle.status() > 0) {
200			if(particle.isVisible()) {
201			if(particle.pT() > 0.5) {
202			pythia.event.append(particle);
203			}
204			}
205			}
206			}
207			puEvent++;
208			}
209
210			// Ultra fast simulation
211
212			if(!ufs.run(pythia.event))
213			{
214			cerr << "Ultra fast simulation failed - aborting" << endl;
215			exit(1);
216			}
217
218	dasu	1.5	// Run analysis
219
220			if(!analysis.run())
221			{
222			cerr << "Analysis failed - aborting" << endl;
223			exit(2);
224			}
225
226	dasu	1.1	// End of event loop. Statistics. Histogram. Done.
227
228			if(!(iEvent % 100)) cout << "Processed event " << iEvent << endl;
229			iEvent++;
230
231			}
232
233			// Save output
234
235			pythia.statistics();
236			pileupPythia.statistics();
237
238			return 0;
239
240			}