Double_t eff_func (Double_t *x, Double_t *par){
    
    // Dummy efficiency function taken from Qiao's talk GlueX-doc-3844 for ease of comparison
    // f = 0.4 + sqrt(Eg);
    
    // constants
    // Double_t const PI = 3.14159;
    Double_t a = par[0];
    Double_t b = par[1];
    
    Double_t MPI =0.139570;
    Double_t PI=3.14159;
    
    Double_t Eg = x[0];
    Double_t f=0;
    
    f = Eg>0? a + b*sqrt(Eg): 0;   // return zero for negative Eg
    if (f > 1) f = 1;             // cap efficiency at 1
    
    // f=1;                          // 100%
    return f;
}
Double_t tanLabCM_func (Double_t *x, Double_t *par){
    
    // Tan (angle_Lab) given the angle_Lab (degrees) in the CM system
    // betastar: Velocity (units of c) of CM relative to Lab (i.e. negative)
    // betaCM: Velocity of particle (units of c) in CM
    
    // constants
    // Double_t const PI = 3.14159;
    Double_t betastar= par[0];
    Double_t betaCM= par[1];
    
    Double_t MPI =0.139570;
    Double_t PI=3.14159;
    
    Double_t thetaCM = x[0]*PI/180.;
    
    Double_t tanLab=0;
    
    if (betastar < -1 || betastar >1) {
        cout << " ***tanLabCM_func: betastar=" << betastar << endl;
        return tanLab;
    }
    Double_t gammastar =  1/sqrt(1-betastar*betastar);
    
    tanLab = (1/gammastar)*(sin(thetaCM)/(cos(thetaCM) + betastar/betaCM));
    Double_t thetaLab = atan(tanLab)*180./PI;
    
    // cout << " betastar=" << betastar << " betaCM=" << betaCM << " thetaCM=" << thetaCM << " tanLab=" << tanLab << " thetalab=" << thetaLab << endl;
    
    return	thetaLab>0? thetaLab : thetaLab + 180.;
}
void AnglesCMLab(void)
{
// Consider transformation of angles from pi0->gg between CM and lab for the study of efficiency calculations in the calorimeters
    // Dec 17, 2018.
//

  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit(111111);
  gStyle->SetPadRightMargin(0.15);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadBottomMargin(0.15);
    
    
    char string[256];
    TString filename = "AnglesCMLab";
    // define function
    
    // Define historgrams
    Int_t nbins=20;
    
    // initialize random number seeds
    
    UInt_t  iseed = 2342411151;
    TRandom1 *r1 = new TRandom1();
    r1->SetSeed(iseed);
    iseed = 2342411155;
    TRandom1 *r2 = new TRandom1();
    r2->SetSeed(iseed);
    iseed = 2342411157;
    TRandom1 *r3 = new TRandom1();
    r3->SetSeed(iseed);
    iseed = 2342411153;
    TRandom1 *r4 = new TRandom1();
    r4->SetSeed(iseed);
    iseed = 2342411353;
    TRandom1 *r5 = new TRandom1();
    r5->SetSeed(iseed);
    
    // kinematics for Ma -> Mb + X (for pi0 Mpi0 -> g g)

    Double_t PI=3.14159;
    Double_t MPI =0.139570;
    Double_t paLab0 =2.;
    Double_t Mg = 0;
    Double_t Ma = MPI;
    Double_t Mb = Mg;
    Double_t EaLab0 = sqrt(paLab0*paLab0 + Ma*Ma);
    Double_t betastar = -paLab0/EaLab0;
    
    Double_t EbCM = Ma/2;
    Double_t pbCM = sqrt(EbCM*EbCM - Mb*Mb);
    Double_t betabCM = pbCM/EbCM;
    Double_t Ethresh = 0.05;      // 50 Mev Threshold
    Double_t BCALa=0.052;
    Double_t BCALb=0.036;
    // Double_t BCALa=0.01;
    // Double_t BCALb=0.01;
    
    // TString tag="palab";
    TString tag="energy";
    
    Double_t sigmaE = EaLab0*sqrt(BCALa*BCALa/EaLab0+BCALb*BCALb);
    Int_t nsigma=2;
    Double_t DeltaE = 2*nsigma*sigmaE;
    // Double_t DeltaE = 2*2*nsigma*sigmaE;
    cout << " Delta E=" << DeltaE << " GeV" << endl;
    
    const Int_t npts=100000;
    // const Int_t npts=10000000;
    
    Double_t xmin=0;
    Double_t xmax=2;
    Int_t npar = 2;
    Double_t aff=0.4;
    Double_t beff=1;
    
    TF1 *eff = new TF1("eff",eff_func,xmin,xmax,npar);
    eff->SetParameters(aff,beff);
    eff->SetParNames("a","b");
    
    Double_t thetaCMmin=0.1;
    Double_t thetaCMmax=180;
    TF1 *tanLabCM = new TF1("tanLabCM",tanLabCM_func,thetaCMmin,thetaCMmax,npar);
    tanLabCM->SetParameters(betastar,betabCM);
    tanLabCM->SetParNames("betastar","betabCM");
 
    TCanvas *c1 = new TCanvas("c1", "c1",200,10,1000,700);
    gPad->SetGridx();
    gPad->SetGridy();
    gPad->SetLogx();
    gPad->SetLogy();
    
    Double_t ymin=0.;
    Double_t ymax=200;
    tanLabCM->SetTitle("#pi^{0} #rightarrow #gamma #gamma");
    // tanLabCM->GetXaxis()->SetRangeUser(xmin,xmax);
    tanLabCM->GetYaxis()->SetRangeUser(ymin,ymax);
    tanLabCM->GetXaxis()->SetTitleSize(0.05);
    tanLabCM->GetYaxis()->SetTitleSize(0.05);
    tanLabCM->GetXaxis()->SetTitle("theta Lab (degrees)");
    tanLabCM->GetYaxis()->SetTitle("theta CM (degrees)");
    tanLabCM->SetLineColor(3);
    tanLabCM->Draw();
    
    TLegend *leg = new TLegend(0.6,0.2,0.8,0.4);
    TString text;
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM,text,"l");
    
    TF1 *tanLabCM2 = tanLabCM->DrawCopy("same");
    betastar = -0.5;
    tanLabCM2->SetParameters(betastar,betabCM);
    tanLabCM2->SetParNames("betastar","betabCM");
    tanLabCM2->SetLineColor(2);
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM2,text,"l");
    
    TF1 *tanLabCM3 = tanLabCM->DrawCopy("same");
    betastar = -0.7;
    tanLabCM3->SetParameters(betastar,betabCM);
    tanLabCM3->SetParNames("betastar","betabCM");
    tanLabCM3->SetLineColor(4);
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM3,text,"l");
    leg->Draw();
    
    TF1 *tanLabCM4 = tanLabCM->DrawCopy("same");
    betastar = -0.8;
    tanLabCM4->SetParameters(betastar,betabCM);
    tanLabCM4->SetParNames("betastar","betabCM");
    tanLabCM4->SetLineColor(1);
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM4,text,"l");
    leg->Draw();
    
    TF1 *tanLabCM5 = tanLabCM->DrawCopy("same");
    betastar = -0.9;
    tanLabCM5->SetParameters(betastar,betabCM);
    tanLabCM5->SetParNames("betastar","betabCM");
    tanLabCM5->SetLineColor(2);
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM5,text,"l");
    leg->Draw();
    
    TF1 *tanLabCM6 = tanLabCM->DrawCopy("same");
    betastar = -0.95;
    tanLabCM6->SetParameters(betastar,betabCM);
    tanLabCM6->SetParNames("betastar","betabCM");
    tanLabCM6->SetLineColor(4);
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM6,text,"l");
    leg->Draw();
    
    TF1 *tanLabCM7 = tanLabCM->DrawCopy("same");
    betastar = -0.99;
    tanLabCM7->SetParameters(betastar,betabCM);
    tanLabCM7->SetParNames("betastar","betabCM");
    tanLabCM7->SetLineColor(1);
    text.Form("#beta=%.7f",abs(betastar));
    leg->AddEntry(tanLabCM7,text,"l");
    leg->Draw();
    
    // Define histograms
    
    nbins=150;
    TH1D *h1_costhestar = new TH1D ("h1_costhestar","Cos(theta_CM)",nbins,-1.5,1.5);
    TH1D *h1_costhestar_eff = new TH1D ("h1_costhestar_eff","Efficiency vs Cos(theta_CM)",nbins,-1.5,1.5);   // must have same binning as previous
    
    TH1D *h1_costhe = new TH1D ("h1_costhe","Cos(theta_Lab)",nbins,-1.5,1.5);
    TH1D *h1_DeltaCos = new TH1D ("h1_DeltaCos","Costhe - CostheMeas",nbins,-1.5,1.5);
    TH1D *h1_costhestar_meas = new TH1D ("h1_costhestar_meas","Cos(theta_CM) Measured",nbins,-1.5,1.5);
    TH2D *h2_Cos_measVSCos = new TH2D("h2_Cos_measVSCos","Cos Meas vs Cos",nbins,-1.5,1.5,nbins,-1.5,1.5);
    TH1D *h1_E1 = new TH1D("h1_E1","E1 (GeV)",100,0,2);
    TH1D *h1_E1meas = new TH1D("h1_E1meas","E1 meas (GeV)",100,0,2);
    TH1D *h1_E2meas = new TH1D("h1_E2meas","E2 meas (GeV)",100,0,2);
    TH1D *h1_Emeas = new TH1D("h1_Emeas","E1+E2 meas (GeV)",200,0,4);
    
    TH2D *h2_costheVScosthestar = new TH2D ("h2_costheVScosthestar","Cos(theta_Lab) vs Cos(theta_CM)",100,-1,1,100,-1,1);
    TH2D *h2_E2VSE1 = new TH2D ("h2_E2VSE1","E2 (GeV) vs E1 (GeV)",100,0,2,100,0,2);
    TH2D *h2_E2measVSE1meas = new TH2D ("h2_E2measVSE1meas","E2meas (GeV) vs E1meas (GeV)",100,0,2,100,0,2);
    // generate photon angular distribution in the CM, then boost to the lab.
    
    
    Double_t Estar = Ma/2;
    Double_t thetamin=0;
    Double_t thetamax=180;
    Double_t cosmax=cos(thetamin*PI/180.);
    Double_t cosmin=cos(thetamax*PI/180.);
    cout << " thetamin=" << thetamin << " thetama=" << thetamax << " cosmin=" << cosmin << " cosmax=" << cosmax << endl;

    
    for (Int_t jj=0; jj<npts; jj++) {
        Double_t EaLab = EaLab0 + (2*r4->Uniform()-1)*DeltaE/2;     // pick random Energy inside uniform bin +/- nsigma x sigma
        Double_t paLab = sqrt(EaLab*EaLab - Ma*Ma);
        TLorentzVector pi0lab (0,0,paLab,EaLab);
        Double_t costhestar = cosmin + (cosmax-cosmin)*r1->Uniform();
        Double_t thetastar = acos(costhestar)*180./PI;
        Double_t phistar = 360*r3->Uniform();
        Double_t sinthestar = sin(thetastar*PI/180.);
        Double_t pxstar = Estar*sinthestar*cos(phistar*PI/180.);
        Double_t pystar = Estar*sinthestar*sin(phistar*PI/180.);
        Double_t pzstar = Estar*costhestar;
        TLorentzVector photonstar1(pxstar,pystar,pzstar,Estar);
        TLorentzVector photonstar2(-pxstar,-pystar,-pzstar,Estar);
        TLorentzRotation Boost(pi0lab.BoostVector() );
        TLorentzVector photon1 = Boost * photonstar1;
        TLorentzVector photon2 = Boost * photonstar2;
        Double_t costhe = cos(photon1.Theta());
        Double_t E1 = photon1.E();
        Double_t E2 = photon2.E();
        
        // check equations:
        // Double_t E2calc = 0.5*(EaLab - paLab*costhestar);
        // cout << " costhestar=" << costhestar << " E2=" << E2 << " E2calc=" << E2calc << endl;
        
        Double_t sigma = E1*sqrt(BCALa*BCALa/E1+BCALb*BCALb);
        Double_t E1meas = r2->Gaus(E1,sigma);
        sigma = E2*sqrt(BCALa*BCALa/E2+BCALb*BCALb);
        Double_t E2meas = r2->Gaus(E2,sigma);
        
        
        if (E1meas < E2meas) {    // always take E1meas to have the higher energy. Switch both gen and meas energies
            Double_t Etemp = E1;
            E1 = E2;
            E2 = Etemp;
            Etemp = E1meas;
            E1meas = E2meas;
            E2meas = Etemp;
            costhestar = -costhestar;
            costhe = cos(photon2.Theta());
            
        }
        
        if (E2meas < Ethresh) continue;     // impose threshold on smallest energy
        if (r4->Uniform() > eff->Eval(E2meas)) continue;   // efficiency for detecting E2. Assume efficiency for E1 is 100%;
        h1_E1->Fill(E1);
        h1_E1meas->Fill(E1meas);
        h1_Emeas->Fill(E1meas+E2meas);
        h2_E2VSE1->Fill(E1,E2);
        h1_costhestar->Fill(costhestar);
        h1_costhe->Fill(costhe);
        h2_costheVScosthestar->Fill(costhestar,costhe);
        h2_E2measVSE1meas->Fill(E1meas,E2meas);
        
        Double_t costhestar2;
        if (tag == "palab") {
    		costhestar2 = (EaLab0 -2*E2meas)/paLab0;
        }
        else if (tag == "energy") {
        	costhestar2 = (E1meas-E2meas)/sqrt( (E1meas+E2meas)*(E1meas+E2meas) - Ma*Ma);
        }
        else {
            cout << "*** Unknown tag=" << tag << endl;
            return;
        }
        if (costhestar2 > 1) {
            cout << "*** costhestar2=" << costhestar2 << " costhestar=" << costhestar << " E1meas=" << E1meas << " E2meas=" << E2meas << endl;
        }
        // if (costhestar2 > 1) costhestar2 = 1;     // impose physicsl limits
        // if (costhestar2 < -1) costhestar2 = -1;
        h1_costhestar_meas->Fill(costhestar2);
        h1_DeltaCos->Fill(costhestar2 - costhestar);
        h2_Cos_measVSCos->Fill(costhestar,costhestar2);
        // cout << " costhestar=" << costhestar << " costhestar2=" << costhestar2 << endl << endl;
        
        // cout << " jj=" << jj << " costhestar=" << costhestar << endl;
    }
    
    
    TCanvas *c2 = new TCanvas("c2", "c2",200,10,1000,1000);
    c2->Divide(2,3);
    
    c2->cd(1);
    ymin=0.1;
    ymax=h1_costhestar->GetMaximum()*1.5;
    gPad->SetLogy();
    text.Form("tag=%s, p_{#pi}=%.2f, E_{thresh}=%.2f, a=%.3f, b=%.3f, #DeltaE=%.3f",tag.Data(),paLab0,Ethresh,BCALa,BCALb,DeltaE);
    h1_costhestar->SetTitle(text);
    // h1_costhestar->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhestar->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhestar->GetXaxis()->SetTitleSize(0.05);
    h1_costhestar->GetYaxis()->SetTitleSize(0.05);
    h1_costhestar->GetXaxis()->SetTitle("Cos(theta_CM)");
    h1_costhestar->GetYaxis()->SetTitle("Events");
    h1_costhestar->SetLineColor(1);
    h1_costhestar->Draw();
    
    ymin=0.1;
    ymax=h1_costhe->GetMaximum()*1.5;
    c2->cd(2);
    h1_costhe->SetTitle(text);
    gPad->SetLogy();
    // h1_costhe->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhe->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhe->GetXaxis()->SetTitleSize(0.05);
    h1_costhe->GetYaxis()->SetTitleSize(0.05);
    h1_costhe->GetXaxis()->SetTitle("Cos(theta_Lab)");
    h1_costhe->GetYaxis()->SetTitle("Events");
    h1_costhe->SetLineColor(1);
    h1_costhe->Draw();
    
    c2->cd(3);
    h2_costheVScosthestar->SetTitle(text);
    // h2_costheVScosthestar->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_costheVScosthestar->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_costheVScosthestar->GetXaxis()->SetTitleSize(0.05);
    h2_costheVScosthestar->GetYaxis()->SetTitleSize(0.05);
    h2_costheVScosthestar->GetXaxis()->SetTitle("Cos(theta_CM)");
    h2_costheVScosthestar->GetYaxis()->SetTitle("Cos(theta_Lab)");
    h2_costheVScosthestar->SetLineColor(1);
    h2_costheVScosthestar->Draw("colz");
    
    c2->cd(4);
    h2_E2VSE1->SetTitle(text);
    // h2_E2VSE1->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_E2VSE1->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_E2VSE1->GetXaxis()->SetTitleSize(0.05);
    h2_E2VSE1->GetYaxis()->SetTitleSize(0.05);
    h2_E2VSE1->GetXaxis()->SetTitle("E1 (GeV)");
    h2_E2VSE1->GetYaxis()->SetTitle("E2 (GeV)");
    h2_E2VSE1->SetLineColor(1);
    h2_E2VSE1->Draw("colz");
    
    c2->cd(5);
    h2_E2measVSE1meas->SetTitle(text);
    // h2_E2measVSE1meas->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_E2measVSE1meas->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_E2measVSE1meas->GetXaxis()->SetTitleSize(0.05);
    h2_E2measVSE1meas->GetYaxis()->SetTitleSize(0.05);
    h2_E2measVSE1meas->GetXaxis()->SetTitle("E1meas (GeV)");
    h2_E2measVSE1meas->GetYaxis()->SetTitle("E2meas (GeV)");
    h2_E2measVSE1meas->SetLineColor(1);
    h2_E2measVSE1meas->Draw("colz");
                          
    c2->cd(6);
    h1_DeltaCos->SetTitle(text);
	// h1_E2measVSE1meas->GetXaxis()->SetRangeUser(xmin,xmax);
    // h1_DeltaCos->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_DeltaCos->GetXaxis()->SetTitleSize(0.05);
    h1_DeltaCos->GetYaxis()->SetTitleSize(0.05);
    h1_DeltaCos->GetXaxis()->SetTitle("CosThe_meas - CosThe");
    h1_DeltaCos->GetYaxis()->SetTitle("Events");
    h1_DeltaCos->SetLineColor(1);
    h1_DeltaCos->Fit("gaus");
    h1_DeltaCos->Draw();
    
    TCanvas *c3 = new TCanvas("c3", "c3",200,10,1000,1000);
    c3->Divide(2,3);
    
    ymin=0;
    ymax=h1_costhestar_meas->GetMaximum()*1.5;
    c3->cd(1);
    h1_costhestar_meas->SetTitle(text);
    // h1_costhestar_meas->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhestar_meas->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhestar_meas->GetXaxis()->SetTitleSize(0.05);
    h1_costhestar_meas->GetYaxis()->SetTitleSize(0.05);
    h1_costhestar_meas->GetXaxis()->SetTitle("Cos(theta_CM) meas");
    h1_costhestar_meas->GetYaxis()->SetTitle("Events");
    h1_costhestar_meas->SetLineColor(1);
    h1_costhestar_meas->Draw();
    
    c3->cd(2);
    h2_Cos_measVSCos->SetTitle(text);
    // h2_Cos_measVSCos->GetXaxis()->SetRangeUser(xmin,xmax);
    // h2_Cos_measVSCos->GetYaxis()->SetRangeUser(ymin,ymax);
    h2_Cos_measVSCos->GetXaxis()->SetTitleSize(0.05);
    h2_Cos_measVSCos->GetYaxis()->SetTitleSize(0.05);
    h2_Cos_measVSCos->GetXaxis()->SetTitle("Cos(theta_CM)");
    h2_Cos_measVSCos->GetYaxis()->SetTitle("Cos(theta_CM) meas");
    h2_Cos_measVSCos->SetLineColor(1);
    h2_Cos_measVSCos->Draw("colz");
    
    ymin=0;
    ymax=h1_E1->GetMaximum()*1.5;
    c3->cd(3);
    h1_E1->SetTitle(text);
    // h1_E1->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_E1->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_E1->GetXaxis()->SetTitleSize(0.05);
    h1_E1->GetYaxis()->SetTitleSize(0.05);
    h1_E1->GetXaxis()->SetTitle("E1 (GeV)");
    h1_E1->GetYaxis()->SetTitle("Events");
    h1_E1->SetLineColor(1);
    h1_E1->Draw();
    
    ymin=0;
    ymax=h1_E1meas->GetMaximum()*1.5;
    c3->cd(4);
    h1_E1meas->SetTitle(text);
    // h1_E1meas->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_E1meas->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_E1meas->GetXaxis()->SetTitleSize(0.05);
    h1_E1meas->GetYaxis()->SetTitleSize(0.05);
    h1_E1meas->GetXaxis()->SetTitle("E1meas (GeV)");
    h1_E1meas->GetYaxis()->SetTitle("Events");
    h1_E1meas->SetLineColor(1);
    h1_E1meas->Draw();
    
    // unpack costhestar distribution and intepret it as an efficiency distribution as a function of E2meas
    
    const Int_t ndim = h1_costhestar_meas->GetNbinsX();
    Double_t xlo = h1_costhestar_meas->GetBinLowEdge(1);
    Double_t width = h1_costhestar_meas->GetBinWidth(1);
    Int_t max_bin = h1_costhestar_meas->GetMaximumBin();
    Double_t max_content = h1_costhestar_meas->GetBinContent(max_bin);     // use as normalization for efficiency
    
    Double_t E2calc[ndim];
    Double_t E2calc_err[ndim];
    Double_t E2calc_eff[ndim];
    Double_t E2calc_efferr[ndim];
    Int_t npoints=0;
    
    for(Int_t j=0;j<ndim;j++) {
        Double_t costhestar_bin = h1_costhestar_meas->GetBinCenter(j+1);
        Double_t content = h1_costhestar_meas->GetBinContent(j+1);
        Double_t efficiency = max_content > 0? content/max_content: 0;
        // cout << " j+1=" << j+1 << " costhestar_bin=" << costhestar_bin << " content=" << content << " maxbin=" << max_bin << " Max=" << max_content << " efficiency=" << efficiency << endl;
        h1_costhestar_eff->SetBinContent(j+1,efficiency);
        // Double_t sigma = sqrt (content*efficiency*(1-efficiency))/max_content;   // assume no uncertainty in maximum  // binomial statistics
        Double_t sigma = sqrt (content)/max_content;   //  just take uncertainty in point, ignore fluctuation in maximum
        h1_costhestar_eff->SetBinError(j+1,sigma);
        // cout << " j=" << j << " EaLab=" << EaLab << " Emeas=" << Emeas << endl;
        
        if (efficiency > 0) {
            E2calc[npoints] = 0.5*(EaLab0 - paLab0*costhestar_bin);      // use nominal pi0 energy in "reconstruction"
            E2calc_err[npoints] = 0.005;
        	E2calc_eff[npoints] = efficiency/0.96;     // normalize to average value of plateau for 100,000 events
            // E2calc_eff[npoints] = efficiency;     // normalize to average value of plateau for 100,000 events
            E2calc_efferr[npoints] = sigma;
            npoints++;
        }
    }
    
    TGraphErrors *gr_E2calc = new TGraphErrors(npoints,E2calc,E2calc_eff,E2calc_err,E2calc_efferr);
    
    ymin=0;
    ymax=h1_costhestar_eff->GetMaximum()*1.5;
    c3->cd(5);
    h1_costhestar_eff->SetTitle(text);
    // h1_costhestar_eff->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_costhestar_eff->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_costhestar_eff->GetXaxis()->SetTitleSize(0.05);
    h1_costhestar_eff->GetYaxis()->SetTitleSize(0.05);
    h1_costhestar_eff->GetXaxis()->SetTitle("Cos(theta_CM)");
    h1_costhestar_eff->GetYaxis()->SetTitle("Efficiency");
    h1_costhestar_eff->SetLineColor(2);
    h1_costhestar_eff->SetMarkerColor(2);
    h1_costhestar_eff->SetMarkerStyle(20);
    h1_costhestar_eff->SetMarkerSize(0.3);
    h1_costhestar_eff->Draw("pE");
    
    ymin=0;
    ymax=h1_Emeas->GetMaximum()*1.5;
    c3->cd(6);
    h1_Emeas->SetTitle(text);
    // h1_Emeas->GetXaxis()->SetRangeUser(xmin,xmax);
    h1_Emeas->GetYaxis()->SetRangeUser(ymin,ymax);
    h1_Emeas->GetXaxis()->SetTitleSize(0.05);
    h1_Emeas->GetYaxis()->SetTitleSize(0.05);
    h1_Emeas->GetXaxis()->SetTitle("Emeas (GeV)");
    h1_Emeas->GetYaxis()->SetTitle("Events");
    h1_Emeas->SetLineColor(1);
    h1_Emeas->Fit("gaus");
    h1_Emeas->Draw();
    
    TCanvas *c4 = new TCanvas("c4", "c4",200,10,1000,700);    //  output summary plot to separate page.
    gPad->SetGridx();
    gPad->SetGridy();
    
    xmin=0.;
    xmax=1.5;
    ymin=0;
    ymax=1.2;
    gr_E2calc->SetTitle(text);
    gr_E2calc->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_E2calc->GetYaxis()->SetRangeUser(ymin,ymax);
    gr_E2calc->GetXaxis()->SetTitleSize(0.05);
    gr_E2calc->GetYaxis()->SetTitleSize(0.05);
    gr_E2calc->GetXaxis()->SetTitle("E2calc (GeV)");
    gr_E2calc->GetYaxis()->SetTitle("Efficiency");
    gr_E2calc->SetLineColor(2);
    gr_E2calc->SetMarkerColor(2);
    gr_E2calc->SetMarkerStyle(20);
    gr_E2calc->SetMarkerSize(0.3);
    gr_E2calc->Draw("Ap");
    eff->Draw("same");
    
    TString outfile;
    outfile.Form("AnglesCMLab_eff_%s_%.0f_%.0f_%.0f_%.0f_%.0f",tag.Data(),1000*paLab0,1000*Ethresh,1000*BCALa,1000*BCALb,1000*DeltaE);
    
    c1->SaveAs(outfile+".pdf(");
    c2->SaveAs(outfile+".pdf");
    c3->SaveAs(outfile+".pdf");
    c4->SaveAs(outfile+".pdf)");
    
}