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 "UFSDataStore.h"

int poisson(double mean, Rndm *rndmPtr) {
  static double oldMean = -1;
  static double g;
  if(mean != oldMean) {
    oldMean = mean;
    if(mean == 0) {
      g = 0;
    }
    else {
      g = exp(-mean);
    }
  }    
  double em = -1;
  double t = 1;
  do {
    em++;
    t *= rndmPtr->flat();
  } while(t > g);
  return em;
}

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);
      }
  }

  // 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;
  char jobName[256];
  sprintf(jobName, "%s-PU%3.3d-%4.4d", argv[1], meanPileupEventCount, runNumber);
  UFSDataStore dataStore(jobName, &ufs);

  // 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 = poisson(meanPileupEventCount, rndmPtr);
    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, rndmPtr))
      {
        cerr << "Ultra fast simulation failed - aborting" << endl;
        exit(1);
      }

    // Store data

    if(!dataStore.run())
      {
        cerr << "Failed to store data" << 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();

  return 0;

}
Revision:	1.16
Committed:	Wed Mar 2 21:57:09 2011 UTC (14 years, 2 months ago) by dasu
Content type:	text/plain
Branch:	MAIN
CVS Tags:	v2
Changes since 1.15:	+23 -2 lines
Log Message:	Added Poisson generator -- for some reason it is missing in Pythia
#	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.13	#include "UFSDataStore.h"
18	dasu	1.1
19	dasu	1.16	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	dasu	1.1	int main(int argc, char **argv) {
41			// Generator. Process selection. LHC initialization. Histogram.
42			Pythia pythia;
43			int nEvents = 0;
44			int runNumber = 0;
45			int meanPileupEventCount = 0;
46			if(argc < 3 \|\| argc > 5)
47			{
48	dasu	1.6	cerr << "Command syntax: " << argv[0] << " ZHmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
49			cerr << "or " << argv[0] << " ZHeebb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
50	dasu	1.10	cerr << "or " << argv[0] << " ZZmmbb[...]" << " <nEvents> [RandomSeed] [meanPileupEventCount]" << endl;
51	dasu	1.6	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	dasu	1.1	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	dasu	1.6	if(strncmp(argv[1], "Zee", 3) == 0)
70	dasu	1.1	{
71			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
72			pythia.readString("PhaseSpace:mHatMin = 60.");
73			pythia.readString("PhaseSpace:mHatMax = -1");
74	dasu	1.6	pythia.readString("23:onMode = 0"); // Z decays to electrons only
75			pythia.readString("23:onIfAny = 11");
76	dasu	1.1	}
77	dasu	1.6	else if(strncmp(argv[1], "Zmm", 3) == 0)
78	dasu	1.3	{
79			pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
80			pythia.readString("PhaseSpace:mHatMin = 60.");
81			pythia.readString("PhaseSpace:mHatMax = -1");
82			pythia.readString("23:onMode = 0"); // Z decays to muons only
83			pythia.readString("23:onIfAny = 13");
84			}
85	dasu	1.6	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	dasu	1.1	{
121			pythia.readString("WeakSingleBoson:ffbar2W = on");
122	dasu	1.6	pythia.readString("24:onMode = 0"); // W decays to muon
123			pythia.readString("24:onIfAny = 13");
124	dasu	1.1	}
125	dasu	1.6	else if(strncmp(argv[1], "Ten", 3) == 0)
126	dasu	1.1	{
127			pythia.readString("Top:all = on");
128	dasu	1.6	pythia.readString("24:onMode = 2"); // W+ decays to electron
129			pythia.readString("24:onIfAny = 11");
130	dasu	1.1	}
131	dasu	1.6	else if(strncmp(argv[1], "Tmn", 3) == 0)
132			{
133			pythia.readString("Top:all = on");
134			pythia.readString("24:onMode = 2"); // W+ decays to muon
135			pythia.readString("24:onIfAny = 13");
136			}
137			else if(strncmp(argv[1], "ZHmmbb", 6) == 0)
138	dasu	1.3	{
139			pythia.readString("HiggsSM:ffbar2HZ = on");
140			pythia.readString("25:m0 = 120");
141			pythia.readString("25:onMode = 0"); // Higgs decays to bbBar only
142			pythia.readString("25:onIfAny = 5");
143			pythia.readString("23:onMode = 0"); // Z decays to mu+,mu- only (for trigger)
144			pythia.readString("23:onIfAny = 13");
145			}
146	dasu	1.6	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	dasu	1.3	else if(strncmp(argv[1], "BB", 2) == 0)
156	dasu	1.1	{
157			pythia.readString("Higgs:useBSM = on");
158	dasu	1.3	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	dasu	1.6	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	dasu	1.1	}
179	dasu	1.9	else if(strncmp(argv[1], "ZZmmbb", 6) == 0)
180			{
181	dasu	1.10	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	dasu	1.9	}
185	dasu	1.1	else
186			{
187			cerr << "Unknown process type " << argv[1] << " - aborting" << endl;
188			exit(1);
189			}
190			}
191
192			// Initialize pythia
193
194			pythia.init( 2212, 2212, 14000.);
195			pythia.particleData.listAll();
196			Rndm * rndmPtr = &pythia.rndm;
197			rndmPtr->init(runNumber);
198
199			// Pileup pythia
200
201			Pythia pileupPythia;
202			pileupPythia.readString("SoftQCD:minBias = on");
203			pileupPythia.init( 2212, 2212, 14000.);
204			pileupPythia.rndm.init(runNumber);
205
206			// Ultra Fast Simulator
207
208	dasu	1.14	UltraFastSim ufs;
209	dasu	1.15	char jobName[256];
210	dasu	1.16	sprintf(jobName, "%s-PU%3.3d-%4.4d", argv[1], meanPileupEventCount, runNumber);
211	dasu	1.15	UFSDataStore dataStore(jobName, &ufs);
212	dasu	1.1
213			// Begin event loop
214			for (int iEvent = 0; iEvent < nEvents; ) {
215
216			// Generate event. Skip if error. List first one.
217			if (!pythia.next()) continue;
218
219	dasu	1.4	if(iEvent < 10) pythia.event.list();
220
221	dasu	1.1	// Add pileup
222
223			int pileupEventCount = 0;
224	dasu	1.16	if(meanPileupEventCount > 0) pileupEventCount = poisson(meanPileupEventCount, rndmPtr);
225	dasu	1.1	for (int puEvent = 0; puEvent < pileupEventCount; ) {
226			if(!pileupPythia.next()) continue;
227	dasu	1.11	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	dasu	1.1	for(int i = 0; i < pileupPythia.event.size(); i++) {
231			Particle& particle = pileupPythia.event[i];
232	dasu	1.11	particle.xProd(vx+particle.xProd());
233			particle.yProd(vy+particle.yProd());
234			particle.zProd(vz+particle.zProd());
235	dasu	1.1	if(particle.status() > 0) {
236			if(particle.isVisible()) {
237			if(particle.pT() > 0.5) {
238			pythia.event.append(particle);
239			}
240			}
241			}
242			}
243			puEvent++;
244			}
245
246			// Ultra fast simulation
247
248	dasu	1.14	if(!ufs.run(pythia.event, rndmPtr))
249	dasu	1.1	{
250			cerr << "Ultra fast simulation failed - aborting" << endl;
251			exit(1);
252			}
253
254	dasu	1.13	// Store data
255	dasu	1.5
256	dasu	1.13	if(!dataStore.run())
257	dasu	1.5	{
258	dasu	1.13	cerr << "Failed to store data" << endl;
259			exit(3);
260	dasu	1.5	}
261
262	dasu	1.1	// End of event loop. Statistics. Histogram. Done.
263
264			if(!(iEvent % 100)) cout << "Processed event " << iEvent << endl;
265			iEvent++;
266
267			}
268
269			// Save output
270
271			pythia.statistics();
272			pileupPythia.statistics();
273
274			return 0;
275
276			}