void can_wbt_ai(void)
{
	auto fil1 = TFile::Open("TMVAOutputBDT_p_can3.root");
	fil1->cd("dataset");
	auto tBDTG_p = (TTree*)gROOT->FindObject("TestTree");

	auto fil2 = TFile::Open("TMVAOutputBDT_theta_can3.root");
	fil2->cd("dataset");
	auto tBDTG_theta = (TTree*)gROOT->FindObject("TestTree");


	auto fil3 = TFile::Open("TMVAOutputBDT_p3.root");
	fil3->cd("dataset");
	auto tBDTG_p3 = (TTree*)gROOT->FindObject("TestTree");

	auto fil4 = TFile::Open("TMVAOutputBDT_theta3.root");
	fil4->cd("dataset");
	auto tBDTG_theta3 = (TTree*)gROOT->FindObject("TestTree");


	TChain *chain = new TChain("t");
	chain->Add("../DATA/hd_root*.root");
	
	TCut mycut1("abs(p-1.0)<0.100");
	TCut mycut2("abs(target-1.0)<0.100");

	auto hpcan = new TH1D("hpcan", "Momentum resolution" , 200, -1.0, 1.0);
	auto hpwbt = new TH1D("hpwbt", "Momentum resolution", 200, -1.0, 1.0);
	auto hpbdt = new TH1D("hpbdt", "Momentum resolution"      , 200, -1.0, 1.0);
	auto hpbdt3 = new TH1D("hpbdt3", "Momentum resolution"    , 200, -1.0, 1.0);
	chain->Project("hpcan", "(p_can-p)/p", mycut1);
	chain->Project("hpwbt", "(p_wbt-p)/p", mycut1);
	tBDTG_p->Project("hpbdt", "(BDTG-target)/target", mycut2);
	tBDTG_p3->Project("hpbdt3", "(BDTG-target)/target", mycut2);

	auto hthetacan = new TH1D("hthetacan", "#theta angle resolution" , 200, -20.0, 40.0);
	auto hthetawbt = new TH1D("hthetawbt", "#theta angle resolution", 200, -20.0, 40.0);
	auto hthetabdt = new TH1D("hthetabdt", "#theta angle resolution"      , 200, -20.0, 40.0);
	auto hthetabdt3 = new TH1D("hthetabdt3", "#theta angle resolution"    , 200, -20.0, 40.0);
	chain->Project("hthetacan", "1000.0*(theta_can-theta)");
	chain->Project("hthetawbt", "1000.0*(theta_wbt-theta)");
	tBDTG_theta->Project("hthetabdt", "1000.0*(BDTG-target)");
	tBDTG_theta3->Project("hthetabdt3", "1000.0*(BDTG-target)");

	// Using Scale method makes histogram draw as points	
	auto scale = hpwbt->Integral()/hpbdt->Integral();
	for(int ibin=1; ibin<=hpwbt->GetNbinsX(); ibin++) hpbdt->SetBinContent(ibin, scale*hpbdt->GetBinContent(ibin) );

	scale = hthetawbt->Integral()/hthetabdt->Integral();
	for(int ibin=1; ibin<=hthetabdt->GetNbinsX(); ibin++) hthetabdt->SetBinContent(ibin, scale*hthetabdt->GetBinContent(ibin) );

	scale = hpwbt->Integral()/hpbdt3->Integral();
	for(int ibin=1; ibin<=hpwbt->GetNbinsX(); ibin++) hpbdt3->SetBinContent(ibin, scale*hpbdt3->GetBinContent(ibin) );

	scale = hthetawbt->Integral()/hthetabdt3->Integral();
	for(int ibin=1; ibin<=hthetabdt->GetNbinsX(); ibin++) hthetabdt3->SetBinContent(ibin, scale*hthetabdt3->GetBinContent(ibin) );


	TCanvas *c1 = new TCanvas("c1", "", 1800, 900);
	c1->Divide(2,1);
	
	c1->cd(1);
	gPad->SetTicks();
	gPad->SetGrid();
	hpwbt->SetLineWidth(3);
	hpwbt->SetLineColor(kBlue);
	hpcan->SetLineWidth(3);
	hpcan->SetLineColor(kRed);
	hpbdt->SetLineWidth(3);
	hpbdt->SetLineColor(kMagenta);
	hpbdt3->SetLineWidth(3);
	hpbdt3->SetLineColor(kBlack);

	hpwbt->SetStats(0);
	hpwbt->SetXTitle("#deltap/p (fraction)");
	hpwbt->Draw();
	hpcan->Draw("same");
	hpbdt->Draw("same");
	hpbdt3->Draw("same");
	
	TLegend *leg = new TLegend(0.6, 0.6, 0.89, 0.89);
	leg->AddEntry(hpwbt, "wire-based");
	leg->AddEntry(hpcan, "candidate");
	leg->AddEntry(hpbdt, "BDTG - can. par.");
	leg->AddEntry(hpbdt3, "BDTG - hits");
	leg->Draw();

	c1->cd(2);
	gPad->SetTicks();
	gPad->SetGrid();
	hthetawbt->SetLineWidth(3);
	hthetawbt->SetLineColor(kBlue);
	hthetacan->SetLineWidth(3);
	hthetacan->SetLineColor(kRed);
	hthetabdt->SetLineWidth(3);
	hthetabdt->SetLineColor(kMagenta);
	hthetabdt3->SetLineWidth(3);
	hthetabdt3->SetLineColor(kBlack);

	hthetawbt->SetStats(0);
	hthetawbt->SetXTitle("#theta (mrad)");
	hthetawbt->Draw();
	hthetacan->Draw("same");
	hthetabdt->Draw("same");
	hthetabdt3->Draw("same");
	
	leg = new TLegend(0.6, 0.6, 0.89, 0.89);
	leg->AddEntry(hthetawbt, "wire-based");
	leg->AddEntry(hthetacan, "candidate");
	leg->AddEntry(hthetabdt, "BDTG - can. par.");
	leg->AddEntry(hthetabdt3, "BDTG - hits");
	leg->Draw();
}