//////////////////////////////////////////////////////////////////////////////////////////
// program to run TPythia6 under ROOT for cc/bb analysis
// This program runs very slow
// -----------------------------------------
// cb: QCD=0, bb =1;cc=2,JPsi=3,DY=4,W=5,Z=6
// -----------------------------------------
// subevt: 1=ele; 2=muon;3=kaon;0=all three
//
//OUTPUT:
// Q_flag: 1=b->x; 2 = c->x; 3 = b->c->x; 4 =b->c->s->x; 6 =s->x ; 0 = from IP
//
// History:
//        01/23/2002  Use cos_70,45,12,35 constants to apply detector acceptance cuts
//        10/21/2002  Include more processes for J/Psi produciton
//        02/14/2003  set single/dilepton stwitch, single = 0 run all; =1 single lepton only 
//
//        02/28/2003  replaced pythia event loop with with a table;
//                    added "forced particle decay" for pi/kaon in volum  z=40cm,r=500cm 
//        03/05/2003  add primary vertex smearing, dx=dy=0, dz =30cm      
//        07/02/2003  add DY, W and Z processes
/////////////////////////////////////////////////////////////////////////////////////////
void run_cc(int cb=2, int subevt=0,int EVT=100000, const char *fname="CC.root")  { 

  gROOT.Reset();

  hfile = new TFile(fname,"recreate","Heavy qurak analysis");

  single= new TNtuple("single","single-particle analysis","etype:q_px:q_py:q_pz:q_id:px:py:pz:E:ptjet:charge:pid:Q_flag:prt_id:prt_px:prt_py:prt_pz:prt_e:prt_m:r:z:vx:vy:vz");

  TNtuple *  mumu  = new TNtuple("mumu","muon-muon correlatoin","etype:px1:py1:pz1:flag1:charge1:pid1:prt_id1:px2:py2:pz2:flag2:charge2:pid2:prt_id2:mass:pt:vx:vy:vz");
  TNtuple * mue   = new TNtuple("mue","muon-electorn correlatoin","etype:px1:py1:pz1:flag1:charge1:pid1:prt_id1:px2:py2:pz2:flag2:charge2:pid2:prt_id2:mass:pt:vx:vy:vz");
  TNtuple * muk   = new TNtuple("muk","muon-kaon correlatoin","etype:px1:py1:pz1:flag1:charge1:pid1:prt_id1:px2:py2:pz2:flag2:charge2:pid2:prt_id2:mass:pt:vx:vy:vz");
  TNtuple * ee    = new TNtuple("ee","ele-ele correlatoin","etype:px1:py1:pz1:flag1:charge1:pid1:prt_id1:px2:py2:pz2:flag2:charge2:pid2:prt_id2:mass:pt:vx:vy:vz");
  TNtuple * ek    = new TNtuple("ek","ele-kaoon correlatoin","etype:px1:py1:pz1:flag1:charge1:pid1:prt_id1:px2:py2:pz2:flag2:charge2:pid2:prt_id2:mass:pt:vx:vy:vz");
  TNtuple * kk    = new TNtuple("kk","kaon-kaoon correlatoin","etype:px1:py1:pz1:flag1:charge1:pid1:prt_id1:px2:py2:pz2:flag2:charge2:pid2:prt_id2:mass:pt:vx:vy:vz");


  TNtuple * evt_g = new TNtuple("evt_g","Global Properties","x1:x2:pid1:pid2:p_px:p_py:p_pz:p_tot:p_pt:p_m:vx:vy:vz");
  TNtuple * evt_q = new TNtuple("evt_q","Global Properties","x1:x2:pid1:pid2:qid1:qid2:q_px1:q_py1:q_pz1:q_m1:q_px2:q_py2:q_pz2:q_m2:vx:vy:vz");

  c1 = new TCanvas("c1","Pythia output: muon/e/k info",10,20,600,500);
  c2 = new TCanvas("c2","Pythia output: parent Q-quark info",100,200,500,400);
  c3 = new TCanvas("c3","Total # of mu/e/k in QQ events",200,400,500,300);

  h1 = new TH1F("h1","mu/e/k ptot",100,0,20);
  h2 = new TH1F("h2","mu/e/k pt",50,0,10);
  h3 = new TH1F("h3","mu/e/k ID",10,0,5);
  h4 = new TH1F("h4","parent ID",600,0,600);
  h5 = new TH1F("h5","mu/e/k Cos_theta",100,-1,1);
  h6 = new TH1F("h6","Q_flag",10,0,5);
  h7 = new TH1F("h7","number of mu/e/k",10,-0.5,9.5);

  h21 = new TH1F("h21","Q(c,b) ptot",50,0,50);
  h22 = new TH1F("h22","Q(c,b) pt",50,0,10);
  h23 = new TH1F("h23","Q(c,b) P_ID",100,-25,25);
  h24 = new TH1F("h24","Q(c,b) Pz",50,-50,50);

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

  TPythia6  pythia;
  //  TPythia  pythia;

  ////////////////////////////////////////////////////////////////////////////////////////////////
  //define user functions
  //For a given child-particle line #, returns Original Qaurk's line # parent_Q 
  //and process type Q_flag=1 for b->l-; 2 for c->l+; 3 for b->c->l+; 0 for others
  ///////////////////////////////////////////////////////////////////////////////////////////////

  int parent_Q(TPythia6 *pythia, int child, int &Q_flag){

    int N_Par,Parent_ID,Parent_KS,Parent_KF,Parent_KF1,origin;
    int abs_parent, c_flag,b_flag,k_flag;

    c_flag=0;
    b_flag=0;
    k_flag=0;
    Q_flag=0;

    N_Par=pythia.GetN();
    Parent_ID = pythia.GetK(child,3);  // parent line #
    Parent_KS=pythia.GetK(Parent_ID,1); // status code KS
    Parent_KF=pythia.GetK(Parent_ID,2); // status code KF
    Parent_KF1=pythia.GetK(Parent_ID,2); // immed parent
    origin =pythia.GetK(Parent_ID,3);   // grand parent line #

    abs_parent=TMath::Abs(Parent_KF1); // immed-parent = mesons and baryons
    //C hadron
    if ((abs_parent>400&&abs_parent<500)||(abs_parent>4000&&abs_parent<5000)){
      c_flag=1;
    }
    //B hadron
    else if((abs_parent>500&&abs_parent<600)||(abs_parent>5000&&abs_parent<6000)){
      b_flag=1;
    }
    while(origin>2){ // find origin by looping parents
      Parent_ID = pythia.GetK(Parent_ID,3); // parent line #
      Parent_KS=pythia.GetK(Parent_ID,1); // status code KS
      Parent_KF=pythia.GetK(Parent_ID,2); // status code KF
      origin =pythia.GetK(Parent_ID,3); // grand parent line #

      abs_parent=abs(Parent_KF); // immed-parent = mesons and baryons
      if ((abs_parent==4)||(abs_parent>400&&abs_parent<500)||(abs_parent>4000&&abs_parent<5000)){
	c_flag=1;
      }
      //B hadron
      else if((abs_parent==5)||(abs_parent>500&&abs_parent<600)||(abs_parent>5000&&abs_parent<6000)){
	b_flag=1;
      }
    }
    if (b_flag==1&&c_flag==0){Q_flag=1;} 
    if (b_flag==0&&c_flag==1){Q_flag=2;} 
    if (b_flag==1&&c_flag==1){Q_flag=3;} 
    if (b_flag==0&&c_flag==0){Q_flag=0;} 
   /*
    cout << "immed_parent = " <<abs_parent <<endl;
    cout << "Parent_ID = " << Parent_ID << "\t Parent_KS =" << Parent_KS; 
    cout << "\t Parent_KF =" << Parent_KF  << "\t Q_flag = " << Q_flag << endl;
    pythia.Pylist(1); // dump out all particles
    */
    return Parent_ID;
  } // end of parent_Q()    
  ////////////////////////////////////////////////////////////////////////////////////////////////

 
  //------set pythia parameters
  if (cb==1){  // bb
  
    //PDFs
    //    pythia.SetMSTP(51,5005);  //1000xNGROUP+NSET
    pythia.SetMSTP(51,4);  //4=GRV94L
    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
  }
  else if (cb==2){  // cc
    //PDFs 
    pythia.SetMSTP(51,4);  //4=GRV94L
    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(4); //1,2=QCD, c-cbar =4; b-bar=5; Z=11(ISUB=1),W=12(ISUB=2)
    //  pythia.SetMSUB(2,1); // for single W,Z production  
  }
  else if (cb==3){ // J/Psi(1S) etc.
    //PDFs 
    //    pythia.SetMSTP(51,5005);  //1000xNGROUP+NSET
    pythia.SetMSTP(51,4);  //4=GRV94L
    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(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)
    //  pythia.SetMSUB(106,1); // gg.J/Psi+gamma 
    //  pythia.SetMSUB(107,1); // ggamma.J/Psi+g 
    //  pythia.SetMSUB(108,1); // gamma gamma.J/Psi+gamma 
  }
  else if (cb==0){ // QCD
    //PDFs 
    //    pythia.SetMSTP(51,5005);  //1000xNGROUP+NSET
    pythia.SetMSTP(51,4);  //4=GRV94L
    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(2); //1= hard scattering; 2=minimum bias
    //    pythia.SetMSUB(86,1); // gg.J/Psi+g 
  }
  else if (cb==4){ // DY process f+fbar.garmma.XXbar

    //pythia.SetMSEL(11); //1,2=QCD, c-cbar =4; b-bar=5; Z=11(ISUB=1),W=12(ISUB=2)
    pythia.SetMSEL(0);
    pythia.SetMSUB(1,1); // for single garmma/Z production
    pythia.SetMSTP(43,1); // 1=garmma only; 2= Z0 only; 3 = Z0 and garmma

  }
  else if (cb==5){ // W+- production
    //    pythia.SetMSEL(12); //1,2=QCD, c-cbar =4; b-bar=5; Z=11(ISUB=?),W=12(ISUB=2)
    pythia.SetMSEL(0);
    pythia.SetMSUB(2,1); //f+f'.W+-
    pythia.SetMSUB(16,1); //f+f'.W+g(jet)
    pythia.SetMSUB(31,1); //f+g.W+f'(jet)
  }
  else if (cb==6){ // Z0 production

    //pythia.SetMSEL(11); //1,2=QCD, c-cbar =4; b-bar=5; Z=11(ISUB=1),W=12(ISUB=2)
    pythia.SetMSEL(0);
    pythia.SetMSUB(1,1); // for single Z production
    pythia.SetMSUB(15,1); //f+f'.Z+g(jet)
    pythia.SetMSUB(30,1); //f+g.Z+f(jet)
    
    pythia.SetMSTP(43,3); // 1=gamma only, 2=Z only, 3(D)=Z/gamma all
  }





  ///////////////////////////////////////////////////
  //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;

  //B0-mixing parameters
  pythia.SetMSTJ(26,2); //
  pythia.SetPARJ(76,0.75); // Bd: Xd=0.75
  pythia.SetPARJ(77,15.0); // Bs: Xs=15


  //primary vertex smearing
  pythia.SetMSTP(151,1); //use vtx smearing
  pythia.SetPARP(151,1); // dx, mm
  pythia.SetPARP(152,1); // dy, mm
  pythia.SetPARP(153,300); // dz, mm?
  pythia.SetPARP(154,0); // t, time?? mm/c?


  //dump some parameters  
  cout << "------- Pythia parameters -----"  << endl;
  //PDF
  cout << "pythia.GetMSTP(51) =" <<pythia.GetMSTP(51) << endl;
  cout << "pythia.GetMSTP(52) =" <<pythia.GetMSTP(52) << endl;
  //K factor
  cout << "pythia.GetMSTP(33) =" <<pythia.GetMSTP(33) << endl;
  cout << "pythia.GetPARP(31) =" <<pythia.GetPARP(31) << endl;
  //process
  cout << "pythia.GetMSEL() =" <<pythia.GetMSEL() << endl;

  // 200GeV p-p run  
  pythia.Initialize("CMS","p","p",200);
  //////////////////////////////////////////////////////////////////////////////////////
  ////// End of PYTHIA setup
  ///////////////////////////////////////////////////////////////////////////////////////


  //  int EVT =1000;
  double mass_e=0.000511;
  double mass_m=0.105658;
  double mass_k=0.493677;

  double mass_e2=mass_e*mass_e; 
  double mass_m2=mass_m*mass_m;
  double mass_k2=mass_k*mass_k;

  int N_Par=0,P_KS,pid1,pid2;
  int particle_ID,particle_ID2,parentI,parentId;
  double px,py,pz,E,mass,ptot,pt,charge,ptjet;
  double vx,vy,vz,vr,prt_id;
  double v0x,v0y,v0z;
  double p1_px,p1_py,p1_pz,p1_E,p1_mass,p1_ptot,p1_pt,p1_id;  // of 1st parent
  double p2_px,p2_py,p2_pz,p2_E,p2_mass,p2_ptot,p2_pt;  // of grandparent(original quark)
  double cosx,abs_cosx,sinx,abs_sinx;
  int  n_good; // number of good lepton/kaons in this events
  int x_pid;
  int Is_muon,Is_ele,Is_kaon,Q_flag;
  int Is_muon2,Is_ele2,Is_kaon2,Q_flag2;
  int MAX_N=100;
  double  varx[100][17]; // [max_particle=MAX_N][property]

  //for particle correlations
  double px_e,py_e,pz_e,E_e,mass_e,ptot_e,pt_e,charge_e,ptjet_e;
  double px_e1,py_e1,pz_e1,E_e1,mass_e1,ptot_e1,pt_e1,charge_e1,ptjet_e1;
  double px_e2,py_e2,pz_e2,E_e2,mass_e2,ptot_e2,pt_e2,charge_e2,ptjet_e2;
  double px_m,py_m,pz_m,E_m,mass_m,ptot_m,pt_m,charge_m,ptjet_m;
  double px_m1,py_m1,pz_m1,E_m1,mass_m1,ptot_m1,pt_m1,charge_m1,ptjet_m1;
  double px_m2,py_m2,pz_m2,E_m2,mass_m2,ptot_m2,pt_m2,charge_m2,ptjet_m2;
  double px_k,py_k,pz_k,E_k,mass_k,ptot_k,pt_k,charge_k,ptjet_k;
  int    qflag_m,qflag_k,qflag_e,qflag_m2,qflag_e2,n_mu,n_e,n_k;
  double flag_m,flag_k,flag_e,flag_m2,flag_e2;
  double Q_m,Q_k,Q_e,Q_m2,Q_e2;
  double mass_em,mass_km,mass_ek,mass_mm,mass_ee;

  double prt_px,prt_py,prt_pz,prt_e,prt_m;

  //global properties
  double x1,x2,p_px,p_py,p_pz,p_m;
  double qid1,q_px1,q_py1,q_pz1,q_m1;
  double qid2,q_px2,q_py2,q_pz2,q_m2;

  float vf_s[24]; // single ntuple
  float vf_g[13];
  float vf_q[17];
  float vf_ll[20];

  ///////////////////////
  //PHENIX Detector acceptance cuts: for electron/muon/kaon
  //////////////////////
  const  double cos_theta_min = 0.9781; // (12 degree), muon acceptance
  const  double cos_theta_max = 0.8192; // (35 degree)
  const  double cos_e_max = 0.3420;  // (+-20 central arms)
  const  double sin_e_max = 0.7071;
  const  double pt_cut_e = 0.5 ; // electron and kaon
  const  double ptot_cut_m = 1.5; // muon = 1.5 
  const  double ptot_cut = 0.5; // for ALL particles = 0.5
  const  double pt_cut = 0.2; // for ALL particles = 0.2
  const  double E_cut = 0.5;  // for all particles = 0.5
  ///////////////////////////////////////////////////////////



  /////////////////////////////////////////////////////////////
  //---------------------------- 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


  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(3); // dump out all particles in details

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

    }

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

    // primary vtx position
    v0x =pythia.GetV(1,1);
    v0y =pythia.GetV(1,2);
    v0z =pythia.GetV(1,3);

    if(i%1000==0){
      cout <<"analyzing event " << i ;
      cout << " \t N_Par =" << N_Par << endl;
    }


    n_good =0; 
    n_mu=n_e=n_k=0;

    //loop all particles in this event
    for (int j=1; j<=N_Par;j++){ // J starts from 1

      //loop only stable and high Energy particles
      //electron, muon, Kaon
      P_KS=pythia.GetK(j,1); // status code KS
      E  =  pythia.GetP(j,4); //(Px, Py,Pz,E,M)
      //      if (P_KS!=1||E<E_cut){  // E>0.5GeV
      if (E<E_cut){  // E>0.5GeV, be careful with pi/kaon decays
	continue;
      }

      //check production vtx position, mainly for decay muons     
      vx =pythia.GetV(j,1); // production position
      vy =pythia.GetV(j,2);
      vz =pythia.GetV(j,3);
      vr=sqrt(vx*vx+vy*vy);

      px =  pythia.GetP(j,1); //(Px, Py,Pz,E,M)
      py =  pythia.GetP(j,2); //(Px, Py,Pz,E,M)
      pz =  pythia.GetP(j,3); //(Px, Py,Pz,E,M)
      E  =  pythia.GetP(j,4); //(Px, Py,Pz,E,M)
      mass  =  pythia.GetP(j,5); //(Px, Py,Pz,E,M)

      pt=sqrt(px*px+py*py);
      ptot = sqrt(px*px+py*py+pz*pz);
      charge=pythia.GetK(j,2)/abs(pythia.GetK(j,2));

      sinx=abs_sinx=1;
      //apply minimum ptot for ALL particles
      if (ptot>ptot_cut && pt>0.1){
	cosx=TMath::Abs(pz/ptot); 
	sinx=py/pt;
	abs_cosx=TMath::Abs(pz/ptot);
	abs_sinx=TMath::Abs(py/pt);
      }
      else {
	continue;
      }

     
      //------- PYTHIA PID----------
      // e =11; muon=13;W=24;Z=23;phonton=22,gluon=21
      // pi+=211;
      //----------------------------
      x_pid=TMath::Abs(pythia.GetK(j,2));

      Is_muon=0;
      Is_ele =0; 
      Is_kaon=0;
      Q_flag=0;
      particle_ID=0;

      // electron acceptance
      if(x_pid==11 && pt>pt_cut_e && abs_cosx<=cos_e_max && abs_sinx<=sin_e_max){ // +-20degree cos70=0.3420; sin(45)=0.707
	Is_ele=1;
	particle_ID=1;
      }   
      //muon acceptance in South arm only
      if(x_pid==13 && ptot>ptot_cut_m && abs_cosx>cos_theta_max && abs_cosx<cos_theta_min){ // cos(12)=0.978; cos(35)=0.819
	Is_muon=1;
	particle_ID=2;}
      // kaon acceptance
      if(x_pid==321 && pt>pt_cut_e && abs_cosx<cos_e_max && sinx<sin_e_max && sinx>0){ // +-20degree i ntheta; 0-45 in phi.TOF
	Is_kaon=1;
	particle_ID=3;
      }  
    
      //select subevent type
      if (subevt==1){ // ele only
	//  Is_ele =0; // run muon only 
	Is_muon=0; // muon only
	Is_kaon=0;
      }
      else if 
	(subevt==2){ // moun only
	Is_ele =0; // run muon only 
	// Is_muon=0; // muon only
	Is_kaon=0;
      }
      else if (subevt==3){ // kaon only
	Is_ele =0; // run muon only 
	Is_muon=0; // muon only
	//Is_kaon=0;
      }
      else if (subevt==0){ // keep all particles 
	//  Is_ele =0; // run muon only 
	//  Is_muon=0; // muon only
	//  Is_kaon=0;
      }

      //if any electron/muon/kaon found
      if (Is_muon||Is_ele||Is_kaon){ // e=11,muons=13,pi+=211,K+=321
		
	//lepton/kaon's immediate parent's information (pi,kaon etc)
	parentI =  pythia.GetK(j,3); // hadron parent particle Line #
	parentId= pythia.GetK(parentI,2);  // parent flavor
	prt_id = (float) abs(parentId);

	prt_px  =  pythia.GetP(parentI,1);   //(Px, Py,Pz,E,M)
	prt_py  =  pythia.GetP(parentI,2);   //(Px, Py,Pz,E,M)
	prt_pz  =  pythia.GetP(parentI,3);   //(Px, Py,Pz,E,M)
	prt_e  =  pythia.GetP(parentI,4);   //(Px, Py,Pz,E,M)
	prt_m  =  pythia.GetP(parentI,5);   //(Px, Py,Pz,E,M)


	//find lepton/kaon's parent's (primary quark) line #, e.g. b in b.c.l process  
	parentI = parent_Q(&pythia,j,Q_flag); 

	p1_px    =  pythia.GetP(parentI,1);   //(Px, Py,Pz,E,M)
	p1_py    =  pythia.GetP(parentI,2);   //(Px, Py,Pz,E,M)
	p1_pz    =  pythia.GetP(parentI,3);   //(Px, Py,Pz,E,M)
	p1_E     =  pythia.GetP(parentI,4);    //(Px, Py,Pz,E,M)
	p1_mass  =  pythia.GetP(parentI,5); //(Px, Py,Pz,E,M)
	p1_id =  pythia.GetK(parentI,2);  // parent flavor

	p1_pt=sqrt(p1_px**2+p1_py**2);
	p1_ptot=sqrt(p1_pt**2+p1_pz**2);
	
	//lepton/kaon's Pt relative the jet axis (parent quark's momentum direction)
	ptjet=ptot*sqrt(1-(p1_px*px+p1_py*py+p1_pz*pz)/ptot/p1_ptot);
	

	//fill varaiables
	if (n_good<=MAX_N){

	  varx[n_good][0]=cb; //event type
	  varx[n_good][1]=p1_px; //primary parent(quark or gluon)'s momentum
	  varx[n_good][2]=p1_py;
	  varx[n_good][3]=p1_pz;
	  varx[n_good][4]=px; //lepton/kaon 's momentum
	  varx[n_good][5]=py;
	  varx[n_good][6]=pz;
	  varx[n_good][7]=E; //
	  varx[n_good][8]=ptjet; //lepton's Pt relative to primary parent's flight directino
	  varx[n_good][9]=charge;
	  varx[n_good][10]=particle_ID; // lepton/kaon'd PID - ele=1;muon=2;kaon=3
	  varx[n_good][11]=Q_flag; // lepton/kaon's weak-decay flavor
	  varx[n_good][12]=prt_id; // lepton/kaon's immed parent particle ID
	  varx[n_good][13]=vr; // production vtx - vr
	  varx[n_good][14]=vz; // production vtz - vz	
	  varx[n_good][15]=mass;  //  

	}
	else {
	  cout << "number of leptons and kaons = " << n_good << " exceeds " << MAX_N << " !!! " << endl;
	}

	//fill ntuple
	//  single= new TNtuple("single","single-particle analysis","etype:q_px:q_py:q_pz:q_id::px:py:pz:E:ptjet:charge:pid:Q_flag:prt_id:prt_px:prt_py:prt_pz:prt_m:r:z:vx:vy:vz");
	vf_s[0]=cb;
	vf_s[1]=p1_px;
	vf_s[2]=p1_py;
	vf_s[3]=p1_pz;
	vf_s[4]=p1_id;
	vf_s[5]=px;
	vf_s[6]=py;
	vf_s[7]=pz;
	vf_s[8]=E;
	vf_s[9]=ptjet;
	vf_s[10]=charge;
	vf_s[11]=particle_ID;
	vf_s[12]=Q_flag;
	vf_s[13]=prt_id;

	vf_s[14]=prt_px;
	vf_s[15]=prt_py;
	vf_s[16]=prt_pz;
	vf_s[17]=prt_e;
	vf_s[18]=prt_m;
	vf_s[19]=vr;
	vf_s[20]=vz;
	vf_s[21]=v0x;
	vf_s[22]=v0y;
	vf_s[23]=v0z;

	single.Fill(vf_s);
	//	single.Fill(cb,p1_px,p1_py,p1_pz,px,py,pz,E,ptjet,charge,particle_ID,Q_flag,prt_id,vr,vz);
	
	//lepton/kaon information
	h1.Fill(ptot);
	h2.Fill(pt);
	h3.Fill(particle_ID);
	h4.Fill(abs(parentId)); // immed parent, may not be the same as the primary parent
	h5.Fill(cosx);
	h6.Fill(Q_flag);
	//primary parent information
	h21.Fill(p1_ptot);
	h22.Fill(p1_pt);
	h23.Fill(p1_id);
	h24.Fill(p1_pz);

 	n_good++; // number of good ele/muons/kaon in this event
      }// if muon,ele,kaon

    } // end of loop event j


    float mass_xx = 0;
    float pt12 = 0;
    float px1,py1,pz1,E1,charge1,m1,flag1;
    float px2,py2,pz2,E2,charge2,m2,flag2;

    int pid1,pid2,prt_id1,prt_id2;

    //two particle correlations
    if (n_good>0){

      //      cout << "n_good = " << n_good << endl;

      for (int ix=0;ix<=n_good;ix++){

	px1 =  varx[ix][4]; 
	py1 =  varx[ix][5]; 
	pz1 =  varx[ix][6]; 
	E1 =   varx[ix][7]; 
	charge1 = varx[ix][9];
	pid1 = (int)varx[ix][10];
	flag1 = varx[ix][11];
	prt_id1 = (int)varx[ix][12];
	m1 =   varx[ix][15];

	for (int jx=ix+1;jx<n_good;jx++){

	  px2 =  varx[jx][4]; 
	  py2 =  varx[jx][5]; 
	  pz2 =  varx[jx][6]; 
	  E2  =   varx[jx][7]; 
	  charge2 = varx[jx][9];
	  pid2 = (int)varx[jx][10];
	  flag2 = varx[jx][11];
	  prt_id2 = (int)varx[jx][12];
	  m2 =   varx[jx][15];


	  mass_xx=sqrt(m1**2 + m2**2 + 2*E1*E2 - 2*(px1*px2 + py1*py2 + pz1*pz2));
	  pt12 = sqrt((px1+px2)**2 + (py1+py2)**2);


	  vf_ll[0]=cb;
	  vf_ll[1]=px1;
	  vf_ll[2]=py1;
	  vf_ll[3]=pz1;
	  vf_ll[4]=flag1;
	  vf_ll[5]=charge1;
	  vf_ll[6]=pid1;
	  vf_ll[7]=prt_id1;

	  vf_ll[8]=px2;
	  vf_ll[9]=py2;
	  vf_ll[10]=pz2;
	  vf_ll[11]=flag2;
	  vf_ll[12]=charge2;
	  vf_ll[13]=pid2;
	  vf_ll[14]=prt_id2;

	  vf_ll[15]=mass_xx;
	  vf_ll[16]=pt12;
	  vf_ll[17]=v0x;
	  vf_ll[18]=v0y;
	  vf_ll[19]=v0z;




	  if (pid1==2&&pid2==2){ // mu-mu
	    mumu.Fill(vf_ll);
	    // mumu.Fill(cb,px1,py1,pz1,flag1,charge1,pid1,px2,py2,pz2,flag2,charge2,pid2,mass_xx,pt12);	  
	  }
	  else if (pid1==2&&pid2==1||pid1==1&&pid2==2){ // mu-e
	    mue.Fill(vf_ll);	    
	    //  mue.Fill(cb,px1,py1,pz1,flag1,charge1,pid1,px2,py2,pz2,flag2,charge2,pid2,mass_xx,pt12);	  
	  }
	  else if (pid1==2&&pid2==3||pid1==3&&pid2==2){ // mu-k
	     muk.Fill(vf_ll);
	     //  muk.Fill(cb,px1,py1,pz1,flag1,charge1,pid1,px2,py2,pz2,flag2,charge2,pid2,mass_xx,pt12);	  
	  }
	  else if (pid1==1&&pid2==1||pid1==1&&pid2==1){ // e-e
	    ee.Fill(vf_ll);	    
	    // ee.Fill(cb,px1,py1,pz1,flag1,charge1,pid1,px2,py2,pz2,flag2,charge2,pid2,mass_xx,pt12);	  
	  }
	  else if (pid1==1&&pid2==3||pid1==3&&pid2==1){ //e-k
	     ek.Fill(vf_ll);
	     // ek.Fill(cb,px1,py1,pz1,flag1,charge1,pid1,px2,py2,pz2,flag2,charge2,pid2,mass_xx,pt12);	  
	  }
	  else if (pid1==3&&pid2==3||pid1==3&&pid2==3){ //k-k
	    kk.Fill(vf_ll);
	    // kk.Fill(cb,px1,py1,pz1,flag1,charge1,pid1,px2,py2,pz2,flag2,charge2,pid2,mass_xx,pt12);	  
	  }

	  else { cout << "ERROR in two particel correlation!!!! \t pid1 =" << pid1 << "pid2 =" << pid2 << endl;}

	}//loop 2nd particle
      }//loop 1st paricles


    } // n_good > 0


    
    //////////////////////////////////////////////////////
    // event global porperties 
    // do more detailed analysis if good particle(s) found
    //////////////////////////////////////////////////////

    if(n_good>0){ // if good particle found

      h7.Fill(n_good); // number of muons/ele/k in this event

      // this is for J/Psi parent information: gg.J/Psi; qq.J/Psi
      for (int j=1; j<=20;j++){ // J starts from 1	
	if((pythia.GetK(j,2)==445)&&(pythia.GetK(j,3)>0)){  // chi_2c ?
	  
	  px =  pythia.GetP(j,1); //(Px, Py,Pz,E,M)
	  py =  pythia.GetP(j,2); //(Px, Py,Pz,E,M)
	  pz =  pythia.GetP(j,3); //(Px, Py,Pz,E,M)
	  E  =  pythia.GetP(j,4); //(Px, Py,Pz,E,M)
	  mass  =  pythia.GetP(j,5); //(Px, Py,Pz,E,M)
	  
	  pt=sqrt(px*px+py*py);
	  ptot = sqrt(px*px+py*py+pz*pz);
	  break;
	}
      }

      x1 =  pythia.GetP(3,3)/100; //(Px, Py,Pz,E,M)
      x2 =  pythia.GetP(4,3)/100; //(Px, Py,Pz,E,M)
      pid1=pythia.GetK(3,2);
      pid2=pythia.GetK(4,2);


      vf_g[0]=x1;
      vf_g[1]=x2;
      vf_g[2]=pid1;
      vf_g[3]=pid2;
      vf_g[4]=px;
      vf_g[5]=py;
      vf_g[6]=pz;
      vf_g[7]=ptot;
      vf_g[8]=pt;
      vf_g[9]=mass;
      vf_g[10]=vx;
      vf_g[11]=vy;
      vf_g[12]=vz;
   

      evt_g.Fill(vf_g);

      //evt_g.Fill(x1,x2,pid1,pid2,px,py,pz,ptot,pt,mass,vx,vy,vz);

      //this is for QQ event analysis
      for (int j=1; j<=20;j++){ // J starts from 1	
	if((abs(pythia.GetK(j,2))==4||abs(pythia.GetK(j,2))==5)&&(pythia.GetK(j,3)==0)){  // cc or bb events
	  
	  //parton-1
	  q_px1 =  pythia.GetP(j,1); //(Px, Py,Pz,E,M)
	  q_py1 =  pythia.GetP(j,2); //(Px, Py,Pz,E,M)
	  q_pz1 =  pythia.GetP(j,3); //(Px, Py,Pz,E,M)
	  q_m1  =  pythia.GetP(j,5); //(Px, Py,Pz,E,M)
	  qid1  =  pythia.GetK(j,2); 
	  //parton-2
	  q_px2 =  pythia.GetP(j+1,1); //(Px, Py,Pz,E,M)
	  q_py2 =  pythia.GetP(j+1,2); //(Px, Py,Pz,E,M)
	  q_pz2 =  pythia.GetP(j+1,3); //(Px, Py,Pz,E,M)
	  q_m2  =  pythia.GetP(j+1,5); //(Px, Py,Pz,E,M)
	  qid2  =  pythia.GetK(j+1,2);
 

	  vf_q[0]=x1;
	  vf_q[1]=x2;
	  vf_q[2]=pid1;
	  vf_q[3]=pid2;
	  vf_q[4]=qid1;
	  vf_q[5]=qid2;

	  vf_q[6]=q_px1;
	  vf_q[7]=q_py1;
	  vf_q[8]=q_pz1;
	  vf_q[9]=q_m1;

	  vf_q[10]=q_px2;
	  vf_q[11]=q_py2;
	  vf_q[12]=q_pz2;
	  vf_q[13]=q_m2;

	  vf_q[14]=vx;
	  vf_q[15]=vy;
	  vf_q[16]=vz;


	  evt_q.Fill(vf_q); 
	  //	  evt_q.Fill(x1,x2,pid1,pid2,qid1,qid2,q_px1,q_py1,q_pz1,q_m1,q_px2,q_py2,q_pz2,q_m2,vx,vy,vz); 
	  break;
	}
      }
    }  // end of global analysis, n_good>0
 
  } // end of event loop  
    
  cout << "------ PYSTAT (1) --------- " << endl;
  pythia.Pystat(1);

  cout << "------ PYSTAT (3) --------- " << endl;
  pythia.Pystat(3);
  
  cout << "------ PYSTAT (4) --------- " << endl;
  pythia.Pystat(4);

  cout << "------ PYSTAT (5) --------- " << endl;
  pythia.Pystat(5);

  cout << "------ PYSTAT (6) --------- " << endl;
  pythia.Pystat(6);



  // draw plots
  c1.Divide(2,3);

  c1.cd(1);
  h1.Draw();
  c1.cd(2);
  h2.Draw();
  c1.cd(3);
  h3.Draw();
  c1.cd(4);
  h4.Draw();
  c1.cd(5);
  h5.Draw();
  c1.cd(6);
  h6.Draw();

  //parent
  c2.Divide(2,2);
  c2.cd(1);
  h21.Draw();
  c2.cd(2);
  h22.Draw();
  c2.cd(3);
  h23.Draw();
  c2.cd(4);
  h24.Draw();

  c3.cd();
  //  c3.SetLogy();
  h7.Draw();

  hfile.Write();
}