// Copyright (C) 2009 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

#include "Pythia.h"

#include <CBook/PBook.h>
#include <CBook/H1D.h>
#include <CBook/H2D.h>


using namespace Pythia8;
int main() {
  // Generator. Process selection. Tevatron initialization. Histogram.
  Pythia pythia;

  int  nEvent = 50000;

//  pythia.readString("WeakSingleBoson:ffbar2gmZ = on");    
//  pythia.readString("HardQCD:all = on");
  pythia.readString("PromptPhoton:all = on");
  pythia.readString("PhaseSpace:mHatMin = 20.");
  pythia.readString("PhaseSpace:mHatMax = 1000.");
  pythia.init( 2212, 2212, 7000.);


  // open histogram file and book histograms
   PBook*  pbook = PBook::open("analysis.pbu");
   H1D* h1 = new H1D("pT(gamma)",140,0,700);
   H1D* h2 = new H1D("eta(gamma)",30,-3,3);



  // Begin event loop. Generate event. Skip if error. List first one.
  for (int n=0; n<nEvent; ++n) {
    if (!pythia.next()) continue;
    if (n < 1) {pythia.info.list(); pythia.event.list();}

       if (n<=10 &&
        n<=100 && (n%10) == 0 ||
        n<=1000 && (n%100) == 0  ||
        n>=1000 && (n%1000) == 0. ) {
        cout << "No events= " << n << " passed"  << endl; };
    // Loop over particles in event. Find last Z0 copy. Fill its pT. 
    int iZ = 0;
    for (int i=0; i < pythia.event.size(); i++) {
      if (!pythia.event[i].isFinal())         continue;
      if (fabs(pythia.event[i].eta()) > 2.37) continue;
      if (pythia.event[i].pT()<25)            continue;
      if (pythia.event[i].idAbs() != 22 )     continue; // fill photons
          h1->fill(pythia.event[i].pT());
          h2->fill(pythia.event[i].y());

    } // end loop oever particles 


  // End of event loop. Statistics. Histogram. Done.
  }
  pythia.statistics();
    // Output histograms
  double sigmapb = pythia.info.sigmaGen() * 1.0E9;
  std::cout << " Total cross section (pb) = " << sigmapb << std::endl;
  std::cout << " Number of events = " << nEvent << std::endl;

  vector<double> stat1;
  stat1.push_back(sigmapb);
  pbook->write("cross",stat1);

  vector<int> stat2;
  stat2.push_back(nEvent);
  pbook->write("events",stat2);


  // check number of entries
   std::cout << " Entries so far = " << pbook->size() << std::endl;

   // print what is written 
   pbook->print();

 // write histograms 
   pbook->write();
   pbook->close();

  return 0;
}