#include <iostream>

/************************************************************
 *
 * This macro takes in the results from plot_waveforms2
 * and plots the histograms.
 *
 ************************************************************/

int plot_results(std::string infilename = "results/all.root"){

  gStyle->SetOptStat(0);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);

  const Int_t colors[] = {kBlack, kRed, kYellow+2, kGreen+2, kBlue, kMagenta, kCyan, 32};

  const Int_t NMAX = 40;
  const Int_t NSECTORS = 32;

  // Get file name only with no directory
  char command[200];
  // cout << infilename << endl;
  sprintf(command,"basename %s | sed 's:results.::' | sed 's:\\.root::' > ___tmp.txt",infilename.c_str());
  // cout << command << endl;
  system(command);
  ifstream IN___tmp("___tmp.txt");
  std::string basename;
  IN___tmp >> basename;
  system("rm -f ___tmp.txt");

  TFile *infile = new TFile(infilename.c_str());

  ////////////////////////////////////////////////////////////
  //                                                        //
  //     Get histograms defined in plot_waveform            //
  //                                                        //
  ////////////////////////////////////////////////////////////
  TH1F *hnvchannels_for_average = (TH1F*)infile->Get("hnvchannels_for_average")->Clone("_hnvchannels_for_average");
  TH2F *hnvchannels_for_average_nwaveforms = (TH2F*)infile->Get("hnvchannels_for_average_nwaveforms")->Clone("_hnvchannels_for_average_nwaveforms");

  TH2F *hthreshold_nvchannels_for_average = (TH2F*)infile->Get("hthreshold_nvchannels_for_average")->Clone("_hthreshold_nvchannels_for_average");

  TH1F *heventnum_nvchannels_for_average_is_32 = (TH1F*)infile->Get("heventnum_nvchannels_for_average_is_32")->Clone("_heventnum_nvchannels_for_average_is_32");
  TH1F *heventnum_nvchannels_for_average_less_than_24 = (TH1F*)infile->Get("heventnum_nvchannels_for_average_less_than_24")->Clone("_heventnum_nvchannels_for_average_less_than_24");
  TH1F *heventnum_nvchannels_for_average_is_0 = (TH1F*)infile->Get("heventnum_nvchannels_for_average_is_0")->Clone("_heventnum_nvchannels_for_average_is_0");
  TH1F *heventnum_nSignals_greater_than_5 = (TH1F*)infile->Get("heventnum_nSignals_greater_than_5")->Clone("_heventnum_nSignals_greater_than_5");
  TH1F *heventnum_nwaveforms_not_32 = (TH1F*)infile->Get("heventnum_nwaveforms_not_32")->Clone("_heventnum_nwaveforms_not_32");

  TH1F *hpedestal[NSECTORS];
  TH1F *hmaxcounts[NSECTORS];
  TH1F *hdiffmaxcounts_pedestal[NSECTORS];
  TH2F *hmaxcounts_pedestal[NSECTORS];

  // TH1F *hpedestal_normalized[NSECTORS];
  TH1F *hwaveform_normalized[NSECTORS];
  TH1F *hdiffmaxcounts_pedestal_normalized[NSECTORS];
  TH1F *hnormalized_integral[NSECTORS];

  TF1 *fdiffmaxcounts_pedestal_normalized[NSECTORS];

  TH1F *hdiffmaxcounts_pedestal_notail[NSECTORS];

  TH1F *hhittime[NSECTORS];
  TH2F *hsector_hittime = (TH2F*)infile->Get("hsector_hittime")->Clone("_hsector_hittime");

  TH1F *hnwaveforms = (TH1F*)infile->Get("hnwaveforms")->Clone("_hnwaveforms");
  TH1F *hnSignals = (TH1F*)infile->Get("hnSignals")->Clone("_hnSignals");
  TH2F *hmultiplicity = (TH2F*)infile->Get("hmultiplicity")->Clone("_hmultiplicity");

  // Histograms from Mike's analysis
  TH1F *hVarX = (TH1F*)infile->Get("hVarX")->Clone("_hVarX");
  TH2F *hSecVsVarX = (TH2F*)infile->Get("hSecVsVarX")->Clone("_hSecVsVarX");
  TH1D *hPrimeTime = (TH1D*)infile->Get("hPrimeTime")->Clone("_hPrimeTime");
  TH1D *hGammaRise = (TH1D*)infile->Get("hGammaRise")->Clone("_hGammaRise");
  TH1D *hGammaFall = (TH1D*)infile->Get("hGammaFall")->Clone("_hGammaFall");
  TH1D *hcVal = (TH1D*)infile->Get("hcVal")->Clone("_hcVal");

  // Correlation between
  // - diffmaxcounts_pedestal
  // - normalized integral
  // - fitted amplitude
  TH2F *hdiffmaxcounts_pedestal_fitAmp[NSECTORS];
  TH2F *hdiffmaxcounts_pedestal_normalized_integral[NSECTORS];
  TH2F *hnormalized_integral_fitAmp[NSECTORS];

  char hname[400];
  char hnewname[400];

  for(Int_t i=0;i<NSECTORS;i++){
    sprintf(hname,"hpedestal%2.2d",i+1);
    sprintf(hnewname,"_hpedestal%2.2d",i+1);
    hpedestal[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hmaxcounts%2.2d",i+1);
    sprintf(hnewname,"_hmaxcounts%2.2d",i+1);
    hmaxcounts[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hdiffmaxcounts_pedestal%2.2d",i+1);
    sprintf(hnewname,"_hdiffmaxcounts_pedestal%2.2d",i+1);
    hdiffmaxcounts_pedestal[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hmaxcounts_pedestal%2.2d",i+1);
    sprintf(hnewname,"_hmaxcounts_pedestal%2.2d",i+1);
    hmaxcounts_pedestal[i] = (TH2F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hwaveform_normalized%2.2d",i+1);
    sprintf(hnewname,"_hwaveform_normalized%2.2d",i+1);
    hwaveform_normalized[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hdiffmaxcounts_pedestal_normalized%2.2d",i+1);
    sprintf(hnewname,"_hdiffmaxcounts_pedestal_normalized%2.2d",i+1);
    hdiffmaxcounts_pedestal_normalized[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hnormalized_integral%2.2d",i+1);
    sprintf(hnewname,"_hnormalized_integral%2.2d",i+1);
    hnormalized_integral[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"fdiffmaxcounts_pedestal_normalized%2.2d",i+1);
    fdiffmaxcounts_pedestal_normalized[i] = new TF1(hname,"[0]*exp(-pow(x-[1],2.)/2./pow([2],2.))",1500,2000);

    sprintf(hname,"hdiffmaxcounts_pedestal_notail%2.2d",i+1);
    sprintf(hnewname,"_hdiffmaxcounts_pedestal_notail%2.2d",i+1);
    hdiffmaxcounts_pedestal_notail[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hhittime%2.2d",i+1);
    sprintf(hnewname,"_hhittime%2.2d",i+1);
    hhittime[i] = (TH1F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hdiffmaxcounts_pedestal_fitAmp%2.2d",i+1);
    sprintf(hnewname,"_hdiffmaxcounts_pedestal_fitAmp%2.2d",i+1);
    hdiffmaxcounts_pedestal_fitAmp[i] = (TH2F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hdiffmaxcounts_pedestal_normalized_integral%2.2d",i+1);
    sprintf(hnewname,"_hdiffmaxcounts_pedestal_normalized_integral%2.2d",i+1);
    hdiffmaxcounts_pedestal_normalized_integral[i] = (TH2F*)infile->Get(hname)->Clone(hnewname);

    sprintf(hname,"hnormalized_integral_fitAmp%2.2d",i+1);
    sprintf(hnewname,"_hnormalized_integral_fitAmp%2.2d",i+1);
    hnormalized_integral_fitAmp[i] = (TH2F*)infile->Get(hname)->Clone(hnewname);
  }

  // Projections
  TH1D *hhittimeEachSector[NSECTORS];
  Double_t MAX_HITTIMEEACHSECTOR = -999;
  for(Int_t i=0;i<NSECTORS;i++){
    sprintf(hname,"hhittimeEachSector%2.2d",i+1);
    hhittimeEachSector[i] = (TH1D*)hsector_hittime->ProjectionX(hname,i+1,i+1);
    if(hhittimeEachSector[i]->GetMaximum() > MAX_HITTIMEEACHSECTOR) MAX_HITTIMEEACHSECTOR = hhittimeEachSector[i]->GetMaximum() * 1.20;
  }
  //______________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////
  //                                                        //
  //     Get maximum counts for each type of histogram      //
  //     so we can draw all sectors at the same scale.      //
  //                                                        //
  ////////////////////////////////////////////////////////////
  Double_t MAX_PEDESTAL = -999;
  Double_t MAX_MAXCOUNTS = -999;
  Double_t MAX_DIFFMAXCOUNTS_PEDESTAL = -999;
  Double_t MAX_WAVEFORM_NORMALIZED = -999;
  Double_t MAX_DIFFMAXCOUNTS_PEDESTAL_NORMALIZED = -999;
  Double_t MAX_DIFFMAXCOUNTS_PEDESTAL_NOTAIL = -999;
  Double_t MAX_NORMALIZED_INTEGRAL = -999;
  Double_t MAX_HITTIME = -999;
  for(Int_t i=0;i<NSECTORS;i++){
    if(hpedestal[i]->GetBinContent(hpedestal[i]->GetMaximumBin()) > MAX_PEDESTAL)
      MAX_PEDESTAL = hpedestal[i]->GetBinContent(hpedestal[i]->GetMaximumBin()) * 1.35;

    if(hmaxcounts[i]->GetBinContent(hmaxcounts[i]->GetMaximumBin()) > MAX_MAXCOUNTS)
      MAX_MAXCOUNTS = hmaxcounts[i]->GetBinContent(hmaxcounts[i]->GetMaximumBin()) * 1.35;

    if(hdiffmaxcounts_pedestal[i]->GetBinContent(hdiffmaxcounts_pedestal[i]->GetMaximumBin()) > MAX_DIFFMAXCOUNTS_PEDESTAL)
      MAX_DIFFMAXCOUNTS_PEDESTAL = hdiffmaxcounts_pedestal[i]->GetBinContent(hdiffmaxcounts_pedestal[i]->GetMaximumBin()) * 1.35;

    if(hwaveform_normalized[i]->GetBinContent(hwaveform_normalized[i]->GetMaximumBin()) > MAX_WAVEFORM_NORMALIZED)
      MAX_WAVEFORM_NORMALIZED = hwaveform_normalized[i]->GetBinContent(hwaveform_normalized[i]->GetMaximumBin()) * 1.35;

    if(hdiffmaxcounts_pedestal_normalized[i]->GetBinContent(hdiffmaxcounts_pedestal_normalized[i]->GetMaximumBin()) > MAX_DIFFMAXCOUNTS_PEDESTAL_NORMALIZED
       )
      MAX_DIFFMAXCOUNTS_PEDESTAL_NORMALIZED = hdiffmaxcounts_pedestal_normalized[i]->GetBinContent(hdiffmaxcounts_pedestal_normalized[i]->GetMaximumBin()) * 1.35;

    if(hdiffmaxcounts_pedestal_notail[i]->GetBinContent(hdiffmaxcounts_pedestal_notail[i]->GetMaximumBin()) > MAX_DIFFMAXCOUNTS_PEDESTAL_NOTAIL
       )
      MAX_DIFFMAXCOUNTS_PEDESTAL_NOTAIL = hdiffmaxcounts_pedestal_notail[i]->GetBinContent(hdiffmaxcounts_pedestal_notail[i]->GetMaximumBin()) * 1.35;

    if(hnormalized_integral[i]->GetBinContent(hnormalized_integral[i]->GetMaximumBin()) > MAX_NORMALIZED_INTEGRAL)
      MAX_NORMALIZED_INTEGRAL = hnormalized_integral[i]->GetBinContent(hnormalized_integral[i]->GetMaximumBin()) * 1.35;

    if(hhittime[i]->GetBinContent(hhittime[i]->GetMaximumBin()) > MAX_HITTIME)
      MAX_HITTIME = hhittime[i]->GetBinContent(hhittime[i]->GetMaximumBin()) * 1.35;
  }
  //___________________________________________________________________________________________________________________________

  TCanvas *canvas = new TCanvas("canvas","canvas",1200,600);
  canvas->Divide(2,2,.002,.002);
  canvas->SetRightMargin(0.10);
  canvas->SetTopMargin(0.02);
  canvas->SetLeftMargin(0.10);
  canvas->SetBottomMargin(0.10);
  canvas->Draw();

  TCanvas *canvasAll = new TCanvas("canvasAll","canvasAll",1200,600);
  canvasAll->SetRightMargin(0.10);
  canvasAll->SetTopMargin(0.02);
  canvasAll->SetLeftMargin(0.10);
  canvasAll->SetBottomMargin(0.10);
  canvasAll->Draw();

  TLatex *latex = new TLatex();
  latex->SetTextColor(kBlack);
  latex->SetTextAlign(12);
  latex->SetTextFont(132);
  latex->SetTextSize(0.060);
  latex->SetNDC(1);

  TLegend *legend = new TLegend(0.15,0.80,0.95,0.93);
  legend->SetTextColor(kBlack);
  legend->SetTextAlign(12);
  legend->SetTextFont(132);
  legend->SetTextSize(0.040);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);
  legend->SetNColumns(6);

  TLegend *legend_sub[4];
  for(Int_t i=0;i<4;i++){
    legend_sub[i] = new TLegend(0.15,0.70,0.95,0.89);
    legend_sub[i]->SetTextColor(kBlack);
    legend_sub[i]->SetTextAlign(12);
    legend_sub[i]->SetTextFont(132);
    legend_sub[i]->SetTextSize(0.050);
    legend_sub[i]->SetBorderSize(0);
    legend_sub[i]->SetFillStyle(0);
    legend_sub[i]->SetNColumns(2);
  }

  char text[200];

  // Counter for plot number
  Int_t nplots = 0;

  // Draw number of waveforms
  canvasAll->cd();
  canvasAll->SetLogy(1);
  hnwaveforms->Draw();
  hnwaveforms->Draw("textsame");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___nwaveforms.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogy(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  canvasAll->SetLogy(1);
  hnvchannels_for_average->Draw();
  hnvchannels_for_average->Draw("textsame");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___nvchannels_for_average.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogy(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  canvasAll->SetLogz(1);
  hnvchannels_for_average_nwaveforms->Draw("colz");
  hnvchannels_for_average_nwaveforms->Draw("textsame");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___nvchannels_for_average_nwaveforms.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogz(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  canvasAll->SetLogz(1);
  hthreshold_nvchannels_for_average->Draw("colz");
  // hthreshold_nvchannels_for_average->Draw("textsame");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___threshold_nvchannels_for_average.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogz(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  legend->Clear();

  heventnum_nvchannels_for_average_is_32->SetLineColor(kRed);
  heventnum_nvchannels_for_average_is_32->Draw();
  legend->AddEntry(heventnum_nvchannels_for_average_is_32,"#channels for pedestal average=32","L");

  heventnum_nvchannels_for_average_less_than_24->SetLineColor(kBlack);
  heventnum_nvchannels_for_average_less_than_24->Draw("same");
  legend->AddEntry(heventnum_nvchannels_for_average_less_than_24,"#channels for pedestal average<24","L");

  heventnum_nvchannels_for_average_is_0->SetLineColor(kGreen+2);
  heventnum_nvchannels_for_average_is_0->Draw("same");
  legend->AddEntry(heventnum_nvchannels_for_average_is_0,"#channels for pedestal average=0","L");

  heventnum_nSignals_greater_than_5->SetLineColor(kBlue);
  heventnum_nSignals_greater_than_5->Draw("same");
  legend->AddEntry(heventnum_nSignals_greater_than_5,"#signals>5","L");

  heventnum_nwaveforms_not_32->SetLineColor(kMagenta);
  heventnum_nwaveforms_not_32->Draw("same");
  legend->AddEntry(heventnum_nwaveforms_not_32,"#waveforms#neq32","L");

  legend->SetNColumns(2);
  legend->Draw("same");

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___eventnums.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogz(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // Which sector had pedestal above readout threshold of 850
  canvasAll->cd();
  hsector_with_pedestal_above_threshold->Draw();
  hsector_with_pedestal_above_threshold_nvchannels_for_average_is_32->SetLineColor(kRed);
  hsector_with_pedestal_above_threshold_nvchannels_for_average_is_32->Draw("same");
  sprintf(text,"figures/%s/summary/%3.3d___sector_above_threshold.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogz(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // Draw nSignals
  canvasAll->cd();
  canvasAll->SetLogy(1);
  hnSignals->SetMinimum(0.5);
  hnSignals->Draw();
  hnSignals->Draw("textsame");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___nSignals.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->SetLogy(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // Draw nSignals
  canvasAll->cd();
  hnvchannels_for_average_nSignals->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___nvchannels_for_average_nSignals.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // Draw sector correlations when nSignals = 2
  canvasAll->cd();
  canvasAll->SetLogz(1);
  hmultiplicity->Draw("colz");
  hmultiplicity->Draw("textsame");

  TLine *line_multiplicity = new TLine(0.5,0.5,32.5,32.5);
  line_multiplicity->SetLineColor(28);
  line_multiplicity->SetLineStyle(2);
  line_multiplicity->SetLineWidth(1);
  line_multiplicity->Draw("same");

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___multiplicity.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  canvasAll->SetLogz(0);
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //                   Draw pedestal distributions                  //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  TF1 *fGaussPedestal[NSECTORS];
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);
    canvas->cd(index+1)->SetLogy(1);

    hpedestal[i]->SetMinimum(0.5);
    hpedestal[i]->SetMaximum(MAX_PEDESTAL);
    // hpedestal[i]->SetAxisRange(900,1200);

    hpedestal[i]->GetXaxis()->SetNdivisions(505);
    hpedestal[i]->GetXaxis()->SetLabelSize(0.050);
    hpedestal[i]->GetYaxis()->SetNdivisions(505);
    hpedestal[i]->GetYaxis()->SetLabelSize(0.050);
    hpedestal[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hpedestal[i]->Draw();
    else
      hpedestal[i]->Draw("same");

    // Fit with Gaussian to get approximate width
    sprintf(text,"fGaussPedestal%2.2d",i+1);
    fGaussPedestal[i] = new TF1(text,"[0] * exp(-pow(x-[1],2.)/2./pow([2],2.))",3200,3600);
    fGaussPedestal[i]->SetParameter(0,hpedestal[i]->GetMaximum());
    fGaussPedestal[i]->SetParameter(1,hpedestal[i]->GetMean());
    fGaussPedestal[i]->SetParameter(2,hpedestal[i]->GetRMS());
    // hpedestal[i]->Fit(fGaussPedestal[i],"EQRN");

    fGaussPedestal[i]->SetLineColor(colors[i % 8]);
    fGaussPedestal[i]->SetLineStyle(2);
    fGaussPedestal[i]->SetLineWidth(1);
    fGaussPedestal[i]->Draw("same");

    cout << 1150. - fGaussPedestal[i]->GetParameter(1) << endl;

    sprintf(text,"sec.%2.2d %6.1f, %4.2f",i+1,fGaussPedestal[i]->GetParameter(1),fGaussPedestal[i]->GetParameter(2)); // fGaussPedestal[i]->GetParError(2)
    legend_sub[index]->AddEntry(hpedestal[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  // canvas->cd(33)->Clear();
  // canvas->cd(33);
  // sprintf(text,"pedestals");
  // latex->DrawLatex(0.20,0.65,text);

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___pedestals.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //                   Draw maxcounts distributions                 //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);

    hmaxcounts[i]->SetMinimum(0.5);
    hmaxcounts[i]->SetMaximum(MAX_MAXCOUNTS);
    // hmaxcounts[i]->SetAxisRange(900,1400);

    hmaxcounts[i]->GetXaxis()->SetNdivisions(505);
    hmaxcounts[i]->GetXaxis()->SetLabelSize(0.050);
    hmaxcounts[i]->GetYaxis()->SetNdivisions(505);
    hmaxcounts[i]->GetYaxis()->SetLabelSize(0.050);
    hmaxcounts[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hmaxcounts[i]->Draw();
    else
      hmaxcounts[i]->Draw("same");

    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hmaxcounts[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___maxcounts.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //              Draw maxcounts - pedestal distributions           //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);
    canvas->cd(index+1)->SetLogy(1);

    hdiffmaxcounts_pedestal[i]->SetMinimum(0.5);
    hdiffmaxcounts_pedestal[i]->SetMaximum(MAX_DIFFMAXCOUNTS_PEDESTAL);
    // hdiffmaxcounts_pedestal[i]->SetAxisRange(-20,300);

    hdiffmaxcounts_pedestal[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hdiffmaxcounts_pedestal[i]->Draw();
    else
      hdiffmaxcounts_pedestal[i]->Draw("same");
    
    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hdiffmaxcounts_pedestal[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;

  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->SetLogy(0);
  }
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //              Draw maxcounts - pedestal distributions           //
  //              (zoomed in, linear scale)                         //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);

    hdiffmaxcounts_pedestal[i]->SetMinimum(0.5);
    hdiffmaxcounts_pedestal[i]->SetMaximum(MAX_DIFFMAXCOUNTS_PEDESTAL * 0.002);
    hdiffmaxcounts_pedestal[i]->SetAxisRange(-20,300);

    hdiffmaxcounts_pedestal[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hdiffmaxcounts_pedestal[i]->Draw();
    else
      hdiffmaxcounts_pedestal[i]->Draw("same");
    
    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hdiffmaxcounts_pedestal[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal_linear.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //                 Draw maxcounts vs pedestal                     //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  canvas->Clear();
  canvas->Divide(6,6,.002,.002);
  for(Int_t i=0;i<NSECTORS;i++){
    canvas->cd(i+1);
    canvas->cd(i+1)->Clear();
    canvas->cd(i+1)->SetLogz(1);

    hmaxcounts_pedestal[i]->GetXaxis()->SetNdivisions(505);
    hmaxcounts_pedestal[i]->GetXaxis()->SetLabelSize(0.050);
    hmaxcounts_pedestal[i]->GetYaxis()->SetNdivisions(505);
    hmaxcounts_pedestal[i]->GetYaxis()->SetLabelOffset(0.001);
    hmaxcounts_pedestal[i]->GetYaxis()->SetLabelSize(0.050);
    // hmaxcounts_pedestal[i]->SetAxisRange(900,1300,"X");
    // hmaxcounts_pedestal[i]->SetAxisRange(900,1400,"Y");
    hmaxcounts_pedestal[i]->Draw("colz");
    
  } // end of loop over sectors

  canvas->cd(33)->Clear();
  canvas->cd(33);
  latex->SetTextSize(0.15);
  sprintf(text,"maxcounts vs");
  latex->DrawLatex(0.20,0.65,text);
  sprintf(text,"pedestals");
  latex->DrawLatex(0.20,0.45,text);
  latex->SetTextSize(0.060);

  sprintf(text,"figures/%s/summary/%3.3d___maxcounts_pedestals.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //            Draw normalized waveform distributions              //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->Clear();
  canvas->Divide(2,2,.002,.002);

  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  TF1 *fGaussWaveform_normalized[NSECTORS];
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);

    hwaveform_normalized[i]->SetMinimum(0);
    hwaveform_normalized[i]->SetMaximum(MAX_WAVEFORM_NORMALIZED);
    hwaveform_normalized[i]->GetXaxis()->SetNdivisions(505);
    hwaveform_normalized[i]->GetXaxis()->SetLabelSize(0.050);
    hwaveform_normalized[i]->GetYaxis()->SetNdivisions(505);
    hwaveform_normalized[i]->GetYaxis()->SetLabelSize(0.050);
    hwaveform_normalized[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hwaveform_normalized[i]->Draw();
    else
      hwaveform_normalized[i]->Draw("same");

    // Fit with Gaussian to get approximate width
    sprintf(text,"fGaussWaveform_normalized%2.2d",i);
    fGaussWaveform_normalized[i] = new TF1(text,"[0] * exp(-pow(x-[1],2.)/2./pow([2],2.))",-20,20);
    fGaussWaveform_normalized[i]->SetParameter(0,hwaveform_normalized[i]->GetMaximum());
    fGaussWaveform_normalized[i]->SetParameter(1,hwaveform_normalized[i]->GetMean());
    fGaussWaveform_normalized[i]->SetParameter(2,hwaveform_normalized[i]->GetRMS());
    hwaveform_normalized[i]->Fit(fGaussWaveform_normalized[i],"EQRN");

    fGaussWaveform_normalized[i]->SetLineColor(colors[i % 8]);
    fGaussWaveform_normalized[i]->SetLineStyle(2);
    fGaussWaveform_normalized[i]->SetLineWidth(1);
    fGaussWaveform_normalized[i]->Draw("same");

    sprintf(text,"sec.%2.2d %4.2f #pm %4.2f",i+1,fGaussWaveform_normalized[i]->GetParameter(2),fGaussWaveform_normalized[i]->GetParError(2));
    legend_sub[index]->AddEntry(hwaveform_normalized[i],text,"L");
    
    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"figures/%s/summary/%3.3d___waveform_normalizeds.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //              Draw maxcounts - pedestal distributions           //
  //              (normalized)                                      //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);
    canvas->cd(index+1)->SetLogy(1);

    hdiffmaxcounts_pedestal_normalized[i]->SetMinimum(0.5);
    hdiffmaxcounts_pedestal_normalized[i]->SetMaximum(MAX_DIFFMAXCOUNTS_PEDESTAL_NORMALIZED);
    hdiffmaxcounts_pedestal_normalized[i]->SetAxisRange(-20,2480);

    hdiffmaxcounts_pedestal_normalized[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_normalized[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_normalized[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_normalized[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_normalized[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0){
      hdiffmaxcounts_pedestal_normalized[i]->Draw("axis");
      fdiffmaxcounts_pedestal_normalized[i]->SetParameter(0,20);
      fdiffmaxcounts_pedestal_normalized[i]->SetParameter(1,1700);
      fdiffmaxcounts_pedestal_normalized[i]->SetParameter(2,100);
      hdiffmaxcounts_pedestal_normalized[i]->Fit(fdiffmaxcounts_pedestal_normalized[i],"EQRN");
      hdiffmaxcounts_pedestal_normalized[i]->Draw();
    }else{
      hdiffmaxcounts_pedestal_normalized[i]->Draw("same");
      fdiffmaxcounts_pedestal_normalized[i]->SetParameter(0,20);
      fdiffmaxcounts_pedestal_normalized[i]->SetParameter(1,1700);
      fdiffmaxcounts_pedestal_normalized[i]->SetParameter(2,100);
      hdiffmaxcounts_pedestal_normalized[i]->Fit(fdiffmaxcounts_pedestal_normalized[i],"EQRN");
    }

    sprintf(text,"sec.%2.2d: %5.3f #pm %5.3f",i+1,fdiffmaxcounts_pedestal_normalized[i]->GetParameter(1),fdiffmaxcounts_pedestal_normalized[i]->GetParError(1));
    legend_sub[index]->AddEntry(hdiffmaxcounts_pedestal_normalized[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal_normalized_normalized.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;

  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->SetLogy(0);
  }
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //              Draw maxcounts - pedestal distributions           //
  //              for normalized waveforms                          //
  //              (zoomed in, linear scale)                         //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);

    hdiffmaxcounts_pedestal_normalized[i]->SetMinimum(0.5);
    hdiffmaxcounts_pedestal_normalized[i]->SetMaximum(MAX_DIFFMAXCOUNTS_PEDESTAL_NORMALIZED * 0.008);
    hdiffmaxcounts_pedestal_normalized[i]->SetAxisRange(-20,300);

    hdiffmaxcounts_pedestal_normalized[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_normalized[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_normalized[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_normalized[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_normalized[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hdiffmaxcounts_pedestal_normalized[i]->Draw();
    else
      hdiffmaxcounts_pedestal_normalized[i]->Draw("same");
    
    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hdiffmaxcounts_pedestal_normalized[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal_normalized_linear.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //                   Draw normalized integrals                    //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);
    canvas->cd(index+1)->SetLogy(1);

    hnormalized_integral[i]->SetMinimum(0.5);
    hnormalized_integral[i]->SetMaximum(MAX_NORMALIZED_INTEGRAL);
    // hnormalized_integral[i]->SetAxisRange(0,400);

    hnormalized_integral[i]->GetXaxis()->SetNdivisions(505);
    hnormalized_integral[i]->GetXaxis()->SetLabelSize(0.050);
    hnormalized_integral[i]->GetYaxis()->SetNdivisions(505);
    hnormalized_integral[i]->GetYaxis()->SetLabelSize(0.050);
    hnormalized_integral[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hnormalized_integral[i]->Draw();
    else
      hnormalized_integral[i]->Draw("same");

    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hnormalized_integral[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___normalized_integral.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;

  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->SetLogy(0);
  }
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //              Draw maxcounts - pedestal distributions           //
  //              (no tail events)                                  //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
    legend_sub[i]->SetHeader("no tail events");
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);
    canvas->cd(index+1)->SetLogy(1);

    hdiffmaxcounts_pedestal_notail[i]->SetMinimum(0.5);
    hdiffmaxcounts_pedestal_notail[i]->SetMaximum(MAX_DIFFMAXCOUNTS_PEDESTAL_NOTAIL);
    // hdiffmaxcounts_pedestal_notail[i]->SetAxisRange(-20,500);

    hdiffmaxcounts_pedestal_notail[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_notail[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_notail[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_notail[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_notail[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hdiffmaxcounts_pedestal_notail[i]->Draw();
    else
      hdiffmaxcounts_pedestal_notail[i]->Draw("same");

    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hdiffmaxcounts_pedestal_notail[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal_notail.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;

  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->SetLogy(0);
  }
  //___________________________________________________________________________________________________________________________

  // Create the sum of all
  // - hdiffmaxcounts_pedestal
  // - hdiffmaxcounts_pedestal_normalized
  // - hnormalized_integral
  // - hittime
  TH1F *hdiffmaxcounts_pedestalAll = (TH1F*)hdiffmaxcounts_pedestal[1]->Clone("hdiffmaxcounts_pedestalAll");
  TH1F *hdiffmaxcounts_pedestal_normalizedAll = (TH1F*)hdiffmaxcounts_pedestal_normalized[1]->Clone("hdiffmaxcounts_pedestal_normalizedAll");
  TH1F *hnormalized_integralAll = (TH1F*)hnormalized_integral[1]->Clone("hnormalized_integralAll");
  TH1F *hhittimeAll = (TH1F*)hhittime[1]->Clone("hhittimeAll");
  for(Int_t i=2;i<NSECTORS;i++){
    hdiffmaxcounts_pedestalAll->Add(hdiffmaxcounts_pedestal[1]);
    hdiffmaxcounts_pedestal_normalizedAll->Add(hdiffmaxcounts_pedestal_normalized[1]);
    hnormalized_integralAll->Add(hnormalized_integral[i]);
    hhittimeAll->Add(hhittime[i]);
  }

  // Draw sum of all diffmaxcounts_pedestal
  canvasAll->cd();
  legend->Clear();
  canvasAll->cd()->SetLogy(1);  
  legend->SetY1NDC(0.80);
  legend->SetY2NDC(0.90);
  hdiffmaxcounts_pedestalAll->SetLineColor(kBlack);
  hdiffmaxcounts_pedestalAll->SetMinimum(0.5);
  hdiffmaxcounts_pedestalAll->GetXaxis()->SetNdivisions(505);
  hdiffmaxcounts_pedestalAll->Draw();
  legend->AddEntry(hdiffmaxcounts_pedestalAll,"sum of maxcounts - pedestal","L");
  hdiffmaxcounts_pedestal_normalizedAll->SetLineColor(kRed);
  hdiffmaxcounts_pedestal_normalizedAll->Draw("same");
  legend->AddEntry(hdiffmaxcounts_pedestal_normalizedAll,"sum of maxcounts - pedestal (normalized)","L");
  legend->Draw("same");

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___sum_diffmaxcounts_pedestal.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->cd()->SetLogy(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // Draw sum of all normalized_integral
  canvasAll->cd();
  canvasAll->cd()->SetLogy(1);  
  hnormalized_integralAll->SetMinimum(0.5);
  hnormalized_integralAll->GetXaxis()->SetNdivisions(505);
  hnormalized_integralAll->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___sum_normalized_integral.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  canvasAll->cd()->SetLogy(0);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //                        Draw hit times                          //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  for(Int_t i=0;i<4;i++){
    canvas->cd(i+1)->Clear();
    legend_sub[i]->Clear();
  }
  for(Int_t i=0;i<NSECTORS;i++){
    // 4 subpads, 8 plots each
    Int_t index = i/(NSECTORS/4);
    canvas->cd(index+1);

    hhittime[i]->SetMinimum(0);
    hhittime[i]->SetMaximum(MAX_HITTIME);
    // hhittime[i]->SetAxisRange(0,400);

    hhittime[i]->GetXaxis()->SetNdivisions(505);
    hhittime[i]->GetXaxis()->SetLabelSize(0.050);
    hhittime[i]->GetYaxis()->SetNdivisions(505);
    hhittime[i]->GetYaxis()->SetLabelSize(0.050);
    hhittime[i]->SetLineColor(colors[i % 8]);

    if(i % (NSECTORS/4) == 0)
      hhittime[i]->Draw();
    else
      hhittime[i]->Draw("same");

    sprintf(text,"sec.%2.2d",i+1);
    legend_sub[index]->AddEntry(hhittime[i],text,"L");

    if(i % (NSECTORS/4) == 7)
      legend_sub[index]->Draw("same");
  } // end of loop over sectors

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___hittime.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // Draw all hit times
  canvasAll->cd();
  canvasAll->cd()->SetRightMargin(0.10);
  hhittimeAll->SetMinimum(0);
  hhittimeAll->GetXaxis()->SetNdivisions(505);
  hhittimeAll->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___hittimeAll.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hsector_hittime->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___sector_hittime.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hVarX->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___VarX.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hSecVsVarX->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___SecVsVarX.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hPrimeTime->SetAxisRange(0,100);
  hPrimeTime->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___PrimeTime.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hSectorVsPrimeTime->SetAxisRange(0,100,"X");
  hSectorVsPrimeTime->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___SectorVsPrimeTime.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  // hfitAmpAll->SetAxisRange(0,100);
  hfitAmpAll->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___fitAmpAll.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  // hSectorVsfitAmp->SetAxisRange(0,100,"X");
  hSectorVsfitAmp->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___SectorVsfitAmp.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hfitFrequency->SetAxisRange(0,100);
  hfitFrequency->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___fitFrequency.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hGammaRise->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___GammaRise.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hSectorVsGammaRise->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___SectorVsGammaRise.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hGammaFall->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___GammaFall.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hcVal->Draw();
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___cVal.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  canvasAll->cd();
  hSectorVscVal->Draw("colz");
  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___SectorVscVal.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //                 Draw maxcounts - pedestal vs fitAmp            //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  canvas->Clear();
  canvas->Divide(6,6,.002,.002);
  for(Int_t i=0;i<NSECTORS;i++){
    canvas->cd(i+1);
    canvas->cd(i+1)->Clear();

    hdiffmaxcounts_pedestal_fitAmp[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_fitAmp[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_fitAmp[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_fitAmp[i]->GetYaxis()->SetLabelOffset(0.001);
    hdiffmaxcounts_pedestal_fitAmp[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_fitAmp[i]->Draw("colz");
    
  } // end of loop over sectors

  canvas->cd(33)->Clear();
  canvas->cd(33);
  latex->SetTextSize(0.15);
  sprintf(text,"maxcounts-pedestal vs");
  latex->DrawLatex(0.20,0.65,text);
  sprintf(text,"fitAmp");
  latex->DrawLatex(0.20,0.45,text);
  latex->SetTextSize(0.060);

  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal_fitAmps.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //      Draw maxcounts - pedestal vs normalized integral          //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  canvas->Clear();
  canvas->Divide(6,6,.002,.002);
  for(Int_t i=0;i<NSECTORS;i++){
    canvas->cd(i+1);
    canvas->cd(i+1)->Clear();
    
    hdiffmaxcounts_pedestal_normalized_integral[i]->GetXaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_normalized_integral[i]->GetXaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_normalized_integral[i]->GetYaxis()->SetNdivisions(505);
    hdiffmaxcounts_pedestal_normalized_integral[i]->GetYaxis()->SetLabelOffset(0.001);
    hdiffmaxcounts_pedestal_normalized_integral[i]->GetYaxis()->SetLabelSize(0.050);
    hdiffmaxcounts_pedestal_normalized_integral[i]->Draw("colz");
    
  } // end of loop over sectors

  canvas->cd(33)->Clear();
  canvas->cd(33);
  latex->SetTextSize(0.15);
  sprintf(text,"maxcounts-pedestal vs");
  latex->DrawLatex(0.20,0.65,text);
  sprintf(text,"norm. integral");
  latex->DrawLatex(0.20,0.45,text);
  latex->SetTextSize(0.060);

  sprintf(text,"figures/%s/summary/%3.3d___diffmaxcounts_pedestal_normalized_integrals.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  ////////////////////////////////////////////////////////////////////
  //                                                                //
  //             Draw normalized integral vs fitAmp                 //
  //                                                                //
  ////////////////////////////////////////////////////////////////////
  canvas->cd();
  canvas->Clear();
  canvas->Divide(6,6,.002,.002);
  for(Int_t i=0;i<NSECTORS;i++){
    canvas->cd(i+1);
    canvas->cd(i+1)->Clear();
    
    hnormalized_integral_fitAmp[i]->GetXaxis()->SetNdivisions(505);
    hnormalized_integral_fitAmp[i]->GetXaxis()->SetLabelSize(0.050);
    hnormalized_integral_fitAmp[i]->GetYaxis()->SetNdivisions(505);
    hnormalized_integral_fitAmp[i]->GetYaxis()->SetLabelOffset(0.001);
    hnormalized_integral_fitAmp[i]->GetYaxis()->SetLabelSize(0.050);
    hnormalized_integral_fitAmp[i]->Draw("colz");
    
  } // end of loop over sectors

  canvas->cd(33)->Clear();
  canvas->cd(33);
  latex->SetTextSize(0.15);
  sprintf(text,"normalized integral vs");
  latex->DrawLatex(0.20,0.65,text);
  sprintf(text,"fitAmp");
  latex->DrawLatex(0.20,0.45,text);
  latex->SetTextSize(0.060);

  sprintf(text,"figures/%s/summary/%3.3d___normalized_integral_fitAmps.pdf",basename.c_str(),nplots);
  canvas->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________

  // TLine *line_mintime = new TLine(

  // Sum all hittime histograms without the noisy sector 21
  /*
  TH1F *hhittimeAllNoSector21 = (TH1F*)hhittimeEachSector[0]->Clone("hhittimeAllNoSector21");
  hhittimeAllNoSector21->Reset();
  for(Int_t i=0;i<NSECTORS;i++){
    if(i!=15) hhittimeAllNoSector21->Add(hhittimeEachSector[i]);
  }
  canvasAll->cd();
  hhittimeAllNoSector21->GetXaxis()->SetNdivisions(505);
  hhittimeAllNoSector21->GetYaxis()->SetNdivisions(505);
  hhittimeAllNoSector21->SetFillColor(kBlue);
  hhittimeAllNoSector21->SetFillStyle(1001);
  hhittimeAllNoSector21->SetBarOffset(-0.10);
  hhittimeAllNoSector21->SetBarWidth(0.25);
  hhittimeAllNoSector21->Draw("BAR");

  latex->DrawLatex(0.50,0.80,"sector 21 excluded");

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___hittimeAllNoSector21.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________
  */

  /*
  // From the above histogram the timing cut is determined to be
  // 85 - 105 ns
  // Make projection of sector # vs hit time with this cut
  TH1D *hsectorDistWithTimingCut = (TH1D*)hsector_hittime->ProjectionY("hsectorDistWithTimingCut");
  hsectorDistWithTimingCut->Reset();
  assert(hsectorDistWithTimingCut->GetNbinsX() == NSECTORS);

  const Double_t minTime = 85.;
  const Double_t maxTime = 105.;

  for(Int_t i=0;i<NSECTORS;i++){

    for(Int_t xbin=0;xbin<=hsector_hittime->GetNbinsX()+1;xbin++){
      if(minTime < hsector_hittime->GetXaxis()->GetBinCenter(xbin) && hsector_hittime->GetXaxis()->GetBinCenter(xbin) < maxTime){
	Double_t previous = hsectorDistWithTimingCut->GetBinContent(i+1,xbin);
	Double_t this_bin  = hsector_hittime->GetBinContent(xbin,i+1);

	// cout << "i = " << i << " xbin = " << xbin << " hsector_hittime->GetXaxis()->GetBinCenter(xbin) = " << hsector_hittime->GetXaxis()->GetBinCenter(xbin)
	//      << " previous = " << hsectorDistWithTimingCut->GetBinContent(i+1,xbin) << " this_bin = " << this_bin << endl;

	hsectorDistWithTimingCut->SetBinContent(i+1,previous + this_bin);
      }
    } // end of loop over time bins
  } // end of loop over sectors

  canvasAll->cd();
  hsectorDistWithTimingCut->Draw();

  // Fit with 1 + lambda cos (2 phi)

  sprintf(text,"mkdir -p figures/%s/summary",basename.c_str());
  system(text);
  sprintf(text,"figures/%s/summary/%3.3d___sectorDistWithTimingCut.pdf",basename.c_str(),nplots);
  canvasAll->SaveAs(text);
  nplots++;
  //___________________________________________________________________________________________________________________________
  */

}