/*Reads CORSIKA longitudinal files*/

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

FILE *inputfile;
FILE *outputfile;		

char line[500];
char ch;
char str_file_in[50];
char str_file_out[50];
char str1[50],str2[50],str3[50],str4[50],str5[50],str6[50],str7[50],str8[50],str9[50],str10[50],str11[50],str12[50];

int ishower=0,nsteps;
float step;

float depth[205],gammas[205],positrons[205],electrons[205],mu_plus[205],mu_minus[205],hadrons[205],charged[205],nuclei[205],cherenkov[250];
float sgammas[205]={205*0.},spositrons[205]={205*0.},selectrons[205]={205*0.},smu_plus[205]={205*0.},smu_minus[205]={205*0.},shadrons[205]={205*0.},scharged[205]={205*0.},snuclei[205]={205*0.},scherenkov[250]={205*0.};
float dgammas[205]={205*0.},dpositrons[205]={205*0.},delectrons[205]={205*0.},dmu_plus[205]={205*0.},dmu_minus[205]={205*0.},dhadrons[205]={205*0.},dcharged[205]={205*0.},dnuclei[205]={205*0.},dcherenkov[250]={205*0.};
float ex[205]={205*0.};

float depth_ed[205],gamma_ed[300],em_ion[300],em_cut[300],mu_ion[300],mu_cut[300],had_ion[300],had_cut[300],neutrino[300],sum[300];

int i;

float p1,p2,p3,p4,p5,p6;
float chi2dof,avdev;

main()
{

printf("Enter DATxxxxxx file code (6 digits)\n");
gets(str1);
strcpy(str_file_in,"DAT");
strcat(str_file_in,str1);
strcat(str_file_in,".long");
strcpy(str_file_out,"OUT");
strcat(str_file_out,str1);
strcat(str_file_out,str1);

inputfile=fopen(str_file_in,"r");
if (inputfile==NULL) {
  (void)printf("Can't open inputfile %s\n",str_file_in);
  return 0;}
 
outputfile=fopen(str_file_out,"w");
if (outputfile==NULL) {
  (void)printf("Can't open outputfile %s\n",str_file_out);
  return 0;}

(void)fprintf(outputfile,"# %s\n",str_file_in);

while (1){
  
  ch=fgetc(inputfile);
  if (ch==EOF){
	(void)printf("%d showers analyzed\n",ishower);  
	(void)fprintf(outputfile,"# EOF\n %d showers analyzed\n",ishower);  
	break;}
  (void)ungetc(ch,inputfile);

  (void)fgets(line,sizeof(line),inputfile);
  (void)sscanf(line,"%s %s %s %s %s %s %s %s %s",&str1,&str2,&str3,&str4,&str5,&str6,&str7,&str8,&str9);

//	if (ishower<=10){(void)printf("%s \n",str1);}


  if ((strcmp(str1,"LONGITUDINAL")==0)&&(strcmp(str2,"DISTRIBUTION")==0)&&(strcmp(str3,"IN")==0)){
	nsteps=atoi(str4);
  	step=atof(str8);
 	//ishower=atoi(str12);
	ishower++;
	//if (ishower<=10){(void)printf("nsteps=%d step=%f ishower=%d \n",nsteps,step,ishower);}
	(void)fprintf(outputfile,"# SHOWER %d\n# Longitudinal particle distribution for %d levels:\n",ishower,nsteps);
     	(void)fgets(line,sizeof(line),inputfile);   
   	for(i=0;i<nsteps;i++) { 
 		(void)fgets(line,sizeof(line),inputfile);   
		(void)sscanf(line,"%f %f %f %f %f %f %f %f %f %f",&depth[i],&gammas[i],&positrons[i],&electrons[i],&mu_plus[i],&mu_minus[i],&hadrons[i],&charged[i],&nuclei[i],&cherenkov[i]);
		(void)fprintf(outputfile,"%.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f\n",depth[i],gammas[i],positrons[i],electrons[i],mu_plus[i],mu_minus[i],hadrons[i],charged[i],nuclei[i],cherenkov[i]);
        sgammas[i]+=gammas[i],spositrons[i]+=positrons[i],selectrons[i]+=electrons[i],smu_plus[i]+=mu_plus[i],smu_minus[i]+=mu_minus[i],shadrons[i]+=hadrons[i],scharged[i]+=charged[i],snuclei[i]+=nuclei[i],scherenkov[i]+=cherenkov[i];
		//if (ishower<=10){(void)printf("%.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f\n",depth[i],gammas[i],positrons[i],electrons[i],mu_plus[i],mu_minus[i],hadrons[i],charged[i],nuclei[i],cherenkov[i]);}			
				}								}

  if ((strcmp(str1,"LONGITUDINAL")==0)&&(strcmp(str2,"ENERGY")==0)&&(strcmp(str3,"DEPOSIT")==0)){
	nsteps=atoi(str5);
  	step=atof(str9);
 	//ishower=atoi(str14);
	//if (ishower<=10){(void)printf("str5=%s, str9=%s, str14=%s,\n nsteps=%d step=%.0f ishower=%d \n",str5,str9,str14,nsteps,step,ishower);}
	(void)fprintf(outputfile,"# Energy deposit distribution for %d levels:\n",nsteps);
     	(void)fgets(line,sizeof(line),inputfile);   
   	for(i=0;i<nsteps;i++) { 
	    	(void)fgets(line,sizeof(line),inputfile);   
		(void)sscanf(line,"%f %f %f %f %f %f %f %f %f %f",&depth_ed[i],&gamma_ed[i],&em_ion[i],&em_cut[i],&mu_ion[i],&mu_cut[i],&had_ion[i],&had_cut[i],&neutrino[i],&sum[i]);
		(void)fprintf(outputfile,"%.0f %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e\n",depth_ed[i],gamma_ed[i],em_ion[i],em_cut[i],mu_ion[i],mu_cut[i],had_ion[i],had_cut[i],neutrino[i],sum[i]);
		//if (ishower<=10){(void)printf("%.0f %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n",depth_ed[i],gamma[i],em_ion[i],em_cut[i],mu_ion[i],mu_cut[i],had_ion[i],had_cut[i],neutrino[i],sum[i]);}
	   			}								}

  if ((strcmp(str1,"FIT")==0)&&(strcmp(str2,"OF")==0)&&(strcmp(str3,"THE")==0)&&(strcmp(str4,"HILLAS")==0)) {
//  if (ishower<=10){(void)printf("str1=%s, str2=%s, str3=%s, str4=%s\n",str1,str2,str3,str4);}  
 	(void)fgets(line,sizeof(line),inputfile);
   	(void)fgets(line,sizeof(line),inputfile);
	(void)sscanf(line,"%s %s %s %s %s %s %s %s",&str1,&str2,&str3,&str4,&str5,&str6,&str7,&str8);
  	if (strcmp(str1,"PARAMETERS")==0) {
		p1=atof(str3);
		p2=atof(str4);
		p3=atof(str5);
		p4=atof(str6);
		p5=atof(str7);
		p6=atof(str8);
		}
   	(void)fgets(line,sizeof(line),inputfile);
	(void)sscanf(line,"%s %s %s",&str1,&str2,&str3);
  	if (strcmp(str1,"CHI**2/DOF")==0) {chi2dof=atof(str3);}
   	(void)fgets(line,sizeof(line),inputfile);
	(void)sscanf(line,"%s %s %s %s %s %s",&str1,&str2,&str3,&str4,&str5,&str6);
  	if (strcmp(str1,"AV.")==0) {avdev=atof(str6);}
	(void)fprintf(outputfile,"# Hillas parameters: %.5e %.5e %.5e %.5e %.5e %.5e\n# chi2/dof= %.5e avdev= %.5e\n# End of shower %d.\n",p1,p2,p3,p4,p5,p6,chi2dof,avdev,ishower);
													} 

/*end while*/}

//average shower
   	for(i=0;i<nsteps;i++) {
        sgammas[i]=sgammas[i]/ishower,spositrons[i]=spositrons[i]/ishower,selectrons[i]=selectrons[i]/ishower,smu_plus[i]=smu_plus[i]/ishower,smu_minus[i]=smu_minus[i]/ishower,shadrons[i]=shadrons[i]/ishower,scharged[i]=scharged[i]/ishower,snuclei[i]=snuclei[i]/ishower,scherenkov[i]=scherenkov[i]/ishower;
                            }

(void)fclose(inputfile);
(void)fclose(outputfile);

//std.dev.
    inputfile=fopen(str_file_in,"r");
    if (inputfile==NULL) {
        (void)printf("Can't open inputfile %s\n",str_file_in);
        return 0;}
    ishower=0;
    while (1){
        
        ch=fgetc(inputfile);
        if (ch==EOF){
            (void)printf("std. dev. of %d showers\n",ishower);
            //(void)fprintf(outputfile,"# EOF\n %d showers analyzed\n",ishower);
            break;}
        (void)ungetc(ch,inputfile);
        
        (void)fgets(line,sizeof(line),inputfile);
        (void)sscanf(line,"%s %s %s %s %s %s %s %s %s",&str1,&str2,&str3,&str4,&str5,&str6,&str7,&str8,&str9);
        
        //	if (ishower<=10){(void)printf("%s \n",str1);}
        
        
        if ((strcmp(str1,"LONGITUDINAL")==0)&&(strcmp(str2,"DISTRIBUTION")==0)&&(strcmp(str3,"IN")==0)){
            nsteps=atoi(str4);
            step=atof(str8);
            //ishower=atoi(str12);
            ishower++;
            //if (ishower<=3){(void)printf("nsteps=%d step=%f ishower=%d \n",nsteps,step,ishower);}
            //(void)fprintf(outputfile,"# SHOWER %d\n# Longitudinal particle distribution for %d levels:\n",ishower,nsteps);
            (void)fgets(line,sizeof(line),inputfile);
            for(i=0;i<nsteps;i++) {
                (void)fgets(line,sizeof(line),inputfile);
                (void)sscanf(line,"%f %f %f %f %f %f %f %f %f %f",&depth[i],&gammas[i],&positrons[i],&electrons[i],&mu_plus[i],&mu_minus[i],&hadrons[i],&charged[i],&nuclei[i],&cherenkov[i]);
                //(void)fprintf(outputfile,"%.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f\n",depth[i],gammas[i],positrons[i],electrons[i],mu_plus[i],mu_minus[i],hadrons[i],charged[i],nuclei[i],cherenkov[i]);
                dgammas[i]+=(sgammas[i]-gammas[i])*(sgammas[i]-gammas[i]);
                dpositrons[i]+=(spositrons[i]-positrons[i])*(spositrons[i]-positrons[i]);
                delectrons[i]+=(selectrons[i]-electrons[i])*(selectrons[i]-electrons[i]);
                dmu_plus[i]+=(smu_plus[i]-mu_plus[i])*(smu_plus[i]-mu_plus[i]);
                dmu_minus[i]+=(smu_minus[i]-mu_minus[i])*(smu_minus[i]-mu_minus[i]);
                dhadrons[i]+=(shadrons[i]-hadrons[i])*(shadrons[i]-hadrons[i]);
                dcharged[i]+=(scharged[i]-charged[i])*(scharged[i]-charged[i]);
                dnuclei[i]+=(snuclei[i]-nuclei[i])*(snuclei[i]-nuclei[i]);
                dcherenkov[i]+=(scherenkov[i]-cherenkov[i])*(scherenkov[i]-cherenkov[i]);
                //if (ishower<=10){(void)printf("%.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f %.0f\n",depth[i],gammas[i],positrons[i],electrons[i],mu_plus[i],mu_minus[i],hadrons[i],charged[i],nuclei[i],cherenkov[i]);}			
            }								}
    }//end of while

    for(i=0;i<nsteps;i++) {
        dgammas[i]=sqrt(dgammas[i]/(ishower-1));
        dpositrons[i]=sqrt(dpositrons[i]/(ishower-1));
        delectrons[i]=sqrt(delectrons[i]/(ishower-1));
        dmu_plus[i]=sqrt(dmu_plus[i]/(ishower-1));
        dmu_minus[i]=sqrt(dmu_minus[i]/(ishower-1));
        dhadrons[i]=sqrt(dhadrons[i]/(ishower-1));
        dcharged[i]=sqrt(dcharged[i]/(ishower-1));
        dnuclei[i]=sqrt(dnuclei[i]/(ishower-1));
        dcherenkov[i]=sqrt(dcherenkov[i]/(ishower-1));
    }

return 0;
}