void reconstruct_D0()
  
{
  TFile *f = new TFile("FileName.root");
  TTreeReader reader("events", f);
  
  TTreeReaderArray<Float_t> pxArray(reader, "ReconstructedParticles.p.x");
  TTreeReaderArray<Float_t> pyArray(reader, "ReconstructedParticles.p.y");
  TTreeReaderArray<Float_t> pzArray(reader, "ReconstructedParticles.p.z");
  TTreeReaderArray<Float_t> energyArray(reader, "ReconstructedParticles.energy");
  TTreeReaderArray<Int_t> pdgArray(reader, "mcparticles2.pdgID");

  /*  TTree *t1 = new TTree("kaons_to_make_D0","tTREE");
  TTree *t2 = new TTree("pions_to_make_D0","tTREE");

  Double_t K_px,K_py,K_pz,K_energy,Pi_px,Pi_py,Pi_pz,Pi_energy;
  
  //Kaons
  TBranch *b_Kpx=t1->Branch("K_px",&K_px,"K_px/F");
  TBranch *b_Kpy=t1->Branch("K_py",&K_py,"K_py/F");
  TBranch *b_Kpz=t1->Branch("K_pz",&K_pz,"K_pz/F");
  TBranch *b_Kenergy=t1->Branch("K_energy",&K_energy,"K_energy/F");

  //Pions
  TBranch *b_Pipx=t2->Branch("Pi_px",&Pi_px,"Pi_px/F");
  TBranch *b_Pipy=t2->Branch("Pi_py",&Pi_py,"Pi_py/F");
  TBranch *b_Pipz=t2->Branch("Pi_pz",&Pi_pz,"Pi_pz/F");
  TBranch *b_Pienergy=t2->Branch("Pi_energy",&Pi_energy,"Pi_energy/F");*/
  
  TH1F *h1 = new TH1F("K_px","",100, -20, 20);
  h1->GetXaxis()->SetTitle("K_px");
  h1->GetXaxis()->SetTitleSize(0.05);
  h1->GetXaxis()->CenterTitle();
  h1->SetLineWidth(2);

  TH1F *h2 = new TH1F("K_py","",100, -20, 20);
  h2->GetXaxis()->SetTitle("K_py");
  h2->GetXaxis()->SetTitleSize(0.05);
  h2->GetXaxis()->CenterTitle();
  h2->SetLineWidth(2);

  TH1F *h3 = new TH1F("K_pz","",100, -20, 20);
  h3->GetXaxis()->SetTitle("K_pz");
  h3->GetXaxis()->SetTitleSize(0.05);
  h3->GetXaxis()->CenterTitle();
  h3->SetLineWidth(2);

  TH1F *h4 = new TH1F("K_energy","",100, 0, 80);
  h4->GetXaxis()->SetTitle("K_energy");
  h4->GetXaxis()->SetTitleSize(0.05);
  h4->GetXaxis()->CenterTitle();
  h4->SetLineWidth(2);
  /////////  

  TH1F *h5 = new TH1F("pi_px","",100, -20, 20);
  h5->GetXaxis()->SetTitle("pi_px");
  h5->GetXaxis()->SetTitleSize(0.05);
  h5->GetXaxis()->CenterTitle();
  h5->SetLineWidth(2);

  TH1F *h6 = new TH1F("pi_py","",100, -20, 20);
  h6->GetXaxis()->SetTitle("pi_py");
  h6->GetXaxis()->SetTitleSize(0.05);
  h6->GetXaxis()->CenterTitle();
  h6->SetLineWidth(2);

  TH1F *h7 = new TH1F("pi_pz","",100, -20, 20);
  h7->GetXaxis()->SetTitle("pi_pz");
  h7->GetXaxis()->SetTitleSize(0.05);
  h7->GetXaxis()->CenterTitle();
  h7->SetLineWidth(2);

  TH1F *h8 = new TH1F("pi_energy","",100, 0, 80);
  h8->GetXaxis()->SetTitle("pi_energy");
  h8->GetXaxis()->SetTitleSize(0.05);
  h8->GetXaxis()->CenterTitle();
  h8->SetLineWidth(2);


  while (reader.Next())
    {
      for (auto pdg: pdgArray)
	{
	  if(pdg==-321)
	    {	  
	      for (auto K_px: pxArray)
		h1->Fill(K_px);
	      for (auto K_py: pyArray)
		h2->Fill(K_py);
	      for (auto K_pz: pzArray)
		h3->Fill(K_pz);
	      for (auto K_energy: energyArray)
		h4->Fill(K_energy);
	    }
	  if(pdg==211)
            {
              for (auto pi_px: pxArray)
                h5->Fill(pi_px);
              for (auto pi_py: pyArray)
                h6->Fill(pi_py);
	      for (auto pi_pz: pzArray)
                h7->Fill(pi_pz);
              for (auto pi_energy: energyArray)
		h8->Fill(pi_energy); 
            }

	}
    }

  TFile * f1 = new TFile("output.root","RECREATE");
  h1->Write();
  h2->Write();
  h3->Write();
  h4->Write();
  h5->Write();
  h6->Write();
  h7->Write();
  h8->Write();

  f1->Close(); 
}