#include "PspGen/EvtGenKine.hh"
#include "EvtVector4R.hh"
#include "EvtRandom.hh"
#include "EvtRandomEngine.hh"
#include "EvtSimpleRandomEngine.hh"
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/WeightContainer.h"
#include "HepMC/Units.h"
#include "HepMC/IO_GenEvent.h"
#include <TROOT.h>
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include <TApplication.h>
#include <iostream>
#include <string>
#include <sstream>
using namespace std;
int main(int argc, char* argv[])
{
int set_limit = 200;
int eventsperset = 500000;
bool silent = true;
for(int set=0; set < set_limit; set++) {
std::stringstream out;
out << set+1001;
string setstring = out.str();
string filename = "phspEvents_500k_set" + setstring + ".out";
cout << "Filename: " << filename << endl;
cout << endl;
HepMC::IO_GenEvent hepMC_out(filename,std::ios::out);
EvtSimpleRandomEngine myRandom(4701+set);
EvtRandom::setRandomEngine(&myRandom);
// variables for the production e+ e- Psi2s
double mass_e=0.511e-3;
double mass_psi2s = 3.68609;
double pe_tot=0.5*sqrt(mass_psi2s*mass_psi2s-2.*mass_e);
// variables for the first decay Psi2s->Chi_c0+gamma
int firstNdaug = 2;
double firstMass[11];
EvtVector4R firstP4[8];
double firstMp = 3.68609;
firstMass[0] = 3.41475;
firstMass[10] = 0.;
// variables for the second decay Chi->c0->K+ K- pi0 pi0
int ndaug = 4;
double mass[30];
EvtVector4R p4[30];
mass[0] = .493677;
mass[1] = .493677;
mass[2] = .1349766;
mass[3] = .1349766;
// 5 events for visual control
if(!silent) {
for (int count = 0; count < 5; count++) {
std::cout << "Event #" << count << std::endl;
EvtGenKine::PhaseSpace(firstNdaug, firstMass, firstP4, firstMp);
for (int t = 0; t < firstNdaug; t++) {
std::cout << firstP4[t] << " m = " << sqrt(firstP4[t]*firstP4[t]) << std::endl;
}
// now use the mass of the Chi_c0 from the first decay as input for the second decay
EvtGenKine::PhaseSpace(ndaug, mass, p4, firstP4[0].mass());
for (int t = 0; t < ndaug; t++) {
std::cout << p4[t] << " m = " << sqrt(p4[t]*p4[t]) << std::endl;
p4[t].applyBoostTo(firstP4[0]);
std::cout << p4[t] << " boostet: m = " << sqrt(p4[t]*p4[t]) << std::endl;
}
}
}
/*
TApplication* rootapp = new TApplication("example",&argc, argv);
TH2F dalitz("Dalitz plot", "Dalitz Plot", 80,0.,5., 80,0.,6.);
TH1F massKpi("invMassKpi", "invMassKpi", 512, 0., 3.);
TH1F massKpipi("invMassKpipi", "invMassKpipi", 512, 0., 3.);
TH1F cosThetapi("cosThetapi", "cosThetapi", 512, -1.2, 1.2);
TH1F cosThetaKpi("cosThetaKpi", "cosThetaKpi", 512, -1.2, 1.2);
*/
for (int count = 0; count < eventsperset; count++) {
if(count % 100000 == 0) cout << "Event " << count << " generated." << endl;
EvtGenKine::PhaseSpace(firstNdaug, firstMass, firstP4, firstMp);
EvtGenKine::PhaseSpace(ndaug, mass, p4, firstP4[0].mass());
for (int t = 0; t < ndaug; t++) {
p4[t].applyBoostTo(firstP4[0]);
}
// dalitz.Fill((p4[0]+p4[1])*(p4[0]+p4[1]), (p4[1]+p4[2])*(p4[1]+p4[2]));
HepMC::GenEvent* evt = new HepMC::GenEvent( 20, count );
evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
// create production vertex
HepMC::GenVertex* vtx_prod = new HepMC::GenVertex();
evt->add_vertex( vtx_prod );
HepMC::GenParticle* eplus_part = new HepMC::GenParticle( HepMC::FourVector(0., 0., pe_tot, sqrt(mass_e*mass_e+pe_tot*pe_tot)), -11, 3 );
HepMC::GenParticle* eminus_part = new HepMC::GenParticle( HepMC::FourVector(0., 0., -pe_tot, sqrt(mass_e*mass_e+pe_tot*pe_tot)), 11, 3 );
vtx_prod->add_particle_in( eplus_part );
vtx_prod->add_particle_in( eminus_part );
HepMC::GenParticle* psi2s_part=new HepMC::GenParticle( HepMC::FourVector(0., 0., 0., mass_psi2s), 20443, 999);
vtx_prod->add_particle_out(psi2s_part);
// create Psi(2S) vertex
HepMC::GenVertex* vtx_psi2S = new HepMC::GenVertex();
evt->add_vertex( vtx_psi2S );
vtx_psi2S->add_particle_in( psi2s_part );
HepMC::GenParticle* chi_c0_part = new HepMC::GenParticle(HepMC::FourVector(firstP4[0].get(1), firstP4[0].get(2), firstP4[0].get(3), firstP4[0].get(0)), 10441, 999 );
vtx_psi2S->add_particle_out(chi_c0_part);
vtx_psi2S->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(firstP4[1].get(1), firstP4[1].get(2), firstP4[1].get(3), firstP4[1].get(0)), 22, 1 ) );
// now do the chi_c0 vertex
HepMC::GenVertex* vtx_chi_c0 = new HepMC::GenVertex();
evt->add_vertex( vtx_chi_c0 );
vtx_chi_c0->add_particle_in( chi_c0_part );
vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[0].get(1),p4[0].get(2),p4[0].get(3), p4[0].get(0)), 321, 1 ));
vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[1].get(1),p4[1].get(2),p4[1].get(3), p4[1].get(0)), -321, 1 ) );
vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[2].get(1),p4[2].get(2),p4[2].get(3), p4[2].get(0)), 111, 1 ) );
vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[3].get(1),p4[3].get(2),p4[3].get(3), p4[3].get(0)), 111, 1 ) );
// evt->print(std::cout);
hepMC_out << evt;
delete evt;
/*
massKpi.Fill( (p4[0]+p4[2]).mass(), 1.);
massKpi.Fill( (p4[0]+p4[3]).mass(), 1.);
massKpi.Fill( (p4[1]+p4[2]).mass(), 1.);
massKpi.Fill( (p4[1]+p4[3]).mass(), 1.);
massKpipi.Fill( (p4[0]+p4[2]+p4[3]).mass(), 1.);
massKpipi.Fill( (p4[1]+p4[2]+p4[3]).mass(), 1.);
cosThetapi.Fill( p4[2].get(3) /
sqrt(p4[2].get(1)*p4[2].get(1) + p4[2].get(2)*p4[2].get(2) + p4[2].get(3)*p4[2].get(3)),
1.);
cosThetaKpi.Fill( (p4[0]+p4[2]).get(3) /
sqrt((p4[0]+p4[2]).get(1)*(p4[0]+p4[2]).get(1) +
(p4[0]+p4[2]).get(2)*(p4[0]+p4[2]).get(2) +
(p4[0]+p4[2]).get(3)*(p4[0]+p4[2]).get(3)),
1.);
*/
}
cout << "Generated PHSP set " << set+1 << ", wrote to " << filename << endl;
}
cout << endl;
cout << "Finished!" << endl;
// massKpipi.Draw();
// rootapp->Run();
return 0;
}