// this is elina's program to make ldf plots at chiquita energies
// it is intended to be run from high_energy/Proton/
// it plots the average density for each type of particle

// change the energy and setangle by hand

{
  char energy[10]; sprintf(energy, "3.16E17"); // 3.16E16, 1.00E17, etc. 
  char energy2[10]; sprintf(energy2, "3.16e17");
  Int_t setangle = 3; // 1 for 0.95, 2 for 0.85, 3 for 0.75
  char angle[10];
  char radians[10];

  if ( setangle == 1) {
    sprintf(angle, "0.95"); sprintf(radians, "0.318");
  }else if( setangle == 2) { 
    sprintf(angle, "0.85"); sprintf(radians, "0.555");
  }else if ( setangle == 3 ) {
    sprintf(angle, "0.75"); sprintf(radians, "0.723");
  }

  Int_t numhists = 10;
  char efine[10]; sprintf(efine, "efine_3");
  char mfine[10]; sprintf(mfine, "mfine_3");
  char gfine[10]; sprintf(gfine, "gfine_3");
  char hnum[10]; // run number
  char filename[100]; 
  char histname[100];

  TH1F* he[numhists]; // input arrays of histograms
  TH1F* hm[numhists];
  TH1F* hg[numhists];


  gStyle->SetOptStat(0);

  //  TH1F* hconst = new TH1F("hconst", "hconst", 100, 0, 1000);
  //  hconst->Sumw2();

  sprintf(histname, "%s_%s", energy, angle);
  TH1F* electrons = new TH1F("electrons", histname, 100, 0, 1000);
  TH1F* muons = new TH1F("muons", histname, 100, 0, 1000);
  TH1F* gammas = new TH1F("gammas", histname, 100, 0, 1000);
  TH1F* ldf = new TH1F("ldf", histname, 100, 0, 1000);
  electrons->Sumw2(); muons->Sumw2(); gammas->Sumw2(); ldf->Sumw2();

  electrons->SetLineWidth(2);
  muons->SetLineWidth(2);
  gammas->SetLineWidth(2);
  ldf->SetLineWidth(2);

  electrons->SetLineColor(kBlue);
  muons->SetLineColor(kRed);
  gammas->SetLineColor(kBlack);
  ldf->SetLineColor(kGreen);

  Int_t i;

  for( i=0; i<numhists; i++ ) {
    sprintf(hnum, "%.3d", i+1 );
    
    // e.g. "3.16E16/0.85/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 j;
  Double_t e_weight[10], m_weight[10], g_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=1; i<=100; i++) {

    e_avg = 0; e_var = 0;
    m_avg = 0; m_var = 0;
    g_avg = 0; g_var = 0;
    l_avg = 0; l_var = 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];
      ldf->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)));
    ldf->SetBinError(i, TMath::Sqrt(l_var/(numhists-1)));

  }    

  TCanvas* c1 = new TCanvas("c1", "Particle Densities", 1);
  c1->SetLogy();

  electrons->Draw();
  for ( i=0; i<10; i++ ) {
    he[i]->SetLineColor(7);
    he[i]->SetLineWidth(1);
    he[i]->Draw("same");
  }
  electrons->Draw("e1same");

  for ( i=0; i<10; i++ ) {
    hm[i]->SetLineColor(9);
    hm[i]->SetLineWidth(1);
    hm[i]->Draw("same");
  }
  muons->Draw("e1same");
 
  for ( i=0; i<10; i++ ) {
    hg[i]->SetLineColor(45);
    hg[i]->SetLineWidth(1);
    hg[i]->Draw("same");
  }
  gammas->Draw("e1same");

  c1->Update();

  //  TLegend* legend1 = new TLegend(0.65,0.75,0.95,0.95);
  //  legend1->AddEntry(gammas, "#gamma","1"); 
  //  legend1->AddEntry(electrons, "e^{+} e^{-}","1");
  //  legend1->AddEntry(muons, "#mu^{+} #mu^{-}","1");
  //  legend1->Draw();

  //  c1->Update();

  TCanvas* c2 = new TCanvas("c2", "Average LDF", 1);
  c2->SetLogy();

  electrons->Draw("e1");
  muons->Draw("e1same");
  ldf->Draw("e1same");

  c2->Update();

  //  TLegend* legend2 = new TLegend(0.65,0.75,0.95,0.95);
  //  legend2->AddEntry(electrons, "e^{+} e^{-}","1");
  //  legend2->AddEntry(muons, "#mu^{+} #mu^{-}","1");
  //  legend2->AddEntry(ldf, "e + #mu","1");
  //  legend2->Draw();

  //  c2->Update();
}