#include<cstdio>
#include<TTree.h>
#include<TSystem.h>
#include<TLorentzVector.h>
#include<TStopwatch.h>
#include<TRandom.h>
#include<TGenPhaseSpace.h>
#include<TApplication.h>
#include<TFile.h>

void nphi_hr()
{

  char tname[4][10]={"c12","fe56","cu63","au197"};
  char mname[5][10]={"d","nb","af","l","b"};
  char mname_full[5][20]={"Direct Production",
			  "No Bound State",
			  "a0/f0 Production",
			  "Lambda Production",
			  "Bound State"};
  //  gSystem->Load("libPhysics");
  
  // Phi-N bound state phase space simulation written by Yi Qiang, May 2008
  //....//
  
  // Defining constants and parameters
  
  // NOTE // The following convention will be used for particles: x_y_z_vw
  // where x is the quantity: ms = mass, mo = momentum, en = energy; 
  // where y is the particle being considered: g = photon from
  // beam, p1 = target proton (in nucleus), f = phi meson, 
  // k1 = K+, and k2 = K-, a = target nucleus, a1,a2 = target child nucleus;
  // where z is the frame being considered: lab or cm.
  // where v indicates whether the quantity is before (bf) or after (af) interaction
  
  // Input offset
  int input_offset = 3;
  // Choose Target
  int target = atoi(gApplication->Argv(input_offset));
  if(target<1||target>4) return;
  
  // Choose Mode
  int mode = atoi(gApplication->Argv(1+input_offset));
  if(mode<1||mode>5) return;

  // 1: direction production: g + p1 -> p1' + k1 + k2
  // 2: no bound state:       g + p1 -> p1' + f -> p1' + k1 + k2
  // 3: a0/f0 production:     g + p1 -> p1' + a/f -> p1' + k1 + k2
  // 4: lambda production:    g + p1 -> l + k1 -> p1' + k2 + k1
  // 5: bound state:          g + p1 + p2 -> p1' + f + p2 
  //                                      -> p1' + n-f ->  p1' + p2' + k1 + k2

  printf("Start simulation of %s on %s target...\n",mname_full[mode-1], tname[target-1]);
  
  float PI = TMath::Pi();		// Pi

  // basic particle mss
  float ms_p = 0.938272; 		// Mass of proton in GeV
  //  float ms_n = 0.939566;		// Mass of neutron in GeV
  float ms_k = 0.493677;		// Mass of K meson in GeV
  // phi
  float ms_f_mean = 1.019413;		// Mean mass of phi meson for breit wigner dist
  float ms_f_width = 0.00443;		// Gamma for breit wigner dist of phi particle
  // lambda
  float ms_l_mean  = 1.5195;            // Mean mass of Lambda 1520
  float ms_l_width = 0.0156;            // Gamma for breit wigner dist of lambda particle
  // a0/f0 
  double ms_af_mean = 0.980;            // Mass of a0/f0
  double ms_af_width = 0.1;             // width of a0/f0

  // Proton beam energy, Kinatic energy threshold on Hydrogen is 2.593 GeV, E = 3.531, P = 3.404
  float t_proton_min = atof(gApplication->Argv(2+input_offset));
  float t_proton_max = atof(gApplication->Argv(3+input_offset));

  float eng_b = 0.0025;                 // Binding energy for Nphi bound state
  float ms_nf_mean = ms_f_mean + ms_p - eng_b;          // Mean mass of Phi-N
  float ms_nf_width = atof(gApplication->Argv(4+input_offset));      // Width of phi-N bound state
  printf("Assuming Phi-N Bound State width: %g\n",ms_nf_width);

  float ms_t_all[4];			// Mass of target nuclei in GeV
  ms_t_all[0] = 11.1779291;             // Carbon 12
  ms_t_all[1] = 52.1030653;             // Iron 56
  ms_t_all[2] = 58.618548;              // Copper 63
  ms_t_all[3] = 183.473194;             // Gold 197
  float ms_t = ms_t_all[target-1];

  // Defining variables
  TStopwatch tictoc;			// Object for tracking time
  gRandom->SetSeed(0);  

  // Input the number of trials
  int num_trials = atoi(gApplication->Argv(5+input_offset));
  
  printf("Total number of trails: %d and T range is %g to %g\n",num_trials, t_proton_min, t_proton_max);

  // Flag of saving a event
  int flag_save_event = 0;

  // Stopwatch to track how long script takes
  tictoc.Start();	

  //----------------------------------------------------------------------------  

  // Defining Variables
  float  ms_f;             // Phi mass
  float  ms_af;            // a0/f0 mass
  float  ms_l;             // lambda mass

  float  mo_g_lab_bf,   en_g_lab_bf;    // photo beam energy
  float  mo_p1_lab_bf,  th_p1_lab_bf,  ph_p1_lab_bf,  en_p1_lab_bf;    // incoming proton 1

  float  mo_f_lab_af,   th_f_lab_af,   ph_f_lab_af,   en_f_lab_af;     // generated phi
  float  mo_af_lab_af,  th_af_lab_af,  ph_af_lab_af,  en_af_lab_af;    // generated a/f
  float  mo_l_lab_af,   th_l_lab_af,   ph_l_lab_af,   en_l_lab_af;     // generated lambda

  float  mo_p1_lab_af,  th_p1_lab_af,  ph_p1_lab_af,  en_p1_lab_af;    // scattered proton 1
  float  mo_a1_lab_af,  th_a1_lab_af,  ph_a1_lab_af,  en_a1_lab_af;    // scattered A-1 system

  float  mo_f_a1_bf,    th_f_a1_bf,    ph_f_a1_bf,    en_f_a1_bf;      // incoming phi in A-1 system
  float  mo_p2_a1_bf,   th_p2_a1_bf,   ph_p2_a1_bf,   en_p2_a1_bf;     // incoming proton 2 in A-1 system
  float  mo_f_p2_bf,    th_f_p2_bf,    ph_f_p2_bf,    en_f_p2_bf;      // incoming phi in p2 system
  float  mo_p2_lab_bf,  th_p2_lab_bf,  ph_p2_lab_bf,  en_p2_lab_bf;    // incoming proton 2 in lab system
  float  mo_nf_lab_af,  th_nf_lab_af,  ph_nf_lab_af,  en_nf_lab_af;    // ph-N bound state
  float  mo_p2_lab_af,  th_p2_lab_af,  ph_p2_lab_af,  en_p2_lab_af;    // scattered proton 2

  float  mo_k1_lab_af,  th_k1_lab_af,  ph_k1_lab_af,  en_k1_lab_af;    // scattered K 1
  float  mo_k2_lab_af,  th_k2_lab_af,  ph_k2_lab_af,  en_k2_lab_af;    // scattered K 2

  float  inv_mass_1;      // invariant mass of K + K + p1'
  float  inv_mass_2;      // invariant mass of K + K + p2'

  float  weight_nf;       // weight of phi-N state
  float  weight_pf;       // weight of X --> p + f
  float  weight_paf;      // weight of X --> p + af
  float  weight_pkk;      // weight of X --> K + K + p
  float  weight_kk;       // weight of X --> K + K
  float  weight_lk;       // weight of X --> l + k
  float  weight_pk;       // weight of X --> p + k
  
  // Defining TTree
  TTree *T = new TTree("T","Tree for Data Storage");

  T->Branch("g.E",&en_g_lab_bf,"g.E/F");
  T->Branch("g.P",&mo_g_lab_bf,"g.P/F");

  T->Branch("p1.bf.P",&mo_p1_lab_bf,"p1.bf.P/F");
  T->Branch("p1.bf.Th",&th_p1_lab_bf,"p1.bf.Th/F");
  T->Branch("p1.bf.Ph",&ph_p1_lab_bf,"p1.bf.Ph/F");
  T->Branch("p1.bf.E",&en_p1_lab_bf,"p1.bf.E/F");

  T->Branch("p1.af.P",&mo_p1_lab_af,"p1.af.P/F");
  T->Branch("p1.af.Th",&th_p1_lab_af,"p1.af.Th/F");
  T->Branch("p1.af.Ph",&ph_p1_lab_af,"p1.af.Ph/F");
  T->Branch("p1.af.E",&en_p1_lab_af,"p1.af.E/F");

  T->Branch("a1.af.P",&mo_a1_lab_af,"a1.af.P/F");
  T->Branch("a1.af.Th",&th_a1_lab_af,"a1.af.Th/F");
  T->Branch("a1.af.Ph",&ph_a1_lab_af,"a1.af.Ph/F");
  T->Branch("a1.af.E",&en_a1_lab_af,"a1.af.E/F");

  T->Branch("k1.af.P",&mo_k1_lab_af,"k1.af.P/F");
  T->Branch("k1.af.Th",&th_k1_lab_af,"k1.af.Th/F");
  T->Branch("k1.af.Ph",&ph_k1_lab_af,"k1.af.Ph/F");
  T->Branch("k1.af.E",&en_k1_lab_af,"k1.af.E/F");

  T->Branch("k2.af.P",&mo_k2_lab_af,"k2.af.P/F");
  T->Branch("k2.af.Th",&th_k2_lab_af,"k2.af.Th/F");
  T->Branch("k2.af.Ph",&ph_k2_lab_af,"k2.af.Ph/F");
  T->Branch("k2.af.E",&en_k2_lab_af,"k2.af.E/F");

  T->Branch("M1",&inv_mass_1,"M1/F");


  if(mode == 2 || mode == 5)
    {
      T->Branch("f.M",&ms_f,"f.M/F");
 
      T->Branch("f.af.P",&mo_f_lab_af,"f.af.P/F");
      T->Branch("f.af.Th",&th_f_lab_af,"f.af.Th/F");
      T->Branch("f.af.Ph",&ph_f_lab_af,"f.af.Ph/F");
      T->Branch("f.af.E",&en_f_lab_af,"f.af.E/F");

      T->Branch("f.bf.a1.P",&mo_f_a1_bf,"f.bf.a1.P/F");
      T->Branch("f.bf.a1.Th",&th_f_a1_bf,"f.bf.a1.Th/F");
      T->Branch("f.bf.a1.Ph",&ph_f_a1_bf,"f.bf.a1.Ph/F");
      T->Branch("f.bf.a1.E",&en_f_a1_bf,"f.bf.a1.E/F");

      T->Branch("w.pf",&weight_pf,"w.pf/F");
  
      if(mode == 5)
	{
	  T->Branch("p2.bf.P",&mo_p2_lab_bf,"p2.bf.P/F");
	  T->Branch("p2.bf.Th",&th_p2_lab_bf,"p2.bf.Th/F");
	  T->Branch("p2.bf.Ph",&ph_p2_lab_bf,"p2.bf.Ph/F");
	  T->Branch("p2.bf.E",&en_p2_lab_bf,"p2.bf.E/F");

	  T->Branch("p2.bf.a1.P",&mo_p2_a1_bf,"p2.bf.a1.P/F");
	  T->Branch("p2.bf.a1.Th",&th_p2_a1_bf,"p2.bf.a1.Th/F");
	  T->Branch("p2.bf.a1.Ph",&ph_p2_a1_bf,"p2.bf.a1.Ph/F");
	  T->Branch("p2.bf.a1.E",&en_p2_a1_bf,"p2.bf.a1.E/F");

	  T->Branch("f.bf.p2.P",&mo_f_p2_bf,"f.bf.p2.P/F");
	  T->Branch("f.bf.p2.Th",&th_f_p2_bf,"f.bf.p2.Th/F");
	  T->Branch("f.bf.p2.Ph",&ph_f_p2_bf,"f.bf.p2.Ph/F");
	  T->Branch("f.bf.p2.E",&en_f_p2_bf,"f.bf.p2.E/F");

	  T->Branch("nf.af.P",&mo_nf_lab_af,"nf.af.P/F");
	  T->Branch("nf.af.Th",&th_nf_lab_af,"nf.af.Th/F");
	  T->Branch("nf.af.Ph",&ph_nf_lab_af,"nf.af.Ph/F");
	  T->Branch("nf.af.E",&en_nf_lab_af,"nf.af.E/F");

	  T->Branch("p2.af.P",&mo_p2_lab_af,"p2.af.P/F");
	  T->Branch("p2.af.Th",&th_p2_lab_af,"p2.af.Th/F");
	  T->Branch("p2.af.Ph",&ph_p2_lab_af,"p2.af.Ph/F");
	  T->Branch("p2.af.E",&en_p2_lab_af,"p2.af.E/F");	  

	  T->Branch("M2",&inv_mass_2,"M2/F");

	  T->Branch("w.nf",&weight_nf,"w.nf/F");
	}
    }

  else if(mode == 3)
    {
      T->Branch("af.M",&ms_af,"af.M/F");
 
      T->Branch("af.af.P",&mo_af_lab_af,"af.af.P/F");
      T->Branch("af.af.Th",&th_af_lab_af,"af.af.Th/F");
      T->Branch("af.af.Ph",&ph_af_lab_af,"af.af.Ph/F");
      T->Branch("af.af.E",&en_af_lab_af,"af.af.E/F");

      T->Branch("w.paf",&weight_paf,"w.paf/F");      
    }

  else if(mode == 4) 
    {
      T->Branch("l.M",&ms_l,"l.M/F");

      T->Branch("l.af.P",&mo_l_lab_af,"l.af.P/F");
      T->Branch("l.af.Th",&th_l_lab_af,"l.af.Th/F");
      T->Branch("l.af.Ph",&ph_l_lab_af,"l.af.Ph/F");
      T->Branch("l.af.E",&en_l_lab_af,"l.af.E/F");

      T->Branch("w.lk",&weight_lk,"w.lk/F");
      T->Branch("w.pk",&weight_pk,"w.pk/F");
    }

  if(mode == 1 || mode == 5) 
    T->Branch("w.pkk",&weight_pkk,"w.pkk/F");

  else if(mode == 2 || mode == 3)  
    T->Branch("w.kk",&weight_kk,"w.kk/F");
  
  //----------------------------------------------------------------------------

  // Fermi distribution
  float sf_t_p_e[49][100];
  char table_file[256];
  
  sprintf(table_file,"./sf_%s.dat",tname[target-1]);

  printf("Reading data file from %s\n",table_file);
  FILE *ff = fopen(table_file,"r");
  if(!ff)
    {
      printf("Cannot find data file: %s, quit...\n",table_file);
      return;
    }

  for(int i=0; i<49; i++)
    for(int j=0; j<100; j++)
      {
	sf_t_p_e[i][j]=0.;
      }

  int nn = 0;
  char line[256];
  float tmp_prob[4900][3];
  float tmp_e;
  float tmp_p;
  float tmpx;
  float sf_prob_sum = 0;
  
  while (fgets(line,256,ff)&&nn<5000)
    {
      int kk = sscanf(line,"%f%f%f",&tmp_p,&tmp_e,&tmpx);
      if(kk<3) continue;
      else
	{
	  tmp_prob[nn][0]=tmp_p;
	  tmp_prob[nn][1]=tmp_e;
	  tmp_prob[nn][2]=tmpx;
	  sf_prob_sum += tmpx;
	  nn++;
	}
    } 
  printf("Total lines read: %d\n",nn);

  // Copying data to a smaller matrix
  for (int i=0; i < nn; i++)
    {
      int tmp_p_indx = int((tmp_prob[i][0]-5.0)/10.);
      if(tmp_p_indx<0||tmp_p_indx>=49) continue;
      int tmp_e_indx = int((tmp_prob[i][1]/2.));
      if(tmp_e_indx<0||tmp_e_indx>=100) continue;
      sf_t_p_e[tmp_p_indx][tmp_e_indx]=tmp_prob[i][2]/sf_prob_sum;
    }

  // Interpolate into a larger matrix
  float zz[3][3];
  float sf_p_e[245][500];   // 5x5 times
  int mag_fac = 5;
  int x0,x1,y0,y1;
  float dx, dy;
  float tmp_zz;
  for (int i=0; i<49; i++)
    for(int j=0; j<100; j++)
      {
	// Read the adjecent grid points: 3x3
	for(int l=-1; l<2; l++)
	  for(int m=-1; m<2; m++)
	    {
	      if(i+l<0||i+l>48||j+m<0||j+m>99) zz[l+1][m+1]=0;
	      else zz[l+1][m+1]=sf_t_p_e[i+l][j+m];
	    }

	for(int l=0; l<mag_fac; l++)
	  for(int m=0; m<mag_fac; m++)
	    {
	      x0 = y0 = 1;
	      if(l-(mag_fac-1.)/2.<-0.25) x1 = 0;
	      else if(l-(mag_fac-1.)/2.>0.25) x1 = 2;
	      else x1 = 1;
	      if(m-(mag_fac-1.)/2.<-0.25) y1 = 0;
	      else if(m-(mag_fac-1.)/2.>0.25) y1 = 2;
	      else y1 = 1;

	      dx = fabs(l-(mag_fac-1.)/2.0);
	      dy = fabs(m-(mag_fac-1.)/2.0);

	      if(x1==1&&y1==1) tmp_zz=zz[1][1];
	      else if(x1==1) tmp_zz=((mag_fac-dy)*zz[1][y0] +
				     dy*zz[1][y1])/mag_fac;
	      else if(y1==1) tmp_zz=((mag_fac-dx)*zz[x0][1] +
				     dx*zz[x1][1])/mag_fac;
  	      else tmp_zz = ((mag_fac-dx)*(mag_fac-dy)*zz[x0][y0] +
 			     (mag_fac-dx)*dy*zz[x0][y1] +
 			     dx*(mag_fac-dy)*zz[x1][y0] +
 			     dx*dy*zz[x1][y1])/pow(mag_fac,2);
	      sf_p_e[i*mag_fac+l][j*mag_fac+m] = tmp_zz/mag_fac/mag_fac;
	    }
      }

  // Making projection Prob(E)

  // Defining some arrys
  float sf_t_proj_ecum[501];   // E miss Accumulated probability
  float sf_t_proj_pcum[500][246];  // P fermi Accumulated probability

  sf_t_proj_ecum[0]=0.;
  for (int i=0; i<100*mag_fac; i++)
    {
      sf_t_proj_pcum[i][0]=0.;
      for (int j=0; j<49*mag_fac; j++)
	{
	  sf_t_proj_pcum[i][j+1] = sf_t_proj_pcum[i][j]+sf_p_e[j][i];
	}
      sf_t_proj_ecum[i+1] = sf_t_proj_ecum[i]+sf_t_proj_pcum[i][49*mag_fac];
    }
  printf("Total probablity is %g\n",sf_t_proj_ecum[100*mag_fac]);

  //----------------------------------------------------------------------------

  // Defining 4-Vectors
  // Primary Frames of Reference
  TLorentzVector lv_g_lab_bf, lv_p1_lab_bf, lv_p1_lab_af;
  TLorentzVector lv_f_lab_af, lv_af_lab_af, lv_l_lab_af;
  TLorentzVector lv_p2_lab_bf, lv_p2_a1_bf, lv_f_a1_bf, lv_f_p2_bf, lv_nf_lab_af;
  TLorentzVector lv_p2_lab_af, lv_k1_lab_af, lv_k2_lab_af;

  // Target Nucleus
  TLorentzVector lv_t_lab_bf, lv_a1_lab_af;
  TLorentzVector lv_tot;

  // Phasespace variables
  TGenPhaseSpace event1, event2;

  // missing mass and missing momentum
  float en_miss_p[2];
  float mo_p[2];
	
  // Attributting 4-Vector for target nuclei
  lv_t_lab_bf.SetXYZM(0,0,0,ms_t);

  ////////////////////////////////////////////
  //                                        //
  // Major loop of events generating        //
  //                                        //
  ////////////////////////////////////////////

  printf("Start simulation...\n");
  // Beginning of major loop where the kinematics are carried out
  for (int i=0; i < num_trials; i++)
    {

      // clear saving event flag
      flag_save_event = 0;

      if((i+1)%500000==0) printf("%d events generated...\n",i+1);
  
      // Generate random missing mass and fermi momentum;
      float tmp_rem=0.;
      float a,b;
      mo_p[0]=mo_p[1]=en_miss_p[0]=en_miss_p[1]=0;

      float n_j = 1;
      if(mode==5) n_j = 2;

      for(int j = 0; j<n_j; j++)
	{
	  a = gRandom->Rndm();
	  for (int k=0; k<100*mag_fac; k++) 
	    { 
	      if ( a <= sf_t_proj_ecum[k+1] )
		{
		  tmp_rem=(a-sf_t_proj_ecum[k])/(sf_t_proj_ecum[k+1]-sf_t_proj_ecum[k]);
		  en_miss_p[j] = 2.0/mag_fac*(k+tmp_rem)/1000.;  // Missing energy in GeV
		  
		  b = gRandom->Uniform(sf_t_proj_pcum[k][49*mag_fac]);
		  for (int l=0; l<49*mag_fac; l++) 
		    {
		      if (b<=sf_t_proj_pcum[k][l+1])
			{
			  tmp_rem = (b-sf_t_proj_pcum[k][l])/(sf_t_proj_pcum[k][l+1]-sf_t_proj_pcum[k][l]);
			  mo_p[j]=(5.+10./mag_fac*(l+tmp_rem))/1000.;  // Fermi momentum in GeV
			  break;
			}
		    }
		  break;
		}
	    }
	}

      // Proton beam energy
      en_g_lab_bf = gRandom->Uniform(t_proton_min,t_proton_max);
      lv_g_lab_bf.SetXYZM(0,0,en_g_lab_bf,0);

      // Clear all variables
      inv_mass_1 = 0;
      inv_mass_2 = 0;
	
      
      weight_nf = 0;
      weight_pf = 0;
      weight_paf = 0;
      weight_kk = 0;
      weight_pkk = 0;
      weight_lk = 0;
      weight_pk = 0;

      ms_f = 0;
      ms_af = 0;
      ms_l = 0;

      lv_p1_lab_af.SetXYZT(0,0,0,0);
      lv_k1_lab_af.SetXYZT(0,0,0,0);
      lv_k2_lab_af.SetXYZT(0,0,0,0);

      lv_p2_lab_bf.SetXYZT(0,0,0,0);
      lv_p2_a1_bf.SetXYZT(0,0,0,0);
      lv_p2_lab_af.SetXYZT(0,0,0,0);

      lv_f_lab_af.SetXYZT(0,0,0,0);
      lv_f_a1_bf.SetXYZT(0,0,0,0);
      lv_f_p2_bf.SetXYZT(0,0,0,0);

      lv_af_lab_af.SetXYZT(0,0,0,0);
      lv_l_lab_af.SetXYZT(0,0,0,0);

      lv_nf_lab_af.SetXYZT(0,0,0,0);

      
      // 4-Vector for residue nucleus: A-1
      lv_a1_lab_af.SetXYZM(0,0,mo_p[0],ms_t+en_miss_p[0]-ms_p);
      lv_a1_lab_af.SetTheta(acos(gRandom->Uniform(-1.0,1.0)));
      lv_a1_lab_af.SetPhi(gRandom->Uniform(-PI,PI));
      
      // 4-Vector for Proton 1 in the nucleus
      lv_p1_lab_bf = lv_t_lab_bf-lv_a1_lab_af;
      lv_tot = lv_p1_lab_bf+lv_g_lab_bf;

      // Direct decay
      if(mode == 1)
	{
	  // g + p1 --> p1' + k1 + k2
	  if(lv_tot.M()>=ms_p+2*ms_k)
	    {
	      double masses1[3] = {ms_p, ms_k, ms_k};
	      event1.SetDecay(lv_tot, 3, masses1);
	      weight_pkk = event1.Generate();
	      lv_p1_lab_af = *(event1.GetDecay(0));
	      lv_k1_lab_af = *(event1.GetDecay(1));
	      lv_k2_lab_af = *(event1.GetDecay(2));
	      
	      flag_save_event = 1;
	    }
	}

      // a0/f0 production
      else if(mode == 3)
	{
	  // Generate random a0/f0 mass
	  while  (!ms_af)
	    {
	      float ms_af_temp = gRandom->BreitWigner(ms_af_mean,ms_af_width);
	      if (ms_af_temp>0) ms_af=ms_af_temp;
	    }

	  // g + p1 -> p1' + af	  
	  if(lv_tot.M()>=ms_p+ms_af)  
	    {
	      double masses1[2] = {ms_p, ms_af};
	      event1.SetDecay(lv_tot, 2, masses1);
	      weight_paf = event1.Generate();
	      lv_p1_lab_af = *(event1.GetDecay(0));
	      lv_af_lab_af = *(event1.GetDecay(1));
	      
	      flag_save_event = 1;

	      // af -> k1 + k2	  
	      if(ms_af>=2*ms_k)
		{
		  double masses2[2] = {ms_k,ms_k};
		  event2.SetDecay(lv_af_lab_af, 2, masses2);
		  weight_kk = event2.Generate();
		  lv_k1_lab_af = *(event2.GetDecay(0));
		  lv_k2_lab_af = *(event2.GetDecay(1));
		}
	    }	
	}  
	  
      // lambda production
      else if(mode == 4)
	{
	  // Generate random lambda mass
	  while  (!ms_l)
	    {
	      float ms_l_temp = gRandom->BreitWigner(ms_l_mean,ms_l_width);
	      if (ms_l_temp>0) ms_l=ms_l_temp;
	    }

	  // g + p1 -> l + k1	  
	  if(lv_tot.M()>=ms_l+ms_k)  
	    {
	      double masses1[2] = {ms_l, ms_k};
	      event1.SetDecay(lv_tot, 2, masses1);
	      weight_lk = event1.Generate();
	      lv_l_lab_af = *(event1.GetDecay(0));
	      lv_k1_lab_af = *(event1.GetDecay(1));
	      
	      flag_save_event = 1;

	      // l -> p1' + k2	  
	      if(ms_l>=ms_p+ms_k)
		{
		  double masses2[2] = {ms_p,ms_k};
		  event2.SetDecay(lv_l_lab_af, 2, masses2);
		  weight_pk = event2.Generate();
		  lv_p1_lab_af = *(event2.GetDecay(0));
		  lv_k2_lab_af = *(event2.GetDecay(1));
		}
	    }	  
	}
      
      // No or Phi-N bound state
      else if(mode == 2 || mode == 5)
	{

	  // Generate random phi mass
	  while  (ms_f <= 0)
	    {
	      ms_f = gRandom->BreitWigner(ms_f_mean,ms_f_width);
	    }

	  // g + p1 -> p1' + f
	  if(lv_tot.M()>=ms_p+ms_f)  
	    {
	      double masses1[2] = {ms_p, ms_f};
	      event1.SetDecay(lv_tot, 2, masses1);
	      weight_pf = event1.Generate();
	      lv_p1_lab_af = *(event1.GetDecay(0));
	      lv_f_lab_af = *(event1.GetDecay(1));
	      
	      lv_f_a1_bf = lv_f_lab_af;
	      lv_f_a1_bf.Boost(-lv_a1_lab_af.BoostVector());
	      

	      // No bound state
	      if(mode == 2)
		{
		  flag_save_event = 1;
		  
		  // f -> k1 + k2	  
		  if(ms_f>=2*ms_k)
		    {
		      double masses2[2] = {ms_k,ms_k};
		      event2.SetDecay(lv_f_lab_af, 2, masses2);
		      weight_kk = event2.Generate();
		      lv_k1_lab_af = *(event2.GetDecay(0));
		      lv_k2_lab_af = *(event2.GetDecay(1));
		    }
		}

	      // Phi-N bound state
	      else if(mode == 5)
		{
		  
		  float ms_nf = 0;
		  // mass of bound state
		  while (ms_nf <= 0)
		    {
		      ms_nf = gRandom->BreitWigner(ms_nf_mean,ms_nf_width);
		    }
		  
		  // decay of nf -> p2' + K + K	      
		  if(ms_nf>=ms_p+2*ms_k)
		    {

		      // f energy in (A-1) rest system
		      float en_f = lv_f_a1_bf.E();
		      float mo_f = lv_f_a1_bf.P();
		      float theta_f = lv_f_a1_bf.Theta();
		      float phi_f = lv_f_a1_bf.Phi();
		      
		      // properties of off-shell proton p2 in A-1 rest system
		      float ms_a1 = lv_a1_lab_af.M();    // Mass for A-1 system
		      float en_miss = en_miss_p[1];      // energy loss for p2 in GeV
		      float mo_p2 = mo_p[1];             // Fermi momentum of p2
		      float ms_a2 = ms_a1+en_miss-ms_p;  // Mass for A-2 nuclear residue
		      float en_p2 = ms_a1-sqrt(mo_p2*mo_p2+ms_a2*ms_a2);  // energy of p2
		      if(en_p2<0) en_p2 = 0;
		      float ms_p2 = 0;
		      if(en_p2>mo_p2) ms_p2 = sqrt(en_p2*en_p2-mo_p2*mo_p2);
		      
		      // cos(theta) of p2 relative to f
		      float cos_theta = (ms_f*ms_f+ms_p2*ms_p2+2*en_f*en_p2-ms_nf*ms_nf)/2/mo_f/mo_p2;
		      // if the angle is real (0<theta<180) and nf bound state is generated
		      if(cos_theta>=-1&&cos_theta<=1)
			{
			  //			  weight_nf=ms_nf/mo_f/mo_p2;
			  weight_nf=sqrt(1-pow(cos_theta,2));
			  
			  lv_p2_a1_bf.SetXYZM(0,0,mo_p2,ms_p2);
			  lv_p2_a1_bf.SetTheta(acos(cos_theta));
			  lv_p2_a1_bf.SetPhi(gRandom->Uniform(-PI,PI));
			  lv_p2_a1_bf.RotateY(theta_f);
			  lv_p2_a1_bf.RotateZ(phi_f);
			  lv_p2_lab_bf=lv_p2_a1_bf;
			  lv_p2_lab_bf.Boost(lv_a1_lab_af.BoostVector());
			  
			  lv_f_p2_bf = lv_f_lab_af;
			  lv_f_p2_bf.Boost(-lv_p2_lab_bf.BoostVector());
			  
			  lv_nf_lab_af = lv_p2_lab_bf+lv_f_lab_af;

			  flag_save_event = 1;
			  
			  // For nf -> p+k+K
			  if(ms_nf>=ms_p+ms_k+ms_k)
			    {
			      double masses2[3] = {ms_p, ms_k, ms_k};
			      event2.SetDecay(lv_nf_lab_af, 3, masses2);
			      weight_pkk = event2.Generate();
			      lv_p2_lab_af = *(event2.GetDecay(0));
			      lv_k1_lab_af = *(event2.GetDecay(1));
			      lv_k2_lab_af = *(event2.GetDecay(2));
			    }
			}
		    }
		}
	    }
	}
	  
	
      // Storing Data
      if(flag_save_event)
	{
	  mo_p1_lab_bf = lv_p1_lab_bf.P();
	  th_p1_lab_bf = lv_p1_lab_bf.Theta();
	  ph_p1_lab_bf = lv_p1_lab_bf.Phi();
	  en_p1_lab_bf = lv_p1_lab_bf.E();

	  mo_p1_lab_af = lv_p1_lab_af.P();
	  th_p1_lab_af = lv_p1_lab_af.Theta();
	  ph_p1_lab_af = lv_p1_lab_af.Phi();
	  en_p1_lab_af = lv_p1_lab_af.E();
	  
	  mo_f_lab_af = lv_f_lab_af.P();
	  th_f_lab_af = lv_f_lab_af.Theta();
	  ph_f_lab_af = lv_f_lab_af.Phi();
	  en_f_lab_af = lv_f_lab_af.E();
	  
	  mo_f_a1_bf = lv_f_a1_bf.P();
	  th_f_a1_bf = lv_f_a1_bf.Theta();
	  ph_f_a1_bf = lv_f_a1_bf.Phi();
	  en_f_a1_bf = lv_f_a1_bf.E();
	  
	  mo_f_p2_bf = lv_f_p2_bf.P();
	  th_f_p2_bf = lv_f_p2_bf.Theta();
	  ph_f_p2_bf = lv_f_p2_bf.Phi();
	  en_f_p2_bf = lv_f_p2_bf.E();

	  mo_af_lab_af = lv_af_lab_af.P();
	  th_af_lab_af = lv_af_lab_af.Theta();
	  ph_af_lab_af = lv_af_lab_af.Phi();
	  en_af_lab_af = lv_af_lab_af.E();

	  mo_l_lab_af = lv_l_lab_af.P();
	  th_l_lab_af = lv_l_lab_af.Theta();
	  ph_l_lab_af = lv_l_lab_af.Phi();
	  en_l_lab_af = lv_l_lab_af.E();
	  
     	  mo_p2_lab_bf = lv_p2_lab_bf.P();
	  th_p2_lab_bf = lv_p2_lab_bf.Theta();
	  ph_p2_lab_bf = lv_p2_lab_bf.Phi();
	  en_p2_lab_bf = lv_p2_lab_bf.E();
	  
	  mo_p2_a1_bf = lv_p2_a1_bf.P();
	  th_p2_a1_bf = lv_p2_a1_bf.Theta();
	  ph_p2_a1_bf = lv_p2_a1_bf.Phi();
	  en_p2_a1_bf = lv_p2_a1_bf.E();

	  mo_p2_lab_af = lv_p2_lab_af.P();
	  th_p2_lab_af = lv_p2_lab_af.Theta();
	  ph_p2_lab_af = lv_p2_lab_af.Phi();
	  en_p2_lab_af = lv_p2_lab_af.E();
	  	  
	  mo_nf_lab_af = lv_nf_lab_af.P();
	  th_nf_lab_af = lv_nf_lab_af.Theta();
	  ph_nf_lab_af = lv_nf_lab_af.Phi();
	  en_nf_lab_af = lv_nf_lab_af.E();
	  
	  mo_k1_lab_af = lv_k1_lab_af.P();
	  th_k1_lab_af = lv_k1_lab_af.Theta();
	  ph_k1_lab_af = lv_k1_lab_af.Phi();
	  en_k1_lab_af = lv_k1_lab_af.E();
	  
	  mo_k2_lab_af = lv_k2_lab_af.P();
	  th_k2_lab_af = lv_k2_lab_af.Theta();
	  ph_k2_lab_af = lv_k2_lab_af.Phi();
	  en_k2_lab_af = lv_k2_lab_af.E();
	  
	  mo_a1_lab_af = lv_a1_lab_af.P();
	  th_a1_lab_af = lv_a1_lab_af.Theta();
	  ph_a1_lab_af = lv_a1_lab_af.Phi();
	  en_a1_lab_af = lv_a1_lab_af.E();
	  
	  inv_mass_1 = (lv_k1_lab_af+lv_k2_lab_af+lv_p1_lab_af).M();	  
	  inv_mass_2 = (lv_k1_lab_af+lv_k2_lab_af+lv_p2_lab_af).M();
	  
	  T->Fill();
	}
    }
  
  // Defining TFile
  sprintf(table_file,"./plot/tree_%s.root",tname[target-1]);
  TFile *tf_data = new TFile(table_file,"update");
  
  // Saving Data
  sprintf(table_file,"T_%s",mname[mode-1]);
  T->Write(table_file,TObject::kOverwrite);
  tf_data->Close();
  tictoc.Stop();
  printf("Time elapsed: %g\n",tictoc.RealTime());
}