{
	TFile* f = new TFile("sim_digi.root");
   	TTree* t = (TTree*)f->Get("er");
   	t->AddFriend("er","reco.root");

	TCanvas *c1 = new TCanvas("c1", "Telescope Coordinate reconstruction", 1400, 1000);
	
	TClonesArray* T1_DoubleSi_points = new TClonesArray("ERQTelescopeSiPoint",100);
	
	TClonesArray* T1_SingleSi_digi = new TClonesArray("ERQTelescopeSiDigi",100);
	TClonesArray* T1_DoubleSi_digiX = new TClonesArray("ERQTelescopeSiDigi",100);
	TClonesArray* T1_DoubleSi_digiY = new TClonesArray("ERQTelescopeSiDigi",100);
	TClonesArray* T1_CsI_digi = new TClonesArray("ERQTelescopeCsIDigi",100);

	TClonesArray* T2_DoubleSi_points = new TClonesArray("ERQTelescopeSiPoint",100);

	TClonesArray* T2_DoubleSi_digiX = new TClonesArray("ERQTelescopeSiDigi",100);
	TClonesArray* T2_DoubleSi_digiY = new TClonesArray("ERQTelescopeSiDigi",100);
	TClonesArray* T2_CsI_digi = new TClonesArray("ERQTelescopeCsIDigi",100);

	
	t->SetBranchAddress("ERQTelescopeSiPoint_T1_DoubleSi_SD2_XY_0_X",&T1_DoubleSi_points);

	t->SetBranchAddress("ERQTelescopeSiDigi_T1_SingleSi_1_X_0",&T1_SingleSi_digi);
	t->SetBranchAddress("ERQTelescopeSiDigi_T1_DoubleSi_SD2_XY_0_X",&T1_DoubleSi_digiX);
	t->SetBranchAddress("ERQTelescopeSiDigi_T1_DoubleSi_SD2_XY_0_Y",&T1_DoubleSi_digiY);
	t->SetBranchAddress("ERQTelescopeCsIDigi_T1_CsI_1_0",&T1_CsI_digi);

	t->SetBranchAddress("ERQTelescopeSiPoint_T2_DoubleSi_SD2_XY_1_X",&T2_DoubleSi_points);

	t->SetBranchAddress("ERQTelescopeSiDigi_T2_DoubleSi_SD2_XY_1_X",&T2_DoubleSi_digiX);
	t->SetBranchAddress("ERQTelescopeSiDigi_T2_DoubleSi_SD2_XY_1_Y",&T2_DoubleSi_digiY);
	t->SetBranchAddress("ERQTelescopeCsIDigi_T2_CsI_1_1",&T2_CsI_digi);

	TH2F* T1deEHe3 = new TH2F("T1deEHe3","T1 de vs E He3",1000,0,0.25,1000,0,0.0014);
	T1deEHe3->GetXaxis()->SetTitle("E");
	T1deEHe3->GetYaxis()->SetTitle("de");
	TH2F* T1deEH3 = new TH2F("T1deEH3","T1 deE H3",1000,0,0.25,1000,0,0.0014);
	TH2F* T1deEOther = new TH2F("T1deEOther","T1 deE Other",1000,0,0.25,1000,0,0.0014);

	TH2F* T2deEHe3 = new TH2F("T2deEHe3","T2 de vs E He3",1000,0,0.25,1000,0,0.04);
	T2deEHe3->GetXaxis()->SetTitle("E");
	T2deEHe3->GetYaxis()->SetTitle("de");
	TH2F* T2deEH3 = new TH2F("T2deEH3","T2 deE H3",1000,0,0.25,1000,0,0.04);
	TH2F* T2deEOther = new TH2F("T2deEOther","T2 deE Other",1000,0,0.25,1000,0,0.04);
	

	TCanvas *c2 = new TCanvas("c2", "Telescope de/E on digi info", 1400, 1000);
	c2->Divide(2,1);

	Long64_t nentries = t->GetEntriesFast();
	//Long64_t nentries = 100;
	for (Long64_t i=0; i<nentries; i++) {
	  t->GetEntry(i);

	  if (T1_DoubleSi_digiX->GetEntriesFast() != 1)
	  	continue;

	  if (T1_SingleSi_digi->GetEntriesFast() != 1)
	  	continue;


	  Double_t de = ((ERQTelescopeSiDigi*)T1_SingleSi_digi->At(0))->GetEdep();
	  Double_t E = ((ERQTelescopeSiDigi*)T1_SingleSi_digi->At(0))->GetEdep() + 
	  			(((ERQTelescopeSiDigi*)T1_DoubleSi_digiX->At(0))->GetEdep()	+
	  			 ((ERQTelescopeSiDigi*)T1_DoubleSi_digiY->At(0))->GetEdep())/2.;
	  
	  for (Int_t iDigi = 0; iDigi < T1_CsI_digi->GetEntriesFast(); iDigi++){
	  	E+= ((ERQTelescopeCsIDigi*)T1_CsI_digi->At(iDigi))->Edep();
	  }

	  Bool_t he3 = kFALSE;
	  Bool_t h3 = kFALSE;
	
	  for (Int_t iPoint = 0; iPoint < T1_DoubleSi_points->GetEntriesFast(); iPoint++){
	  	ERQTelescopeSiPoint* p = (ERQTelescopeSiPoint*)T1_DoubleSi_points->At(iPoint);
	  	if (p->GetPDG() == 1000020030)
	  		he3 = kTRUE;
	  	if (p->GetPDG() == 1000010030)
	  		h3 = kTRUE;
	  }

	  if (he3)
	  	T1deEHe3->Fill(E,de);
	  else if (h3)
	  	T1deEH3->Fill(E,de);
	  else
	  	T1deEOther->Fill(E,de);
	}

	c2->cd(1);
	T1deEHe3->SetMarkerColor(kBlue);
	T1deEHe3->SetMarkerStyle(21);
	T1deEHe3->Draw();
	T1deEH3->SetMarkerColor(kRed);
	T1deEH3->SetMarkerStyle(21);
	T1deEH3->Draw("same");
	T1deEOther->SetMarkerColor(kGreen);
	T1deEOther->SetMarkerStyle(21);
	T1deEOther->Draw("same");

	auto legend = new TLegend(0.7,0.5,0.95,0.7);
	legend->SetHeader("Ions");
	legend->AddEntry(T1deEHe3,"He3","p");
	legend->AddEntry(T1deEH3,"H3","p");
	legend->AddEntry(T1deEOther,"Others","p");
	legend->Draw();

	for (Long64_t i=0; i<nentries; i++) {
	  t->GetEntry(i);

	  if (T2_DoubleSi_digiX->GetEntriesFast() != 1)
	  	continue;

	  Double_t de = (((ERQTelescopeSiDigi*)T2_DoubleSi_digiX->At(0))->GetEdep()	+
	  			 ((ERQTelescopeSiDigi*)T2_DoubleSi_digiY->At(0))->GetEdep())/2.;
	  Double_t E = de;
	  
	  for (Int_t iDigi = 0; iDigi < T2_CsI_digi->GetEntriesFast(); iDigi++){
	  	E+= ((ERQTelescopeCsIDigi*)T2_CsI_digi->At(iDigi))->Edep();
	  }

	  Bool_t he3 = kFALSE;
	  Bool_t h3 = kFALSE;
	
	  for (Int_t iPoint = 0; iPoint < T2_DoubleSi_points->GetEntriesFast(); iPoint++){
	  	ERQTelescopeSiPoint* p = (ERQTelescopeSiPoint*)T2_DoubleSi_points->At(iPoint);
	  	if (p->GetPDG() == 1000020030)
	  		he3 = kTRUE;
	  	if (p->GetPDG() == 1000010030)
	  		h3 = kTRUE;
	  }

	  if (he3)
	  	T2deEHe3->Fill(E,de);
	  else if (h3)
	  	T2deEH3->Fill(E,de);
	  else
	  	T2deEOther->Fill(E,de);
	}

	c2->cd(2);
	T2deEHe3->SetMarkerColor(kBlue);
	T2deEHe3->SetMarkerStyle(21);
	T2deEHe3->Draw();
	T2deEH3->SetMarkerColor(kRed);
	T2deEH3->SetMarkerStyle(21);
	T2deEH3->Draw("same");
	T2deEOther->SetMarkerColor(kGreen);
	T2deEOther->SetMarkerStyle(21);
	T2deEOther->Draw("same");

	legend = new TLegend(0.7,0.5,0.95,0.7);
	legend->SetHeader("Ions");
	legend->AddEntry(T2deEHe3,"He3","p");
	legend->AddEntry(T2deEH3,"H3","p");
	legend->AddEntry(T2deEOther,"Others","p");
	legend->Draw();
}