// main01.cc is a part of the PYTHIA event generator.
// Copyright (C) 2013 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 fits on one slide in a talk. 
// It studies the charged multiplicity distribution at the LHC.

#include "Pythia.h"
using namespace Pythia8; 
int main() {
  // Generator. Process selection. LHC initialization. Histogram.
  Pythia pythia;

  // pythia.readString("Main:timesAllowErrors = 500");
  // pythia.readString("6:m0 = 172.5");
  // pythia.readString("23:m0 = 91.1876");
  // pythia.readString("23:mWidth = 2.4952");
  // pythia.readString("24:m0 = 80.399");
  // pythia.readString("24:mWidth = 2.085");
  // pythia.readString("StandardModel:sin2thetaW = 0.23113");
  // pythia.readString("StandardModel:sin2thetaWbar = 0.23146");
  // pythia.readString("ParticleDecays:limitTau0 = on");
  // pythia.readString("ParticleDecays:tau0Max = 10.0");
    
  // pythia.readString("Tune:pp = 5");
  // //  pythia.readString("PDF:useLHAPDF = on");
  // //  pythia.readString("PDF:LHAPDFset = CT10.LHgrid");
  // pythia.readString("MultipartonInteractions:bProfile = 4");
  // pythia.readString("MultipartonInteractions:a1 = 0.10");
  // pythia.readString("MultipartonInteractions:pT0Ref = 1.70");
  // pythia.readString("MultipartonInteractions:ecmPow = 0.16");
  // pythia.readString("BeamRemnants:reconnectRange = 4.67");
  // pythia.readString("SpaceShower:rapidityOrder=0");

  pythia.readString("Beams:eCM = 14000.");    
  pythia.readString("Beams:idA = 2212");
  pythia.readString("Beams:idB = 2212");

  
  pythia.readString("25:m0 = 125");
  pythia.readString("36:m0 = 500");

  pythia.readString("Higgs:useBSM = on");
  pythia.readString(" HiggsBSM:gg2A3 = on");
  pythia.readString("PhaseSpace:mHatMin = 20.");
  pythia.readString("25:onMode = off");
  pythia.readString("25:onIfMatch = 5 5");
  pythia.readString("36:onMode = off");
  pythia.readString("36:onIfMatch = 23 25");
  pythia.readString("23:onMode = off");
  pythia.readString("23:onIfMatch = 13 13");
  

  pythia.particleData.listChanged();

  // Create an LHAup object that can access relevant information in pythia.
  LHAupFromPYTHIA8 myLHA(&pythia.process, &pythia.info);

  // Open a file on which LHEF events should be stored, and write header.
  myLHA.openLHEF("signal_A_Zh.lhe"); 



  pythia.init();
  //pythia.particleData.listChanged();
  //return 0;
  //  Hist mult("charged multiplicity", 100, -0.5, 799.5);
  // Begin event loop. Generate event. Skip if error. List first one.

  // Store initialization info in the LHAup object. 
  myLHA.setInit();

  // Write out this initialization info on the file.
  myLHA.initLHEF();



  for (int iEvent = 0; iEvent < 10; ++iEvent) {
    if (!pythia.next()) continue;
    // Find number of all final charged particles and fill histogram.
    //int nCharged = 0;
    cout << "Event number: " << iEvent << endl;
    for (int i = 0; i < pythia.event.size(); ++i) {
      //      if (pythia.event[i].isFinal() && abs(pythia.event[i].id()) == 7 ) {
      if (abs(pythia.event[i].id()) == 7 and 
	  pythia.event[i].id() != pythia.event[ pythia.event[i].daughter1()  ].id()) {
	cout << "Found " << pythia.event[i].id() << " status: " 
	     << pythia.event[i].status()
	     << " with daughter: "
	     << pythia.event[ pythia.event[i].daughter1() ].id() << " and "
	     << pythia.event[ pythia.event[i].daughter2() ].id() 
	     << endl;
      }
      if (abs(pythia.event[i].id()) == 6 and 
	  pythia.event[i].id() != pythia.event[ pythia.event[i].daughter1()  ].id()) {
	cout << "Found " << pythia.event[i].id() << " status: " 
	     << pythia.event[i].status()
	     << " with daughter: "
	     << pythia.event[ pythia.event[i].daughter1() ].id() << " and "
	     << pythia.event[ pythia.event[i].daughter2() ].id() 
	     << endl;
      }

    }

      // Store event info in the LHAup object. 
      myLHA.setEvent();

      // Write out this event info on the file. 
      // With optional argument (verbose =) false the file is smaller.
      myLHA.eventLHEF();

    //    mult.fill( nCharged );
  // End of event loop. Statistics. Histogram. Done.
  }
  pythia.stat();

  // Update the cross section info based on Monte Carlo integration during run.
  myLHA.updateSigma();

  // Write endtag. Overwrite initialization info with new cross sections.
  myLHA.closeLHEF(true);


  //cout << mult; 
  return 0;
}