#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#define GetFakeRate_Kaon_C
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooHistPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"
#include <string>
#include <cmath>
#include <TH2.h>
#include <TStyle.h>
#include "GetFakeRate_Kaon.h"
using namespace RooFit ;

void GetFakeRate_Kaon::Loop()
{
	bool isMC = false; 
	 TFile file("hists_data17.root","RECREATE");

////Create and Fill Histograms

  TH1D *kmass = new TH1D("kmass","kaon mass with fit",100 , 400, 600);
  kmass->Sumw2();

  TH1D *npix_MCreal = new TH1D("npix_MCreal","truth matched number of real pixel hits",10,-0.5,9.5);
  npix_MCreal->Sumw2();
  TH1D *npix_MCfake = new TH1D("npix_MCfake","truth matched number of fake pixel hits",10,-0.5,9.5);
  npix_MCfake->Sumw2();
  TH1D *npix_side = new TH1D("npix_side","number of pixel hits in  sideband region",10,-0.5,9.5);
  npix_side->Sumw2();
  TH1D *npix_sig = new TH1D("npix_sig","number of pixel hits in signal region",10,-0.5,9.5);
  npix_sig->Sumw2();
  TH1D *npix_sig_test = new TH1D("npix_sig_test","number of pixel hits in signal test region",10,-0.5,9.5);
  npix_sig_test->Sumw2();
  TH1D *npix_side_test = new TH1D("npix_side_test","number of pixel hits in sideband test region",10,-0.5,9.5);
  npix_side_test->Sumw2();

  TH1D *nSCT_MCreal = new TH1D("nSCT_MCreal","truth matched number of real SCT hits",18,-0.5,17.5);
  nSCT_MCreal->Sumw2();
  TH1D *nSCT_MCfake = new TH1D("nSCT_MCfake","truth matched number of fake SCT hits",18,-0.5,17.5);
  nSCT_MCfake->Sumw2();
  TH1D *nSCT_side = new TH1D("nSCT_side","number of SCT hits in  sideband region",18,-0.5,17.5);
  nSCT_side->Sumw2();
  TH1D *nSCT_sig = new TH1D("nSCT_sig","number of SCT hits in signal region",18,-0.5,17.5);
  nSCT_sig->Sumw2();
  TH1D *nSCT_sig_test = new TH1D("nSCT_sig_test","number of SCT hits in signal test region",18,-0.5,17.5);
  nSCT_sig_test->Sumw2();
  TH1D *nSCT_side_test = new TH1D("nSCT_side_test","number of SCT hits in sideband test region",18,-0.5,17.5);
  nSCT_side_test->Sumw2();

  TH1D *nTRT_MCreal = new TH1D("nTRT_MCreal","truth matched number of real TRT hits",70,-0.5,69.5);
  nTRT_MCreal->Sumw2();
  TH1D *nTRT_MCfake = new TH1D("nTRT_MCfake","truth matched number of fake TRT hits",70,-0.5,69.5);
  nTRT_MCfake->Sumw2();
  TH1D *nTRT_side = new TH1D("nTRT_side","number of TRT hits in  sideband region",70,-0.5,69.5);
  nTRT_side->Sumw2();
  TH1D *nTRT_sig = new TH1D("nTRT_sig","number of TRT hits in signal region",70,-0.5,69.5);
  nTRT_sig->Sumw2();
  TH1D *nTRT_sig_test = new TH1D("nTRT_sig_test","number of TRT hits in signal test region",70,-0.5,69.5);
  nTRT_sig_test->Sumw2();
  TH1D *nTRT_side_test = new TH1D("nTRT_side_test","number of TRT hits in sideband test region",70,-0.5,69.5);
  nTRT_side_test->Sumw2();

  TH1D *prob = new TH1D("prob","truth match probability",100,0,1);
  prob->Sumw2();
  TH1D *prob_sig = new TH1D("prob_sig","truth match probability - signal region",100,0,1);
  prob_sig->Sumw2();
  TH1D *prob_sideband = new TH1D("prob_sideband","truth match probability - sideband region",100,0,1);
  prob_sideband->Sumw2();

///////
//CUTS
///////
double_t DeltaRmax = 0.4; //cut out above this
double_t d0 = .15; //cut out below this
double_t xylengthmin = 2.0;
double_t xylengthmax = 3.0;
double_t xylengthmax_endcap = 30.0;
double_t cosThstcutmax = 0.92; //cut out over this
double_t kpTcutmax = 20000; //cut out over this
double_t kpTcutmin = 3000; //cut out below this
double_t eta_barrel_max = 2.2;
double_t eta_endcap_min = 1.6;
double_t eta_endcap_max = 2.7;

// Signal and sideband region definitions:
double_t mean_guess = 499; //from fit
double_t sigma_guess = 11; //from fit

//Definition of real and fake truth match criteria:
double_t cutoff_real = 0.5;
double_t cutoff_fake = 0.5;

int rrand;

   Long64_t nentries = fChain->GetEntriesFast();

    fChain->SetBranchStatus("*",0);  // disable all branches
    fChain->SetBranchStatus("kaon_pi2Npixelhits",1);  // activate branchname
    fChain->SetBranchStatus("kaon_pi2TruthMatchProb",1);  // activate branchname
    fChain->SetBranchStatus("kaon_pi1Npixelhits",1);  // activate branchname
    fChain->SetBranchStatus("kaon_pi1TruthMatchProb",1);  // activate branchname
    fChain->SetBranchStatus("Ntracks",1);  // activate branchname
    fChain->SetBranchStatus("kaon_mpipi",1);  // activate branchname
    fChain->SetBranchStatus("kaon_pipiDeltaR",1);
    fChain->SetBranchStatus("kaon_pi1d0",1);
    fChain->SetBranchStatus("kaon_pi2d0",1);
    fChain->SetBranchStatus("kaon_xylength",1);
    fChain->SetBranchStatus("kaon_cosThetaStar",1);
    fChain->SetBranchStatus("kaon_pipivertex_x",1);
    fChain->SetBranchStatus("kaon_pipivertex_y",1);
    fChain->SetBranchStatus("kaon_pipivertex_z",1);
    fChain->SetBranchStatus("kaon_pipipT",1);
    fChain->SetBranchStatus("kaon_pi1pT",1);
    fChain->SetBranchStatus("kaon_pi2pT",1);
    fChain->SetBranchStatus("kaon_pi1eta",1);
    fChain->SetBranchStatus("kaon_pi2eta",1);
    fChain->SetBranchStatus("kaon_pi1NSCThits",1);
    fChain->SetBranchStatus("kaon_pi2NSCThits",1);
    fChain->SetBranchStatus("kaon_pi1NTRThits",1);
    fChain->SetBranchStatus("kaon_pi2NTRThits",1);

   Long64_t nbytes = 0, nb = 0;

	for (Long64_t jentry=0; jentry<nentries;jentry++) {
	Long64_t ientry = LoadTree(jentry);
	if (ientry < 0) break;
	nb = fChain->GetEntry(jentry);   nbytes += nb;

	b_kaon_pi2Npixelhits->GetEntry(ientry);
        b_kaon_pi2TruthMatchProb->GetEntry(ientry);
        b_kaon_pi1Npixelhits->GetEntry(ientry);
        b_kaon_pi1TruthMatchProb->GetEntry(ientry);
        b_Ntracks->GetEntry(ientry);
        b_kaon_mpipi->GetEntry(ientry);
	b_kaon_pipiDeltaR->GetEntry(ientry);
	b_kaon_pi1d0->GetEntry(ientry);
	b_kaon_pi2d0->GetEntry(ientry);
	b_kaon_cosThetaStar->GetEntry(ientry);
	b_kaon_pipivertex_x->GetEntry(ientry);
	b_kaon_pipivertex_y->GetEntry(ientry);
     	b_kaon_pipipT->GetEntry(ientry);
      	b_kaon_pi1pT->GetEntry(ientry);
      	b_kaon_pi2pT->GetEntry(ientry);
      	b_kaon_pi1eta->GetEntry(ientry);
      	b_kaon_pi2eta->GetEntry(ientry);
        b_kaon_pi1NSCThits->GetEntry(ientry);
        b_kaon_pi2NSCThits->GetEntry(ientry);
        b_kaon_pi1NTRThits->GetEntry(ientry);
        b_kaon_pi2NTRThits->GetEntry(ientry);

		for (int i=0; i<Ntracks; i++){
                rrand = rand() % 2;

	        double_t kpT = (*kaon_pipipT)[i];
	        double_t pi1pT = (*kaon_pi1pT)[i];
	        double_t pi2pT = (*kaon_pi2pT)[i];
	        double_t pi1d0 = (*kaon_pi1d0)[i];
	        double_t pi2d0 = (*kaon_pi2d0)[i];
	        double_t xylength = (*kaon_xylength)[i];
	        double_t cosThetaStar = (*kaon_cosThetaStar)[i];
	        double_t pipiDeltaR = (*kaon_pipiDeltaR)[i];
	        double_t pi1eta = (*kaon_pi1eta)[i];
	        double_t pi2eta = (*kaon_pi2eta)[i];
	        double_t kaonmass = (*kaon_mpipi)[i];
	        double_t npix1 = (*kaon_pi1Npixelhits)[i];
	        double_t npix2 = (*kaon_pi2Npixelhits)[i];
	        double_t truth1 = (*kaon_pi1TruthMatchProb)[i];
	        double_t truth2 = (*kaon_pi2TruthMatchProb)[i];
	        double_t nSCT1 = (*kaon_pi1NSCThits)[i];
	        double_t nSCT2 = (*kaon_pi2NSCThits)[i];
	        double_t nTRT1 = (*kaon_pi1NTRThits)[i];
	        double_t nTRT2 = (*kaon_pi2NTRThits)[i];

	        bool passThetaStarcut = true;
	        bool passDeltaRcut = true;
	        bool passpi1d0cut = true;
	        bool passpi2d0cut = true;
	        bool passxylengthcut = false;
	        bool passkpTcut = true;
	        bool pi1barrel = true;
	        bool pi2barrel = true;

        if (abs(pi1eta)>2.2) pi1barrel=false;
        if (abs(pi2eta)>2.2) pi2barrel=false;
        if (abs(cosThetaStar) > cosThstcutmax) passThetaStarcut = false; //Cos(ThetaStar) is the angle btwn vertex and (pion BOOSTED in k0 frame) so backgrounds will be close to 1
        if (pipiDeltaR > DeltaRmax) passDeltaRcut = false;
        if (abs(pi1d0)<d0) passpi1d0cut = false;
        if (abs(pi2d0)<d0) passpi2d0cut = false;
        if (pi1barrel && pi2barrel && xylength>xylengthmin && xylength<xylengthmax) passxylengthcut = true;
        if (kpT > kpTcutmax || kpT <  kpTcutmin) passkpTcut = false;
        if (passThetaStarcut && passDeltaRcut && passpi1d0cut && passpi2d0cut && passxylengthcut && passkpTcut) {
        if(kaonmass<400 || kaonmass>600) continue;

		kmass->Fill(kaonmass);

	if(rrand==0){ // fill build hists

		if(truth1<cutoff_fake) {
			npix_MCfake->Fill(npix1);
                        npix_MCfake->Fill(npix2);
                        nSCT_MCfake->Fill(nSCT1);
                        nSCT_MCfake->Fill(nSCT2);
                        nTRT_MCfake->Fill(nTRT1);
                        nTRT_MCfake->Fill(nTRT2);
			}
                if(truth1>cutoff_real) {
			npix_MCreal->Fill(npix1);
                        npix_MCreal->Fill(npix2);
                        nSCT_MCreal->Fill(nSCT1);
                        nSCT_MCreal->Fill(nSCT2);
                        nTRT_MCreal->Fill(nTRT1);
                        nTRT_MCreal->Fill(nTRT2);
			}
		if(kaonmass<(mean_guess-3*sigma_guess) || kaonmass>(mean_guess+3*sigma_guess)) {
			npix_side->Fill(npix1);
                        npix_side->Fill(npix2);
                        nSCT_side->Fill(nSCT1);
                        nSCT_side->Fill(nSCT2);
                        nTRT_side->Fill(nTRT1);
                        nTRT_side->Fill(nTRT2);
			}
		if(kaonmass>(mean_guess-2*sigma_guess) && kaonmass<(mean_guess+2*sigma_guess)) {
			npix_sig->Fill(npix1);
                        npix_sig->Fill(npix2);
                        nSCT_sig->Fill(nSCT1);
                        nSCT_sig->Fill(nSCT2);
                        nTRT_sig->Fill(nTRT1);
                        nTRT_sig->Fill(nTRT2);
			}

		}//fill build hists

	if(rrand==1){ //fill test hists

		if(kaonmass<(mean_guess-3*sigma_guess) || kaonmass>(mean_guess+3*sigma_guess)) {
                        npix_side_test->Fill(npix1);
                        npix_side_test->Fill(npix2);
                        nSCT_side_test->Fill(nSCT1);
                        nSCT_side_test->Fill(nSCT2);
                        nTRT_side_test->Fill(nTRT1);
                        nTRT_side_test->Fill(nTRT2);
                        }
                if(kaonmass>(mean_guess-2*sigma_guess) && kaonmass<(mean_guess+2*sigma_guess)) {
                        npix_sig_test->Fill(npix1);
                        npix_sig_test->Fill(npix2);
                        nSCT_sig_test->Fill(nSCT1);
                        nSCT_sig_test->Fill(nSCT2);
                        nTRT_sig_test->Fill(nTRT1);
                        nTRT_sig_test->Fill(nTRT2);
                        }
		}//fill test hists

        prob->Fill(truth1);
        prob->Fill(truth2);
        if(kaonmass>(mean_guess-2*sigma_guess) && kaonmass<(mean_guess+2*sigma_guess)) prob_sig->Fill(truth1);
        if(kaonmass>(mean_guess-2*sigma_guess) && kaonmass<(mean_guess+2*sigma_guess)) prob_sig->Fill(truth2);
        if(kaonmass<(mean_guess-3*sigma_guess) || kaonmass>(mean_guess+3*sigma_guess)) prob_sideband->Fill(truth1);
        if(kaonmass<(mean_guess-3*sigma_guess) || kaonmass>(mean_guess+3*sigma_guess)) prob_sideband->Fill(truth2);

		} //if pass cuts statement
		} //loop over i
	} //loop over jentry

	double_t nevents_signal_region = npix_sig->Integral();
	double_t nevents_sideband_region = npix_side->Integral();
        double_t nevents_signal_region_test = npix_sig_test->Integral();
        double_t nevents_sideband_region_test = npix_side_test->Integral();
cout << "nevents_signal_region " << nevents_signal_region << endl;
cout << "nevents_sideband_region " << nevents_sideband_region << endl;
cout << "nevents_signal_region_test " << nevents_signal_region_test << endl;
cout << "nevents_sideband_region_test " << nevents_sideband_region_test << endl;

////Create Kaon Mass Plot and Fit

  //create roofit variable
        RooRealVar x("x","x",400,600) ;
  //create roofit gaussian
        RooRealVar mean("mean","mean",498,490,510) ;
        RooRealVar sigma("sigma","sigma",9,7,11) ;
        RooRealVar sig_yield("sig_yield","signal yield",10000,0,3000000);
        RooGaussian gauss("gauss","signal p.d.f.",x,mean,sigma) ;
  //create roofit quadratic
        RooRealVar c1("c1","coefficient of x term",-10, -400, -1) ;
        RooRealVar c2("c2","coefficient of x^2 term",10, 1, 30) ;
        RooRealVar c3("c3","coefficient of x^3 term",-0.01, -1, -0.0001) ;
        RooRealVar bkg_yield("bkg_yield","background yield",500000,0,4000000);
        RooPolynomial bkg("bkg","background", x, RooArgList(c1,c2,c3)) ;
  //create composite model model(x) = sig_yield*gauss(x) + bkg_yield*bkg(x)
        RooAddPdf model("model","model",RooArgList(gauss,bkg),RooArgList(sig_yield,bkg_yield)) ;
  //imposrt data from hist
        RooDataHist data("data","dataset with K0 mass",x,kmass) ;
  //fit model to data
        RooFitResult* r_model = model.fitTo(data,Extended(),Save());
  //roofit's version of a canvas:
        RooPlot* xframe = x.frame(Title("Kaon Mass with Gaussian Fit and Quadradic Background")) ;
        xframe->SetName("kaonmass_fit");
        xframe->SetXTitle("Kaon Mass (GeV)");
        xframe->SetYTitle("Number of Events");
  //plot
        data.plotOn(xframe) ;
        model.plotOn(xframe) ;
        model.plotOn(xframe,Components(bkg),LineColor(kGreen));

 //create integrals
        x.setRange("signal",mean.getVal()-2*sigma.getVal(),mean.getVal()+2*sigma.getVal());
        RooAbsReal* fsigregion_model = model.createIntegral(x,NormSet(x),Range("signal")); //this is the fraction of total events that are in the signal region
        RooAbsReal* fsigregion_bkg = bkg.createIntegral(x,NormSet(x),Range("signal")); //this is the fraction of background events that are in the signal region
	//the number of event in the signal region that are actually signal (not background):
	Double_t nsigevents = fsigregion_model->getVal()*(sig_yield.getVal()+bkg_yield.getVal())-fsigregion_bkg->getVal()*bkg_yield.getVal();
	//calc error:
	Double_t error_nsigevents = sqrt(
                                                pow(fsigregion_model->getPropagatedError(*r_model)/fsigregion_model->getVal(),2)
                                                +(pow(sig_yield.getError(),2)+pow(bkg_yield.getError(),2))/pow(sig_yield.getVal()+bkg_yield.getVal(),2) 
                                                + pow(fsigregion_bkg->getPropagatedError(*r_model)/fsigregion_bkg->getVal(),2) 
                                                + pow(bkg_yield.getError()/bkg_yield.getVal(),2)
                                         	)*nsigevents;

	//the fraction of events in the signal region that are actually signal (not background):
	Double_t fsig = nsigevents/(fsigregion_model->getVal()*(sig_yield.getVal()+bkg_yield.getVal()));
        //calc error:
	Double_t error_fsig = sqrt(
                                        pow(error_nsigevents/nsigevents,2)
                                        +pow(fsigregion_model->getPropagatedError(*r_model)/fsigregion_model->getVal(),2)
                                        +(pow(sig_yield.getError(),2)+pow(bkg_yield.getError(),2))/pow(sig_yield.getVal()+bkg_yield.getVal(),2)
                                        )*fsig;

cout << "fsigregion_model = " << fsigregion_model->getVal() << " +/- " << fsigregion_model->getPropagatedError(*r_model) << endl;
cout << "fsigregion_bkg = " << fsigregion_bkg->getVal() << " +/- " << fsigregion_bkg->getPropagatedError(*r_model) << endl;
cout << "nsigevents = " << nsigevents << " +/- " << error_nsigevents << endl;
cout << "fsig = " << fsig << " +/- " << error_fsig << endl;

//make ratio plot
        TCanvas *cratio = new TCanvas("cratio","cratio",800,800);
        TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
        pad1->SetBottomMargin(0); // Upper and lower plot are joined
        pad1->SetGridx();         // Vertical grid
        pad1->Draw();             // Draw the upper pad: pad1
        pad1->cd();               // pad1 becomes the current pad
        xframe->Draw();
        TLegend *l = new TLegend(0.1,0.6,0.4,0.9);
        std::string str2 = "fraction of in signal region = " +  std::to_string(fsig) + " +/- " + std::to_string(error_fsig);
        const char *chr2 = str2.c_str();
        l->AddEntry((TObject*)0,chr2,"");
        std::string str3 = "mean = " +  std::to_string(mean.getVal()) + " +/- " + std::to_string(mean.getError());
        const char *chr3 = str3.c_str();
        l->AddEntry((TObject*)0,chr3,"");
        std::string str4 = "sigma = " +  std::to_string(sigma.getVal()) + " +/- " + std::to_string(sigma.getError());
        const char *chr4 = str4.c_str();
        l->AddEntry((TObject*)0,chr4,"");
        std::string str5 = "num signal events = " +  std::to_string(sig_yield.getVal()) + " +/- " + std::to_string(sig_yield.getError());
        const char *chr5 = str5.c_str();
        l->AddEntry((TObject*)0,chr5,"");
        std::string str6 = "num background events = " +  std::to_string(bkg_yield.getVal()) + " +/- " + std::to_string(bkg_yield.getError());
        const char *chr6 = str6.c_str();
        l->AddEntry((TObject*)0,chr6,"");
        l->Draw();
        cratio->cd();          // Go back to the main canvas before defining pad2
        TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
        pad2->SetTopMargin(0);
        pad2->SetBottomMargin(0.2);
        pad2->SetGridx(); // vertical grid
        pad2->Draw();
        pad2->cd();       // pad2 becomes the current pad
        TH1F *fitratio = new TH1F("fitratio","fitratio",100,400,600);
        RooDataHist* data_from_model = model.generateBinned(x,0,kTRUE);
        TH1F *kmass_fromroomodel = (TH1F*) data_from_model->createHistogram("kmass_fromroomodel",x,Binning(100,400,600));
        TH1F *kmass_fromroodata = (TH1F*) data.createHistogram("kmass_fromroodata",x,Binning(100,400,600));
        kmass_fromroomodel->Sumw2();
        kmass_fromroodata->Sumw2();
        fitratio->Sumw2();
        fitratio->Divide(kmass_fromroomodel,kmass_fromroodata,1,1);
        fitratio->GetYaxis()->SetTitle("ratio fit/data");
        fitratio->SetStats(kFALSE);
        fitratio->SetTitle("");
        fitratio->GetYaxis()->SetTitleSize(20);
        fitratio->GetYaxis()->SetTitleFont(43);
        fitratio->GetYaxis()->SetTitleOffset(1.55);
        fitratio->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
        fitratio->GetYaxis()->SetLabelSize(15);
        fitratio->GetXaxis()->SetTitleSize(20);
        fitratio->GetXaxis()->SetTitleFont(43);
        fitratio->GetXaxis()->SetTitleOffset(1.55);
        fitratio->GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
        fitratio->GetXaxis()->SetLabelSize(15);
        fitratio->Draw();
        cratio->Print("PLOTS/kmass_ratio.pdf");


////Now for the analysis

//Create the real track distributions
        TH1D *npix_real_norm = new TH1D("npix_real_norm","number of real pixel hits (no truth info) in signal region - normalized",10,-0.5,9.5);
        npix_real_norm->Sumw2();
        npix_real_norm->Add(npix_sig,npix_side,1/nevents_signal_region,-1*(1-fsig)/nevents_sideband_region);
                for(int i=1;i<11;i++){
                if(npix_real_norm->GetBinContent(i)<0) npix_real_norm->SetBinContent(i,0);
                }
        double_t scale_npix_real_norm = npix_real_norm->Integral();
        npix_real_norm->Scale(1/scale_npix_real_norm);

        TH1D *nSCT_real_norm = new TH1D("nSCT_real_norm","number of real SCT hits (no truth info) in signal region - normalized",18,-0.5,17.5);
        nSCT_real_norm->Sumw2();
        nSCT_real_norm->Add(nSCT_sig,nSCT_side,1/nevents_signal_region,-1*(1-fsig)/nevents_sideband_region);
                for(int i=1;i<19;i++){
                if(nSCT_real_norm->GetBinContent(i)<0) nSCT_real_norm->SetBinContent(i,0);
                }
        double_t scale_nSCT_real_norm = nSCT_real_norm->Integral();
        nSCT_real_norm->Scale(1/scale_nSCT_real_norm);

        TH1D *nTRT_real_norm = new TH1D("nTRT_real_norm","number of real TRT hits (no truth info) in signal region - normalized",70,-0.5,69.5);
        nTRT_real_norm->Sumw2();
        nTRT_real_norm->Add(nTRT_sig,nTRT_side,1/nevents_signal_region,-1*(1-fsig)/nevents_sideband_region);
                for(int i=1;i<71;i++){
                if(nTRT_real_norm->GetBinContent(i)<0) nTRT_real_norm->SetBinContent(i,0);
                }
        double_t scale_nTRT_real_norm = nTRT_real_norm->Integral();
        nTRT_real_norm->Scale(1/scale_nTRT_real_norm);

//Create the fake track distributions, calculate and save chi-square values for each of the initial guesses of fake fraction

        TH1D *npix_fake_norm = new TH1D("npix_fake_norm","pixel hits data-driven fake track distribution - normalized",10,-0.5,9.5);
        npix_fake_norm->Sumw2();
        TH1D *npix_sig_pred = new TH1D("npix_sig_pred","npix_sig_pred",10,-0.5,9.5);
        npix_sig_pred->Sumw2();

        TH1D *nSCT_fake_norm = new TH1D("nSCT_fake_norm","SCT hits data-driven fake track distribution - normalized",18,-0.5,17.5);
        nSCT_fake_norm->Sumw2();
        TH1D *nSCT_sig_pred = new TH1D("nSCT_sig_pred","nSCT_sig_pred",18,-0.5,17.5);
        nSCT_sig_pred->Sumw2();

        TH1D *nTRT_fake_norm = new TH1D("nTRT_fake_norm","TRT hits data-driven fake track distribution - normalized",70,-0.5,69.5);
        nTRT_fake_norm->Sumw2();
        TH1D *nTRT_sig_pred = new TH1D("nTRT_sig_pred","nTRT_sig_pred",70,-0.5,69.5);
        nTRT_sig_pred->Sumw2();


	//Calculate vector magnitudes for BDM
	double_t mag_npix_sig_test = 0;
	for(int i=1;i<11;i++){
	mag_npix_sig_test = mag_npix_sig_test + pow(npix_sig_test->GetBinContent(i),2);
	}
	mag_npix_sig_test = sqrt(mag_npix_sig_test) ;

cout << "test unit vector " << endl;
	double_t shouldbe1 = 0;
	for(int i=1;i<11;i++){
	shouldbe1 = shouldbe1 + pow(npix_sig_test->GetBinContent(i)/mag_npix_sig_test,2);
	}
cout << "should be 1 => " << shouldbe1 << endl;

        double_t mag_nSCT_sig_test = 0;
        for(int i=1;i<19;i++){
        mag_nSCT_sig_test = mag_nSCT_sig_test + pow(nSCT_sig_test->GetBinContent(i),2);
        }
        mag_nSCT_sig_test = sqrt(mag_nSCT_sig_test) ;

        double_t mag_nTRT_sig_test = 0;
        for(int i=1;i<71;i++){
        mag_nTRT_sig_test = mag_nTRT_sig_test + pow(nTRT_sig_test->GetBinContent(i),2);
        }
        mag_nTRT_sig_test = sqrt(mag_nTRT_sig_test) ;



        double_t chi2_pix[500] = { };
        double_t chi2_SCT[500] = { };
        double_t chi2_TRT[500] = { };
        double_t delta2_pix[500] = { };
        double_t sigma2_pix[500] = { };
        double_t delta2_SCT[500] = { };
        double_t sigma2_SCT[500] = { };
        double_t delta2_TRT[500] = { };
        double_t sigma2_TRT[500] = { };
        double_t bdm_pix[500] = { };
        double_t bdm_SCT[500] = { };
        double_t bdm_TRT[500] = { };
        double_t mag_npix_sig_pred = 0;
        double_t mag_nSCT_sig_pred = 0;
        double_t mag_nTRT_sig_pred = 0;
        double_t cumulative1 = 0;
        double_t cumulative2 = 0;
        double_t f[500] = { };
        double_t ffake;
        double_t delta;
        double_t sigma_square;
        double_t nevents_sig_prediction_region;

    for(int a=1;a<501;a++){
        ffake = .002*a; //initial guess at fake fraction
        f[a-1] = ffake;

        cumulative1 = 0;
        cumulative2 = 0;

	//fake templates:
        npix_fake_norm->Add(npix_side,npix_real_norm,1/nevents_sideband_region,-1*(1-ffake));
                for(int i=1;i<11;i++){
                if(npix_fake_norm->GetBinContent(i)<0) npix_fake_norm->SetBinContent(i,0);
                }
        double_t scale_npix_fake_norm = npix_fake_norm->Integral();
        npix_fake_norm->Scale(1/scale_npix_fake_norm);

        nSCT_fake_norm->Add(nSCT_side,nSCT_real_norm,1/nevents_sideband_region,-1*(1-ffake));
                for(int i=1;i<19;i++){
                if(nSCT_fake_norm->GetBinContent(i)<0) nSCT_fake_norm->SetBinContent(i,0);
                }
        double_t scale_nSCT_fake_norm = nSCT_fake_norm->Integral();
        nSCT_fake_norm->Scale(1/scale_nSCT_fake_norm);

        nTRT_fake_norm->Add(nTRT_side,nTRT_real_norm,1/nevents_sideband_region,-1*(1-ffake));
                for(int i=1;i<71;i++){
                if(nTRT_fake_norm->GetBinContent(i)<0) nTRT_fake_norm->SetBinContent(i,0);
                }
        double_t scale_nTRT_fake_norm = nTRT_fake_norm->Integral();
        nTRT_fake_norm->Scale(1/scale_nTRT_fake_norm);

	//predictions for signal region:
	npix_sig_pred->Add(npix_real_norm,npix_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<11;i++){
                if(npix_sig_pred->GetBinContent(i)<0) npix_sig_pred->SetBinContent(i,0);
                }
        npix_sig_pred->Scale(nevents_signal_region/npix_sig_pred->Integral());
	
        nSCT_sig_pred->Add(nSCT_real_norm,nSCT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<19;i++){
                if(nSCT_sig_pred->GetBinContent(i)<0) nSCT_sig_pred->SetBinContent(i,0);
                }
        nSCT_sig_pred->Scale(nevents_signal_region/nSCT_sig_pred->Integral());

        nTRT_sig_pred->Add(nTRT_real_norm,nTRT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<71;i++){
                if(nTRT_sig_pred->GetBinContent(i)<0) nTRT_sig_pred->SetBinContent(i,0);
                }
        nTRT_sig_pred->Scale(nevents_signal_region/nTRT_sig_pred->Integral());

	//Calculate Chi^2
	for(int i=1;i<11;i++){ 
	if(npix_sig_test->GetBinContent(i)==0 || npix_sig_pred->GetBinContent(i)==0) continue;
	delta = npix_sig_test->GetBinContent(i)/nevents_signal_region_test - npix_sig_pred->GetBinContent(i)/nevents_signal_region ;
        delta2_pix[a-1] = delta2_pix[a-1] + pow(delta,2);
	sigma_square = npix_sig_test->GetBinContent(i)/nevents_signal_region_test + npix_sig_pred->GetBinContent(i)/nevents_signal_region ;
	sigma2_pix[a-1] = sigma2_pix[a-1] + sigma_square;
	chi2_pix[a-1] = chi2_pix[a-1] + pow(delta,2)/sigma_square ;
	}

	for(int i=1;i<19;i++){
	if(nSCT_sig_test->GetBinContent(i)==0 || nSCT_sig_pred->GetBinContent(i)==0) continue;
        delta = nSCT_sig_test->GetBinContent(i)/nevents_signal_region_test - nSCT_sig_pred->GetBinContent(i)/nevents_signal_region ;
        delta2_SCT[a-1] = delta2_SCT[a-1] + pow(delta,2);
        sigma_square = nSCT_sig_test->GetBinContent(i)/nevents_signal_region_test + nSCT_sig_pred->GetBinContent(i)/nevents_signal_region ;
        sigma2_SCT[a-1] = sigma2_SCT[a-1] + sigma_square;       
        chi2_SCT[a-1] = chi2_SCT[a-1] + pow(delta,2)/sigma_square ;
        }

	for(int i=1;i<71;i++){
	if(nTRT_sig_test->GetBinContent(i)==0 || nTRT_sig_pred->GetBinContent(i)==0) continue;
        delta = nTRT_sig_test->GetBinContent(i)/nevents_signal_region_test - nTRT_sig_pred->GetBinContent(i)/nevents_signal_region ;
        delta2_TRT[a-1] = delta2_TRT[a-1] + pow(delta,2);
        sigma_square = nTRT_sig_test->GetBinContent(i)/nevents_signal_region_test + nTRT_sig_pred->GetBinContent(i)/nevents_signal_region ;
        sigma2_TRT[a-1] = sigma2_TRT[a-1] + sigma_square;       
        chi2_TRT[a-1] = chi2_TRT[a-1] + pow(delta,2)/sigma_square ;
        }

	//Calculate BDM

	// get vector magnitudes
	mag_npix_sig_pred = 0;
        for(int i=1;i<11;i++){
        mag_npix_sig_pred = mag_npix_sig_pred + pow(npix_sig_pred->GetBinContent(i),2);
        }
        mag_npix_sig_pred = sqrt(mag_npix_sig_pred) ;

	mag_nSCT_sig_pred = 0;
        for(int i=1;i<19;i++){
        mag_nSCT_sig_pred = mag_nSCT_sig_pred + pow(nSCT_sig_pred->GetBinContent(i),2);
        }
        mag_nSCT_sig_pred = sqrt(mag_nSCT_sig_pred) ;

	mag_nSCT_sig_pred = 0;
        for(int i=1;i<71;i++){
        mag_nTRT_sig_pred = mag_nTRT_sig_pred + pow(nTRT_sig_pred->GetBinContent(i),2);
        }
        mag_nTRT_sig_pred = sqrt(mag_nTRT_sig_pred) ;

	//bdm is basically just the cos of the angle in between the n-dim vecotors - for exact match bdm = 1
	for(int i=1;i<11;i++){
	if(npix_sig_test->GetBinContent(i)==0 || npix_sig_pred->GetBinContent(i)==0) continue;
	bdm_pix[a-1] = bdm_pix[a-1] + sqrt((npix_sig_test->GetBinContent(i)/nevents_signal_region_test)*(npix_sig_pred->GetBinContent(i)/nevents_signal_region)) ;
	}

	for(int i=1;i<19;i++){
	if(nSCT_sig_test->GetBinContent(i)==0 || nSCT_sig_pred->GetBinContent(i)==0) continue;
        bdm_SCT[a-1] = bdm_SCT[a-1] + sqrt((nSCT_sig_test->GetBinContent(i)/nevents_signal_region_test)*(nSCT_sig_pred->GetBinContent(i)/nevents_signal_region)) ;
        }

        for(int i=1;i<71;i++){
	if(nTRT_sig_test->GetBinContent(i)==0 || nTRT_sig_pred->GetBinContent(i)==0) continue;
        bdm_TRT[a-1] = bdm_TRT[a-1] + sqrt((nTRT_sig_test->GetBinContent(i)/nevents_signal_region_test)*(nTRT_sig_pred->GetBinContent(i)/nevents_signal_region)) ;
        }

    } //end loop over fake fractions

//Make profile histograms for test values

        TProfile *bdmvsf_pix_pfx = new TProfile("bdmvsf_pix_pfx","bdmvsf_pix_pfx",500,0.001,1.001);
        bdmvsf_pix_pfx->Sumw2();
        TProfile *bdmvsf_SCT_pfx = new TProfile("bdmvsf_SCT_pfx","bdmvsf_SCT_pfx",500,0.001,1.001);
        bdmvsf_SCT_pfx->Sumw2();
        TProfile *bdmvsf_TRT_pfx = new TProfile("bdmvsf_TRT_pfx","bdmvsf_TRT_pfx",500,0.001,1.001);
        bdmvsf_TRT_pfx->Sumw2();

        TProfile *chi2vsf_pix_pfx = new TProfile("chi2vsf_pix_pfx","chi2vsf_pix_pfx",500,0.001,1.001);
        chi2vsf_pix_pfx->Sumw2();
        TProfile *chi2vsf_SCT_pfx = new TProfile("chi2vsf_SCT_pfx","chi2vsf_SCT_pfx",500,0.001,1.001);
        chi2vsf_SCT_pfx->Sumw2();
        TProfile *chi2vsf_TRT_pfx = new TProfile("chi2vsf_TRT_pfx","chi2vsf_TRT_pfx",500,0.001,1.001);
        chi2vsf_TRT_pfx->Sumw2();

        TProfile *delta2vf_pix_pfx = new TProfile("delta2vf_pix_pfx","delta2vf_pix_pfx",500,0.001,1.001);
        delta2vf_pix_pfx->Sumw2();
        TProfile *sigma2vf_pix_pfx = new TProfile("sigma2vf_pix_pfx","sigma2vf_pix_pfx",500,0.001,1.001);
        sigma2vf_pix_pfx->Sumw2();
        TProfile *delta2vf_SCT_pfx = new TProfile("delta2vf_SCT_pfx","delta2vf_SCT_pfx",500,0.001,1.001);
        delta2vf_SCT_pfx->Sumw2();
        TProfile *sigma2vf_SCT_pfx = new TProfile("sigma2vf_SCT_pfx","sigma2vf_SCT_pfx",500,0.001,1.001);
        sigma2vf_SCT_pfx->Sumw2();
        TProfile *delta2vf_TRT_pfx = new TProfile("delta2vf_TRT_pfx","delta2vf_TRT_pfx",500,0.001,1.001);
        delta2vf_TRT_pfx->Sumw2();
        TProfile *sigma2vf_TRT_pfx = new TProfile("sigma2vf_TRT_pfx","sigma2vf_TRT_pfx",500,0.001,1.001);
        sigma2vf_TRT_pfx->Sumw2();

        //Fill hists:
        for(int i=0;i<500;i++){
        bdmvsf_pix_pfx->Fill(f[i],bdm_pix[i]);
        bdmvsf_SCT_pfx->Fill(f[i],bdm_SCT[i]);
        bdmvsf_TRT_pfx->Fill(f[i],bdm_TRT[i]);

        delta2vf_pix_pfx->Fill(f[i],delta2_pix[i]);
        sigma2vf_pix_pfx->Fill(f[i],sigma2_pix[i]);
        delta2vf_SCT_pfx->Fill(f[i],delta2_SCT[i]);
        sigma2vf_SCT_pfx->Fill(f[i],sigma2_SCT[i]);
        delta2vf_TRT_pfx->Fill(f[i],delta2_TRT[i]);
        sigma2vf_TRT_pfx->Fill(f[i],sigma2_TRT[i]);

        chi2vsf_pix_pfx->Fill(f[i],chi2_pix[i]);
        chi2vsf_SCT_pfx->Fill(f[i],chi2_SCT[i]);
        chi2vsf_TRT_pfx->Fill(f[i],chi2_TRT[i]);
        }

////PLOTS

        //// For BDM
        TCanvas *c_bdmvsf_pix = new TCanvas("c_bdmvsf_pix","c_bdmvsf_pix",1);
        bdmvsf_pix_pfx->SetStats(kFALSE);
        bdmvsf_pix_pfx->GetXaxis()->SetTitle("fraction of fake events");
        bdmvsf_pix_pfx->GetYaxis()->SetTitle("BDM");
        	if(isMC) bdmvsf_pix_pfx->SetTitle("BDM shape comparison test - pixel - MC");
                if(!isMC) bdmvsf_pix_pfx->SetTitle("BDM shape comparison test - pixel - data");
        bdmvsf_pix_pfx->Draw();
        c_bdmvsf_pix->Print("PLOTS/CTIDE_08Aug2018/c_bdmvsf_pix.pdf");

        TCanvas *c_bdmvsf_SCT = new TCanvas("c_bdmvsf_SCT","c_bdmvsf_SCT",1);
        bdmvsf_SCT_pfx->SetStats(kFALSE);
        bdmvsf_SCT_pfx->GetXaxis()->SetTitle("fraction of fake events");
        bdmvsf_SCT_pfx->GetYaxis()->SetTitle("BDM");
	        if(isMC) bdmvsf_SCT_pfx->SetTitle("BDM shape comparison test - SCT - MC");
                if(!isMC) bdmvsf_SCT_pfx->SetTitle("BDM shape comparison test - SCT - data");
        bdmvsf_SCT_pfx->Draw();
        c_bdmvsf_SCT->Print("PLOTS/CTIDE_08Aug2018/c_bdmvsf_SCT.pdf");

        TCanvas *c_bdmvsf_TRT = new TCanvas("c_bdmvsf_TRT","c_bdmvsf_TRT",1);
        bdmvsf_TRT_pfx->SetStats(kFALSE);
        bdmvsf_TRT_pfx->GetXaxis()->SetTitle("fraction of fake events");
        bdmvsf_TRT_pfx->GetYaxis()->SetTitle("BDM");
	        if(isMC) bdmvsf_TRT_pfx->SetTitle("BDM shape comparison test - TRT - MC");
                if(!isMC) bdmvsf_TRT_pfx->SetTitle("BDM shape comparison test - TRT - data");
        bdmvsf_TRT_pfx->Draw();
        c_bdmvsf_TRT->Print("PLOTS/CTIDE_08Aug2018/c_bdmvsf_TRT.pdf");

	////For Chi^2
        TCanvas *cpix_chi2 = new TCanvas("cpix_chi2","cpix_chi2",1);
        chi2vsf_pix_pfx->SetStats(kFALSE);
        chi2vsf_pix_pfx->GetXaxis()->SetTitle("fraction of fake tracks");
        chi2vsf_pix_pfx->GetYaxis()->SetTitle("#chi^{2}");
        	if(isMC) chi2vsf_pix_pfx->SetTitle("#chi^{2} shape comparison test - pixel hit - MC");
                if(!isMC) chi2vsf_pix_pfx->SetTitle("#chi^{2} shape comparison test - pixel hit - data");
        chi2vsf_pix_pfx->Draw();
        cpix_chi2->Print("PLOTS/CTIDE_08Aug2018/cpix_chi2.pdf");

        TCanvas *cSCT_chi2 = new TCanvas("cSCT_chi2","cSCT_chi2",1);
        chi2vsf_SCT_pfx->SetStats(kFALSE);
        chi2vsf_SCT_pfx->GetXaxis()->SetTitle("fraction of fake tracks");
        chi2vsf_SCT_pfx->GetYaxis()->SetTitle("#chi^{2}");
        	if(isMC) chi2vsf_SCT_pfx->SetTitle("#chi^{2} shape comparison test - SCT hit - MC");
                if(!isMC) chi2vsf_SCT_pfx->SetTitle("#chi^{2} shape comparison test - SCT hit - data");
        chi2vsf_SCT_pfx->Draw();
        cSCT_chi2->Print("PLOTS/CTIDE_08Aug2018/cSCT_chi2.pdf");

        TCanvas *cTRT_chi2 = new TCanvas("cTRT_chi2","cTRT_chi2",1);
        chi2vsf_TRT_pfx->SetStats(kFALSE);
        chi2vsf_TRT_pfx->GetXaxis()->SetTitle("fraction of fake tracks");
        chi2vsf_TRT_pfx->GetYaxis()->SetTitle("#chi^{2}");
        	if(isMC) chi2vsf_TRT_pfx->SetTitle("#chi^{2} shape comparison test - TRT hit - MC");
                if(!isMC) chi2vsf_TRT_pfx->SetTitle("#chi^{2} shape comparison test - TRT hit - data");
        chi2vsf_TRT_pfx->Draw();
        cTRT_chi2->Print("PLOTS/CTIDE_08Aug2018/cTRT_chi2.pdf");

/*        TCanvas *cdelta2_pix = new TCanvas("cdelta2_pix","cdelta2_pix",1);
        delta2vf_pix_pfx->SetStats(kFALSE);
        delta2vf_pix_pfx->Draw();
        cdelta2_pix->Print("PLOTS/CTIDE_08Aug2018/delta2_pix.pdf");

        TCanvas *csigma2_pix = new TCanvas("csigma2_pix","csigma2_pix",1);
        sigma2vf_pix_pfx->SetStats(kFALSE);
        sigma2vf_pix_pfx->Draw();
        csigma2_pix->Print("PLOTS/CTIDE_08Aug2018/sigma2_pix.pdf");

        TCanvas *cdelta2_SCT = new TCanvas("cdelta2_SCT","cdelta2_SCT",1);
        delta2vf_SCT_pfx->SetStats(kFALSE);
        delta2vf_SCT_pfx->Draw();
        cdelta2_SCT->Print("PLOTS/CTIDE_08Aug2018/delta2_SCT.pdf");

        TCanvas *csigma2_SCT = new TCanvas("csigma2_SCT","csigma2_SCT",1);
        sigma2vf_SCT_pfx->SetStats(kFALSE);
        sigma2vf_SCT_pfx->Draw();
        csigma2_SCT->Print("PLOTS/CTIDE_08Aug2018/sigma2_SCT.pdf");

        TCanvas *cdelta2_TRT = new TCanvas("cdelta2_TRT","cdelta2_TRT",1);
        delta2vf_TRT_pfx->SetStats(kFALSE);
        delta2vf_TRT_pfx->Draw();
        cdelta2_TRT->Print("PLOTS/CTIDE_08Aug2018/delta2_TRT.pdf");

        TCanvas *csigma2_TRT = new TCanvas("csigma2_TRT","csigma2_TRT",1);
        sigma2vf_TRT_pfx->SetStats(kFALSE);
        sigma2vf_TRT_pfx->Draw();
        csigma2_TRT->Print("PLOTS/CTIDE_08Aug2018/sigma2_TRT.pdf");
*/
	//comparing test and build distributions: 
	//signal region
/*	TCanvas *cpix_sig_testvsbuild = new TCanvas("cpix_sig_testvsbuild","cpix_sig_testvsbuild",1);
	npix_sig->SetStats(kFALSE);
	npix_sig_test->SetStats(kFALSE);
	npix_sig->GetXaxis()->SetTitle("number of pixel hits");
	if(isMC) npix_sig->SetTitle("signal region - pixel hits - MC");
	if(!isMC) npix_sig->SetTitle("signal region - pixel hits - data");
	npix_sig->SetFillStyle(3004);
	npix_sig->SetFillColor(kBlue);
	npix_sig_test->SetFillStyle(3005);
	npix_sig_test->SetFillColor(kRed);
	npix_sig->Draw("bar");
	npix_sig_test->Draw("same bar");
	TLegend *lpix_sig_testvsbuild = new TLegend(.6,.6,.9,.9);
	lpix_sig_testvsbuild->AddEntry("npix_sig","build","f");
	lpix_sig_testvsbuild->AddEntry("npix_sig_test","test","f");
	lpix_sig_testvsbuild->Draw();
	cpix_sig_testvsbuild->Print("PLOTS/CTIDE_08Aug2018/cpix_sig_testvsbuild.pdf");

        TCanvas *cSCT_sig_testvsbuild = new TCanvas("cSCT_sig_testvsbuild","cSCT_sig_testvsbuild",1);
	nSCT_sig->SetStats(kFALSE);
	nSCT_sig_test->SetStats(kFALSE);
        nSCT_sig->GetXaxis()->SetTitle("number of SCT hits");
        if(isMC) nSCT_sig->SetTitle("signal region - SCT hits - MC");
        if(!isMC) nSCT_sig->SetTitle("signal region - SCT hits - data");
        nSCT_sig->SetFillStyle(3004);
        nSCT_sig->SetFillColor(kBlue);
        nSCT_sig_test->SetFillStyle(3005);
        nSCT_sig_test->SetFillColor(kRed);
        nSCT_sig->Draw("bar");
        nSCT_sig_test->Draw("same bar");
        TLegend *lSCT_sig_testvsbuild = new TLegend(.6,.6,.9,.9);
        lSCT_sig_testvsbuild->AddEntry("nSCT_sig","build","f");
        lSCT_sig_testvsbuild->AddEntry("nSCT_sig_test","test","f");
        lSCT_sig_testvsbuild->Draw();
        cSCT_sig_testvsbuild->Print("PLOTS/CTIDE_08Aug2018/cSCT_sig_testvsbuild.pdf");

        TCanvas *cTRT_sig_testvsbuild = new TCanvas("cTRT_sig_testvsbuild","cTRT_sig_testvsbuild",1);
	nTRT_sig->SetStats(kFALSE);
	nTRT_sig_test->SetStats(kFALSE);
        nTRT_sig->GetXaxis()->SetTitle("number of TRT hits");
        if(isMC) nTRT_sig->SetTitle("signal region - TRT hits - MC");
        if(!isMC) nTRT_sig->SetTitle("signal region - TRT hits - data");
        nTRT_sig->SetFillStyle(3004);
        nTRT_sig->SetFillColor(kBlue);
        nTRT_sig_test->SetFillStyle(3005);
        nTRT_sig_test->SetFillColor(kRed);
        nTRT_sig->Draw("bar");
        nTRT_sig_test->Draw("same bar");
        TLegend *lTRT_sig_testvsbuild = new TLegend(.6,.6,.9,.9);
        lTRT_sig_testvsbuild->AddEntry("nTRT_sig","build","f");
        lTRT_sig_testvsbuild->AddEntry("nTRT_sig_test","test","f");
        lTRT_sig_testvsbuild->Draw();
        cTRT_sig_testvsbuild->Print("PLOTS/CTIDE_08Aug2018/cTRT_sig_testvsbuild.pdf");

	//sideband region
        TCanvas *cpix_side_testvsbuild = new TCanvas("cpix_side_testvsbuild","cpix_side_testvsbuild",1);
	npix_side->SetStats(kFALSE);
	npix_side_test->SetStats(kFALSE);
        npix_side->GetXaxis()->SetTitle("number of pixel hits");
        if(isMC) npix_side->SetTitle("sideband region - pixel hits - MC");
        if(!isMC) npix_side->SetTitle("sideband region - pixel hits - data");
        npix_side->SetFillStyle(3004);
        npix_side->SetFillColor(kBlue);
        npix_side_test->SetFillStyle(3005);
        npix_side_test->SetFillColor(kRed);
        npix_side->Draw("bar");
        npix_side_test->Draw("same bar");
        TLegend *lpix_side_testvsbuild = new TLegend(.6,.6,.9,.9);
        lpix_side_testvsbuild->AddEntry("npix_side","build","f");
        lpix_side_testvsbuild->AddEntry("npix_side_test","test","f");
        lpix_side_testvsbuild->Draw();
        cpix_side_testvsbuild->Print("PLOTS/CTIDE_08Aug2018/cpix_side_testvsbuild.pdf");

        TCanvas *cSCT_side_testvsbuild = new TCanvas("cSCT_side_testvsbuild","cSCT_side_testvsbuild",1);
	nSCT_side->SetStats(kFALSE);
        nSCT_side_test->SetStats(kFALSE);
        nSCT_side->GetXaxis()->SetTitle("number of SCT hits");
        if(isMC) nSCT_side->SetTitle("sideband region - SCT hits - MC");
        if(!isMC) nSCT_side->SetTitle("sideband region - SCT hits - data");
        nSCT_side->SetFillStyle(3004);
        nSCT_side->SetFillColor(kBlue);
        nSCT_side_test->SetFillStyle(3005);
        nSCT_side_test->SetFillColor(kRed);
        nSCT_side->Draw("bar");
        nSCT_side_test->Draw("same bar");
        TLegend *lSCT_side_testvsbuild = new TLegend(.6,.6,.9,.9);
        lSCT_side_testvsbuild->AddEntry("nSCT_side","build","f");
        lSCT_side_testvsbuild->AddEntry("nSCT_side_test","test","f");
        lSCT_side_testvsbuild->Draw();
        cSCT_side_testvsbuild->Print("PLOTS/CTIDE_08Aug2018/cSCT_side_testvsbuild.pdf");

        TCanvas *cTRT_side_testvsbuild = new TCanvas("cTRT_side_testvsbuild","cTRT_side_testvsbuild",1);
	nTRT_side->SetStats(kFALSE);
        nTRT_side_test->SetStats(kFALSE);
        nTRT_side->GetXaxis()->SetTitle("number of TRT hits");
        if(isMC) nTRT_side->SetTitle("sideband region - TRT hits - MC");
        if(!isMC) nTRT_side->SetTitle("sideband region - TRT hits - data");
        nTRT_side->SetFillStyle(3004);
        nTRT_side->SetFillColor(kBlue);
        nTRT_side_test->SetFillStyle(3005);
        nTRT_side_test->SetFillColor(kRed);
        nTRT_side->Draw("bar");
        nTRT_side_test->Draw("same bar");
        TLegend *lTRT_side_testvsbuild = new TLegend(.6,.6,.9,.9);
        lTRT_side_testvsbuild->AddEntry("nTRT_side","build","f");
        lTRT_side_testvsbuild->AddEntry("nTRT_side_test","test","f");
        lTRT_side_testvsbuild->Draw();
        cTRT_side_testvsbuild->Print("PLOTS/CTIDE_08Aug2018/cTRT_side_testvsbuild.pdf");

	//predictions for signal region for different values of ffake:
	ffake = 0.1;
        npix_sig_pred->Add(npix_real_norm,npix_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<11;i++){
                if(npix_sig_pred->GetBinContent(i)<0) npix_sig_pred->SetBinContent(i,0);
                }
        npix_sig_pred->Scale(nevents_signal_region/npix_sig_pred->Integral());

        nSCT_sig_pred->Add(nSCT_real_norm,nSCT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<19;i++){
                if(nSCT_sig_pred->GetBinContent(i)<0) nSCT_sig_pred->SetBinContent(i,0);
                }
        nSCT_sig_pred->Scale(nevents_signal_region/nSCT_sig_pred->Integral());

        nTRT_sig_pred->Add(nTRT_real_norm,nTRT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<71;i++){
                if(nTRT_sig_pred->GetBinContent(i)<0) nTRT_sig_pred->SetBinContent(i,0);
                }
        nTRT_sig_pred->Scale(nevents_signal_region/nTRT_sig_pred->Integral());	

		//plots prediction and signal
		TCanvas *cpix_sig_testvspred_01 = new TCanvas("cpix_sig_testvspred_01","cpix_sig_testvspred_01",1);
		npix_sig_pred->SetStats(kFALSE);
		npix_sig_pred->GetXaxis()->SetTitle("number of pixel hits");
		if(isMC) npix_sig_pred->SetTitle("signal region - pixel hits - MC");
		if(!isMC) npix_sig_pred->SetTitle("signal region - pixel hits - data");
	        npix_sig_pred->SetFillStyle(3004);
	        npix_sig_pred->SetFillColor(kBlue);
	        npix_sig_test->SetFillStyle(3005);
	        npix_sig_test->SetFillColor(kRed);
		npix_sig_pred->Draw("bar");
		npix_sig_test->Draw("same bar");
		TLegend *lpix_sig_testvspred_01 = new TLegend(.6,.6,.9,.9);
                lpix_sig_testvspred_01->AddEntry((TObject*)0,"f_{fake} = 0.1","");
		lpix_sig_testvspred_01->AddEntry("npix_sig_pred","prediction","f");
		lpix_sig_testvspred_01->AddEntry("npix_sig_test","test","f");
		lpix_sig_testvspred_01->Draw();
		cpix_sig_testvspred_01->Print("PLOTS/CTIDE_08Aug2018/cpix_sig_testvspred_01.pdf");

                TCanvas *cSCT_sig_testvspred_01 = new TCanvas("cSCT_sig_testvspred_01","cSCT_sig_testvspred_01",1);
                nSCT_sig_pred->SetStats(kFALSE);
                nSCT_sig_pred->GetXaxis()->SetTitle("number of SCT hits");
                if(isMC) nSCT_sig_pred->SetTitle("signal region - SCT hits - MC");
                if(!isMC) nSCT_sig_pred->SetTitle("signal region - SCT hits - data");
                nSCT_sig_pred->SetFillStyle(3004);
                nSCT_sig_pred->SetFillColor(kBlue);
                nSCT_sig_test->SetFillStyle(3005);
                nSCT_sig_test->SetFillColor(kRed);
                nSCT_sig_pred->Draw("bar");
                nSCT_sig_test->Draw("same bar");
                TLegend *lSCT_sig_testvspred_01 = new TLegend(.6,.6,.9,.9);
		lSCT_sig_testvspred_01->AddEntry((TObject*)0,"f_{fake} = 0.1","");
                lSCT_sig_testvspred_01->AddEntry("nSCT_sig_pred","prediction","f");
                lSCT_sig_testvspred_01->AddEntry("nSCT_sig_test","test","f");
                lSCT_sig_testvspred_01->Draw();
                cSCT_sig_testvspred_01->Print("PLOTS/CTIDE_08Aug2018/cSCT_sig_testvspred_01.pdf");

                TCanvas *cTRT_sig_testvspred_01 = new TCanvas("cTRT_sig_testvspred_01","cTRT_sig_testvspred_01",1);
                nTRT_sig_pred->SetStats(kFALSE);
                nTRT_sig_pred->GetXaxis()->SetTitle("number of TRT hits");
                if(isMC) nTRT_sig_pred->SetTitle("signal region - TRT hits - MC");
                if(!isMC) nTRT_sig_pred->SetTitle("signal region - TRT hits - data");
                nTRT_sig_pred->SetFillStyle(3004);
                nTRT_sig_pred->SetFillColor(kBlue);
                nTRT_sig_test->SetFillStyle(3005);
                nTRT_sig_test->SetFillColor(kRed);
                nTRT_sig_pred->Draw("bar");
                nTRT_sig_test->Draw("same bar");
                TLegend *lTRT_sig_testvspred_01 = new TLegend(.6,.6,.9,.9);
		lTRT_sig_testvspred_01->AddEntry((TObject*)0,"f_{fake} = 0.1","");
                lTRT_sig_testvspred_01->AddEntry("nTRT_sig_pred","prediction","f");
                lTRT_sig_testvspred_01->AddEntry("nTRT_sig_test","test","f");
                lTRT_sig_testvspred_01->Draw();
                cTRT_sig_testvspred_01->Print("PLOTS/CTIDE_08Aug2018/cTRT_sig_testvspred_01.pdf");

        ffake = 0.2;
        npix_sig_pred->Add(npix_real_norm,npix_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<11;i++){
                if(npix_sig_pred->GetBinContent(i)<0) npix_sig_pred->SetBinContent(i,0);
                }
        npix_sig_pred->Scale(nevents_signal_region/npix_sig_pred->Integral());

        nSCT_sig_pred->Add(nSCT_real_norm,nSCT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<19;i++){
                if(nSCT_sig_pred->GetBinContent(i)<0) nSCT_sig_pred->SetBinContent(i,0);
                }
        nSCT_sig_pred->Scale(nevents_signal_region/nSCT_sig_pred->Integral());

        nTRT_sig_pred->Add(nTRT_real_norm,nTRT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<71;i++){
                if(nTRT_sig_pred->GetBinContent(i)<0) nTRT_sig_pred->SetBinContent(i,0);
                }
        nTRT_sig_pred->Scale(nevents_signal_region/nTRT_sig_pred->Integral());

                //plots prediction and signal
                TCanvas *cpix_sig_testvspred_02 = new TCanvas("cpix_sig_testvspred_02","cpix_sig_testvspred_02",1);
                npix_sig_pred->GetXaxis()->SetTitle("number of pixel hits");
                if(isMC) npix_sig_pred->SetTitle("signal region - pixel hits - MC");
                if(!isMC) npix_sig_pred->SetTitle("signal region - pixel hits - data");
                npix_sig_pred->Draw("bar");
                npix_sig_test->Draw("same bar");
                TLegend *lpix_sig_testvspred_02 = new TLegend(.6,.6,.9,.9);
		lpix_sig_testvspred_02->AddEntry((TObject*)0,"f_{fake} = 0.2","");
                lpix_sig_testvspred_02->AddEntry("npix_sig_pred","prediction","f");
                lpix_sig_testvspred_02->AddEntry("npix_sig_test","test","f");
                lpix_sig_testvspred_02->Draw();
                cpix_sig_testvspred_02->Print("PLOTS/CTIDE_08Aug2018/cpix_sig_testvspred_02.pdf");

                TCanvas *cSCT_sig_testvspred_02 = new TCanvas("cSCT_sig_testvspred_02","cSCT_sig_testvspred_02",1);
                nSCT_sig_pred->GetXaxis()->SetTitle("number of SCT hits");
                if(isMC) nSCT_sig_pred->SetTitle("signal region - SCT hits - MC");
                if(!isMC) nSCT_sig_pred->SetTitle("signal region - SCT hits - data");
                nSCT_sig_pred->Draw("bar");
                nSCT_sig_test->Draw("same bar");
                TLegend *lSCT_sig_testvspred_02 = new TLegend(.6,.6,.9,.9);
		lSCT_sig_testvspred_02->AddEntry((TObject*)0,"f_{fake} = 0.2","");
                lSCT_sig_testvspred_02->AddEntry("nSCT_sig_pred","prediction","f");
                lSCT_sig_testvspred_02->AddEntry("nSCT_sig_test","test","f");
                lSCT_sig_testvspred_02->Draw();
                cSCT_sig_testvspred_02->Print("PLOTS/CTIDE_08Aug2018/cSCT_sig_testvspred_02.pdf");

                TCanvas *cTRT_sig_testvspred_02 = new TCanvas("cTRT_sig_testvspred_02","cTRT_sig_testvspred_02",1);
                nTRT_sig_pred->GetXaxis()->SetTitle("number of TRT hits");
                if(isMC) nTRT_sig_pred->SetTitle("signal region - TRT hits - MC");
                if(!isMC) nTRT_sig_pred->SetTitle("signal region - TRT hits - data");
                nTRT_sig_pred->Draw("bar");
                nTRT_sig_test->Draw("same bar");
                TLegend *lTRT_sig_testvspred_02 = new TLegend(.6,.6,.9,.9);
		lTRT_sig_testvspred_02->AddEntry((TObject*)0,"f_{fake} = 0.2","");
                lTRT_sig_testvspred_02->AddEntry("nTRT_sig_pred","prediction","f");
                lTRT_sig_testvspred_02->AddEntry("nTRT_sig_test","test","f");
                lTRT_sig_testvspred_02->Draw();
                cTRT_sig_testvspred_02->Print("PLOTS/CTIDE_08Aug2018/cTRT_sig_testvspred_02.pdf");

        ffake = 0.3;
        npix_sig_pred->Add(npix_real_norm,npix_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<11;i++){
                if(npix_sig_pred->GetBinContent(i)<0) npix_sig_pred->SetBinContent(i,0);
                }
        npix_sig_pred->Scale(nevents_signal_region/npix_sig_pred->Integral());

        nSCT_sig_pred->Add(nSCT_real_norm,nSCT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<19;i++){
                if(nSCT_sig_pred->GetBinContent(i)<0) nSCT_sig_pred->SetBinContent(i,0);
                }
        nSCT_sig_pred->Scale(nevents_signal_region/nSCT_sig_pred->Integral());

        nTRT_sig_pred->Add(nTRT_real_norm,nTRT_fake_norm,(fsig+(1-fsig)*(1-ffake)),ffake*(1-fsig));
                for(int i=1;i<71;i++){
                if(nTRT_sig_pred->GetBinContent(i)<0) nTRT_sig_pred->SetBinContent(i,0);
                }
        nTRT_sig_pred->Scale(nevents_signal_region/nTRT_sig_pred->Integral());

                //plots prediction and signal
                TCanvas *cpix_sig_testvspred_03 = new TCanvas("cpix_sig_testvspred_03","cpix_sig_testvspred_03",1);
                npix_sig_pred->GetXaxis()->SetTitle("number of pixel hits");
                if(isMC) npix_sig_pred->SetTitle("signal region - pixel hits - MC");
                if(!isMC) npix_sig_pred->SetTitle("signal region - pixel hits - data");
                npix_sig_pred->Draw("bar");
                npix_sig_test->Draw("same bar");
                TLegend *lpix_sig_testvspred_03 = new TLegend(.6,.6,.9,.9);
                lpix_sig_testvspred_03->AddEntry((TObject*)0,"f_{fake} = 0.3","");
                lpix_sig_testvspred_03->AddEntry("npix_sig_pred","prediction","f");
                lpix_sig_testvspred_03->AddEntry("npix_sig_test","test","f");
                lpix_sig_testvspred_03->Draw();
                cpix_sig_testvspred_03->Print("PLOTS/CTIDE_08Aug2018/cpix_sig_testvspred_03.pdf");

                TCanvas *cSCT_sig_testvspred_03 = new TCanvas("cSCT_sig_testvspred_03","cSCT_sig_testvspred_03",1);
                nSCT_sig_pred->GetXaxis()->SetTitle("number of SCT hits");
                if(isMC) nSCT_sig_pred->SetTitle("signal region - SCT hits - MC");
                if(!isMC) nSCT_sig_pred->SetTitle("signal region - SCT hits - data");
                nSCT_sig_pred->Draw("bar");
                nSCT_sig_test->Draw("same bar");
                TLegend *lSCT_sig_testvspred_03 = new TLegend(.6,.6,.9,.9);
                lSCT_sig_testvspred_03->AddEntry((TObject*)0,"f_{fake} = 0.3","");
                lSCT_sig_testvspred_03->AddEntry("nSCT_sig_pred","prediction","f");
                lSCT_sig_testvspred_03->AddEntry("nSCT_sig_test","test","f");
                lSCT_sig_testvspred_03->Draw();
                cSCT_sig_testvspred_03->Print("PLOTS/CTIDE_08Aug2018/cSCT_sig_testvspred_03.pdf");

                TCanvas *cTRT_sig_testvspred_03 = new TCanvas("cTRT_sig_testvspred_03","cTRT_sig_testvspred_03",1);
                nTRT_sig_pred->GetXaxis()->SetTitle("number of TRT hits");
                if(isMC) nTRT_sig_pred->SetTitle("signal region - TRT hits - MC");
                if(!isMC) nTRT_sig_pred->SetTitle("signal region - TRT hits - data");
                nTRT_sig_pred->Draw("bar");
                nTRT_sig_test->Draw("same bar");
                TLegend *lTRT_sig_testvspred_03 = new TLegend(.6,.6,.9,.9);
                lTRT_sig_testvspred_03->AddEntry((TObject*)0,"f_{fake} = 0.3","");
                lTRT_sig_testvspred_03->AddEntry("nTRT_sig_pred","prediction","f");
                lTRT_sig_testvspred_03->AddEntry("nTRT_sig_test","test","f");
                lTRT_sig_testvspred_03->Draw();
                cTRT_sig_testvspred_03->Print("PLOTS/CTIDE_08Aug2018/cTRT_sig_testvspred_03.pdf");

*/
file.Write();

}//very end