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"
#include "bbHAnalysis.cc"
//#include "ZHAnalysis.cc"

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