int status1 = gSystem->Load("/star/data01/pwg/pusheax/lib/libLHAPDF.so");
int status = gSystem->Load("/star/data01/pwg/pusheax/simulation/pythia6/libPythia6.so");
/*********************************************************
 *                                                       *
 * Macros for p+p 200 GeV simulation on Pythia 6.4.28.   *
 *                                                       *
 * first argument is an output file                      *
 * second argument is  a tune id see (arxiv:1005.3457v5) *
 *  if mTuneId = -1 then vanilla pythia 6.4.28 are used  *
 *                                                       *
 * Some useful tune id's:                                *
 *    TUNE                 |  ID                         *
 *  --------------------------------                     *
 *   Perugia 2011          |  350                        *
 *   Perugia 2011 noCR     |  354                        *
 *   Perugia 0             |  320                        *
 *   Tune A (used by Zibi) |  100                        *
 *                                                       *
 *********************************************************/
#include "TROOT.h"
#include "TSystem.h"
#include "TChain.h"
#include "TFile.h"
#include "TString.h"
#include "TPythia6.h"
#include "TRandom3.h"
#include <iostream>

//
// Uncomment line below to turn on verbose output
//
// #define MY_DEBUG

const Char_t *defaultOutFile = "oFemtoDst_pp200_pythia6.root";


void FemtoDstMaker_mesons_pp200_pythia6(const Char_t* outFileName = defaultOutFile,
                                        Int_t mTuneId = -1)
{
  std::cout << "**************************" << std::endl
            << "*    FemtoDstMaker       *" << std::endl
            << "*        PYTHIA6         *" << std::endl
            << "*   p+p 200 GeV mesons   *" << std::endl
            << "*         Start          *" << std::endl
            << "**************************" << std::endl << endl;

  //Constants, cut and initial parameters
  Int_t    mNEvents2Gen = 50000;
  Double_t mCollisionEnergy = 200.;

  Bool_t  mIsFullField = false;

  Float_t cParticleMom[2] = {0.15, 1.55};
  Float_t cTrackDcaGlobal[2] = {-0.01, 2.};
  Float_t cTrackEta[2] = {-1., 1.};
  Bool_t  cRemoveProtons = true;

  //
  // Load libraries
  //
  std::cout << "Loading libraries..." << std::endl;
  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();
  gSystem->Load("libMinuit");
  gSystem->Load("StRefMultCorr");
  gSystem->Load("StFlowMaker");   // should be placed before HBT
  gSystem->Load("StHbtMaker");
  gSystem->Load("StarClassLibrary");
  gSystem->Load("libgsl");
  gSystem->Load("libgslcblas");
  gSystem->Load("StPicoDstMakerRun11");
  gSystem->Load("StFemtoDstMaker");


  std::cout << "Libraries have been successfully loaded" << std::endl;
  //
  // Create chain
  //
  StChain *mChain = new StChain("StChain");
  mChain->SetDebug(0);
  StMuDebug::setLevel(0);
  //
  // Random initialization
  //
  TRandom3 randy(0);
  Int_t mSeed = floor(randy.Uniform(1., 100000.));
  std::cout << "Random seed: " << mSeed << std::endl;
  //
  // TPythia6 initialization
  //
  TPythia6* mTPythia = new TPythia6();
  mTPythia->SetMSEL(1);                //mbias
  mTPythia->SetMRPY(1, mSeed);
  if (mTuneId != -1)   mTPythia->SetMSTP(5, mTuneId);
  mTPythia->Initialize("CMS", "p", "p", mCollisionEnergy);
  mTPythia->Pyexec();
  //
  // StFemtoDstPythia6Maker initialization
  //
  StFemtoDstPythia6Maker *mFemtoDstMaker =
      new StFemtoDstPythia6Maker(mTPythia, outFileName);

  mFemtoDstMaker->SetIsFullField(mIsFullField);
  mFemtoDstMaker->SetParticleMomentum(cParticleMom[0], cParticleMom[1]);
  mFemtoDstMaker->SetTrackEta(cTrackEta[0], cTrackEta[1]);
  mFemtoDstMaker->SetTrackDca(cTrackDcaGlobal[0], cTrackDcaGlobal[1]);
  mFemtoDstMaker->SetTrackDcaGlobal(cTrackDcaGlobal[0], cTrackDcaGlobal[1]);
  mFemtoDstMaker->SetRemoveProtons(cRemoveProtons);
  //
  // Chain initialization and loop over events
  //
  mChain->Init();
  Int_t iReturn = 0;
  Int_t mNEventsProcessed = 0;
  Float_t mPercentCounter = 0.05;
  Float_t mProgress = 0.;

  unsigned int mSubCounter = 1;              // time
  time_t mStartTime, mStopTime, mDiffTime;   // time
  float mFrac;                               // time
  mStartTime = time(0);                      // time

  for(Int_t iEvent = 0; iEvent < mNEvents2Gen; iEvent++)
  {
    mProgress = (Float_t)iEvent/(Float_t)mNEvents2Gen;
    mNEventsProcessed++;

    if(mProgress >= mPercentCounter)
    {
      mPercentCounter += 0.05;
      mStopTime = time(0);
      mDiffTime = difftime(mStopTime, mStartTime);
      mFrac = (float)mDiffTime*(float)(mNEvents2Gen - iEvent)/(float)iEvent;
      mSubCounter++;
      std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
                << std::endl;
    }

    mChain->Clear();
    iReturn = mChain->Make(iEvent);
#ifdef MY_DEBUG
    mTPythia->Pylist(2);
#endif
    if(iReturn != 0)
    {
      std::cout << "Error has been occured. Event processing has been stopped"
          << std::endl;
      break;
    }
  }

  mChain->Finish();

  delete mFemtoDstMaker;
  delete mChain;

  std::cout << "***********************" << std::endl
            << "*       Finish        *" << std::endl
            << "***********************" << std::endl << std::endl;
}