// this is elina's program to make ldf plots at chiquita energies
// it is intended to be run from high_energy/Proton/
// it creates a plot of the ldf of each energy for a given angle

{
  char angle[10]; sprintf(angle, "0.95"); // 0.95, 0.85, 0.75
  char radians[10]; sprintf(radians, "0.318"); // 0.318, 0.555, 0.723
  char efine[10]; sprintf(efine, "efine_3");
  char mfine[10]; sprintf(mfine, "mfine_3");
  char gfine[10]; sprintf(gfine, "gfine_3");
  char hnum[10]; char filename[100]; char histname[100];
  Int_t numhists=10;

  // Input histograms
  TH1F* he[numhists]; TH1F* hm[numhists]; TH1F* hg[numhists];
  TH1F* hldf[4];

  gStyle->SetOptStat(0);
  sprintf(histname, "LDF at cos %s", angle);
  hldf[0] = new TH1F("ldf1e16", histname, 100, 0, 1000);
  hldf[1] = new TH1F("ldf3e16", histname, 100, 0, 1000);
  hldf[2] = new TH1F("ldf1e17", histname, 100, 0, 1000);
  hldf[3] = new TH1F("ldf3e17", histname, 100, 0, 1000);

  Int_t i;

  for( i=0; i<4; i++ ) {
    hldf[i]->Sumw2(); hldf[i]->SetLineWidth(2); 
  }
  
  hldf[0]->SetLineColor(kRed);
  hldf[1]->SetLineColor(kBlue);
  hldf[2]->SetLineColor(kGreen);
  hldf[3]->SetLineColor(kBlack);

  Int_t k;
  Int_t j;
  Double_t e_weight[10], m_weight[10], l_weight[10];
  Double_t e_var, m_var, g_var, l_var;
  Double_t e_avg, m_avg, g_avg, l_avg;
  
  for( k=0; k<4; k++ ){
    char energy[10]; char energy2[10];
    if (k==0) { sprintf(energy, "1.00E16"); sprintf(energy2, "1.00e16"); 
    }else if (k==1) { sprintf(energy, "3.16E16"); sprintf(energy2, "3.16e16");
    }else if (k==2) { sprintf(energy, "1.00E17"); sprintf(energy2, "1.00e17");
    }else if (k==3) { sprintf(energy, "3.16E17"); sprintf(energy2, "3.16e17");
    }

    for( i=0; i<numhists; i++ ) {
      sprintf(hnum, "%.3d", i+1 );
      
      // eg: "3.16E17/0.85/007/proton_3.16e16-ev_0.555-rad.ev007_3hist.root"
      sprintf(filename, "%s/%s/%s/proton_%s-ev_%s-rad.ev%s_3hist.root", 
	      energy, angle, hnum, energy2, radians, hnum);
      
      TFile* fin = new TFile(filename);
      
      he[i] = (TH1F*)gDirectory->Get(efine);
      hm[i] = (TH1F*)gDirectory->Get(mfine);
      //    hg[i] = (TH1F*)gDirectory->Get(gfine);

    }

  for (i=1; i<=100; i++) {

    e_var = 0; e_avg = 0;
    m_var = 0; m_avg = 0;
    g_var = 0; g_avg = 0;
    l_var = 0; l_avg = 0;

    for (j=0; j<numhists; j++ ) {

      //      e_weight[j]=(he[j]->GetBinContent(i))/numhists; e_avg += e_weight[j];
      //      electrons->AddBinContent(i, e_weight[j]);
      
      //      m_weight[j]=(hm[j]->GetBinContent(i))/numhists; m_avg += m_weight[j];
      //     muons->AddBinContent(i, m_weight[j]);
      
      //      g_weight[j]=(hg[j]->GetBinContent(i))/numhists; 
      //      g_avg += g_weight[j];
      //      gammas->AddBinContent(i, g_weight[j]);
      
      l_weight[j]=(he[j]->GetBinContent(i)+hm[j]->GetBinContent(i))/numhists;
      l_avg += l_weight[j];
      hldf[k]->AddBinContent(i, l_weight[j]);

    }

    for ( j=0; j<numhists; j++ ) {
      //      e_var += (numhists*e_weight[j]-e_avg)*(numhists*e_weight[j]-e_avg);
      //      m_var += (numhists*m_weight[j]-m_avg)*(numhists*m_weight[j]-m_avg);
      //      g_var += (numhists*g_weight[j]-g_avg)*(numhists*g_weight[j]-g_avg);
      l_var += (numhists*l_weight[j]-l_avg)*(numhists*l_weight[j]-l_avg);
    }
    
    //    electrons->SetBinError(i, TMath::Sqrt(e_var/(numhists-1)));
    //    muons->SetBinError(i, TMath::Sqrt(m_var/(numhists-1)));
    //    gammas->SetBinError(i, TMath::Sqrt(g_var/(numhists-1)));
    hldf[k]->SetBinError(i, TMath::Sqrt(l_var/(numhists-1)));
    
  }    
  
  TCanvas* c1 = new TCanvas("c1", "Average LDF", 1);
  
  c1->SetLogy();
  
  hldf[0]->Draw("e1");
  hldf[1]->Draw("e1same");
  hldf[2]->Draw("e1same");
  hldf[3]->Draw("e1same");
  c1->Update();
  
  }