{  
  char infile[128], datafile[128];
  sprintf(infile,"/home/elina/chicos/sim/iron/1MeV/energies.txt");
  FILE *fenergy = fopen(infile, "r");
  sprintf(datafile,"/home/elina/chicos/sim/iron/1MeV/sim.txt");
  FILE *fdata = fopen(datafile, "r");
  
  Int_t ndata=4000; // n times 2
  Double_t nenergies1=678472;
  char fileline[10];
  TNtupleD* simin1 = new TNtupleD("simin1","simin","i:LogE:x:y:z:phi:theta:chisq:nsites");
  TNtupleD* simout1 = new TNtupleD("simout1","simout","i:LogE:x:y:z:phi:theta:chisq:nsites");
  
  TH1F* inputspec1 = new TH1F("inputspec1","Input Spectrum (Full Area)",15,16,19);
  inputspec1->SetLineWidth(2);
  inputspec1->SetLineColor(kRed);
  inputspec1->SetMarkerStyle(20);
  inputspec1->SetMarkerColor(kRed);
  inputspec1->Sumw2();
  inputspec1->SetXTitle("Log(E)");
  inputspec1->SetYTitle("Number of Showers");
  
  TH1F* datain1 = new TH1F("datain1","Showers which pass the cut (Full Area)",15,16,19);
  datain1->SetLineWidth(2);
  datain1->SetLineColor(kRed);
  datain1->SetMarkerStyle(20);
  datain1->SetMarkerColor(kRed);
  datain1->Sumw2();
  datain1->SetXTitle("Log(E)");
  datain1->SetYTitle("Number of Showers");
  
  TH1F* dataout1 = new TH1F("dataout1","Reconstructed showers (Full Area)",15,16,19);
  dataout1->SetLineWidth(2);
  dataout1->SetLineColor(kRed);
  dataout1->SetMarkerStyle(20);
  dataout1->SetMarkerColor(kRed);
  dataout1->Sumw2();
  dataout1->SetXTitle("Log(E)");
  dataout1->SetYTitle("Number of Showers");
  
  Int_t i=0,k=0;
  Double_t LogEnergy1[nenergies1],LogE;
  
  while( fgets(fileline,7,fenergy) ) {
    sscanf(fileline,"%lf",LogEnergy1+i);
    ++i;
  }
  cout<<i<<endl;
  fclose(fenergy);
  
  for(i=0;i<nenergies1;i++){
    LogE = LogEnergy1[i];
    inputspec1->Fill(LogE);
  }
  
  char dataline[500]; 
  i=0;k=0;
  Double_t energy[ndata], x[ndata], y[ndata], z[ndata], theta[ndata], phi[ndata], chisq[ndata], nsites[ndata];
  while( fgets(dataline,499,fdata) && i<ndata ) {
    sscanf(dataline,"%lf %lf %lf %lf %lf %lf %lf %lf",
           energy+i,x+i,y+i,z+i,theta+i,phi+i,chisq+i,nsites+i);
    if(energy[i]>0){
      ++i;
    }
    if(i==k+2){ fgets(dataline,499,fdata); k=i; }
  }
  cout<<i<<endl;
  fclose(fdata);
  
  Double_t j=1;
  for(i=0;i<ndata;i+=2){
    simin1->Fill((Double_t)j, TMath::Log10(energy[i]), 
                 x[i], y[i], z[i], theta[i], phi[i], chisq[i], nsites[i]);
    datain1->Fill(TMath::Log10(energy[i]));
    simout1->Fill((Double_t)j, TMath::Log10(energy[i+1]), 
                  x[i+1], y[i+1], z[i+1], theta[i+1], phi[i+1], chisq[i+1], nsites[i]);
    dataout1->Fill(TMath::Log10(energy[i+1]));
    j++;
  }  

  TH1F* aperture1 = new TH1F("aperture1","Aperture (Full Area)",15,16,19);
  aperture1->SetLineWidth(2);
  aperture1->SetLineColor(kRed);
  aperture1->SetMarkerStyle(20);
  aperture1->SetMarkerColor(kRed);
  aperture1->Sumw2();
  aperture1->SetXTitle("Log(E)");
  aperture1->SetYTitle("N(out)/N(in)");
  
  aperture1->Divide(dataout1,inputspec1);
  
  TH1F* LogA1 = new TH1F("LogA1","LogA (Full Area)",15,16,19);
  LogA1->SetLineWidth(2);
  LogA1->SetLineColor(kRed);
  LogA1->SetMarkerStyle(20);
  LogA1->SetMarkerColor(kRed);
  LogA1->Sumw2();
  LogA1->SetXTitle("Log(E)");
  LogA1->SetYTitle("Aperture [m^2 sr]");
      for(i=1;i<16;i++){
    if(aperture1->GetBinContent(i)>0){
      LogA1->SetBinContent(i,TMath::Log10(aperture1->GetBinContent(i))+7.099);
    }else{
      LogA1->SetBinContent(i,0);
    }
  }
        
  //  TF1* specfn = new TF1("specfn","[0]+[1]*(x-17)",16,18);
}