{
  gStyle->SetPalette(1);
  // TFile *file= new TFile("../../Outputs/Analysis/PhysicsTree.root");
  // TFile *file= new TFile("../../Outputs/Analysis/configJuly2021_ana_1MeV_perfect_Decay.root");
  // TFile* file = new TFile("../../Outputs/Analysis/configJuly2021_ana_1MeV.root");
  // TFile* file = new TFile("../../Outputs/Analysis/configJune2022_ana_0MeV_real_100MeVc_50um.root");
  // TFile* file = new TFile("../../Outputs/Analysis/configJune2022_ana_0MeV_real_100MeVc_I200um_O300um.root");
  // TFile* file = new TFile("../../Outputs/Analysis/configJune2022_ana_0MeV_real_100MeVc_I200um_O300um_newtar.root");
  TFile* file = new TFile("../../Outputs/Analysis/test_newtar_ana.root");
  // TFile* file = new TFile("../../Outputs/Analysis/test_oldtar_ana.root");

  TTree* tree = (TTree*)file->Get("PhysicsTree");

  TCanvas* c1 = new TCanvas("c1", "c1", 1000, 1000);
  c1->Divide(3, 4);
  c1->cd(1);
  tree->Draw("Ex>>h1(400,-10,10)", "", "");
  h1->SetXTitle("Ex_{missingmass}");
  h1->SetYTitle("counts / 50 keV");
  c1->cd(2);
  tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(400,-150,250,200,-100,100)", "", "colz");
  TBox* box = new TBox(-75, -15, 75, 15);
  box->SetFillStyle(0);
  box->SetLineColor(kRed);
  box->Draw("same");
  c1->cd(3);
  tree->Draw("VertexX:VertexZ>>h4(400,-150,250,200,-100,100)", "", "colz");
  box->Draw("same");
  double xmin = -75;
  int xmax = 75;
  double ymin = -15;
  int ymax = 15;
  int binx1 = h4->GetXaxis()->FindBin(xmin);
  int binx2 = h4->GetXaxis()->FindBin(xmax);
  int biny1 = h4->GetYaxis()->FindBin(ymin);
  int biny2 = h4->GetYaxis()->FindBin(ymax);
  double integralall = h3->Integral(binx1, binx2, biny1, biny2);
  double integralmult2 = h4->Integral();
  double integralrec = h4->Integral(binx1, binx2, biny1, biny2);
  double integralEx = h1->Integral();
  cout << "Simulated = " << integralall << endl;
  cout << "Mult2 from vertex = " << integralmult2 << " (" << integralmult2 / integralall * 100 << "%)" << endl;
  cout << "Rec. within target = " << integralrec << endl;
  cout << "Rec. within target(%) = " << integralrec / integralall * 100 << endl;
  cout << "Rec. Ex (%) = " << integralEx / integralall * 100 << endl;
  c1->cd(4);
  tree->Draw("fRC_Vertex_Position_X-VertexX>>h5(100,-2,2)", "", "");
  c1->cd(5);
  tree->Draw("fRC_Vertex_Position_Y-VertexY>>h6(100,-2,2)", "", "");
  c1->cd(6);
  tree->Draw("fRC_Vertex_Position_Z-VertexZ>>h7(100,-2,2)", "", "");
  c1->cd(7);
  // tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[0]>>h8(100,-10,10)","Catana.GetEventMultiplicity()==2","");
  tree->Draw("fRC_Kinetic_Energy[1]-E1>>h8(100,-50,50)", "", "");
  c1->cd(8);
  // tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[1]>>h9(100,-10,10)","Catana.GetEventMultiplicity()==2","");
  tree->Draw("fRC_Kinetic_Energy[0]-E2>>h9(100,-50,50)", "", "");
  c1->cd(9);
  tree->Draw("fRC_Kinetic_Energy[0]:E2>>h10(300,0,300,300,0,300)", "", "colz");
  c1->cd(10);
  tree->Draw("Strasse.EventMultiplicity>>h11(10,0,10)", "", "colz");
  double integralM2 = h11->Integral(3, 3);
  double integralM2bis = h11->Integral(2, 3);
  cout << "Mult2 Strasse.EventMult = " << integralM2 << " (" << integralM2 / integralall * 100 << "%)" << endl;
  cout << "Mult1and2 Strasse.EventMult = " << integralM2bis << " (" << integralM2bis / integralall * 100 << "%)"
       << endl;
}