void test_p(int EVT=5, const char *fname="test_p.root")  { 

  gROOT.Reset();

  gSystem.Load("libEG.so");
  gSystem.Load("libPythia6.so");
  gSystem.Load("libEGPythia6");

  TPythia6  pythia;
  //

  //PDFs
  //    pythia.SetMSTP(51,5005);  //1000xNGROUP+NSET
  pythia.SetMSTP(51,7);  //4=GRV94L; 7=CTEQ5L
  pythia.SetMSTP(52,1);  // 1 =pythia internal; 2 =external PDFLIB
  //K factor
  pythia.SetMSTP(33,1); // use commom K factor 
  pythia.SetPARP(31,1); // K=1.0

  //Process
  //  pythia.SetMSEL(5); //1,2=QCD, c-cbar =4; b-bar=5; Z=11(ISUB=?),W=12(ISUB=2)
  //  pythia.SetMSUB(2,1); // for single W,Z production

  pythia->SetMSEL(0); //
  pythia->SetMSUB(86,1); // gg->J/Psi+g 
  //  pythia->SetMSUB(87,1); // gg->chi0c+g (=10441)
  //  pythia->SetMSUB(88,1); // gg->chi1c+g (=20443)
  //  pythia->SetMSUB(89,1); // gg->chi2c+g (=445)

  ///////////////////////////////////////////////////
  //Global propety of event
  //Kinematic cuts
  ///////////////////////////////////////////////////
  pythia.SetCKIN(3,0.5); // pt_min=0.5Gev
  pythia.SetCKIN(4,-1); // pt_max=none

  //Kt factor 
  pythia.SetPARP(62,1.0); // Q=Kt=1.0GeV
  pythia.SetPARP(64,1.0); // Kt scale factor
  pythia.SetMSTP(64,2); // Kt**2=PARP(64)*(1-z)Q**2

  //Particle decay
  int mu_KC=pythia.Pycomp(13); // get KC code from KF
  int e_KC=pythia.Pycomp(11); // get KC code from KF
  int pi_KC=pythia.Pycomp(211); // get KC code from KF
  int kp_KC=pythia.Pycomp(321); // K+-
  int ks0_KC=pythia.Pycomp(310); // K_S0
  int kl0_KC=pythia.Pycomp(130); // K_L0

  pythia.SetMSTJ(22,4); //  decay vertex is within cylindrical volum
  pythia.SetPARJ(73,5000.0); // r=1000mm; DC_r=2.0m; EMCal=5.0m
  pythia.SetPARJ(74,400.0);  // z=+-400mm; 

  //
  

 //allow pi/kaon/muon decay in the volume
  cout <<"----- Particle decays -------" << endl;
  pythia.SetMDCY(mu_KC,1,1); 
  cout <<"muon MDCY(1) =" <<  pythia.GetMDCY(mu_KC,1) << endl;
  pythia.SetMDCY(pi_KC,1,1); 
  cout <<"Pi MDCY(1) =" <<  pythia.GetMDCY(pi_KC,1) << endl;
  pythia.SetMDCY(kp_KC,1,1); 
  cout <<"kaon MDCY(1) =" <<  pythia.GetMDCY(kp_KC,1) << endl;
  pythia.SetMDCY(ks0_KC,1,1); 
  cout <<"kaon_S0 MDCY(1) =" <<  pythia.GetMDCY(ks0_KC,1) << endl;
  pythia.SetMDCY(kl0_KC,1,1); 
  cout <<"kaon_L0 MDCY(1) =" <<  pythia.GetMDCY(kl0_KC,1) << endl;

  //check some partricles lifetime
  cout << "------ Particle Lifetime: (cm) " << endl; 
  cout <<"Muon c*t = " <<  pythia.GetPMAS(mu_KC,4)/10.0 << endl;
  cout <<"pi c*t = " <<  pythia.GetPMAS(pi_KC,4)/10 << endl;
  cout <<"K+- c*t = " <<  pythia.GetPMAS(kp_KC,4)/10 << endl;
  cout <<"K_S0 c*t = " <<  pythia.GetPMAS(ks0_KC,4)/10 << endl;
  cout <<"K_L0 c*t = " <<  pythia.GetPMAS(kl0_KC,4)/10 << endl;

  cout << "Particle decay volume R (cm) = " << pythia.GetPARJ(73)/10 <<endl;
  cout << "Particle decay volume Z (cm) = " << pythia.GetPARJ(74)/10 <<endl;


  // 200GeV p-p run  
  pythia.Initialize("CMS","p","p",200);

/////////////////////////////////////////////////////////////
  //---------------------------- event loop -------------------
  /////////////////////////////////////////////////////////////

  //generate a random number 
  //  gRandom.SetSeed(1);//default = 65539
  
  //define my own random generator
  //  delete gRandom;
 
 //  gRandom = new TRandom(0); // old one
  gRandom = new TRandom3(0); // new one, better one, use time as the seed

  Float_t randx = 0; // gRandom.Rndm(1);

  randx =  gRandom.Rndm();

  Int_t  Iseed =  (int) 1000000*randx;
  cout << "Iseed = " << Iseed << "\t";

  pythia.SetMRPY(1,Iseed); // default=19780503
  cout << "pythia.GetMRPY(1) = " << pythia.GetMRPY(1) << endl;
  pythia.SetMRPY(2,0); // default=0

  // --------- 
  int j=0, x_pid1=0, x_pid2=0;

  for (int i=0;i<EVT;i++){
    
    pythia.GenerateEvent();   
    pythia.ImportParticles();   

    //   pythia.Pytest(0); //  test program
    if (i<3){
      cout <<"----- Dump EVT = " << i << endl;   
      //    cout << " ---- PYLIST(1) ----" << endl;
      //      pythia.Pylist(1); // dump out all particles

      cout << " ---- PYLIST(3) ----" << endl;
      pythia.Pylist(1); // dump out all particles in details ==3;

      //     cout << " ---- PYLIST(11) ----" << endl;
      //  pythia.Pylist(11); // dump out all particles

    }

    //check output
    //kinematic analysis
    N_Par=pythia.GetN(); // number of particles generated in this event

    cout << "Total # of particle in this event =" << N_Par << endl;

    int LinePart2 = pythia.GetMSTI(7); // line numbers for documentation of outgoing partons/particles
    int LinePart1 = pythia.GetMSTI(8); // from hard scattering for 2->2 processes (else = 0)

    int flav1 = pythia.GetMSTI(21);
    int flav2 = pythia.GetMSTI(22);

    flaot x1=0, x2=0; //calculated from line3,4

    cout << "Outgoing paton lines = " << LinePart1 << "\t" << LinePart2 << endl;

    j = LinePart1;
    //    x_pid1 = TMath::Abs(pythia.GetK(j,2));
    x_pid1 = pythia.GetK(j,2);

    j = LinePart2;
    x_pid2 = pythia.GetK(j,2);

    cout << " >>>>> parton1 = " << x_pid1 << " \t parton2 =" << x_pid2 << endl;
    cout << " >>>>> parton1_f = " << flav1  << " \t parton2_f =" << flav2 << endl;

    cout << "PAR(33)x1 = " << pythia.GetPARI(33) << "\tPAR(34)x2 = " << pythia.GetPARI(34) << endl;

    cout << "Cal::  x1 = " <<  pythia.GetP(3,3)/99.9956 << "\t Cal:: x2 = " << pythia.GetP(4,3)/99.9956 << endl; ; //(Px, Py,Pz,E,M)
    cout << "Cal::  x1 = " <<  pythia.GetP(5,4)/99.9956 << "\t Cal:: x2 = " << pythia.GetP(6,4)/99.9956 << endl; ; //(Px, Py,Pz,E,M)

  }// event loop

}