True = new TH1D ("True", "True", Nbins, Bin_min, Bin_max);
Reco = new TH1D ("Reco", "Reco", Nbins, Bin_min, Bin_max);

          TFile *t=new TFile("Input.root");
	  
	  TTree *truthTuple = (TTree*)t->Get("truthTuple");
	  TTree *recoTuple  = (TTree*)t->Get("recoTuple");

	  // Projection of Real Statistics Files/* IT HAS TO BE DONE WITH THE REAL STATISTICS FILE*/ 
	  truthTuple->Project("True", "variable", "eppid>3");
          recoTuple->Project("Reco", "(sqrt((feppe+fempe)**2-(feppx+fempx)**2-(feppy+fempy)**2-(feppz+fempz)**2))**2", "eppid>3");
	  

int TotalEntries=0;
for (int i=0; i <Nbins;i++ )
{
TotalEntries=TotalEntries+Reco->GetBinContent(i+1);// histograms start on bin=1 not bin=0
}
cout << "Total Reco Entries: "<< TotalEntries<< endl;

	  // Variables to be read and filled recoTuple
	  float feppx;
	  float feppy;
	  float feppz;
	  float feppe;
	  float fepp3;
	  float fepcosth;
	  int eppid, ntrk;
	      
	  recoTuple->SetBranchAddress("feppx", &feppx);
	  recoTuple->SetBranchAddress("feppy", &feppy);
	  recoTuple->SetBranchAddress("feppz", &feppz);
	  recoTuple->SetBranchAddress("feppe", &feppe);
	  recoTuple->SetBranchAddress("fepp3", &fepp3);
	  recoTuple->SetBranchAddress("epcosth", &fepcosth);
	  recoTuple->SetBranchAddress("eppid", &eppid);
	  recoTuple->SetBranchAddress("ntrk", &ntrk);
	     
	  float fempx;
	  float fempy;
	  float fempz;
	  float fempe;
	  float femp3;
	  float femcosth;
	  int empid;
	      
	  recoTuple->SetBranchAddress("fempx", &fempx);
	  recoTuple->SetBranchAddress("fempy", &fempy);
	  recoTuple->SetBranchAddress("fempz", &fempz);
	  recoTuple->SetBranchAddress("fempe", &fempe);
	  recoTuple->SetBranchAddress("femp3", &femp3);
	  recoTuple->SetBranchAddress("emcosth", &femcosth);
	  recoTuple->SetBranchAddress("empid", &empid);


	  float fpi0px;
	  float fpi0py;
	  float fpi0pz;
	  float fpi0pe;
	    	      
	  recoTuple->SetBranchAddress("fpi0px", &fpi0px);
	  recoTuple->SetBranchAddress("fpi0py", &fpi0py);
	  recoTuple->SetBranchAddress("fpi0pz", &fpi0pz);
	  recoTuple->SetBranchAddress("fpi0pe", &fpi0pe);

	  
	  TH1D * RecoFill;
	  RecoFill = new TH1D("RecoFill",  "RecoFill", Nbins, Bin_max, Bin_max);
	
	  long NEntriesReco=(long)recoTuple->GetEntries();
	  int kkk=0;

	  for (int k=0; k<NEntriesReco; k++)
	    {
		  
	      recoTuple->GetEntry(k);
	

	      if( eppid>3)
		{
		  
	      	  if (kkk % 10000 == 0 && kkk != 0) 
	      	    {      	
	      	      cout<<"*** FILLING *** "<< kkk << " : " <<endl;
	      	    }
		  		  
	      	  kkk++;
	
	      	  RecoFill->Fill((sqrt((feppe+fempe)**2-(feppx+fempx)**2-(feppy+fempy)**2-(feppz+fempz)**2))**2);
	      	}
	    }

	  cout<< "RecoFill: "<<	  RecoFill->GetEntries()<<endl;