#include <string.h>
#include <iostream>
#include <vector>
#include "TFile.h"
#include "TLorentzVector.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TRandom3.h"
#include "TVector3.h"                                                                           
#include <fstream>

using namespace std;

struct myGen_t {                                                     
   int runNumber;
   int eventNumber;
   double weight;
   vector<int> pidVecGeant;
   vector<int> pidVecPDG;
   vector<int> chargeVec;
   vector<TLorentzVector> mo4Vec;
   vector<TVector3> vertex3Vec;
};


//FUNCTION PROTOTYPES
void printUsage(); // print usage
void writeEventLund(myGen_t myGen,std::ofstream *lundFile);                          
void writeEventGenr8(myGen_t myGen,std::ofstream *genr8File);



void OneBodyDecay(TRandom3 *gRandom,double m_mom, double m_d1, double m_d2, TLorentzVector *D1, TLorentzVector *D2)
{
  double pi = 3.14159;
  double energy_D1 = (pow(m_mom,2)+pow(m_d1,2)-pow(m_d2,2))/(2*m_mom);
  double p_D1      = sqrt(pow(energy_D1,2)-pow(m_d1,2));
  
  double energy_D2 = m_mom - energy_D1; 

  
  //gRandom =  new TRandom3();
  //gRandom->SetSeed(0);

  double phi      = gRandom->Uniform(-pi,pi);    
  double cosTheta = gRandom->Uniform(-1.0,1.0);
  
  double p_D1_ZCM           = p_D1*cosTheta;
  double p_D1_Transverse    = sqrt(pow(p_D1,2) - pow(p_D1_ZCM,2));
  double p_D1_XCM           = p_D1_Transverse*cos(phi);
  double p_D1_YCM           = p_D1_Transverse*sin(phi);

  double p_D2_ZCM           = -p_D1_ZCM;
  double p_D2_XCM           = -p_D1_XCM;
  double p_D2_YCM           = -p_D1_YCM;           

  *D1 =  TLorentzVector(p_D1_XCM,p_D1_YCM,p_D1_ZCM,energy_D1);
  *D2 =  TLorentzVector(p_D2_XCM,p_D2_YCM,p_D2_ZCM,energy_D2);

}



//MAIN
int main(int argc, char **argv) {
  
  char  inFileName[120];
  char  outFileName[120];
  char  outFileNameLund[120];                                                                   
  char  outFileNameGenr8[120];
  char *argptr;
  //  int nToRun = 100000;
  int nToRun = 250000;
  int rSeed = 0;

  //double pi = 3.14159;
  //  double mHyperon      = 2.153;//2.96;
  //double mHyperonCenter= 2.153;//2.96;
  //double gHyperon      = 0.107;//0.8458;
  
  double mHyperon      = 2.96749;// will change it to 
  double mHyperonCenter= 2.96749;
  double gHyperon      = 0.807657;

  double mKplus        = 0.4937; 
  double mXiStarCenter = 1.535;
  double mXiStar       = mXiStarCenter;
  double gXiStar       = 9.9/1000;
  double mpi0          = 0.1349768;
  double mXiGS         = 1.3217; 
  double mPro          = .938;
  double mpim          = 0.1395701;
  double mLambda       = 1.11568;



  TLorentzVector target(0,0,0,mPro);

  TLorentzVector x4Vec(0,0,0,0);



  //Set the default output file
  sprintf(outFileName,"./outFile.root");
  sprintf(outFileNameLund,"./myOutFile.lund");                                                  
  sprintf(outFileNameGenr8,"./myOutFile.genr8");

  int outType  = 1;                                                                            
  int runNumber = 99999;

  //if (argc == 1) printUsage();
  for (int i=1; i<argc; i++) {
    argptr = argv[i];
    if (*argptr == '-') {
      argptr++;
      switch (*argptr) {
      case 'h':
        printUsage();
        break;
      case 'n':
        nToRun = atoi(++argptr);
        break;
      case 's':
        rSeed = atoi(++argptr);
        break;
      case 'O':
	strcpy(outFileName,++argptr);
	break;
      case 'o':
	strcpy(outFileNameLund,++argptr);
	strcpy(outFileNameGenr8,argptr);                                          
	break;
      case 't':                                                                                 
	outType = atoi(++argptr);                                                              
	break;                                                                     
      case 'r':                                                                               
	runNumber = atoi(++argptr);                                                             
	break;     
      default:
        printUsage();
        break;
      }
    } else {
      strcpy(inFileName,argptr);
    }
  }

  //Create histograms here
  TH1D *hDummy1d     = new TH1D("hDummy1d","",100,0.0,1.0);
  TH2D *hDummy2d     = new TH2D("hDummy2d","",100,0.0,1.0,100,0.0,1.0);

  TH1D * hEGamma       = new TH1D("hEGamma","",400,7.0,11.0);  
  TH1D * hTVal            = new TH1D("hTVal","",200,0.0,20);  
  TH2D * hMoVsThetaY      = new TH2D("hMoVsThetaY","",90,0.0,90.0,100,1.5,12.0);  
  TH2D * hMoVsThetaK1     = new TH2D("hMoVsThetaK1","",90,0.0,90.0,100,0.0,10.0);  
  TH2D * hMoVsThetaK2     = new TH2D("hMoVsThetaK2","",90,0.0,90.0,100,0.0,10.0);  
  //TH2D * hMoVsThetaXiStar = new TH2D("hMoVsThetaXiStar","",90,0.0,90.0,100,0.0,10.0); 
  TH1D * hMassHyperon     = new TH1D("hMassHyperon","",2000,2.0,4.0);
  TH1D * hMassHyperon_V2  = new TH1D("hMassHyperon_V2","",2000,2.0,4.0);
  TH1D * hMassK2          = new TH1D("hMassK2","",500,0.0,2.0);   
  TH2D * hMoVsThetaKLow   = new TH2D("hMoVsThetaKLow","",900,0.0,90.0,1000,0.0,10.0);  
  TH2D * hMoVsThetaKHigh  = new TH2D("hMoVsThetaKHigh","",900,0.0,90.0,1000,0.0,10.0);  
  //TH1D * hXiStarRecon     = new TH1D("hXiStarRecon","",500,1.4,2.8);  
  //TH1D * hpXiStarRecon    = new TH1D("hpXiStarRecon","",1000,0.0,1.5);  
  TH1D * hEXiStarRecon        = new TH1D("hEXiStarRecon","",1000,1.4,2.4);  
  TH1D * hXiStarRecon         = new TH1D("hXiStarRecon","",1000,1.4,2.4);  
  TH2D * hXiStarvsCosTheta    = new TH2D("hXiStarvsCosTheta","",500,1.46,1.62,200,-1.0,1.0);  
  TH2D * hXiStarvsCosThetaCM  = new TH2D("hXiStarvsCosThetaCM","",500,1.46,1.62,200,-1.0,1.0);  
  TH2D * hMassXiMinusVsP      = new TH2D("hMassXiMinusVsP","",500,1.0,2.0,800,0.0,8.0);   
  TH2D * hMassPi0VsP          = new TH2D("hMassPi0VsP","",500,0.0,0.5,800,0.0,8.0);   
  TH2D * hXiStarvsP           = new TH2D("hXiStarvsP","",500,1.45,1.63,900,1.0,10.0);  
  TH1D * hXiStar_P            = new TH1D("hXiStar_P","",1000,0.0,10.0);  
  TH2D * hXiStarvsTheta       = new TH2D("hXiStarvsTheta","",500,1.45,1.63,500,-1.0,1.0);  
  TH1D * hYStar               = new TH1D("hYStar","",1000,0.0,5.0);  
  TH1D * k1P                  = new TH1D("k1P","",1000,0.0,5.0);  
  TH1D * hYStarLow            = new TH1D("hYStarLow","",3000,1.5,4.5);  
  TH1D * hYStarHigh           = new TH1D("hYStarHigh","",3000,1.5,4.5);  
  TH2D * hMLambdaVsP          = new TH2D("hLambdaVsP","",300,1.0,1.3,1000,0.0,10.0);  
  TH2D * hProtonMvsP          = new TH2D("hProtonMvsP","",300,0.7,1.0,1000,0.0,10.0);  
  TH2D * hPiMinus1MvsP        = new TH2D("hPiMinus1MvsP","",400,0.0,0.2,1000,0.0,10.0);  
  TH2D * hPiMinus2MvsP        = new TH2D("hPiMinus2MvsP","",400,0.0,0.2,1000,0.0,10.0);  

  TH1D * hDeltaPx                  = new TH1D("hDeltaPx","",101,-1.0,1.0);  
  TH1D * hDeltaPy                  = new TH1D("hDeltaPy","",101,-1.0,1.0);  
  TH1D * hDeltaPz                  = new TH1D("hDeltaPz","",101,-1.0,1.0);  
  TH1D * hDeltaE                   = new TH1D("hDeltaE","",101,-1.0,1.0);  

  //  TH1D * hXpx               = new TH1D("hXpx","",1000,5.0,5.0);  
  //TH1D * hXpy               = new TH1D("hXpy","",1000,5.0,5.0);  
  //TH1D * hXpz               = new TH1D("hXpz","",1000,5.0,5.0);  
  //TH1D * hXe               = new TH1D("hXe","",1000,5.0,5.0);  

  //  double eGammaLow_V1  = 3;
  double eGammaLow_V1  = 7.00;
  double eGammaHigh_V1 = 11.00;

  TH1D * hEGamma_V1    = new TH1D("hEGamma_V1","",500,eGammaLow_V1,eGammaHigh_V1);

  TFile *eFile=new TFile("eDist.root");
  TH1D* hEDist=(TH1D*)eFile->Get("hEDist"); 
  /*  
  TFile *inCoBrem=new TFile("cobrems.root"); //Get the file that has histo to MC       
  TH1D* hGvsE=(TH1D*)inCoBrem->Get("cobrem_vs_E"); //Get the histogram to MC    
  hGvsE->GetXaxis()->SetRangeUser(eGammaLow_V1,eGammaHigh_V1);//Set Range to be used for MC 
  double gMax = hGvsE->GetMaximum();//Find maximum                              
 
  //Setup the output file                                                       
  // TFile *outFile = new TFile(outFileName,"RECREATE");
  */
  //Setup the random generator
  TRandom3 *gRandom =  new TRandom3();
  gRandom->SetSeed(rSeed);
 
  //  double eff;
  double eGamma; 
  /*
  for(Int_t i=0;i<nToRun;i++){
    //    if (i%10000 == 0 && i > 0) cout<<"Events processed = "<<i<<endl;
    //cout<<i<<endl;
    double yMax = gMax*1.02;
    double yVal = 0.0;
    double testValY = yMax + 10.0;
    eGamma = -1;
    int eBin;
    // cout<<"hi_1"<<endl;
    //test val is testValY = gMax*1.02 + 10.0                                   
    while(testValY>yVal){//Monte Carlo the event to get brem spectrum  
      
	  eGamma = gRandom->Uniform(eGammaLow_V1,eGammaHigh_V1);//Grab a photon energy
      testValY = gRandom->Uniform(0.0,yMax); //Grab a test value    
      eBin = hGvsE->GetXaxis()->FindBin(eGamma); //Find bin                 
      yVal = hGvsE->GetBinContent(eBin); //Use bin to find value   
    }

    hEGamma_V1->Fill(eGamma);
    }
  }
  */
  
  //eGamma = gRandom->Uniform(eGammaLow_V1,eGammaHigh_V1);//Grab a photon energy

  



  std::ofstream lundFile;
  std::ofstream genr8File;
  if (outType == 2) lundFile.open(outFileNameLund);
  if (outType == 3) genr8File.open(outFileNameGenr8);
  
  
  for(Int_t i=0;i<nToRun;i++){
    if (i%10000 == 0 && i > 0) cout<<"Events processed = "<<i<<endl;
    
    bool keep = false;
    // double eGamma = 0;
    double testVal = 0;
    //eGammaLow_V1  = 7.0;//was 8.0//was 3 prior
    //eGammaHigh_V1 = 11.0;//was 9.0
    /*
    while (keep == false)
      {
      eGamma = gRandom->Uniform(eGammaLow,eGammaHigh);
      testVal = gRandom->Uniform(0.0,1.0/eGammaLow);
      if (testVal < 1.0/eGamma) keep = true; 
      }
    *///This will shape the photon distriibution 
    eGamma = gRandom->Uniform(eGammaLow_V1,eGammaHigh_V1);

    //int eBin = hEDist->FindBin(eGamma);
    //cout<<"e:e = "<<eGamma<<" : "<<hEDist->GetBinContent(eBin) <<endl;
    //eGamma = hEDist->GetBinContent(eBin);

    hEGamma->Fill(eGamma);
    hEGamma_V1->Fill(eGamma);    
    TLorentzVector beam(0,0,eGamma,eGamma);
    TLorentzVector pTot4Vec = beam + target;
    //double sVal = pTot4Vec.Mag2();
    double wVal = pTot4Vec.Mag();  //Invariant mass of initial state
    
    //
    //beam and target to cm frame
    //
    
    TLorentzVector beamCM = beam;
    beamCM.Boost(-pTot4Vec.BoostVector());
    TLorentzVector targetCM = target;
    targetCM.Boost(-pTot4Vec.BoostVector());
    
    //    mXiStar = gRandom->Uniform(mXiStarCenter-gXiStar,mXiStarCenter+gXiStar);   
    keep = false;
    int iTmp = 0;
    while (keep == false){
      //mHyperon = gRandom->BreitWigner(mHyperonCenter,gHyperon);  
      mHyperon = gRandom->Uniform(2.16,3.26);//3.78);//(2.0,3.1);//was 2.04-2.26  
      //     mHyperon = gRandom->Uniform(2.046,2.6);  
      hMassHyperon->Fill(mHyperon);    
      if (mHyperon >= mKplus+mXiStarCenter-gXiStar && wVal >= (mHyperon+mKplus)) keep = true;//I took the 3* off the gXistar
      iTmp++;
      if (iTmp > 1000){
	cout<<"Oh nooo!!!!! : "<<endl;
        //cout<<"test1: "<<mHyperon<<" : "<<mKplus+mXiStar<<endl;
	//cout<<"test2: "<<wVal<<" : "<<mHyperon+mKplus<<endl;
	//cout<<"wVal - mKplus : mXiStar + mKplus = "<<wVal - mKplus<<" : "<<mXiStar + mKplus<<endl;
	//cout<<"mXiStar = "<<mXiStar<<endl;
      }
    } 

    //You are going to try and change the code to such that the hyperon is the only dynamic variable 
    hMassHyperon_V2->Fill(mHyperon);    
    
    keep = false;
    while (keep == false){
      mXiStar = gRandom->BreitWigner(mXiStarCenter,gXiStar);  
      //cout<<"mXiStarCenter+gXiStar"<<mXiStarCenter+gXiStar<<endl;
      //cout<<"mXiStarCenter-gXiStar"<<mXiStarCenter-gXiStar<<endl;
      //      mXiStar = gRandom->Uniform(1.5251,1.5449);  
      if (mXiStar >= mpi0+mXiGS && mHyperon >= (mKplus+mXiStar)) keep = true;
      iTmp++;
      if (iTmp > 1000){
	cout<<"Oh nooo222222222!!!!! : "<<endl;
	cout<<"test1: "<<mXiStar<<" : "<<mpi0+mXiGS<<endl;
	cout<<"test2: "<< mHyperon<<" : "<<mKplus+mXiStar<<endl;
      }
    }
    

    TLorentzVector HyperonCM4Vec;
    TLorentzVector Kaon1CM4Vec;

   
    OneBodyDecay(gRandom,wVal,mHyperon,mKplus,&HyperonCM4Vec,&Kaon1CM4Vec);
    //Since wVal is invariant, it is the same in every frame
    //The center of mass frame is especially nice. Let us "decay"
    //the initial state into the final state = proton and rho meson
    //

    
    /*
      cout<<"Momentum Along X, for Y* and K1 is \n"<<HyperonCM4Vec.Px()<<" , "<<Kaon1CM4Vec.Px()<<endl;
    cout<<"Momentum Along Y, for Y* and K1 is \n"<<HyperonCM4Vec.Py()<<" , "<<Kaon1CM4Vec.Py()<<endl;
    cout<<"Momentum Along Z, for Y* and K1 is \n"<<HyperonCM4Vec.Pz()<<" , "<<Kaon1CM4Vec.Pz()<<endl;
    cout<<"The Mass for Y* and K1 is \n"<<HyperonCM4Vec.M()<<" , "<<Kaon1CM4Vec.M()<<endl;
    */
    
    
    keep = false;
    while (keep == false){      
      //DOING THE SAME FOR EXCITED HYPERON
      //Find extreme values of t (theta = 0 and pi)
      
      //SAME BUT FOR THE HYPERON
      TLorentzVector pKplus1CM4VecExtA = TLorentzVector(0,0, Kaon1CM4Vec.P(),Kaon1CM4Vec.E());
      TLorentzVector pKplus1CM4VecExtB = TLorentzVector(0,0,-Kaon1CM4Vec.P(),Kaon1CM4Vec.E());
      TLorentzVector p1MinusP3ExtA_Kaon = beamCM - pKplus1CM4VecExtA;
      TLorentzVector p1MinusP3ExtB_Kaon = beamCM - pKplus1CM4VecExtB;
      
      //      double tValAbsExtA_Kaon = abs(p1MinusP3ExtA_Kaon.Mag2());
      double tValAbsExtB_Kaon = abs(p1MinusP3ExtB_Kaon.Mag2());
      
      TLorentzVector p1MinusP3_Kaon = targetCM - Kaon1CM4Vec;
      
      double tValAbs_Kaon = abs(p1MinusP3_Kaon.Mag2());
      //double bVal = -0.5;
      double bVal_Kaon = -1.05;// This some value for runs until 240, then was changed to                                    //from run 270-300 the tslope is -1.2   
                              //from run 300+ to XYZ the tslope is -1.1
                              //from runs 312-323 tslope of 1.02
                              //from runs 350->361 so tslope of 1.12
	                      //Runs 400 to XXX have a tslope of 1.05 
      //  double expTest         = exp(bVal*tValAbs);
      //double expTestMax      = exp(bVal*tValAbsE to xtB);     
      double expTest_Kaon    = exp(bVal_Kaon*tValAbs_Kaon);
      double expTestMax_Kaon = exp(bVal_Kaon*tValAbsExtB_Kaon);
      //     double randTest = gRandom->Uniform(0.0,expTestMax);     
      //  if (randTest < expTest) keep = true;
 
      double randTest = gRandom->Uniform(0.0,expTestMax_Kaon);     
      if (randTest < expTest_Kaon) keep = true;
      
    }
    //K+Y*->K+(K+Xi*-)
    //Boost back to lab of K+K+Xi*-
    
  
    TLorentzVector Hyperon4Vec = HyperonCM4Vec;
    Hyperon4Vec.Boost(pTot4Vec.BoostVector());
    hYStar->Fill(Hyperon4Vec.M());

    TLorentzVector Kaon1_4Vec = Kaon1CM4Vec;
    Kaon1_4Vec.Boost(pTot4Vec.BoostVector());
    k1P->Fill(Kaon1_4Vec.P());

    TLorentzVector missEP = pTot4Vec-Hyperon4Vec-Kaon1_4Vec;
    //hMassHyperon->Fill(Hyperon4Vec.M());    
    
    TLorentzVector p1MinusP3_Kaon = beam - Kaon1_4Vec;
    double tValAbs_Kaon = abs(p1MinusP3_Kaon.Mag2());
    /////////////Hyperon Decay YRF////////////////

    TLorentzVector XiStarYRF4Vec;
    TLorentzVector Kaon2YRF4Vec;
    OneBodyDecay(gRandom,mHyperon,mXiStar,mKplus,&XiStarYRF4Vec,&Kaon2YRF4Vec);
    
    /*
    cout<<"Momentum Along X, for Y* and K1 is \n"<<XiStarYRF4Vec.Px()<<" , "<<Kaon2YRF4Vec.Px()<<endl;
    cout<<"Momentum Along Y, for Y* and K1 is \n"<<XiStarYRF4Vec.Py()<<" , "<<Kaon2YRF4Vec.Py()<<endl;
    cout<<"Momentum Along Z, for Y* and K1 is \n"<<XiStarYRF4Vec.Pz()<<" , "<<Kaon2YRF4Vec.Pz()<<endl;
    cout<<"The Mass for Y* and K1 is \n"<<XiStarYRF4Vec.M()<<" , "<<Kaon2YRF4Vec.M()<<endl;
    */
    
    TLorentzVector XiStar4Vec = XiStarYRF4Vec;
    XiStar4Vec.Boost(Hyperon4Vec.BoostVector());

    TLorentzVector Kaon2_4Vec = Kaon2YRF4Vec;
    Kaon2_4Vec.Boost(Hyperon4Vec.BoostVector());
    
    TLorentzVector XiMinus4Vec_CM;
    TLorentzVector Pi04Vec_CM;
    OneBodyDecay(gRandom,XiStar4Vec.M(),mXiGS,mpi0,&XiMinus4Vec_CM,&Pi04Vec_CM);

    TLorentzVector XiMinus4Vec = XiMinus4Vec_CM;
    XiMinus4Vec.Boost(XiStar4Vec.BoostVector());

    TLorentzVector Pi04Vec = Pi04Vec_CM;
    Pi04Vec.Boost(XiStar4Vec.BoostVector());

    hMassXiMinusVsP->Fill(XiMinus4Vec.M(),XiMinus4Vec.P());
    hMassPi0VsP->Fill(Pi04Vec.M(),Pi04Vec.P());

    TLorentzVector KLow;
    TLorentzVector KHigh;
    
    double theta_1 = Kaon1_4Vec.Theta()*180.0/3.14159;
    double theta_2 = Kaon2_4Vec.Theta()*180.0/3.14159;

    if(theta_1<theta_2)
      {
	KLow  = Kaon1_4Vec;
	KHigh = Kaon2_4Vec;
      }
    else
      {
	KLow  = Kaon2_4Vec;
	KHigh = Kaon1_4Vec;
      }

    TLorentzVector HyperonLow  = XiStar4Vec + KLow ;
    TLorentzVector HyperonHigh = XiStar4Vec + KHigh;
    
    hYStarLow->Fill(HyperonLow.M());
    hYStarHigh->Fill(HyperonHigh.M());
    
    //Decay the ground state cascade

    TLorentzVector Lambda4Vec_CM;
    TLorentzVector PiMinus4Vec_CM;
    OneBodyDecay(gRandom,XiMinus4Vec.M(),mLambda,mpim,&Lambda4Vec_CM,&PiMinus4Vec_CM);

    TLorentzVector Lambda4Vec = Lambda4Vec_CM;
    Lambda4Vec.Boost(XiMinus4Vec.BoostVector());

    TLorentzVector PiMinus4Vec = PiMinus4Vec_CM;
    PiMinus4Vec.Boost(XiMinus4Vec.BoostVector());

    hMLambdaVsP->Fill(Lambda4Vec.M(),Lambda4Vec.P());
    hPiMinus1MvsP->Fill(PiMinus4Vec.M(),PiMinus4Vec.P());
    TLorentzVector Proton4Vec_CM;
    TLorentzVector PiMinus_2_4Vec_CM;
    OneBodyDecay(gRandom,Lambda4Vec.M(),mPro,mpim,&Proton4Vec_CM,&PiMinus_2_4Vec_CM);
    TLorentzVector Proton4Vec = Proton4Vec_CM;
    Proton4Vec.Boost(Lambda4Vec.BoostVector());

    TLorentzVector PiMinus_2_4Vec = PiMinus_2_4Vec_CM;
    PiMinus_2_4Vec.Boost(Lambda4Vec.BoostVector());

    hProtonMvsP->Fill(Proton4Vec.M(),Proton4Vec.P());
    hPiMinus2MvsP->Fill(PiMinus_2_4Vec.M(),PiMinus_2_4Vec.P());

    //Filling Structure

//Make a vertex vector for the event                                                  

    // TLorentzVector pX4Vec = pTot4Vec - Kaon1_4Vec - Kaon2_4Vec - XiMinus4Vec - Pi04Vec;
    //hXpx->Fill(pX4Vec.Px());
    //hXpy->Fill(pX4Vec.Py());
    //hXpz->Fill(pX4Vec.Pz());
    //hXe->Fill(pX4Vec.E());

    TVector3 vertVec(0,0,0);
   //Load the myGen structure with final-state particles                                 
   myGen_t myGen;
   myGen.weight      = 1.0;
   myGen.runNumber   = runNumber;
   myGen.eventNumber = i+1;
   myGen.vertex3Vec.push_back(vertVec);
   //K+ 1                                                                              
   myGen.mo4Vec.push_back(Kaon1_4Vec);
   myGen.pidVecPDG.push_back(321);
   myGen.pidVecGeant.push_back(11);
   myGen.chargeVec.push_back(1);
  //K+ 2                                                                              
   myGen.mo4Vec.push_back(Kaon2_4Vec);
   myGen.pidVecPDG.push_back(321);
   myGen.pidVecGeant.push_back(11);
   myGen.chargeVec.push_back(1);
   //Xi -                                                                               
   //myGen.mo4Vec.push_back(XiMinus4Vec);
   //myGen.pidVecPDG.push_back(3312);
   //myGen.pidVecGeant.push_back(23);
   //myGen.chargeVec.push_back(-1);
   //Pi0                                                                                 
   myGen.mo4Vec.push_back(Pi04Vec);
   myGen.pidVecPDG.push_back(111);
   myGen.pidVecGeant.push_back(7);
   myGen.chargeVec.push_back(0);
   //Lambda
   //myGen.mo4Vec.push_back(Lambda4Vec);
   //myGen.pidVecPDG.push_back(3122);
   //myGen.pidVecGeant.push_back(18);
   //myGen.chargeVec.push_back(0);
   //Pi-Minus
   myGen.mo4Vec.push_back(PiMinus4Vec);
   myGen.pidVecPDG.push_back(-211);
   myGen.pidVecGeant.push_back(9);
   myGen.chargeVec.push_back(-1);
   //Proton
   myGen.mo4Vec.push_back(Proton4Vec);
   myGen.pidVecPDG.push_back(2212);
   myGen.pidVecGeant.push_back(14);
   myGen.chargeVec.push_back(1);
   //PiMinus2
   myGen.mo4Vec.push_back(PiMinus_2_4Vec);
   myGen.pidVecPDG.push_back(-211);
   myGen.pidVecGeant.push_back(9);
   myGen.chargeVec.push_back(-1);
   
   
   
   if (outType == 2) writeEventLund(myGen, &lundFile);                           
   if (outType == 3) writeEventGenr8(myGen, &genr8File);
 
   //START TEST
   //cout<<"---------------"<<endl;
   x4Vec = beam + target;
   for (int i = 0; i<(int)myGen.pidVecGeant.size(); i++){
     x4Vec -= myGen.mo4Vec[i];
   }
   hDeltaPx->Fill(x4Vec.Px());
   hDeltaPy->Fill(x4Vec.Py());
   hDeltaPz->Fill(x4Vec.Pz());
   hDeltaE->Fill(x4Vec.E());
   //END TEST

    hMassK2->Fill(Kaon2_4Vec.M());
    //hTVal->Fill(tValAbs);
   
    hTVal->Fill(tValAbs_Kaon);
    hMoVsThetaY->Fill(Hyperon4Vec.Theta()*180.0/3.14159,Hyperon4Vec.P());      
    hMoVsThetaK1->Fill(Kaon1_4Vec.Theta()*180.0/3.14159,Kaon1_4Vec.P());      
    hMoVsThetaK2->Fill(Kaon2_4Vec.Theta()*180.0/3.14159,Kaon2_4Vec.P());      

    hMoVsThetaKLow->Fill(KLow.Theta()*180.0/3.14159,KLow.P());      
    hMoVsThetaKHigh->Fill(KHigh.Theta()*180.0/3.14159,KHigh.P());      
    hEXiStarRecon->Fill(XiStar4Vec.E());
    hXiStarRecon->Fill(XiStar4Vec.M());
    hXiStarvsCosTheta->Fill(XiStar4Vec.M(),XiStar4Vec.CosTheta());
    hXiStarvsCosThetaCM->Fill(XiStarYRF4Vec.M(),XiStarYRF4Vec.CosTheta());
    hXiStarvsP->Fill(XiStar4Vec.M(),XiStar4Vec.P());
    hXiStar_P->Fill(XiStar4Vec.P());
    hXiStarvsTheta->Fill(XiStar4Vec.M(),XiStar4Vec.Theta());
  }


if (outType == 1 || outType == 3){                                                   
   TFile *outFile = new TFile(outFileName,"RECREATE");
   outFile->cd();

  //outFile->cd();
  //Write histograms to file here

   hDeltaPx->Write();
   hDeltaPy->Write();
   hDeltaPz->Write();
   hDeltaE->Write();

  hDummy1d->Write();
  hDummy2d->Write();
  hEGamma_V1->Write();
  hEGamma->Write();
  hMassK2->Write();
  hTVal->Write();

  //hXpx->Write();
  //hXpy->Write();
  //hXpz->Write();
  //hXe->Write();
  
  hMoVsThetaY->Write();
  hMoVsThetaK1->Write();
  hMoVsThetaK2->Write();
  
  hMoVsThetaKLow->Write();
  hMoVsThetaKHigh->Write();
  hEXiStarRecon->Write();

  hXiStarvsCosTheta->Write();
  hXiStarvsCosThetaCM->Write();
  hMassXiMinusVsP->Write();
  hMassPi0VsP->Write();
  hXiStarRecon->Write();
  hXiStarvsP->Write();
  hXiStar_P->Write();
  hXiStarvsTheta->Write();
  
  //outFile->Write();

  hMassHyperon->Write();
  k1P->Write();
  hYStarLow->Write();
  hYStarHigh->Write();
  hMLambdaVsP->Write();
  hProtonMvsP->Write();
  hPiMinus1MvsP->Write();
  hPiMinus2MvsP->Write();
  hMassHyperon_V2->Write();    



  outFile->Write();


}


if (outType == 2) lundFile.close();                                               
if (outType == 3) genr8File.close();



return 0;

  

}

void printUsage(){
  fprintf(stderr,"\nUsage:");
  fprintf(stderr,"\nmyEventGen [-switches]\n");
  fprintf(stderr,"\nSWITCHES:\n");
  fprintf(stderr,"-h\tPrint this message\n");
  fprintf(stderr,"-s<arg>\tRandom seed. Default is 0\n");
  fprintf(stderr,"-o<arg>\tOutFileName. Default is myOutFile.txt\n\n");
  fprintf(stderr,"-r<arg>\tRun number. Default is 99999\n\n");
  fprintf(stderr,"-t<arg>\tOutFile type: 1->Root, 2->Lund, 3->Genr8\n\n");

  cout<<"The current default operation is equivalent to the command:"<<endl;
  //  cout<<"myEventGen -s0 -n100000 -omyOutFile.root\n"<<endl;
  cout<<"myEventGen -s0 -t1 -r99999 -n100000 -omyOutFile.root\n"<<endl;
  exit(0);
}

void writeEventLund(myGen_t myGen, std::ofstream *lundFile){
 //See https://gemc.jlab.org/gemc/html/documentation/generator/lund.html                 
 //LUND HEADER                                                                           
 *lundFile<<myGen.pidVecPDG.size();//  1) Number of particles                            
 *lundFile<<"\t0.938";             //  2) Mass of target:          NOT REQUIRED          
 *lundFile<<"\t1";                 //  3) Atomic number of target: NOT REQUIRED          
 *lundFile<<"\t0";                 //  4) Target polarization:     NOT REQUIRED          
 *lundFile<<"\t0";                 //  5) Beam polarization                              
 *lundFile<<"\t11";                //  6) Beam type:               NOT REQUIRED          
 *lundFile<<"\t10";                //  7) Beam energy:             NOT REQUIRED          
 *lundFile<<"\t2212";              //  8) Interacted nucleon ID:   NOT REQUIRED          
 *lundFile<<"\t0";                 //  9) Process ID:              NOT REQUIRED          
 *lundFile<<"\t"<<myGen.weight;    // 10) Event weight:            NOT REQUIRED          
 endl(*lundFile); //Write the header            
 //LUND particle                                                                         
 for (int i = 0; i<(int)myGen.pidVecPDG.size(); i++){
   *lundFile<<i+1;                          //  1) Particle index                        
   *lundFile<<"\t0";                        //  2) Lifetime:            NOT REQUIRED     
   *lundFile<<"\t1";                        //  3) Type (1 is active)                    
   *lundFile<<"\t"<<myGen.pidVecPDG[i];     //  4) Particle ID                           
   *lundFile<<"\t0";                        //  5) Index of parent:     NOT REQUIRED     
   *lundFile<<"\t-1";                       //  6) Index of daughter:   NOT REQUIRED     
   *lundFile<<"\t"<<myGen.mo4Vec[i].Px();   //  7) Momentum in x (GeV)                   
   *lundFile<<"\t"<<myGen.mo4Vec[i].Py();   //  8) Momentum in y (GeV)                   
   *lundFile<<"\t"<<myGen.mo4Vec[i].Pz();   //  9) Momentum in z (GeV)                   
   *lundFile<<"\t"<<myGen.mo4Vec[i].E();    // 10) Energy (GeV):        NOT REQUIRED     
   *lundFile<<"\t"<<myGen.mo4Vec[i].Mag();  // 11) Mass (GeV):          NOT REQUIRED     
   *lundFile<<"\t"<<myGen.vertex3Vec[i].X();// 12) Vertex x (cm)                         
   *lundFile<<"\t"<<myGen.vertex3Vec[i].Y();// 12) Vertex x (cm)                         
   *lundFile<<"\t"<<myGen.vertex3Vec[i].Z();// 12) Vertex x (cm)                         
   endl(*lundFile); //Write the particle                                                 
 }
}

void writeEventGenr8(myGen_t myGen, std::ofstream *genr8File){
 //See page 5 of https://halldweb.jlab.org/DocDB/0000/000011/001/genr8_v_1_0.pdf         
 //Genr8 HEADER                                                                          
 *genr8File<<myGen.runNumber;                // 1) Run number                            
 *genr8File<<"\t"<<myGen.eventNumber;        // 2) Event number                          
 *genr8File<<"\t"<<myGen.pidVecGeant.size(); // 3) Number of particles                   
 endl(*genr8File); //Write the header                                                    
 //Genr8 particle                                                                        
 double pxTot = 0;
 double pyTot = 0;
 double pzTot = 0;
 double eTot = 0;

 for (int i = 0; i<(int)myGen.pidVecGeant.size(); i++){
   *genr8File<<i+1;                          // 1) Particle index                        
   *genr8File<<"\t"<<myGen.pidVecGeant[i];   // 2) Particle ID                           
   *genr8File<<"\t"<<myGen.mo4Vec[i].Mag();  // 3) Mass                                  
   endl(*genr8File);                         //Write line                                
   *genr8File<<"\t"<<myGen.chargeVec[i];     // 1) Charge                                
   *genr8File<<"\t"<<myGen.mo4Vec[i].Px();   // 2) Momentum in x                         
   *genr8File<<"\t"<<myGen.mo4Vec[i].Py();   // 3) Momentum in y                         
   *genr8File<<"\t"<<myGen.mo4Vec[i].Pz();   // 4) Momentum in z                         
   *genr8File<<"\t"<<myGen.mo4Vec[i].E();    // 5) Energy                                

   pxTot+= myGen.mo4Vec[i].Px();
   pyTot+= myGen.mo4Vec[i].Py();
   pzTot+= myGen.mo4Vec[i].Pz();
   eTot+=  myGen.mo4Vec[i].E();

   endl(*genr8File); //Write the particle                                                

   //cout<<"p4v = "<<pxTot<<" , "<<pyTot<<" , "<<pzTot<<" , "<<eTot<<endl;
 }
}