#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(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 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 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;
 
  TLorentzVector target(0,0,0,mPro);

  //Set the default output file
  sprintf(outFileName,"./myOutFile.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);
	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_V1       = new TH1D("hEGamma_V1","",400,3.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","",1500,1.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);  

  //  double eGammaLow_V1  = 3;
  double eGammaLow_V1  = 7;
  double eGammaHigh_V1 = 11;

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

  double bin_1    =15.6;
  double bin_1_er =3.95;
  double bin_2 = 7.2;
  double bin_2_er= 2.68; 
  double bin_3=39.6;
  double bin_3_er=6.3;
  double bin_4 = 56.7;
  double bin_4_er=7.53;
  double bin_5=16;
  double bin_5_er=4;
  double bin_6=35.33;
  double bin_6_er=5.94;
  double bin_7=26.65;
  double bin_7_er=4.86;
  double bin_8=19.16;
  double bin_8_er=4.4;
  double bin_9=20.52; 
  double bin_9_er=4.53;
  double bin_10=22.76;
  double bin_10_er=4.77;


  heff->SetBinContent(1,bin_1/3.97843e+12);
  heff->SetBinContent(2,bin_2/7.88475e+12);
  heff->SetBinContent(3,bin_3/1.38321e+13);
  heff->SetBinContent(4,bin_4/1.91518e+13);
  heff->SetBinContent(5,bin_5/1.57437e+13);
  heff->SetBinContent(6,bin_6/5.62737e+12);
  heff->SetBinContent(7,bin_7/6.66557e+12);
  heff->SetBinContent(8,bin_8/7.26285e+12);
  heff->SetBinContent(9,bin_9/6.35828e+12);
  heff->SetBinContent(10,bin_10/7.97716e+12);
  //  heff->Draw();
 
  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
  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   
    }

    
    //cout<<"hi_2"<<endl;
    //cout<<"E"<<endl;
    hEGamma_V1->Fill(eGamma);
    eff=heff->FindBin(eGamma);
    //cout<<"The efficiency for this bin is: "<<eff<<endl;
    hEg_eff->Fill(eGamma/eff);
  }
  
  
  
  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;
    double eGammaLow  = 3.0;//was 8.0
    double eGammaHigh = 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; 
      }*/
    //hEGamma->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());
    
    
    //cout<<"hi:1"<<endl;
    keep = false;
    int iTmp = 0;
    while (keep == false){
      //      mHyperon = gRandom->BreitWigner(mHyperonCenter,gHyperon);  
      mHyperon = gRandom->Uniform(2.03,2.27);  
      if (mHyperon >= mKplus+mXiStarCenter-3*gXiStar && wVal >= (mHyperon+mKplus)) keep = true;
      //if (wVal >= (mHyperon+mKplus)) keep = true;
      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;
      }
    } 
    //cout<<"hi:2"<<endl;
    hMassHyperon->Fill(mHyperon);

    //You are going to try and change the code to such that the hyperon is the only dynamic variable 

    mXiStar = gRandom->Uniform(mXiStarCenter-gXiStar,mXiStarCenter+gXiStar);  
    /*keep = false;
    while (keep == false){
      //      mXiStar = gRandom->BreitWigner(mXiStarCenter,gXiStar);  
      mXiStar = gRandom->Uniform(mXiStarCenter-gXiStar,mXiStarCenter+gXiStar);  
      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(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());
  
    TLorentzVector Kaon1_4Vec = Kaon1CM4Vec;
    Kaon1_4Vec.Boost(pTot4Vec.BoostVector());
    
    TLorentzVector p1MinusP3_Kaon = beam - Kaon1_4Vec;
    double tValAbs_Kaon = abs(p1MinusP3_Kaon.Mag2());
    /////////////Hyperon Decay YRF////////////////

    TLorentzVector XiStarYRF4Vec;
    TLorentzVector Kaon2YRF4Vec;
    OneBodyDecay(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(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;
      }


    //Filling Structure

//Make a vertex vector for the event                                                  
    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);

 if (outType == 2) writeEventLund(myGen, &lundFile);                           
 if (outType == 3) writeEventGenr8(myGen, &genr8File);


    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){                                                   
   TFile *outFile = new TFile(outFileName,"RECREATE");
   outFile->cd();

  //outFile->cd();
  //Write histograms to file here
  hDummy1d->Write();
  hDummy2d->Write();
  hEGamma_V1->Write();
  hMassK2->Write();
  hTVal->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();
  hEg_eff->Write();
  //outFile->Write();

  hMassHyperon->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                                                                        
 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                                
   endl(*genr8File); //Write the particle                                                
 }
}