double m23SqMax(double m1, double m2, double m3, double m12, double M123){
  double m1Sq = pow(m1,2);
  double m2Sq = pow(m2,2);
  double m3Sq = pow(m3,2);
  double m12Sq = pow(m12,2);
  double M123Sq = pow(M123,2);
  double e2Star = (m12Sq - m1Sq + m2Sq)/(2*m12);
  double e3Star = (M123Sq - m12Sq - m3Sq)/(2*m12);
  double e2StarSq = pow(e2Star,2);
  double e3StarSq = pow(e3Star,2);
  double m23SqMaxVal = pow(e2Star+e3Star,2) - pow(sqrt(e2StarSq-m2Sq)-sqrt(e3StarSq-m3Sq),2);
  return m23SqMaxVal;
}

double m23SqMin(double m1, double m2, double m3, double m12, double M123){
  double m1Sq = pow(m1,2);
  double m2Sq = pow(m2,2);
  double m3Sq = pow(m3,2);
  double m12Sq = pow(m12,2);
  double M123Sq = pow(M123,2);
  double e2Star = (m12Sq - m1Sq + m2Sq)/(2*m12);
  double e3Star = (M123Sq - m12Sq - m3Sq)/(2*m12);
  double e2StarSq = pow(e2Star,2);
  double e3StarSq = pow(e3Star,2);
  double m23SqMinVal = pow(e2Star+e3Star,2) - pow(sqrt(e2StarSq-m2Sq)+sqrt(e3StarSq-m3Sq),2);
  return m23SqMinVal;
}

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

  TH1D * hEGamma = new TH1D("hEGamma","",200,7.5,9.5);  
  TH2D * hDalitz = new TH2D("hDalitz","",120,0.0,12.0,200,0.0,20.0);  

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

  double pi = 3.14159;

  double mPro       = 0.938;
  double mPip       = 0.13957;
    
  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 = 8.5;
    double testVal = 0;
    if (iShow != 1) {
      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
    /////////
    double M123 = wVal;
    double m1 = mPip;
    double m2 = mPip;
    double m3 = mPro;
    double m12SqLimitLow = pow(m1+m2,2);
    double m12SqLimitHigh = pow(M123-m3,2);

    double m23SqLimitLow = pow(m2+m3,2);
    double m23SqLimitHigh = pow(M123-m1,2);
    double m12Sq=0;
    double m23Sq=0;
    
    keep = false;
    while (keep == false){
      m12Sq = gRandom->Uniform(m12SqLimitLow,m12SqLimitHigh);
      m23Sq = gRandom->Uniform(m23SqLimitLow,m23SqLimitHigh);
      double m12 = sqrt(m12Sq); 
      double m23SqMinVal = m23SqMin(m1,m2,m3,m12,M123);
      double m23SqMaxVal = m23SqMax(m1,m2,m3,m12,M123);
      if (m23Sq > m23SqMinVal && m23Sq < m23SqMaxVal) keep = true; 
    }    
    hDalitz->Fill(m12Sq,m23Sq);

  }//End event generation
  if (iShow == 1) hDalitz->Draw("colz"); 
  if (iShow == 2) hDalitz->Draw("colz"); 
}
void mc4a(int iShow){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  TH1D * hEGamma = new TH1D("hEGamma","",200,7.5,9.5);  
  TH1D * hTVal = new TH1D("hTVal","",200,0.0,20);  
  TH2D * hMoVsThetaPro = new TH2D("hMoVsThetaPro","",90,0.0,90.0,100,0.0,10.0);  
  TH2D * hMoVsThetaRho = new TH2D("hMoVsThetaRho","",90,0.0,90.0,100,0.0,10.0);  
  TH1D * hMassRho = new TH1D("hMassRho","",200,0.0,2.0);  
  TH1D * hMassPipPim = new TH1D("hMassPipPim","",200,0.0,2.0);  
  TH2D * hMoVsThetaPim = new TH2D("hMoVsThetaPim","",90,0.0,90.0,100,0.0,10.0);  
  TH2D * hMoVsThetaPip = new TH2D("hMoVsThetaPip","",90,0.0,90.0,100,0.0,10.0);  
  TH2D * hDalitz = new TH2D("hDalitz","",120,0.0,12.0,200,0.0,20.0);

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

  double pi = 3.14159;

  double mPro       = 0.938;
  double mPip       = 0.13957;
  double mRhoCenter = 0.7755;
  double rhoWidth   = 0.1464;
  double mRho       = mRhoCenter;
    
  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 = 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());
    
    keep = false;
    while (keep == false){
      mRho = gRandom->BreitWigner(mRhoCenter,rhoWidth);  
      if (mRho >= 2*mPip && wVal >= (mRho+mPro)) keep = true;
    }
    hMassRho->Fill(mRho);

    //TLorentzVector pRhoCM4Vec = TLorentzVector(pRhoXCM,pRhoYCM,pRhoZCM,eRhoCM);
    //TLorentzVector pProCM4Vec = TLorentzVector(pProXCM,pProYCM,pProZCM,eProCM);
    TLorentzVector pRhoCM4Vec;
    TLorentzVector pProCM4Vec;


    //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;
    keep = false;
    while (keep == false){      
      //
      //Now phase-space generate the angles
      double phiRho      = gRandom->Uniform(-pi,pi);    
      double cosThetaRho = gRandom->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
      pRhoCM4Vec = TLorentzVector(pRhoXCM,pRhoYCM,pRhoZCM,eRhoCM);
      pProCM4Vec = TLorentzVector(pProXCM,pProYCM,pProZCM,eProCM);
      //
      //Find extreme values of t (theta = 0 and pi)
      TLorentzVector pProCM4VecExtA = TLorentzVector(0,0,pProCM,eProCM);
      TLorentzVector pProCM4VecExtB = TLorentzVector(0,0,-pProCM,eProCM);
      TLorentzVector p1MinusP3ExtA = targetCM - pProCM4VecExtA;
      TLorentzVector p1MinusP3ExtB = targetCM - pProCM4VecExtB;
      double tValAbsExtA = abs(p1MinusP3ExtA.Mag2());
      double tValAbsExtB = abs(p1MinusP3ExtB.Mag2());
      //cout<<"ExtA, ExtB = "<<p1MinusP3ExtA.Mag2()<<" , "<<p1MinusP3ExtB.Mag2()<<endl;
      //keep = true;
      TLorentzVector p1MinusP3 = targetCM - pProCM4Vec;
      double tValAbs = abs(p1MinusP3.Mag2());
      double bVal = -0.5;
      double expTest = exp(bVal*tValAbs);
      double expTestMax = exp(bVal*tValAbsExtB);
      double randTest = gRandom->Uniform(0.0,expTestMax);     
      if (randTest < expTest) keep = true;
      
    }

    //
    //Boost back to lab
    //
    TLorentzVector pPro4Vec = pProCM4Vec;
    pPro4Vec.Boost(pTot4Vec.BoostVector());
  
    TLorentzVector pRho4Vec = pRhoCM4Vec;
    pRho4Vec.Boost(pTot4Vec.BoostVector());
    
    TLorentzVector p1MinusP3 = target - pPro4Vec;
    double tValAbs = abs(p1MinusP3.Mag2());
    hTVal->Fill(tValAbs);
    hMoVsThetaPro->Fill(pPro4Vec.Theta()*180.0/3.14159,pPro4Vec.P());      
    hMoVsThetaRho->Fill(pRho4Vec.Theta()*180.0/3.14159,pRho4Vec.P());      

    //////Rho decay////////////
    //Let RFR = Rest frame of Rho
    double ePipRFR = mRho/2.0; 
    double pPipRFR = sqrt(pow(ePipRFR,2) - pow(mPip,2));    
    double ePimRFR = ePipRFR;
    double pPimRFR = pPipRFR;
    double phiPipRFR      = gRandom->Uniform(-pi,pi);    
    double cosThetaPipRFR = gRandom->Uniform(-1.0,1.0);
    //Obtain the components of the pip 4-vector
    double pPipZRFR          = pPipRFR*cosThetaPipRFR;
    double pPipTransverseRFR = sqrt(pow(pPipRFR,2) - pow(pPipZRFR,2));    
    double pPipXRFR          = pPipTransverseRFR*cos(phiPipRFR);
    double pPipYRFR          = pPipTransverseRFR*sin(phiPipRFR);
    //Load the RFR 4-vectors for pip and pim
    TLorentzVector pPipRFR4Vec = TLorentzVector(pPipXRFR,pPipYRFR,pPipZRFR,ePipRFR);
    TLorentzVector pPimRFR4Vec = TLorentzVector(-pPipXRFR,-pPipYRFR,-pPipZRFR,ePipRFR);
    ///////////////////////////
    //Boost the pions to the lab frame from the Rest-Frame of the Rho (RFR)
    //
    TLorentzVector pPip4Vec = pPipRFR4Vec;
    pPip4Vec.Boost(pRho4Vec.BoostVector());    

    TLorentzVector pPim4Vec = pPimRFR4Vec;
    pPim4Vec.Boost(pRho4Vec.BoostVector());    

    TLorentzVector pPipPim4Vec = pPip4Vec + pPim4Vec;
    hMassPipPim->Fill(pPipPim4Vec.Mag());      

    hMoVsThetaPip->Fill(pPip4Vec.Theta()*180.0/3.14159,pPip4Vec.P());      
    hMoVsThetaPim->Fill(pPim4Vec.Theta()*180.0/3.14159,pPim4Vec.P());      

    TLorentzVector pProPip4Vec = pPro4Vec + pPim4Vec;

    hDalitz->Fill(pPipPim4Vec.Mag2(),pProPip4Vec.Mag2());

  }
  if (iShow == 1){
    hEGamma->Draw();
  }
  if (iShow == 2){
    hMoVsThetaPro->Draw("colz");      
  } 
  if (iShow == 3){
    hMassRho->Draw();      
  } 
  if (iShow == 4){
    hMassPipPim->Draw();      
  } 
  if (iShow == 5){
    hMoVsThetaPip->Draw("colz");      
  } 
  if (iShow == 6){
    hMoVsThetaPim->Draw("colz");      
  } 
  if (iShow == 7){
    hTVal->Draw();      
  } 
  if (iShow == 8){
    hMoVsThetaRho->Draw("colz");      
  } 
  if (iShow == 9){
    hDalitz->Draw("colz");      
  } 
}