Skip to content
Snippets Groups Projects
Commit b868b0b1 authored by Jan Schulze's avatar Jan Schulze
Browse files

PspGen: added set functionality

parent 90eb82c2
No related branches found
No related tags found
No related merge requests found
......@@ -13,3 +13,5 @@ lib Psi2STo2K2PiGam :
exe ParserTestApp : ParserTestApp.cc Psi2STo2K2PiGam : ;
exe Mpsi2STo2K2PiGamTestApp : Mpsi2STo2K2PiGamTestApp.cc Psi2STo2K2PiGam : ;
exe qaMixedSampleApp : qaMixedSampleApp.cc Psi2STo2K2PiGam : ;
exe qaMixedSampleCacheApp : qaMixedSampleCacheApp.cc Psi2STo2K2PiGam : ;
......@@ -192,9 +192,10 @@ void Psi2STo2K2PiGamHist::initRootStuff(const std::string& fileName)
initHistMap();
_dataTuple = new TNtuple("dataTuple","dataTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipi:phipipi:costhetapiViapipi:phipiViapipi:weight");
_dataTuple = new TNtuple("dataTuple","dataTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipiViaK1pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype");
_mcTuple = new TNtuple("mcTuple","mcTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipiViaK1pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype");
_mcTuple = new TNtuple("mcTuple","mcTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipi:phipipi:costhetapiViapipi:phipiViapipi:weight");
}
......@@ -298,6 +299,8 @@ void Psi2STo2K2PiGamHist::fillHistos(const std::string& theKey, const Psi2STo2K2
void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEvtData* theData, double weight)
{
float datatype = 1.; // 1 = real data, 2 = pwamc
float Phichic0_HeliPsi2S = theData->chic0_HeliPsi2S_4V.Phi();
float CosThetachic0_HeliPsi2S = theData->chic0_HeliPsi2S_4V.CosTheta();
float PhiK1pi1pi2 = theData->KpPiPi_HeliChic0_4V.Phi();
......@@ -320,10 +323,10 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv
float mK2pi1pi2 = theData->KmPiPi_HeliChic0_4V.M();
float mK1pi1 = theData->KpPi0_HeliChic0_4V.M();
float mK1pi2 = theData->KpPi1_HeliChic0_4V.M();
float MK2pi1 = theData->KmPi0_HeliChic0_4V.M();
float MK2pi2 = theData->KmPi1_HeliChic0_4V.M();
float mK2pi1 = theData->KmPi0_HeliChic0_4V.M();
float mK2pi2 = theData->KmPi1_HeliChic0_4V.M();
float Mpipi = theData->PiPi_HeliChic0_4V.M();
float mpipi = theData->PiPi_HeliChic0_4V.M();
float CosThetaPiPiFromK1PiPi = theData->PiPi_HeliKpPi0Pi0_4V.CosTheta();
float PhiPiPiFromK1PiPi = theData->PiPi_HeliKpPi0Pi0_4V.Phi();
......@@ -356,23 +359,23 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv
fourVectors[19] = mK2pi1pi2;
fourVectors[20] = mK1pi1;
fourVectors[21] = mK1pi2;
fourVectors[22] = MK2pi1;
fourVectors[23] = MK2pi2;
fourVectors[22] = mK2pi1;
fourVectors[23] = mK2pi2;
fourVectors[24] = Mpipi;
fourVectors[24] = mpipi;
fourVectors[25] = CosThetaPiPiFromK1PiPi;
fourVectors[26] = PhiPiPiFromK1PiPi;
fourVectors[27] = CosThetaPiFromPiPi;
fourVectors[28] = PhiPiFromPiPi;
fourVectors[29] = evtweight;
fourVectors[30] = datatype;
// cout << evtweight << endl;
theTuple->Fill(fourVectors);
}
void Psi2STo2K2PiGamHist::plotCosPsi(const std::string& theKey, const Psi2STo2K2PiGamEvtData* theData, double weight){
std::string cosPsiString="cosPsi"+theKey;
TH1F* theHisto=_hist1DMap[cosPsiString];
......
This diff is collapsed.
This diff is collapsed.
......@@ -19,133 +19,159 @@
#include <iostream>
#include <string>
#include <sstream>
using namespace std;
int main(int argc, char* argv[])
{
HepMC::IO_GenEvent hepMC_out("phspEvents_1mio_RndSeed4720.out",std::ios::out);
EvtSimpleRandomEngine myRandom(4720);
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[8];
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
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;
int set_limit = 200;
int eventsperset = 500000;
bool silent = true;
for(int set=100; 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[8];
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;
}
}
}
// 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 < 1000000; 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]);
/*
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.);
*/
}
// 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();
// massKpipi.Draw();
// rootapp->Run();
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment