void mc1a(int iShow){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  TRandom3 * randGen =  new TRandom3();

  TH2D * h2dThrown   = new TH2D("h2dThrown","",100,-1.0,1.0,100,0.0,1.0);
  TH2D * h2dAccepted = new TH2D("h2dAccepted","",100,-1.0,1.0,100,0.0,1.0);
  TH1D * h1dAccepted = new TH1D("h1dAccepted","",100,-1.0,1.0);

  int nToGen = 1000000;
  for (int i = 0; i<nToGen; i++) {
    double xVal = randGen->Uniform(-1,1);
    double yVal = randGen->Uniform(0,1);
    h2dThrown->Fill(xVal,yVal);
    
    double functionVal = pow(xVal,2);
    if (yVal < functionVal) {
      h2dAccepted->Fill(xVal,yVal);
      h1dAccepted->Fill(xVal);
      
    }

  }
  
  h2dThrown->SetMinimum(0);
  if (iShow == 1) h2dThrown->Draw("colz");
  if (iShow == 2) h2dAccepted->Draw("colz");
  if (iShow == 3) h1dAccepted->Draw();
}

void mc1b(void){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  TRandom3 * randGen =  new TRandom3();
  TH1D * h1dAccepted = new TH1D("h1dAccepted","",100,-1.0,1.0);

  int nToKeep = 1000000;
  for (int i = 0; i<nToKeep; i++) {
    bool keep = false;
    while (keep == false) {
      double xVal = randGen->Uniform(-1,1);
      double yVal = randGen->Uniform(0,1);
      
      double functionVal = pow(xVal,2);
      if (yVal < functionVal) {
	h1dAccepted->Fill(xVal);
	keep = true;
      }
    }
  }
  h1dAccepted->Draw();
}

void mc1c(int iShow){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  TH1D * hEGamma = new TH1D("hEGamma","",200,7.5,9.5);  
  TH2D * hMoVsThetaPro = new TH2D("hMoVsThetaPro","",90,0.0,90.0,100,0.0,10.0);  

  TRandom3 * randGen =  new TRandom3();

  double mPro = 0.938;
  double mRho = 0.770;

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

  double eGammaLow  = 8.0;
  double eGammaHigh = 9.0;

  int nToGen = 1000000;
  for (int i = 0; i<nToGen; i++) {

    bool keep = false;
    double eGamma = 0;
    double testVal = 0;
    while (keep == false){
      eGamma = randGen->Uniform(eGammaLow,eGammaHigh);
      testVal = randGen->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
    //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
    //
    double eRhoCM = (sVal + pow(mRho,2) - pow(mPro,2))/(2*wVal); 
    double pRhoCM = sqrt(pow(eRhoCM,2) - pow(mRho,2));
    double eProCM = wVal - eRhoCM;
    double pProCM = pRhoCM;
  
    //
    //Now phase-space generate the angles
    double pi = 3.14159;
    double phiRho      = randGen->Uniform(-pi,pi);    
    //double cosThetaRho = randGen->Uniform(-1.0,1.0);
    double cosThetaRho = randGen->Uniform(-1.0,1.0);
    //Obtain the components of the rho 4-vector
    double pRhoZCM          = pRhoCM*cosThetaRho;
    double pRhoTransverseCM = sqrt(pow(pRhoCM,2) - pow(pRhoZCM,2));
    double pRhoXCM          = pRhoTransverseCM*cos(phiRho);
    double pRhoYCM          = pRhoTransverseCM*sin(phiRho);
    //Obtain the components of the proton 4-vector
    double pProXCM          = -pRhoXCM;
    double pProYCM          = -pRhoYCM;
    double pProZCM          = -pRhoZCM;
    //Create the rho CM 4-vector and the proton CM 4-vector
    TLorentzVector pRhoCM4Vec = TLorentzVector(pRhoXCM,pRhoYCM,pRhoZCM,eRhoCM);
    TLorentzVector pProCM4Vec = TLorentzVector(pProXCM,pProYCM,pProZCM,eProCM);
    //
    //Boost back to lab frame. First find beta and gamma of CM frame
    //Treat the entire system as a single particle of mass W = sqrt(s).
    //This means that the (total momentum) = gamma betaCM W
    //and that the (total energy) = gamma W
    // => gammaCM = (total energy)/ W
    //and betaCM = (total momentum)/(total energy)
    double betaCM  = pTot4Vec.P()/pTot4Vec.E();
    double gammaCM = pTot4Vec.E()/pTot4Vec.Mag();
    //
    //Boost back to lab
    //
    double pProX = pProXCM;
    double pProY = pProYCM;
    double pProZ = pProZCM*gammaCM + eProCM*gammaCM*betaCM;
    double ePro  = eProCM*gammaCM  + pProZCM*gammaCM*betaCM;
    TLorentzVector pPro4Vec = TLorentzVector(pProX,pProY,pProZ,ePro);
    //
    double pRhoX = pRhoXCM;
    double pRhoY = pRhoYCM;
    double pRhoZ = pRhoZCM*gammaCM + eRhoCM*gammaCM*betaCM;
    double eRho  = eRhoCM*gammaCM  + pRhoZCM*gammaCM*betaCM;
    TLorentzVector pRho4Vec = TLorentzVector(pRhoX,pRhoY,pRhoZ,eRho);    
    //
    //Test using a different (and easier) way to boost
    //
    //Boost from Lab to CM frame and test
    TLorentzVector pProCM4Vec_Test = pPro4Vec;
    pProCM4Vec_Test.Boost(-pTot4Vec.BoostVector());//Notice negative sign
    if (iShow == 2){ 
      cout<<"pProCM4Vec_Test/pProCM4Vec = "<<
	pProCM4Vec_Test.Px()/pProCM4Vec.Px()<<" x, "<<
	pProCM4Vec_Test.Py()/pProCM4Vec.Py()<<" y, "<<
	pProCM4Vec_Test.Pz()/pProCM4Vec.Pz()<<" z, "<<
	pProCM4Vec_Test.E()/pProCM4Vec.E()<<" e, "<<endl;
    }
    //Boost from CM to Lab frame and test
    TLorentzVector pPro4Vec_Test = pProCM4Vec;
    pPro4Vec_Test.Boost(pTot4Vec.BoostVector());//Use positive sign this time
    if (iShow == 3){ 
      cout<<"pPro4Vec_Test/pPro4Vec = "<<
	pPro4Vec_Test.Px()/pPro4Vec.Px()<<" x, "<<
	pPro4Vec_Test.Py()/pPro4Vec.Py()<<" y, "<<
	pPro4Vec_Test.Pz()/pPro4Vec.Pz()<<" z, "<<
	pPro4Vec_Test.E()/pPro4Vec.E()<<" e, "<<endl;
    }
    
    hMoVsThetaPro->Fill(pPro4Vec.Theta()*180.0/3.14159,pPro4Vec.P());      
  }
  if (iShow == 1){
    hEGamma->Draw();
  }
  if (iShow == 4){
    hMoVsThetaPro->Draw("colz");      
  } 
}