#include "veto_psim.h"

#define RAD2DEG (180.0/3.14159)
//
// Histogram filling routine. This gets called from the  main loop for every event.                           
//
error_t hist_fill(Banks_t *banks)
{
  int i, Ngam, Ngam_det, Ngam_veto;
  float enr, px, py, pz, mass, mom;
  float theta, phi, x, y;
  float Mpi0, Ppi0;

  TVector3 V3[2];
  TLorentzVector Gamma[2], pi0;

  primMCPART_t *mcpart = banks->MCPART;

  Ngam = 0;
  Ngam_det = 0;
  Ngam_veto = 0;

  Npart->Fill(mcpart->bank.nrow);
  if(mcpart->bank.nrow){
    for(i=0;i<mcpart->bank.nrow;i++){
      MCstate->Fill(mcpart->mcpart[i].state);
      MCtype->Fill(mcpart->mcpart[i].part.type);
      enr = mcpart->mcpart[i].part.v.E;
      px = mcpart->mcpart[i].part.v.px;
      py = mcpart->mcpart[i].part.v.py;
      pz = mcpart->mcpart[i].part.v.pz;
      mass = mcpart->mcpart[i].part.m;
      mom = mcpart->mcpart[i].part.p;
      theta = mcpart->mcpart[i].part.theta;
      phi = mcpart->mcpart[i].part.phi;

      cout<<i<<" "<<mcpart->mcpart[i].state;
      cout<<" "<<mcpart->mcpart[i].part.type;
      cout<<" "<<phi<<" "<<theta<<" "<<mass<<" "<<mom;
      cout<<" "<<enr<<" "<<px<<" "<<py<<" "<<pz<<endl;

      if(mcpart->mcpart[i].part.type==1){
	MCe->Fill(enr);
	MCpx->Fill(px);
	MCpy->Fill(py);
	MCpz->Fill(pz);

	x = Z_TGT2VETO*tan(theta)*cos(phi);
	y = Z_TGT2VETO*tan(theta)*sin(phi);
	V3[Ngam].SetXYZ(x,y,Z_TGT2VETO);
	hit_xy[0]->Fill(x,y);
	if(fabs(x)<=60.0 && fabs(y)<=60.0){
	  hit_xy[1]->Fill(x,y);
	  Ngam_det++;
	  if(drand48()<=0.01) Ngam_veto++;
	}

	Gamma[Ngam].SetPxPyPzE(px,py,pz,enr);
	Ngam++;
	//       	if(Ngam==2) break;
      }
    }
  }
  //  cout<<"Number of photons "<<Ngam<<endl;

  if(Ngam==2){
    pi0 = Gamma[0]+Gamma[1];
    Mpi0 = pi0.M();
    Ppi0 = pi0.P();

    //    cout<<"Pi0 "<<Mpi0<<" "<<Ppi0<<endl;
    pi0Mass->Fill(Mpi0);
    pi0Mom->Fill(Ppi0);

    Ngen++;
    if(Ngam_det==2){
      Ndet++;
      if(Ngam_veto==0) N0veto++;
      if(Ngam_veto==1) N1veto++;
      if(Ngam_veto==2) N2veto++;
    }
  }
  return NOERROR;
}