// 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 angle for a given energy

{
  char energy[10]; sprintf(energy, "1.00E16");
  char energy2[10]; sprintf(energy2, "1.00e16");
  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 %s", energy);
  hldf[0] = new TH1F("ldf95", histname, 100, 0, 1000);
  hldf[1] = new TH1F("ldf85", histname, 100, 0, 1000);
  hldf[2] = new TH1F("ldf75", histname, 100, 0, 1000);

  int i;

  for( i=0; i<3; i++ ) {
    hldf[i]->Sumw2(); hldf[i]->SetLineWidth(1); 
  }
  
  hldf[0]->SetLineColor(kRed);
  hldf[1]->SetLineColor(kBlue);
  hldf[2]->SetLineColor(kGreen);
  
  Int_t k;
  char angle[10]; char radians[10];

  for( k=0; k<3; k++ ){
    if (k==0) { sprintf(angle, "0.95"); sprintf(radians, "0.318"); 
    }else if (k==1) { sprintf(angle, "0.85"); sprintf(radians, "0.555");
    }else if (k==2) { sprintf(angle, "0.75"); sprintf(radians, "0.723");
    }
    
    for( i=0; i<10; 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);
      
    }
    
    Int_t i, 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( 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");
  c1->Update();
  
}