#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.14159265358979323844
#define nimt 50    // number of mt bins
#define delmt 20.0  // bin size for mt
#define dely 1.0   // make spectra from y=0 +/- dely/2
#define delycorr 10.0// Correlation file only includes pions in this rap. range
//#define rap0 2.92  // Rapidity to be uses as p=0 in correlation program
int main(){
  double dndp_pi[nimt],dndp_nucl[nimt],dndp_k[nimt];
  double px,py,pz,mt,x,y,z,t,mass,rapidity,e;
  double eperpbar_pi,eperpbar_nucl,eperpbar_k,teffective;
  double bimp,phi;
  char cdummy[200];
  int imt,id,i,ipart,npart,nevents,ievent,n_pi,n_nucl,n_k;
  FILE *fptrin,*fptrcorr,*fptrspectra;
  scanf("%d",&nevents);

  for(imt=0;imt<nimt;imt++) {
    dndp_pi[imt]=dndp_nucl[imt]=dndp_k[imt]=0.0;
  }
  eperpbar_pi=eperpbar_k=eperpbar_nucl=0.0;
  n_pi=n_nucl=n_k=0;

  fptrin=fopen("results/phasespace.dat","r");
  /* read in header, and ignore it */
  for(i=0;i<3;i++){
    fgets(cdummy,200,fptrin);
  }

  for(ievent=1;ievent<=nevents;ievent++){
    fscanf(fptrin,"%d %d %lf %lf",&i,&npart,&bimp,&phi);

    for(ipart=1;ipart<=npart;ipart++){
      fscanf(fptrin,"%d %d %lf %lf %lf %lf",&i,&id,&px,&py,&pz,&e);
      fscanf(fptrin," %lf %lf %lf %lf %lf",&mass,&x,&y,&z,&t);
      px=px*1000.0;py=py*1000.0,pz=pz*1000.0;e=e*1000.0;mass=mass*1000.0;
      mt=sqrt(mass*mass+px*px+py*py);
      imt=(int)((mt-mass)/delmt);
      rapidity=.5*log((e+pz)/(e-pz));
      if(fabs(rapidity)<0.5){
	if(id==211||id==-211||id==111){
	  n_pi=n_pi+1;
	  eperpbar_pi=eperpbar_pi+mt;
	  if(imt<nimt&&fabs(rapidity)<0.5*dely) dndp_pi[imt]=dndp_pi[imt]+1.0;
	}
	if(id==2112||id==2212){
	  n_nucl=n_nucl+1;
	  eperpbar_nucl=eperpbar_nucl+mt;
	  if(imt<nimt&&fabs(rapidity)<0.5*dely) 
	    dndp_nucl[imt]=dndp_nucl[imt]+1.0;
	}
	if(id==311||id==-311||id==321||id==-321){
	  n_k=n_k+1;
	  eperpbar_k=eperpbar_k+mt;
	  if(imt<nimt&&fabs(rapidity)<0.5*dely) dndp_k[imt]=dndp_k[imt]+1.0;
	}
      }
    }
  }
  fclose(fptrin);
  printf("%d pions %d nucleons %d kaons\n",n_pi,n_nucl,n_k);

  /* write spectra for pions */
  fptrspectra=fopen("results/spectra.pi.dat","w+");
  eperpbar_pi=eperpbar_pi/(double)n_pi;
  mass=137.0;
  printf("For pions, <et>-m=%g\n",eperpbar_pi-mass);
  teffective=0.25*eperpbar_pi-0.5*mass
    +0.25*sqrt(eperpbar_pi*eperpbar_pi+4.0*mass*eperpbar_pi-4.0*mass*mass);
  printf("T_eff for pions  = %g\n",teffective);
  for(imt=0;imt<nimt;imt++){
    mt=mass+0.5*delmt+imt*delmt;
    fprintf(fptrspectra,"%lf   %lf\n",mt-mass,
	    1.0e6*dndp_pi[imt]/(3.0*nevents*2.0*pi*mt*delmt*dely));
  } 
  fclose(fptrspectra);

  /* write spectra for nucleons */
  fptrspectra=fopen("results/spectra.nucl.dat","w+");
  eperpbar_nucl=eperpbar_nucl/(double)n_nucl;
  mass=939.0;
  printf("For nucleons, <et>-m=%g\n",eperpbar_nucl-mass);
  teffective=0.25*eperpbar_nucl-0.5*mass
    +0.25*sqrt(eperpbar_nucl*eperpbar_nucl+4.0*mass*eperpbar_nucl-4.0*mass*mass);
  printf("T_eff for nucleons  = %g\n",teffective);
  for(imt=0;imt<nimt;imt++){
    mt=mass+0.5*delmt+imt*delmt;
    fprintf(fptrspectra,"%lf   %lf\n",mt-mass,
	    1.0e6*dndp_nucl[imt]/(3.0*nevents*2.0*pi*mt*delmt*dely));
  } 
  fclose(fptrspectra);

  /* write spectra for kaons */
  fptrspectra=fopen("results/spectra.k.dat","w+");
  eperpbar_k=eperpbar_k/(double)n_k;
  mass=497.0;
  printf("For kaons, <et>-m=%g\n",eperpbar_k-mass);
  teffective=0.25*eperpbar_k-0.5*mass
    +0.25*sqrt(eperpbar_k*eperpbar_k+4.0*mass*eperpbar_k-4.0*mass*mass);
  printf("T_eff for kaons  = %g\n",teffective);
  for(imt=0;imt<nimt;imt++){
    mt=mass+0.5*delmt+imt*delmt;
    fprintf(fptrspectra,"%lf   %lf\n",mt-mass,
	    1.0e6*dndp_k[imt]/(3.0*nevents*2.0*pi*mt*delmt*dely));
  } 
  fclose(fptrspectra);

  return 0;
}