// main06.cc is a part of the PYTHIA event generator.
// Copyright (C) 2014 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This is a simple test program.
// It studies event properties of LEP1 events.

#include "Pythia8/Pythia.h"
#include "TFile.h"
#include "TTree.h"


using namespace Pythia8;
using namespace std;
int main(int argc, char *argv[])
{
  // need seed
  if(argc !=3)
    {
      cout<<"Wrong number of arguments"<<endl;
      return -1;
    }
  int seed=atoi(argv[1]);
  int number_sim=atoi(argv[2]);
  string seed_s=argv[1];
  TFile *f_out = new TFile(("out_"+seed_s + ".root").c_str(), "RECREATE");
  TTree *t= new TTree("t", "t");
  
  int particles_max=200;
  Double_t PX[particles_max];
  Double_t PY[particles_max];
  Double_t PZ[particles_max];
  Double_t PE[particles_max];
  Int_t    PDG[particles_max];

  t->Branch("PX", &PX, "PX[200]/D");
  t->Branch("PY", &PY, "PY[200]/D");
  t->Branch("PZ", &PZ, "PZ[200]/D");
  t->Branch("PZ", &PZ, "PZ[200]/D");
  t->Branch("PDG", &PDG, "PDG[200]/I"); 
  
  // Generator.
  Pythia pythia;
  // set deed
  pythia.readString("Random:seed =" + seed_s);
  
  // Allow no substructure in e+- beams: normal for corrected LEP data.
  pythia.readString("PDF:lepton = off");
  // Process selection.
  pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
  // Switch off all Z0 decays and then switch back on those to quarks.
  pythia.readString("23:onMode = off");
  pythia.readString("23:onIfAny = 1 2 3 4 5");
  
  // initialization at Z0 mass.
  pythia.readString("Beams:idA =  11");
  pythia.readString("Beams:idB = -11");
  double mZ = pythia.particleData.m0(23);
  pythia.settings.parm("Beams:eCM", mZ);
  pythia.init();


  // Begin event loop. Generate event. Skip if error. List first few.
  for (int iEvent = 0; iEvent < number_sim; ++iEvent) {
    if (!pythia.next()) continue;
    // Find and histogram charged multiplicity.
    if( pythia.event.size()>200) continue;
    
    for (int i = 0; i < pythia.event.size(); ++i)
      {
        if(pythia.event[i].isFinal() )
          {
            PX[i]= pythia.event[i].px();
            PY[i]=pythia.event[i].py(); 
            PZ[i]=pythia.event[i].pz(); 
            PE[i]=pythia.event[i].e();
            //cout<<pythia.event[i].id()<<endl;
            PDG[i]=pythia.event[i].id();
            
          }
        
      }//particle loop in event
    t->Fill();
  }//event loop
  t->Write();
  f_out->Close();

  // Done.
  return 0;
}