// ----------------------------------------------------------------------------
// $Id: root_GoGenMC.C,v 1.13 2001/10/10 16:52:57 toshi Exp $
// ----------------------------------------------------------------------------
//
//  Example file for event generation interactively with Root.
//  
//  In this file, 1) event generation (e+e- -> ZH) with pandora_pythia
//                2) out put a root file
//        
//                               Apr. 26 2001  Toshi
//
//   :: How to run ::
// 
//   1) type root
//   2)
//      root [0]  .x root_GoGenMC.C
//
// For example, if you like to generate 20 events, type like
//      root [0]  .x root_GoGenMC.C(20)

//__________________________________________________________________________
Int_t root_GoGenMC(Int_t   nEvent  = 100){

  gROOT->Reset();

  // In order to refer to shareable libraries in this simple form,
  // include their paths in a .rootrc file.

  gSystem->Load("libLCDEvent");
  gSystem->Load("libPandora");
  gSystem->Load("libPandoraSMProc");
  if (strcmp(gSystem->GetName(),"WinNT")) {
    gSystem->Load("libPythia6");
    gSystem->Load("libTauola");
  }
  gSystem->Load("libEG");
  gSystem->Load("libEGPythia6");
  gSystem->Load("libPandPyth");

  //------------ Out put file (Root)
  Char_t* outputfilename="test.root";
  TFile* ofile = new TFile(outputfilename,"RECREATE","LCD Event ROOT file");

  LCDEvent* event=new LCDEvent();
  Int_t comp    = 2;       // by default file is compressed
  ofile->SetCompressionLevel(comp);
  TTree* tree = new TTree("T","LCD ROOT tree");
  tree->SetAutoSave(1000000000);  // autosave when 1 Gbyte written      
  tree->Branch("LCDEvent","LCDEvent", &event);

  //Pyjet -> McPart converter
  // LCDPyj2McPart py2mc;
  LCDPandPyth2McPart py2mc;


  //------------< Setup Generator (Pandora-Pythia) >------------------------
  printf("Start Generator Setup! \n");

#include "pandpythconsts.h"

  Int_t iseed_pan = 7;  // default pandora initial seed

  /*    define the pandora event selection in this space     */
  
  double ECM  = 500.;
  double PolE = -0.8; // Polarization Electron
  double PolP =  0.0; // Polarization Positron

  ebeam b1(ECM/2.0, PolE, electron, electron);
  b1.setup(NLC500);
  b1.ISRoff();            // Turn off initial-state radiation
  b1.beamstrahlungoff();  // Turn off beamstrahlung

  ebeam b2(ECM/2.0, PolP, positron, positron);
  b2.setup(NLC500);
  b2.ISRoff();
  b2.beamstrahlungoff();

  double Hmass = 120.; // Higgs Mass
//double Hmass = 140.; // Higgs Mass
  eetoZHiggs pr(Hmass);

  //For example:: Z0 decays to anything. H0 decays to bb-bar only.
  pr.onlyDecay(allStates,bOnly); 

  pandora P(b1,b2,pr);
  
  pandorarun  PR(&P, nEvent, iseed_pan);
  
  //----- For PYTHIA ------
  //Or you can also change PYTHIA seed here.
  int iseed_pyth = 22846; PR.newseed(iseed_pyth); // Change Seed

  //PR.partonshoweroff();   // Turn off parton shower
  //PR.fragmentationoff();  // Turn off fragmentation
  //PR.FSRoff();            // Turn off final-state QCD QED radiation

  PR.StdHepIOoff();  // Don't use stdhep IO

  //       end of definition

  PR.initialize(0);
  
  // If you want/don't want to decay some specific particle(s),
  // you can call PartDecayOn/PartDecayOff here.   
  // For example: Do not allow to decay K0S and Lambda
  //int ipart1 = 310;  // K0S
  //int ipart2 = 3122; // Lambda
  //PR.PartDecayOff(ipart1);
  //PR.PartDecayOff(ipart2);

  printf("End of Generator Setup!\n");
  //------------< End Generator (Pandora-Pythia) setup>-------------------

  gBenchmark->Start("gen_eetoZH");

  Int_t iEvent;
  for (iEvent=0 ; iEvent < nEvent ; iEvent++) {
    PR.getevent();//Event Generation

    event->Clean();

    // Make McPart
    //py2mc.MakeMcPart(event->MCparticles());
    py2mc.MakeMcPart(&PR, event->MCparticles());

    // Store Event Header
    event->EventHeader()->SetEvtNum(iEvent);
    event->EventHeader()->SetProcessName(pr->name);

    // Store beam information
    event->BEAM()->SetEcm(ECM);
    event->BEAM()->SetPolElec(PolE);
    event->BEAM()->SetPolPosi(PolP);

    tree->Fill();      // Fill information into output root file
  }

  gBenchmark->Stop("gen_eetoZH");
  gBenchmark->Print("gen_eetoZH");

  PR.terminate();

  ofile->Write();
  tree->Print();
  ofile->Close();

  return iEvent;
}