/*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[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;

read_long()
{

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[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[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));
    }
    
//Root
    TCanvas *c1 = new TCanvas("long","Longitudinal distribution",10,10,500,500);

    const int n=nsteps;
    
    gr_gam=new TGraph(n,depth,gammas);
    gr_ele=new TGraph(n,depth,electrons);
    gr_pos=new TGraph(n,depth,positrons);
    gr_mup=new TGraph(n,depth,mu_plus);
    gr_mun=new TGraph(n,depth,mu_minus);
    gr_had=new TGraph(n,depth,hadrons);
    gr_cha=new TGraph(n,depth,charged);
    gr_nuc=new TGraph(n,depth,nuclei);
    gr_che=new TGraph(n,depth,cherenkov);

    gr_mgam=new TGraphErrors(n,depth,sgammas,ex,dgammas);
    gr_mele=new TGraphErrors(n,depth,selectrons,ex,delectrons);
    gr_mpos=new TGraphErrors(n,depth,spositrons,ex,dpositrons);
    gr_mmup=new TGraphErrors(n,depth,smu_plus,ex,dmu_plus);
    gr_mmun=new TGraphErrors(n,depth,smu_minus,ex,dmu_minus);
    gr_mhad=new TGraphErrors(n,depth,shadrons,ex,dhadrons);
    gr_mcha=new TGraphErrors(n,depth,scharged,ex,dcharged);
    gr_mnuc=new TGraphErrors(n,depth,snuclei,ex,dnuclei);
    gr_mche=new TGraphErrors(n,depth,scherenkov,ex,dcherenkov);
    
    gr_gam->SetTitle("n_{#\gamma} x depth");
    gr_ele->SetTitle("n_{e^{-}} x depth");
    gr_pos->SetTitle("n_{e^{+} } x depth");
    gr_mup->SetTitle("n_{#\mu^{+} } x depth");
    gr_mun->SetTitle("n_{#\mu^{-}} x depth");
    gr_had->SetTitle("n_{hadrons} x depth");
    gr_cha->SetTitle("n_{charged} x depth");
    gr_nuc->SetTitle("n_{nuclei} x depth");
    gr_che->SetTitle("n_{cherenkov} x depth");

    gr_mgam->SetTitle("n_{#\gamma}^{av} x depth");
    gr_mele->SetTitle("n_{e^{-}}^{av} x depth");
    gr_mpos->SetTitle("n_{e^{+}}^{av} x depth");
    gr_mmup->SetTitle("n_{#\mu^{+}}^{av} x depth");
    gr_mmun->SetTitle("n_{#\mu^{-}}^{av} x depth");
    gr_mhad->SetTitle("n_{hadrons}^{av} x depth");
    gr_mcha->SetTitle("n_{charged}^{av} x depth");
    gr_mnuc->SetTitle("n_{nuclei}^{av} x depth");
    gr_mche->SetTitle("n_{cherenkov}^{av} x depth");
    
    gr_gam->SetMarkerStyle(7);
    gr_ele->SetMarkerStyle(7);
    gr_pos->SetMarkerStyle(7);
    gr_mup->SetMarkerStyle(7);
    gr_mun->SetMarkerStyle(7);
    gr_had->SetMarkerStyle(7);
    gr_cha->SetMarkerStyle(7);
    gr_nuc->SetMarkerStyle(7);
    gr_che->SetMarkerStyle(7);
    
    gr_mgam->SetMarkerStyle(7);
    gr_mele->SetMarkerStyle(7);
    gr_mpos->SetMarkerStyle(7);
    gr_mmup->SetMarkerStyle(7);
    gr_mmun->SetMarkerStyle(7);
    gr_mhad->SetMarkerStyle(7);
    gr_mcha->SetMarkerStyle(7);
    gr_mnuc->SetMarkerStyle(7);
    gr_mche->SetMarkerStyle(7);
    
    gr_gam->GetXaxis()->SetTitle("g/cm^{2}");
    gr_ele->GetXaxis()->SetTitle("g/cm^{2}");
    gr_pos->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mup->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mun->GetXaxis()->SetTitle("g/cm^{2}");
    gr_had->GetXaxis()->SetTitle("g/cm^{2}");
    gr_cha->GetXaxis()->SetTitle("g/cm^{2}");
    gr_nuc->GetXaxis()->SetTitle("g/cm^{2}");
    gr_che->GetXaxis()->SetTitle("g/cm^{2}");
    
    gr_mgam->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mele->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mpos->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mmup->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mmun->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mhad->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mcha->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mnuc->GetXaxis()->SetTitle("g/cm^{2}");
    gr_mche->GetXaxis()->SetTitle("g/cm^{2}");
    
    gr_gam->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("gammas.pdf");
    
    gr_ele->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("electrons.pdf");

    gr_pos->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("positrons.pdf");
    
    gr_mup->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("muons_plus.pdf");
    
    gr_mun->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("muons_minus.pdf");

    gr_had->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("hadrons.pdf");

    gr_cha->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("charged.pdf");

    gr_nuc->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("nuclei.pdf");
    
/*    gr_che->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("cherenkov.pdf");
*/


    gr_mgam->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("gammas_av.pdf");
    
    gr_mele->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("electrons_av.pdf");
    
    gr_mpos->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("positrons_av.pdf");
    
    gr_mmup->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("muons_plus_av.pdf");
    
    gr_mmun->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("muons_minus_av.pdf");
    
    gr_mhad->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("hadrons_av.pdf");
    
    gr_mcha->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("charged_av.pdf");
    
    gr_mnuc->Draw("AP");
    //TLatex text(23,1000,"primary: proton");
    //text.DrawClone();
    c1->Update();
    c1.SaveAs("nuclei_av.pdf");
    
    /*    gr_mche->Draw("AP");
     //TLatex text(23,1000,"primary: proton");
     //text.DrawClone();
     c1->Update();
     c1.SaveAs("cherenkov_av.pdf");
     */

    
return 0;
}