#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define pi 3.14159265358979323844
#define ntaumax 100

int main(){
  double x,y,z,t,tau,px,py,pz;
  int itau,count_proton[ntaumax]={0},count_pion[ntaumax]={0},dummy;
  FILE *fptrproton,*fptrpion,*fptrlastint;

  fptrproton=fopen("proton.dat","r");
  while(feof(fptrproton)==0){
    fscanf(fptrproton,"%lf %d %lf %lf %lf %lf %lf %lf",
	   &t,&dummy,&px,&py,&pz,&x,&y,&z);
    if(feof(fptrproton)==0){
      if(fabs(t)<fabs(z)){
	printf("tau is imaginary, t=%lf z=%lf\n",t,z);
	exit(1);
      }
      tau=sqrt(t*t-z*z);
      itau=floor(tau);
      if(itau<ntaumax) count_proton[itau]=count_proton[itau]+1;
    }
  }	
  fclose(fptrproton);
  
  fptrpion=fopen("pion.dat","r");
  while(feof(fptrpion)==0){
    fscanf(fptrpion,"%lf %d %lf %lf %lf %lf %lf %lf",
	   &t,&dummy,&px,&py,&pz,&x,&y,&z);
    if(fabs(t)<fabs(z)){
      printf("tau is imaginary, t=%lf z=%lf\n",t,z);
      exit(1);
    }
    tau=sqrt(t*t-z*z);
    itau=floor(tau);
    if(itau<ntaumax&&feof(fptrpion)==0) count_pion[itau]=count_pion[itau]+1;
  }	
  fclose(fptrpion);
  
  fptrlastint=fopen("last_int.dat","w");
  for(itau=0;itau<ntaumax;itau++){
    printf("%9d %9d %9d\n",
	   itau,count_proton[itau],count_pion[itau]);
    fprintf(fptrlastint,"%9d %9d %9d\n",
	    itau,count_proton[itau],count_pion[itau]);
  }
  fclose(fptrlastint);
  
  return(0);
}