#ifndef __CINT__
 #include "TPythia6.h"
 #include "TTree.h"
 #include "TClonesArray.h"
 #include "TH1.h"
 #include "TStyle.h"
 #include "TCanvas.h"
 #include "TMCParticle6.h"
 using namespace std;
 #endif

 void loadLibraries()
 {
 #ifdef __CINT__
  
  gSystem->Load("libEG");
  gSystem->Load("/home/brun/pythia6/libPythia6");
  gSystem->Load("libEGPythia6");
 #endif
 }

 void pythiatest2() {
 loadLibraries();

 TPythia6* pythia=TPythia6::Instance();

 pythia->SetupTest();
 pythia->Initialize("CMS","p","p", 900);

 const int nEvents=10;

// TClonesArray *fpart = new TClonesArray("TMCParticle");
 //TMCParticle *prt = new TMCParticle();
 
 TTree* part = new TTree("part","A tree of particles");
 TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
 part->Branch("particles", &particles);
 
 for (Int_t i = 0; i < nEvents; i++) {
   //fpart->Clear();
   pythia->GenerateEvent();
   
   //TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
   
   for (Int_t j=0; j<particles->GetEntries(); j++) {
     prt=(TMCParticle*)particles->At(j);
     
     Int_t fKF = 0, Int_t fKS = 0 ;
     fKF=prt->GetKF();
     fKS=prt->GetKS();
     
     Int_t finalstate;
     //the next line is wrong
     //finalstate = TDatabasePDG::Instance()->GetParticle(fKS);
     //it should be
     TParticlePDG *pdg = TDatabasePDG::Instance()->GetParticle(fKS);
     
     Int_t charge;
     charge  = TDatabasePDG::Instance()->GetParticle(fKF)->Charge();
     printf("  prt %d\n, charge=%d\n", j,charge);
         
     if (charge == 0) continue;
     //if (finalstate != 1) continue;
     
   }
  part->Fill();
 }
 part->Print();
 part->StartViewer();    
 
  
  TH1D* hist1 = new TH1D("for 10 events", "charged particles d#sigma/dp_{#perp} spectrum for 14000GeV", 50, 0, 5);
  hist1->SetXTitle("p_{#perp} [GeV]");
  hist1->SetYTitle("d#sigma/dp_{#perp}   (#mub/GeV)");
  part->Draw("sqrt(pow(particles.fPx,2)+pow(particles.fPy,2))>>for 10 events");  
  hist1->Draw();
  
}