#include <string>
#include <vector>

using namespace std;

void run_ascii(
    const string& asciiFile = "/lustre/nyx/cbm/users/slebedev/hades/data/pluto.ascii.0.evt",
    int nEvents = 100000)
{
    double minMomentum = 0.0001; //[GeV/c]
    double maxMomentum = 0.9; //[GeV/c]
    //double minVertexZ = -70.;
    //double maxVertexZ = 0.;
    // IMPORTANT!!!
    // the values should correspond to values in  HRich700DigiParCreator::initZVertex()
    vector<double> zVec;
    zVec.push_back(0.);
    zVec.push_back(-25.);
    zVec.push_back(-50.);
    zVec.push_back(-75.);
    int nParticleEvent = 1;
    double ebeam = 1.58; //[GeV]
    double impactParameter = 0.; // [fm]
    double minPhi = -1.*TMath::Pi();
    double maxPhi = TMath::Pi();
    double minTheta = 0.;
    double maxTheta = 90.*TMath::DegToRad();
    int plutoFlag = 3;
    TRandom* random = new TRandom();
    random->SetSeed(0);

    ofstream file(asciiFile.c_str());
    int pdg = 2;
    for(Int_t iev = 0; iev < nEvents; iev++) {
	//double vz = random->Uniform(minVertexZ, maxVertexZ);
	double vz = zVec[iev % zVec.size()];
	double vx = random->Gaus(0., 0.1);
	double vy = random->Gaus(0., 0.1);

	file << std::scientific;
	file << (iev+1) << " " << nParticleEvent << " " << ebeam << " " << impactParameter << " " << plutoFlag << endl;

	// Select electron or positron, switch after all zVec
	if (iev % zVec.size() == 0) {
	    if (pdg == 2) pdg = 3;
	    else pdg = 2;
	}
	double p = random->Uniform(minMomentum, maxMomentum);
	double theta = random->Uniform(minTheta, maxTheta);
	double phi = random->Uniform(minPhi, maxPhi);
	double pt = p*TMath::Sin(theta);
	double px = pt*TMath::Cos(phi);
	double py = pt*TMath::Sin(phi);
	double pz = p*TMath::Cos(theta);

	file << p << " " << px << " " << py << " " << pz << " "<< vx << " " << vy << " " << vz <<" " << pdg << " " << 0 << " " << 0 << " " << 1.0 <<endl;
    }// events

    file.close();
    cout << "Done. Pluto ascii file:" << asciiFile << endl;
}