RecoLocalMuon/CSCValidation/macros/myFunctions.C - master - githubmirror/cmssw

RecoLocalMuon/CSCValidation/macros/myFunctions.C
2 553 строки80 KB
Tim Cox
apply code-checks and code-format
15 июл 2022, 20:45
15 июл 2022, 20:4581ecb03
// myFunctions.C contains various functions required by makePlots.C
// for making plots of root file output from CSCValidation
// Original author: Andy Kubik (NWU)
// - Updated by many people over the past 15 years
// - This version from Sicheng Wang (UCSB) Jul-2022
// - Please contact CSC DPG for current status - 12.07.2022

void drawColoredChamberLines(int station, int nchamber1[4][36]);

TFile *OpenFiles(std::string path) {
  TFile *f;
  f = new TFile(path.c_str(), "READ");
  return f;
}

void printEmptyChambers(std::string histoname, std::string oname, TFile *f) {
  TH2I *plot = (TH2I *)f->Get(histoname.c_str());
  std::string endcap, chamber;
  int limitr, limitc;
  std::vector<string> deadchambers;

  for (int e = 0; e < 2; e++) {
    for (int s = 0; s < 4; s++) {
      if (s == 0)
        limitr = 4;
      if (s == 1 || s == 2)
        limitr = 2;
      if (s == 3)
        limitr = 1;
      for (int r = 0; r < limitr; r++) {
        if (s == 0)
          limitc = 36;
        if (s != 0 && r == 0)
          limitc = 18;
        if (s != 0 && r == 1)
          limitc = 36;
        for (int c = 0; c < limitc; c++) {
          int type = 0;
          if (s == 0 && r == 0)
            type = 2;
          else if (s == 0 && r == 1)
            type = 3;
          else if (s == 0 && r == 2)
            type = 4;
          else if (s == 0 && r == 3)
            type = 1;
          else
            type = (s + 1) * 2 + (r + 1);
          if (e == 0)
            type = type + 10;
          if (e == 1)
            type = 11 - type;
          int bin = plot->GetBin((c + 1), type);
          float content = plot->GetBinContent(bin);
          std::ostringstream oss;
          if (e == 0)
            endcap = "+";
          if (e == 1)
            endcap = "-";
          oss << "ME " << endcap << (s + 1) << "/" << (r + 1) << "/" << (c + 1);
          chamber = oss.str();
          if (content == 0) {
            if (oname == "wire digis" && (s == 0 && r == 3))
              continue;
            else
              deadchambers.push_back(chamber);
          }
        }
      }
    }
  }

  int n_dc = deadchambers.size();
  ofstream file;
  file.open("deadchamberlist.txt", ios::app);
  file << "Chambers with missing " << oname << "...\n" << endl;
  if (n_dc > 0) {
    for (int n = 0; n < n_dc; n++) {
      file << deadchambers[n] << endl;
    }
  }
  file << "\n\n\n\n";
  file.close();
}

void make1DPlot(std::string histoname, TFile *f1, std::string histotitle, int statoption, std::string savename) {
  TH1F *h1 = (TH1F *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.7);
    h1->GetXaxis()->SetTitleSize(0.06);
    h1->GetYaxis()->SetTitleSize(0.06);
    h1->GetYaxis()->CenterTitle(kTRUE);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
}

void make1DPlot(std::string histoname,
                TFile *f1,
                std::string histotitle,
                std::string xtitle,
                std::string ytitle,
                int statoption,
                std::string savename) {
  TH1F *h1 = (TH1F *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitle(xtitle.c_str());
    h1->GetYaxis()->SetTitle(ytitle.c_str());
    h1->GetXaxis()->SetTitleOffset(0.85);
    h1->GetYaxis()->SetTitleOffset(0.85);
    h1->GetXaxis()->SetTitleSize(0.05);
    h1->GetYaxis()->SetTitleSize(0.05);
    h1->GetYaxis()->CenterTitle(kTRUE);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
}

void make2DPlot(std::string histoname, TFile *f1, std::string histotitle, int statoption, std::string savename) {
  TH2F *h1 = (TH2F *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.7);
    h1->GetXaxis()->SetTitleSize(0.06);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
}

void make2DPlot(std::string histoname,
                TFile *f1,
                std::string histotitle,
                std::string xtitle,
                std::string ytitle,
                int statoption,
                std::string savename) {
  TH2F *h1 = (TH2F *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetTitle(xtitle.c_str());
    h1->GetYaxis()->SetTitle(ytitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.85);
    h1->GetYaxis()->SetTitleOffset(0.85);
    h1->GetXaxis()->SetTitleSize(0.05);
    h1->GetYaxis()->SetTitleSize(0.05);
    h1->GetYaxis()->CenterTitle(kTRUE);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
}

void makeProfile(std::string histoname, TFile *f1, std::string histotitle, int statoption, std::string savename) {
  TProfile *h1 = (TProfile *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->GetYaxis()->SetRangeUser(h1->GetYmin(), h1->GetYmax());
    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.7);
    h1->GetXaxis()->SetTitleSize(0.06);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
}
void makeProfile(std::string histoname,
                 TFile *f1,
                 std::string histotitle,
                 std::string xtitle,
                 std::string ytitle,
                 int statoption,
                 std::string savename) {
  TProfile *h1 = (TProfile *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->GetYaxis()->SetRangeUser(h1->GetYmin(), h1->GetYmax());
    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetTitle(xtitle.c_str());
    h1->GetYaxis()->SetTitle(ytitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.85);
    h1->GetYaxis()->SetTitleOffset(0.85);
    h1->GetXaxis()->SetTitleSize(0.05);
    h1->GetYaxis()->SetTitleSize(0.05);
    h1->GetYaxis()->CenterTitle(kTRUE);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
}

void makeCSCOccupancy(std::string histoname, TFile *f1, std::string histotitle, std::string savename) {
  TH1F *h1 = (TH1F *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1) {
    gStyle->SetOptStat(kFALSE);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    h1->UseCurrentStyle();

    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.7);
    h1->GetXaxis()->SetTitleSize(0.06);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    h1->GetXaxis()->SetBinLabel(2, "Total Events");
    h1->GetXaxis()->SetBinLabel(4, "# Events with Wires");
    h1->GetXaxis()->SetBinLabel(6, "# Events with Strips");
    h1->GetXaxis()->SetBinLabel(8, "# Events with Wires&Strips");
    h1->GetXaxis()->SetBinLabel(10, "# Events with Rechits");
    h1->GetXaxis()->SetBinLabel(12, "# Events with Segments");
    h1->GetXaxis()->SetBinLabel(14, "Events Rejected");

    h1->Draw();

    c->Update();
    c->Print(savename.c_str(), "png");
  }

  delete c;
}

void make1DPlot2(std::string histoname1,
                 std::string histoname2,
                 int statoption,
                 TFile *f1,
                 std::string t1,
                 std::string t2,
                 std::string savename) {
  // use this if you want two plots on one canvas

  TH1F *a1 = (TH1F *)f1->Get(histoname1.c_str());
  TH1F *b1 = (TH1F *)f1->Get(histoname2.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);
  c->Divide(1, 2);

  if (a1) {
    c->cd(1);
    if (statoption == 0)
      gStyle->SetOptStat(kFALSE);
    else
      gStyle->SetOptStat(statoption);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.4);
    gStyle->SetTitleH(0.09);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    a1->UseCurrentStyle();
    a1->SetTitle(t1.c_str());
    a1->GetXaxis()->SetLabelSize(0.06);
    a1->GetYaxis()->SetLabelSize(0.06);
    a1->GetXaxis()->SetTitleOffset(0.7);
    a1->GetXaxis()->SetTitleSize(0.06);
    a1->GetXaxis()->SetNdivisions(208, kTRUE);
    a1->Draw();
  }

  if (b1) {
    gStyle->SetHistFillColor(72);
    b1->UseCurrentStyle();

    //t2 = t2 + " (run " + run + ")";
    b1->SetTitle(t2.c_str());
    b1->GetXaxis()->SetLabelSize(0.06);
    b1->GetYaxis()->SetLabelSize(0.06);
    b1->GetXaxis()->SetTitleOffset(0.7);
    b1->GetXaxis()->SetTitleSize(0.06);
    b1->GetXaxis()->SetNdivisions(508, kTRUE);
    c->cd(2);
    b1->Draw();
  }

  c->Update();
  c->Print(savename.c_str(), "png");
  delete c;
}

void makeEffGif(std::string histoname, TFile *f1, std::string histotitle, std::string savename) {
  TH1F *ho = (TH1F *)f1->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  TH1F *hn = new TH1F("tmp", histotitle.c_str(), 20, 0.5, 20.5);

  if (ho) {
    float Num = 1;
    float Den = 1;
    for (int i = 0; i < 20; i++) {
      Num = ho->GetBinContent(i + 1);
      Den = ho->GetBinContent(i + 21);
      //getEfficiency(Num, Den, eff);
      float Eff = 0.;
      float EffE = 0.;
      if (fabs(Den) > 0.000000001) {
        Eff = Num / Den;
        if (Num < Den) {
          EffE = sqrt((1. - Eff) * Eff / Den);
        }
      }
      hn->SetBinContent(i + 1, Eff);
      hn->SetBinError(i + 1, EffE);
    }

    gStyle->SetOptStat(kFALSE);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    c->SetFillStyle(4000);
    hn->UseCurrentStyle();

    hn->SetTitle(histotitle.c_str());
    hn->GetXaxis()->SetLabelSize(0.04);
    hn->GetYaxis()->SetLabelSize(0.04);
    hn->GetXaxis()->SetTitleOffset(0.7);
    hn->GetXaxis()->SetTitleSize(0.06);
    hn->GetXaxis()->SetNdivisions(208, kTRUE);
    hn->GetYaxis()->SetRangeUser(0.5, 1.1);
    hn->SetMarkerStyle(6);
    hn->GetXaxis()->SetBinLabel(1, "ME +1/1b");
    hn->GetXaxis()->SetBinLabel(2, "ME +1/2");
    hn->GetXaxis()->SetBinLabel(3, "ME +1/3");
    hn->GetXaxis()->SetBinLabel(4, "ME +1/1a");
    hn->GetXaxis()->SetBinLabel(5, "ME +2/1");
    hn->GetXaxis()->SetBinLabel(6, "ME +2/2");
    hn->GetXaxis()->SetBinLabel(7, "ME +3/1");
    hn->GetXaxis()->SetBinLabel(8, "ME +3/2");
    hn->GetXaxis()->SetBinLabel(9, "ME +4/1");
    hn->GetXaxis()->SetBinLabel(10, "ME +4/2");
    hn->GetXaxis()->SetBinLabel(11, "ME -1/1b");
    hn->GetXaxis()->SetBinLabel(12, "ME -1/2");
    hn->GetXaxis()->SetBinLabel(13, "ME -1/3");
    hn->GetXaxis()->SetBinLabel(14, "ME -1/1a");
    hn->GetXaxis()->SetBinLabel(15, "ME -2/1");
    hn->GetXaxis()->SetBinLabel(16, "ME -2/2");
    hn->GetXaxis()->SetBinLabel(17, "ME -3/1");
    hn->GetXaxis()->SetBinLabel(18, "ME -3/2");
    hn->GetXaxis()->SetBinLabel(19, "ME -4/1");
    hn->GetXaxis()->SetBinLabel(20, "ME -4/2");
    hn->Draw();
    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
  delete hn;
}

void Draw2DProfile(std::string histoname, TFile *f1, std::string title, std::string option, std::string savename) {
  TProfile2D *test = (TProfile2D *)f1->Get(histoname.c_str());
  TH2D *plot = test->ProjectionXY("test2", option.c_str());

  if (plot) {
    TCanvas *c = new TCanvas("c", "my canvas", 1);
    gStyle->SetPalette(1, 0);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    plot->SetStats(kFALSE);
    plot->GetYaxis()->SetBinLabel(1, "ME- 4/1");
    plot->GetYaxis()->SetBinLabel(2, "ME- 3/2");
    plot->GetYaxis()->SetBinLabel(3, "ME- 3/1");
    plot->GetYaxis()->SetBinLabel(4, "ME- 2/2");
    plot->GetYaxis()->SetBinLabel(5, "ME- 2/1");
    plot->GetYaxis()->SetBinLabel(6, "ME- 1/3");
    plot->GetYaxis()->SetBinLabel(7, "ME- 1/2");
    plot->GetYaxis()->SetBinLabel(8, "ME- 1/1b");
    plot->GetYaxis()->SetBinLabel(9, "ME- 1/1a");
    plot->GetYaxis()->SetBinLabel(10, "ME+ 1/1a");
    plot->GetYaxis()->SetBinLabel(11, "ME+ 1/1b");
    plot->GetYaxis()->SetBinLabel(12, "ME+ 1/2");
    plot->GetYaxis()->SetBinLabel(13, "ME+ 1/3");
    plot->GetYaxis()->SetBinLabel(14, "ME+ 2/1");
    plot->GetYaxis()->SetBinLabel(15, "ME+ 2/2");
    plot->GetYaxis()->SetBinLabel(16, "ME+ 3/1");
    plot->GetYaxis()->SetBinLabel(17, "ME+ 3/2");
    plot->GetYaxis()->SetBinLabel(18, "ME+ 4/1");

    plot->SetTitle(title.c_str());

    for (int i = 1; i < 37; i++) {
      ostringstream oss1;
      oss1 << i;
      string ch = oss1.str();
      plot->GetXaxis()->SetBinLabel(i, ch.c_str());
    }

    c->SetRightMargin(0.12);
    plot->GetYaxis()->SetNdivisions(20, kFALSE);
    plot->GetXaxis()->SetNdivisions(36, kFALSE);
    plot->GetXaxis()->SetTitle("Chamber #");
    c->SetGrid();

    plot->Draw("colz");
    c->Update();
    c->Print(savename.c_str(), "png");
    delete c;
  }
}

void Draw2DEfficiency(std::string histo, TFile *f1, std::string title, std::string savename) {
  TCanvas *c = new TCanvas("c", "my canvas", 1);
  gStyle->SetPalette(1, 0);
  gPad->SetFillColor(4000);
  c->SetFillStyle(4000);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetOptStat(10);

  TH2F *num = (TH2F *)f1->Get(histo.c_str());
  TH2F *denom = (TH2F *)f1->Get("Efficiency/hEffDenominator");
  if (histo.find("Tight") != std::string::npos)
    denom = (TH2F *)f1->Get("Efficiency/hEffDenominatorTight");

  TH2F *plot = new TH2F("plot", title.c_str(), 36, 0.5, 36.5, 20, 0.5, 20.5);

  if (num && denom)
    plot->Divide(num, denom, 1., 1., "B");

  //plot->SetStats(kFALSE);

  plot->GetYaxis()->SetBinLabel(1, "ME- 4/2");  //jun28
  plot->GetYaxis()->SetBinLabel(2, "ME- 4/1");
  plot->GetYaxis()->SetBinLabel(3, "ME- 3/2");
  plot->GetYaxis()->SetBinLabel(4, "ME- 3/1");
  plot->GetYaxis()->SetBinLabel(5, "ME- 2/2");
  plot->GetYaxis()->SetBinLabel(6, "ME- 2/1");
  plot->GetYaxis()->SetBinLabel(10, "ME- 1/1a");
  plot->GetYaxis()->SetBinLabel(7, "ME- 1/3");
  plot->GetYaxis()->SetBinLabel(8, "ME- 1/2");
  plot->GetYaxis()->SetBinLabel(9, "ME- 1/1b");
  plot->GetYaxis()->SetBinLabel(12, "ME+ 1/1b");
  plot->GetYaxis()->SetBinLabel(13, "ME+ 1/2");
  plot->GetYaxis()->SetBinLabel(14, "ME+ 1/3");
  plot->GetYaxis()->SetBinLabel(11, "ME+ 1/1a");
  plot->GetYaxis()->SetBinLabel(15, "ME+ 2/1");
  plot->GetYaxis()->SetBinLabel(16, "ME+ 2/2");
  plot->GetYaxis()->SetBinLabel(17, "ME+ 3/1");
  plot->GetYaxis()->SetBinLabel(18, "ME+ 3/2");
  plot->GetYaxis()->SetBinLabel(19, "ME+ 4/1");
  plot->GetYaxis()->SetBinLabel(20, "ME+ 4/2");  //jun28

  for (int i = 1; i < 37; i++) {
    ostringstream oss1;
    oss1 << i;
    string ch = oss1.str();
    plot->GetXaxis()->SetBinLabel(i, ch.c_str());
  }

  c->SetRightMargin(0.12);

  plot->GetYaxis()->SetNdivisions(20, kFALSE);
  plot->GetXaxis()->SetNdivisions(36, kFALSE);

  plot->GetXaxis()->SetTitle("Chamber #");

  c->SetGrid();

  plot->Draw("COLZ");

  c->Update();
  c->Print(savename.c_str(), "png");
  delete c;
  delete plot;
}

void Draw2DTempPlot(std::string histo, TFile *f1, bool includeME11, std::string savename, bool hasLabels = false) {
  TCanvas *c = new TCanvas("c", "my canvas", 1);
  gStyle->SetPalette(1, 0);
  gPad->SetFillColor(4000);
  c->SetFillStyle(4000);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetOptStat(10);

  TH2I *plot = (TH2I *)f1->Get(histo.c_str());

  //plot->SetStats(kFALSE);

  if (!hasLabels) {
    if (includeME11) {
      plot->GetYaxis()->SetBinLabel(1, "ME- 4/2");
      plot->GetYaxis()->SetBinLabel(2, "ME- 4/1");
      plot->GetYaxis()->SetBinLabel(3, "ME- 3/2");
      plot->GetYaxis()->SetBinLabel(4, "ME- 3/1");
      plot->GetYaxis()->SetBinLabel(5, "ME- 2/2");
      plot->GetYaxis()->SetBinLabel(6, "ME- 2/1");
      plot->GetYaxis()->SetBinLabel(10, "ME- 1/1a");
      plot->GetYaxis()->SetBinLabel(7, "ME- 1/3");
      plot->GetYaxis()->SetBinLabel(8, "ME- 1/2");
      plot->GetYaxis()->SetBinLabel(9, "ME- 1/1b");
      plot->GetYaxis()->SetBinLabel(12, "ME+ 1/1b");
      plot->GetYaxis()->SetBinLabel(13, "ME+ 1/2");
      plot->GetYaxis()->SetBinLabel(14, "ME+ 1/3");
      plot->GetYaxis()->SetBinLabel(11, "ME+ 1/1a");
      plot->GetYaxis()->SetBinLabel(15, "ME+ 2/1");
      plot->GetYaxis()->SetBinLabel(16, "ME+ 2/2");
      plot->GetYaxis()->SetBinLabel(17, "ME+ 3/1");
      plot->GetYaxis()->SetBinLabel(18, "ME+ 3/2");
      plot->GetYaxis()->SetBinLabel(19, "ME+ 4/1");
      plot->GetYaxis()->SetBinLabel(20, "ME+ 4/2");
    } else {
      plot->GetYaxis()->SetBinLabel(1, "ME- 4/1");
      plot->GetYaxis()->SetBinLabel(2, "ME- 3/2");
      plot->GetYaxis()->SetBinLabel(3, "ME- 3/1");
      plot->GetYaxis()->SetBinLabel(4, "ME- 2/2");
      plot->GetYaxis()->SetBinLabel(5, "ME- 2/1");
      plot->GetYaxis()->SetBinLabel(6, "ME- 1/3");
      plot->GetYaxis()->SetBinLabel(7, "ME- 1/2");
      plot->GetYaxis()->SetBinLabel(8, "ME- 1/1b");
      plot->GetYaxis()->SetBinLabel(9, "ME- 1/1a");
      plot->GetYaxis()->SetBinLabel(10, "ME+ 1/1a");
      plot->GetYaxis()->SetBinLabel(11, "ME+ 1/1b");
      plot->GetYaxis()->SetBinLabel(12, "ME+ 1/2");
      plot->GetYaxis()->SetBinLabel(13, "ME+ 1/3");
      plot->GetYaxis()->SetBinLabel(14, "ME+ 2/1");
      plot->GetYaxis()->SetBinLabel(15, "ME+ 2/2");
      plot->GetYaxis()->SetBinLabel(16, "ME+ 3/1");
      plot->GetYaxis()->SetBinLabel(17, "ME+ 3/2");
      plot->GetYaxis()->SetBinLabel(18, "ME+ 4/1");
    }

    for (int i = 1; i < 37; i++) {
      ostringstream oss1;
      oss1 << i;
      string ch = oss1.str();
      plot->GetXaxis()->SetBinLabel(i, ch.c_str());
    }

    plot->GetYaxis()->SetNdivisions(20, kFALSE);
    plot->GetXaxis()->SetNdivisions(36, kFALSE);

    plot->GetXaxis()->SetTitle("Chamber #");

    c->SetGrid();
  }

  c->SetRightMargin(0.12);
  plot->Draw("COLZ");

  c->Update();
  c->Print(savename.c_str(), "png");
  delete c;
}

void Draw2DTempPlot2(std::string histo, TFile *f1, bool includeME11, std::string savename) {
  TCanvas *c = new TCanvas("c", "my canvas", 1);
  gStyle->SetPalette(1, 0);
  gPad->SetFillColor(4000);
  c->SetFillStyle(4000);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);

  TH2F *plot = (TH2F *)f1->Get(histo.c_str());

  plot->SetStats(kFALSE);

  if (includeME11) {
    plot->GetYaxis()->SetBinLabel(1, "ME- 4/2");
    plot->GetYaxis()->SetBinLabel(2, "ME- 4/1");
    plot->GetYaxis()->SetBinLabel(3, "ME- 3/2");
    plot->GetYaxis()->SetBinLabel(4, "ME- 3/1");
    plot->GetYaxis()->SetBinLabel(5, "ME- 2/2");
    plot->GetYaxis()->SetBinLabel(6, "ME- 2/1");
    plot->GetYaxis()->SetBinLabel(10, "ME- 1/1a");
    plot->GetYaxis()->SetBinLabel(7, "ME- 1/3");
    plot->GetYaxis()->SetBinLabel(8, "ME- 1/2");
    plot->GetYaxis()->SetBinLabel(9, "ME- 1/1b");
    plot->GetYaxis()->SetBinLabel(12, "ME+ 1/1b");
    plot->GetYaxis()->SetBinLabel(13, "ME+ 1/2");
    plot->GetYaxis()->SetBinLabel(14, "ME+ 1/3");
    plot->GetYaxis()->SetBinLabel(11, "ME+ 1/1a");
    plot->GetYaxis()->SetBinLabel(15, "ME+ 2/1");
    plot->GetYaxis()->SetBinLabel(16, "ME+ 2/2");
    plot->GetYaxis()->SetBinLabel(17, "ME+ 3/1");
    plot->GetYaxis()->SetBinLabel(18, "ME+ 3/2");
    plot->GetYaxis()->SetBinLabel(19, "ME+ 4/1");
    plot->GetYaxis()->SetBinLabel(20, "ME+ 4/2");
  } else {
    plot->GetYaxis()->SetBinLabel(1, "ME- 4/1");
    plot->GetYaxis()->SetBinLabel(2, "ME- 3/2");
    plot->GetYaxis()->SetBinLabel(3, "ME- 3/1");
    plot->GetYaxis()->SetBinLabel(4, "ME- 2/2");
    plot->GetYaxis()->SetBinLabel(5, "ME- 2/1");
    plot->GetYaxis()->SetBinLabel(6, "ME- 1/3");
    plot->GetYaxis()->SetBinLabel(7, "ME- 1/2");
    plot->GetYaxis()->SetBinLabel(8, "ME- 1/1b");
    plot->GetYaxis()->SetBinLabel(9, "ME- 1/1a");
    plot->GetYaxis()->SetBinLabel(10, "ME+ 1/1a");
    plot->GetYaxis()->SetBinLabel(11, "ME+ 1/1b");
    plot->GetYaxis()->SetBinLabel(12, "ME+ 1/2");
    plot->GetYaxis()->SetBinLabel(13, "ME+ 1/3");
    plot->GetYaxis()->SetBinLabel(14, "ME+ 2/1");
    plot->GetYaxis()->SetBinLabel(15, "ME+ 2/2");
    plot->GetYaxis()->SetBinLabel(16, "ME+ 3/1");
    plot->GetYaxis()->SetBinLabel(17, "ME+ 3/2");
    plot->GetYaxis()->SetBinLabel(18, "ME+ 4/1");
  }

  for (int i = 1; i < 37; i++) {
    ostringstream oss1;
    oss1 << i;
    string ch = oss1.str();
    plot->GetXaxis()->SetBinLabel(i, ch.c_str());
  }

  c->SetRightMargin(0.12);

  plot->GetYaxis()->SetNdivisions(20, kFALSE);
  plot->GetXaxis()->SetNdivisions(36, kFALSE);

  plot->GetXaxis()->SetTitle("Chamber #");

  c->SetGrid();

  plot->Draw("COLZ");

  c->Update();
  c->Print(savename.c_str(), "png");
  delete c;
}

void GlobalPosfromTree(
    std::string graphname, TFile *f1, int endcap, int station, std::string type, std::string savename) {
  TTree *t1;
  TBranch *b1;
  struct posRecord {
    int endcap;
    int station;
    int ring;
    int chamber;
    int layer;
    float localx;
    float localy;
    float globalx;
    float globaly;
  } points;

  if (type == "rechit") {
    t1 = (TTree *)f1->Get("recHits/rHPositions");
    b1 = t1->GetBranch("rHpos");
    b1->SetAddress(&points);
  }

  if (type == "segment") {
    t1 = (TTree *)f1->Get("Segments/segPositions");
    b1 = t1->GetBranch("segpos");
    b1->SetAddress(&points);
  }

  int n1 = (int)t1->GetEntries();
  const int nevents1 = n1;
  TVectorF globx1(nevents1);
  TVectorF globy1(nevents1);
  int nstation1 = 0;
  const int num_of_rings = 4;
  const int num_of_chambers = 36;
  int nchamber1[num_of_rings][num_of_chambers];
  for (int i = 0; i < num_of_rings; i++) {
    for (int j = 0; j < num_of_chambers; j++) {
      nchamber1[i][j] = 0;
    }
  }

  for (int i = 0; i < nevents1; i++) {
    b1->GetEntry(i);
    if (points.station == station && points.endcap == endcap) {
      globx1[nstation1] = points.globalx;
      globy1[nstation1] = points.globaly;
      nstation1++;
      nchamber1[points.ring - 1][points.chamber - 1]++;
    }
  }

  TCanvas *c = new TCanvas("c", "my canvas", 1);
  c->SetCanvasSize(700, 700);
  //TGraph *graph1 = new TGraph(nstation1,globx1,globy1);
  TGraph *graph1 = new TGraph(globx1, globy1);

  std::string name1 = graphname;

  gStyle->SetPalette(1, 0);
  gStyle->SetTitleW(0.9);
  gStyle->SetTitleH(0.1);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);
  gPad->SetFillColor(4000);
  c->SetFillStyle(4000);

  //asking to show the total number of entries for rechit and seg plots
  //i.e. SetOptStat(10), returns 0 entries for some reason
  gStyle->SetOptStat(0);

  graph1->GetXaxis()->SetLimits(-720, 720);
  graph1->GetYaxis()->SetLimits(-720, 720);
  graph1->GetXaxis()->SetRangeUser(-720, 720);
  graph1->GetYaxis()->SetRangeUser(-720, 720);

  graph1->SetTitle(name1.c_str());
  graph1->UseCurrentStyle();
  graph1->Draw("AP");

  //  drawChamberLines(station);
  drawColoredChamberLines(station, nchamber1);

  c->Print(savename.c_str(), "png");
  delete c;

}  // end GlobalPosfromTree

void GlobalPosfromChain(std::string graphname,
                        TChain *ch1,
                        int endcap,
                        int station,
                        std::string type,
                        std::string savename,
                        const int maxNumPoints = 10'000'000) {
  struct posRecord {
    int endcap;
    int station;
    int ring;
    int chamber;
    int layer;
    float localx;
    float localy;
    float globalx;
    float globaly;
  } points;

  if (type == "rechit")
    ch1->SetBranchAddress("rHpos", &points);
  else if (type == "segment")
    ch1->SetBranchAddress("segpos", &points);

  int n1 = ch1->GetEntries();
  const int nevents1 = TMath::Min(n1, maxNumPoints);
  TVectorF globx1(nevents1);
  TVectorF globy1(nevents1);
  int nstation1 = 0;
  const int num_of_rings = 4;
  const int num_of_chambers = 36;
  int nchamber1[num_of_rings][num_of_chambers];
  for (int i = 0; i < num_of_rings; i++) {
    for (int j = 0; j < num_of_chambers; j++) {
      nchamber1[i][j] = 0;
    }
  }

  for (int i = 0; i < nevents1; i++) {
    ch1->GetEntry(i);
    if (points.station == station && points.endcap == endcap) {
      globx1[nstation1] = points.globalx;
      globy1[nstation1] = points.globaly;
      nstation1++;
      nchamber1[points.ring - 1][points.chamber - 1]++;
    }
  }

  TCanvas *c = new TCanvas("c", "my canvas", 1);
  c->SetCanvasSize(700, 700);
  gStyle->SetPalette(1, 0);
  gStyle->SetTitleW(0.9);
  gStyle->SetTitleH(0.1);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);
  gPad->SetFillColor(4000);
  c->SetFillStyle(4000);

  //asking to show the total number of entries for rechit and seg plots
  //i.e. SetOptStat(10), returns 0 entries for some reason
  gStyle->SetOptStat(0);

  TGraph *graph1 = new TGraph(globx1, globy1);
  std::string name1 = graphname;
  graph1->GetXaxis()->SetLimits(-720, 720);
  graph1->GetYaxis()->SetLimits(-720, 720);
  graph1->GetXaxis()->SetRangeUser(-720, 720);
  graph1->GetYaxis()->SetRangeUser(-720, 720);

  graph1->SetTitle(name1.c_str());
  graph1->UseCurrentStyle();
  graph1->Draw("AP");

  drawColoredChamberLines(station, nchamber1);

  c->Print(savename.c_str(), "png");

}  // end GlobalPosfromChain

void drawChamberLines(int station) {
  gStyle->SetLineWidth(2);
  float pi = 3.14159;
  TVector3 x(0, 0, 1);
  int linecolor = 1;
  //for alternating colors, set 2 diff colors here
  int lc1 = 1;
  int lc2 = 1;

  if (station == 1) {
    TVector3 p1(101, 9.361, 0);
    TVector3 p2(101, -9.361, 0);
    TVector3 p3(260, -22.353, 0);
    TVector3 p4(260, 22.353, 0);

    TLine *line1;
    TLine *line2;
    TLine *line3;
    TLine *line4;

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc1)
        linecolor = lc2;
      else
        linecolor = lc1;

      line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
      line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
      line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
      line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      p1.Rotate(pi / 18, x);
      p2.Rotate(pi / 18, x);
      p3.Rotate(pi / 18, x);
      p4.Rotate(pi / 18, x);
    }

    TVector3 q1(281.49, 25.5, 0);
    TVector3 q2(281.49, -25.5, 0);
    TVector3 q3(455.99, -41.87, 0);
    TVector3 q4(455.99, 41.87, 0);

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc2)
        linecolor = lc1;
      else
        linecolor = lc2;

      line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
      line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
      line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
      line4 = new TLine(q4(0), q4(1), q1(0), q1(1));

      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      q1.Rotate(pi / 18, x);
      q2.Rotate(pi / 18, x);
      q3.Rotate(pi / 18, x);
      q4.Rotate(pi / 18, x);
    }

    TVector3 r1(511.99, 31.7, 0);
    TVector3 r2(511.99, -31.7, 0);
    TVector3 r3(676.15, -46.05, 0);
    TVector3 r4(676.15, 46.05, 0);

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc1)
        linecolor = lc2;
      else
        linecolor = lc1;

      line1 = new TLine(r1(0), r1(1), r2(0), r2(1));
      line2 = new TLine(r2(0), r2(1), r3(0), r3(1));
      line3 = new TLine(r3(0), r3(1), r4(0), r4(1));
      line4 = new TLine(r4(0), r4(1), r1(0), r1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      r1.Rotate(pi / 18, x);
      r2.Rotate(pi / 18, x);
      r3.Rotate(pi / 18, x);
      r4.Rotate(pi / 18, x);
    }
  }

  if (station == 2) {
    TVector3 p1(146.9, 27.0, 0);
    TVector3 p2(146.9, -27.0, 0);
    TVector3 p3(336.56, -62.855, 0);
    TVector3 p4(336.56, 62.855, 0);

    p1.Rotate(pi / 36, x);
    p2.Rotate(pi / 36, x);
    p3.Rotate(pi / 36, x);
    p4.Rotate(pi / 36, x);

    TLine *line1;
    TLine *line2;
    TLine *line3;
    TLine *line4;

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc1)
        linecolor = lc2;
      else
        linecolor = lc1;

      line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
      line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
      line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
      line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      p1.Rotate(pi / 9, x);
      p2.Rotate(pi / 9, x);
      p3.Rotate(pi / 9, x);
      p4.Rotate(pi / 9, x);
    }

    TVector3 q1(364.02, 33.23, 0);
    TVector3 q2(364.02, -33.23, 0);
    TVector3 q3(687.08, -63.575, 0);
    TVector3 q4(687.08, 63.575, 0);

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc2)
        linecolor = lc1;
      else
        linecolor = lc2;

      line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
      line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
      line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
      line4 = new TLine(q4(0), q4(1), q1(0), q1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      q1.Rotate(pi / 18, x);
      q2.Rotate(pi / 18, x);
      q3.Rotate(pi / 18, x);
      q4.Rotate(pi / 18, x);
    }
  }

  if (station == 3) {
    TVector3 p1(166.89, 30.7, 0);
    TVector3 p2(166.89, -30.7, 0);
    TVector3 p3(336.59, -62.855, 0);
    TVector3 p4(336.59, 62.855, 0);

    p1.Rotate(pi / 36, x);
    p2.Rotate(pi / 36, x);
    p3.Rotate(pi / 36, x);
    p4.Rotate(pi / 36, x);

    TLine *line1;
    TLine *line2;
    TLine *line3;
    TLine *line4;

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc1)
        linecolor = lc2;
      else
        linecolor = lc1;

      line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
      line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
      line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
      line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      p1.Rotate(pi / 9, x);
      p2.Rotate(pi / 9, x);
      p3.Rotate(pi / 9, x);
      p4.Rotate(pi / 9, x);
    }

    TVector3 q1(364.02, 33.23, 0);
    TVector3 q2(364.02, -33.23, 0);
    TVector3 q3(687.08, -63.575, 0);
    TVector3 q4(687.08, 63.575, 0);

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc2)
        linecolor = lc1;
      else
        linecolor = lc2;

      line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
      line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
      line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
      line4 = new TLine(q4(0), q4(1), q1(0), q1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      q1.Rotate(pi / 18, x);
      q2.Rotate(pi / 18, x);
      q3.Rotate(pi / 18, x);
      q4.Rotate(pi / 18, x);
    }
  }

  if (station == 4) {
    TVector3 p1(186.99, 34.505, 0);
    TVector3 p2(186.99, -34.505, 0);
    TVector3 p3(336.41, -62.825, 0);
    TVector3 p4(336.41, 62.825, 0);

    p1.Rotate(pi / 36, x);
    p2.Rotate(pi / 36, x);
    p3.Rotate(pi / 36, x);
    p4.Rotate(pi / 36, x);

    TLine *line1;
    TLine *line2;
    TLine *line3;
    TLine *line4;

    for (int i = 0; i < 36; i++) {
      if (linecolor == lc1)
        linecolor = lc2;
      else
        linecolor = lc1;

      line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
      line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
      line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
      line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
      line1->SetLineColor(linecolor);
      line2->SetLineColor(linecolor);
      line3->SetLineColor(linecolor);
      line4->SetLineColor(linecolor);

      line1->Draw();
      line2->Draw();
      line3->Draw();
      line4->Draw();

      p1.Rotate(pi / 9, x);
      p2.Rotate(pi / 9, x);
      p3.Rotate(pi / 9, x);
      p4.Rotate(pi / 9, x);
    }
  }
}

void compare1DPlot(
    std::string histoname, TFile *f1, TFile *f2, std::string histotitle, int statoption, std::string savename) {
  // used to compare two of the same histogram from different releases/runs/etc

  TH1F *h2 = (TH1F *)f1->Get(histoname.c_str());
  TH1F *h1 = (TH1F *)f2->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  if (h1 && h2) {
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.5);
    gStyle->SetTitleH(0.07);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    gPad->SetFillColor(4000);
    c->SetFillStyle(4000);
    gStyle->SetOptStat(statoption);
    h1->UseCurrentStyle();
    h2->UseCurrentStyle();
    h2->SetFillColor(52);

    h1->SetTitle(histotitle.c_str());
    h1->GetXaxis()->SetLabelSize(0.04);
    h1->GetYaxis()->SetLabelSize(0.04);
    h1->GetXaxis()->SetTitleOffset(0.7);
    h1->GetXaxis()->SetTitleSize(0.06);
    h1->GetXaxis()->SetNdivisions(208, kTRUE);

    TLegend *leg = new TLegend(0.79, 0.74, 0.89, 0.84);
    leg->AddEntry(h1, "ref", "f");
    leg->AddEntry(h2, "new", "l");

    h1->Draw();
    h2->Draw("same e");
    leg->Draw();
  }

  c->Update();
  c->Print(savename.c_str(), "png");
  delete c;
}

void compareEffGif(std::string histoname, TFile *f1, TFile *f2, std::string histotitle, std::string savename) {
  TH1F *h1 = (TH1F *)f1->Get(histoname.c_str());
  TH1F *h2 = (TH1F *)f2->Get(histoname.c_str());

  TCanvas *c = new TCanvas("c", "my canvas", 1);

  TH1F *hn1 = new TH1F("tmp1", histotitle.c_str(), 20, 0.5, 20.5);
  TH1F *hn2 = new TH1F("tmp2", histotitle.c_str(), 20, 0.5, 20.5);

  if (h1 && h2) {
    float Num = 1;
    float Den = 1;
    for (int i = 0; i < 20; i++) {
      Num = h1->GetBinContent(i + 1);
      Den = h1->GetBinContent(i + 21);
      //getEfficiency(Num, Den, eff);
      float Eff = 0.;
      float EffE = 0.;
      if (fabs(Den) > 0.000000001) {
        Eff = Num / Den;
        if (Num < Den) {
          EffE = sqrt((1. - Eff) * Eff / Den);
        }
      }
      hn1->SetBinContent(i + 1, Eff);
      hn1->SetBinError(i + 1, EffE);
    }

    Num = 1;
    Den = 1;
    for (int i = 0; i < 20; i++) {
      Num = h2->GetBinContent(i + 1);
      Den = h2->GetBinContent(i + 21);
      //getEfficiency(Num, Den, eff);
      float Eff = 0.;
      float EffE = 0.;
      if (fabs(Den) > 0.000000001) {
        Eff = Num / Den;
        if (Num < Den) {
          EffE = sqrt((1. - Eff) * Eff / Den);
        }
      }
      hn2->SetBinContent(i + 1, Eff);
      hn2->SetBinError(i + 1, EffE);
    }

    gStyle->SetOptStat(kFALSE);
    gStyle->SetHistFillColor(92);
    gStyle->SetFrameFillColor(4000);
    gStyle->SetTitleW(0.7);
    gStyle->SetTitleH(0.07);
    gPad->SetFillColor(4000);
    gStyle->SetStatColor(0);
    gStyle->SetTitleFillColor(0);
    c->SetFillStyle(4000);
    hn1->UseCurrentStyle();

    hn1->SetTitle(histotitle.c_str());
    hn1->GetXaxis()->SetLabelSize(0.04);
    hn1->GetYaxis()->SetLabelSize(0.04);
    hn1->GetXaxis()->SetTitleOffset(0.7);
    hn1->GetXaxis()->SetTitleSize(0.06);
    hn1->GetXaxis()->SetNdivisions(208, kTRUE);
    hn1->GetYaxis()->SetRangeUser(0.5, 1.1);
    hn1->SetMarkerStyle(6);
    hn1->SetMarkerColor(kBlue);
    hn2->SetMarkerStyle(6);
    hn2->SetMarkerColor(kRed);
    hn1->GetXaxis()->SetBinLabel(1, "ME +1/1b");
    hn1->GetXaxis()->SetBinLabel(2, "ME +1/2");
    hn1->GetXaxis()->SetBinLabel(3, "ME +1/3");
    hn1->GetXaxis()->SetBinLabel(4, "ME +1/1a");
    hn1->GetXaxis()->SetBinLabel(5, "ME +2/1");
    hn1->GetXaxis()->SetBinLabel(6, "ME +2/2");
    hn1->GetXaxis()->SetBinLabel(7, "ME +3/1");
    hn1->GetXaxis()->SetBinLabel(8, "ME +3/2");
    hn1->GetXaxis()->SetBinLabel(9, "ME +4/1");
    hn1->GetXaxis()->SetBinLabel(10, "ME +4/2");
    hn1->GetXaxis()->SetBinLabel(11, "ME -1/1b");
    hn1->GetXaxis()->SetBinLabel(12, "ME -1/2");
    hn1->GetXaxis()->SetBinLabel(13, "ME -1/3");
    hn1->GetXaxis()->SetBinLabel(14, "ME -1/1a");
    hn1->GetXaxis()->SetBinLabel(15, "ME -2/1");
    hn1->GetXaxis()->SetBinLabel(16, "ME -2/2");
    hn1->GetXaxis()->SetBinLabel(17, "ME -3/1");
    hn1->GetXaxis()->SetBinLabel(18, "ME -3/2");
    hn1->GetXaxis()->SetBinLabel(19, "ME -4/1");
    hn1->GetXaxis()->SetBinLabel(20, "ME -4/2");

    TLegend *leg = new TLegend(0.79, 0.79, 0.89, 0.89);
    leg->AddEntry(hn1, "new", "p");
    leg->AddEntry(hn2, "ref", "p");

    hn1->Draw();
    hn2->Draw("same");
    leg->Draw();
    c->Update();
    c->Print(savename.c_str(), "png");
  }
  delete c;
  delete hn1;
  delete hn2;
}

void GlobalPosfromTreeCompare(
    std::string graphname, TFile *f1, TFile *f2, int endcap, int station, std::string type, std::string savename) {
  struct posRecord {
    int endcap;
    int station;
    int ring;
    int chamber;
    int layer;
    float localx;
    float localy;
    float globalx;
    float globaly;
  } points1, points2;

  TTree *t1;
  TTree *t2;
  TBranch *b1;
  TBranch *b2;

  if (type == "rechit") {
    t1 = (TTree *)f1->Get("recHits/rHPositions");
    t2 = (TTree *)f2->Get("recHits/rHPositions");
    b1 = t1->GetBranch("rHpos");
    b2 = t2->GetBranch("rHpos");
    b1->SetAddress(&points1);
    b2->SetAddress(&points2);
  }

  if (type == "segment") {
    t1 = (TTree *)f1->Get("Segments/segPositions");
    t2 = (TTree *)f2->Get("Segments/segPositions");
    b1 = t1->GetBranch("segpos");
    b2 = t2->GetBranch("segpos");
    b1->SetAddress(&points1);
    b2->SetAddress(&points2);
  }

  int n1 = (int)t1->GetEntries();
  int n2 = (int)t2->GetEntries();

  const int nevents1 = n1;
  const int nevents2 = n2;

  float globx1[nevents1];
  float globy1[nevents1];
  float globx2[nevents2];
  float globy2[nevents2];
  int nstation1 = 0;
  int nstation2 = 0;

  for (int i = 0; i < nevents1; i++) {
    b1->GetEntry(i);
    if (points1.station == station && points1.endcap == endcap) {
      globx1[nstation1] = points1.globalx;
      globy1[nstation1] = points1.globaly;
      nstation1++;
    }
  }
  for (int i = 0; i < nevents2; i++) {
    b2->GetEntry(i);
    if (points2.station == station && points2.endcap == endcap) {
      globx2[nstation2] = points2.globalx;
      globy2[nstation2] = points2.globaly;
      nstation2++;
    }
  }

  std::string name1 = graphname + " (New)";
  std::string name2 = graphname + " (Ref)";
  TCanvas *c = new TCanvas("c", "my canvas", 1);
  c->SetCanvasSize(1300, 700);
  c->Divide(2, 1);
  TGraph *graph1 = new TGraph(nstation1, globx1, globy1);
  TGraph *graph2 = new TGraph(nstation2, globx2, globy2);

  gStyle->SetTitleW(0.6);
  gStyle->SetTitleH(0.1);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetOptStat(10);

  c->cd(1);
  graph1->SetTitle(name1.c_str());
  graph1->UseCurrentStyle();
  drawChamberLines(station);
  graph1->Draw("AP");
  c->cd(2);
  graph2->SetTitle(name2.c_str());
  graph2->UseCurrentStyle();
  drawChamberLines(station);
  graph2->Draw("AP");

  //c->Update();
  c->Print(savename.c_str(), "png");
  delete c;
}

void NikolaiPlots(TFile *f_in, int flag) {
  gROOT->SetStyle("Plain");  // to get rid of gray color of pad and have it white
  gStyle->SetPalette(1, 0);  //

  std::ostringstream ss, ss1;

  std::string folder;
  std::string input_histName, input_title_X, input_title_Y, slice_title_X;
  Int_t ny;
  Float_t ylow, yhigh;
  std::string result_histName, result_histTitle, result_title_Y, result_histNameEntries, result_histTitleEntries;

  if (flag == 1) {  // gas gain results
    folder = "GasGain/";

    input_histName = "gas_gain_rechit_adc_3_3_sum_location_ME_";
    input_title_X = "Location=(layer-1)*nsegm+segm";
    input_title_Y = "3X3 ADC Sum";

    slice_title_X = "3X3 ADC Sum Location";

    ny = 30;
    ylow = 1.0;
    yhigh = 31.0;
    result_histName = "mean_gas_gain_vs_location_csc_ME_";
    result_histTitle = "Mean 3X3 ADC Sum";
    result_title_Y = "Location=(layer-1)*nsegm+segm";

    result_histNameEntries = "entries_gas_gain_vs_location_csc_ME_";
    result_histTitleEntries = "Entries 3X3 ADC Sum";
  }

  if (flag == 2) {  // AFEB timing results
    folder = "AFEBTiming/";

    input_histName = "afeb_time_bin_vs_afeb_occupancy_ME_";
    input_title_X = "AFEB";
    input_title_Y = "Time Bin";

    slice_title_X = "AFEB";

    ny = 42;
    ylow = 1.0;
    yhigh = 42.0;
    result_histName = "mean_afeb_time_bin_vs_afeb_csc_ME_";
    result_histTitle = "AFEB Mean Time Bin";
    result_title_Y = "AFEB";

    result_histNameEntries = "entries_afeb_time_bin_vs_afeb_csc_ME_";
    result_histTitleEntries = "Entries AFEB Time Bin";
  }

  if (flag == 3) {  // Comparator timing results
    folder = "CompTiming/";

    input_histName = "comp_time_bin_vs_cfeb_occupancy_ME_";
    input_title_X = "CFEB";
    input_title_Y = "Time Bin";

    slice_title_X = "CFEB";

    ny = 5;
    ylow = 1.0;
    yhigh = 6.0;
    result_histName = "mean_comp_time_bin_vs_cfeb_csc_ME_";
    result_histTitle = "Comparator Mean Time Bin";
    result_title_Y = "CFEB";

    result_histNameEntries = "entries_comp_time_bin_vs_cfeb_csc_ME_";
    result_histTitleEntries = "Entries Comparator Time Bin";
  }

  if (flag == 4) {  // Strip ADC timing results
    folder = "ADCTiming/";

    input_histName = "adc_3_3_weight_time_bin_vs_cfeb_occupancy_ME_";
    input_title_X = "CFEB";
    input_title_Y = "Time Bin";

    slice_title_X = "CFEB";

    ny = 5;
    ylow = 1.0;
    yhigh = 6.0;
    result_histName = "mean_adc_time_bin_vs_cfeb_csc_ME_";
    result_histTitle = "ADC 3X3 Mean Time Bin";
    result_title_Y = "CFEB";

    result_histNameEntries = "entries_adc_time_bin_vs_cfeb_csc_ME_";
    result_histTitleEntries = "Entries ADC 3X3 Time Bin";
  }

  std::vector<std::string> xTitle;
  xTitle.push_back("ME+1/1 CSC");
  xTitle.push_back("ME+1/2 CSC");
  xTitle.push_back("ME+1/3 CSC");
  xTitle.push_back("ME+2/1 CSC");
  xTitle.push_back("ME+2/2 CSC");
  xTitle.push_back("ME+3/1 CSC");
  xTitle.push_back("ME+3/2 CSC");
  xTitle.push_back("ME+4/1 CSC");
  xTitle.push_back("ME+4/2 CSC");
  xTitle.push_back("ME-1/1 CSC");
  xTitle.push_back("ME-1/2 CSC");
  xTitle.push_back("ME-1/3 CSC");
  xTitle.push_back("ME-2/1 CSC");
  xTitle.push_back("ME-2/2 CSC");
  xTitle.push_back("ME-3/1 CSC");
  xTitle.push_back("ME-3/2 CSC");
  xTitle.push_back("ME-4/1 CSC");
  xTitle.push_back("ME-4/2 CSC");

  TH2F *h2[500];
  TH2F *h;
  Int_t esr[18] = {111, 112, 113, 121, 122, 131, 132, 141, 142, 211, 212, 213, 221, 222, 231, 232, 241, 242};
  Int_t entries[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  Int_t k = 0;
  TCanvas *c1 = new TCanvas("c1", "canvas");
  c1->cd();

  //if(flag==2) { // adding special case for AFEB timing
  ss.str("");
  ss << "mean_afeb_time_bin_vs_csc_ME";
  ss1.str("");
  ss1 << "Mean AFEB time bin vs CSC and ME";
  gStyle->SetOptStat(0);
  TH2F *hb = new TH2F(ss.str().c_str(), ss1.str().c_str(), 36, 1.0, 37.0, 18, 1.0, 19.0);
  hb->SetStats(kFALSE);
  hb->GetXaxis()->SetTitle("CSC #");
  hb->GetZaxis()->SetLabelSize(0.03);
  hb->SetOption("COLZ");

  hb->GetYaxis()->SetBinLabel(1, "ME- 4/2");
  hb->GetYaxis()->SetBinLabel(2, "ME- 4/1");
  hb->GetYaxis()->SetBinLabel(3, "ME- 3/2");
  hb->GetYaxis()->SetBinLabel(4, "ME- 3/1");
  hb->GetYaxis()->SetBinLabel(5, "ME- 2/2");
  hb->GetYaxis()->SetBinLabel(6, "ME- 2/1");
  hb->GetYaxis()->SetBinLabel(7, "ME- 1/3");
  hb->GetYaxis()->SetBinLabel(8, "ME- 1/2");
  hb->GetYaxis()->SetBinLabel(9, "ME- 1/1");
  hb->GetYaxis()->SetBinLabel(10, "ME+ 1/1");
  hb->GetYaxis()->SetBinLabel(11, "ME+ 1/2");
  hb->GetYaxis()->SetBinLabel(12, "ME+ 1/3");
  hb->GetYaxis()->SetBinLabel(13, "ME+ 2/1");
  hb->GetYaxis()->SetBinLabel(14, "ME+ 2/2");
  hb->GetYaxis()->SetBinLabel(15, "ME+ 3/1");
  hb->GetYaxis()->SetBinLabel(16, "ME+ 3/2");
  hb->GetYaxis()->SetBinLabel(17, "ME+ 4/1");
  hb->GetYaxis()->SetBinLabel(18, "ME+ 4/2");
  //}

  for (Int_t jesr = 0; jesr < 18; jesr++) {
    ss.str("");
    ss << result_histName.c_str() << esr[jesr];
    ss1.str("");
    ss1 << result_histTitle;
    TH2F *h = new TH2F(ss.str().c_str(), ss1.str().c_str(), 40, 0.0, 40.0, ny, ylow, yhigh);
    h->SetStats(kFALSE);
    h->GetXaxis()->SetTitle(xTitle[jesr].c_str());
    h->GetYaxis()->SetTitle(result_title_Y.c_str());
    h->GetZaxis()->SetLabelSize(0.03);
    h->SetOption("COLZ");

    ss.str("");
    ss << result_histNameEntries.c_str() << esr[jesr];
    ss1.str("");
    ss1 << result_histTitleEntries;
    TH2F *hentr = new TH2F(ss.str().c_str(), ss1.str().c_str(), 40, 0.0, 40.0, ny, ylow, yhigh);
    hentr->SetStats(kFALSE);
    hentr->GetXaxis()->SetTitle(xTitle[jesr].c_str());
    hentr->GetYaxis()->SetTitle(result_title_Y.c_str());
    hentr->GetZaxis()->SetLabelSize(0.03);
    hentr->SetOption("COLZ");

    //TH2F* ha;

    if (flag == 2) {  // adding special cases for AFEB timing
      ss.str("");
      ss << "normal_afeb_time_bin_vs_csc_ME_" << esr[jesr];
      ss1.str("");
      ss1 << "Normalized AFEB time bin, %";
      //TH2F *ha=new TH2F(ss.str().c_str(),ss1.str().c_str(),40,0.0,40.0,16,0.0,16.0);
      //ha->SetStats(kFALSE);
      //ha->GetXaxis()->SetTitle(xTitle[jesr].c_str());
      //ha->GetYaxis()->SetTitle("Time Bin");
      //ha->GetZaxis()->SetLabelSize(0.03);
      //ha->SetOption("COLZ");
    }

    for (Int_t csc = 1; csc < 37; csc++) {
      Int_t idchamber = esr[jesr] * 100 + csc;
      ss.str("");
      ss << folder.c_str() << input_histName.c_str() << idchamber;
      f_in->cd();
      TH2F *h2[1];
      h2[k] = (TH2F *)f_in->Get(ss.str().c_str());
      if (h2[k] != NULL) {
        // saving original, adding X,Y titles, color and "BOX" option
        h2[k]->GetXaxis()->SetTitle(input_title_X.c_str());
        h2[k]->GetYaxis()->SetTitle(input_title_Y.c_str());
        h2[k]->GetYaxis()->SetTitleOffset(1.2);
        h2[k]->SetFillColor(4);
        h2[k]->SetOption("BOX");
        gStyle->SetOptStat(1001111);

        // saving Y projection of the whole 2D hist for given chamber
        ss.str("");
        ss << input_histName.c_str() << idchamber << "_Y_all";
        TH1D *h1d = h2[k]->ProjectionY(ss.str().c_str(), 1, h2[k]->GetNbinsX(), "");
        h1d->GetYaxis()->SetTitle("Entries");
        h1d->GetYaxis()->SetTitleOffset(1.2);
        gStyle->SetOptStat(1001111);

        if (flag == 2 && h1d->GetEntries() > 0) {  // adding spec. case for afeb timing
          Float_t entr = h1d->GetEntries();
          for (Int_t m = 1; m < h1d->GetNbinsX(); m++) {
            Float_t w = h1d->GetBinContent(m);
            w = 100.0 * w / entr;
            //ha->SetBinContent(csc+1,m,w);
          }
          Float_t mean = h1d->GetMean();
          Int_t me;
          if (jesr < 9)
            me = 10 + jesr;
          if (jesr > 8)
            me = 18 - jesr;
          hb->SetBinContent(csc, me, mean);
        }
        delete h1d;

        // saving slices, finding MEAN in each slice, fill 2D hist
        for (Int_t j = 1; j <= h2[k]->GetNbinsX(); j++) {
          Int_t n = j;
          ss.str("");
          ss << input_histName.c_str() << idchamber << "_Y_" << n;
          TH1D *h1d = h2[k]->ProjectionY(ss.str().c_str(), j, j, "");
          if (h1d->GetEntries() > 0) {
            Float_t mean = h1d->GetMean();
            Float_t entr = h1d->GetEntries();
            entries[jesr] = entries[jesr] + 1;
            h->SetBinContent(csc + 1, j, mean);
            hentr->SetBinContent(csc + 1, j, entr);
            ss.str("");
            ss << slice_title_X << " " << n;
            h1d->GetXaxis()->SetTitle(ss.str().c_str());
            h1d->GetYaxis()->SetTitle("Entries");
            h1d->GetYaxis()->SetTitleOffset(1.2);
            gStyle->SetOptStat(1001111);
          }
          delete h1d;
        }
      }
    }
    if (entries[jesr] > 0) {
      h->SetStats(kFALSE);
      hentr->SetStats(kFALSE);
      c1->Update();

      // printing

      h->Draw();
      ss.str("");
      ss << result_histName.c_str() << esr[jesr] << ".png";
      c1->Print(ss.str().c_str(), "png");

      hentr->Draw();
      ss.str("");
      ss << result_histNameEntries.c_str() << esr[jesr] << ".png";
      c1->Print(ss.str().c_str(), "png");
    }
    delete h;
    delete hentr;
    //if(flag==2) delete ha;
  }
  if (flag == 2) {
    hb->Draw();
    ss.str("");
    ss << "mean_afeb_time_bin_vs_csc_ME"
       << ".png";
    c1->Print(ss.str().c_str(), "png");

    c1->Update();
    //delete hb;
  }
  delete hb;
  delete c1;
}

void drawColoredChamberLines(int station, int nchamber1[4][36]) {
  // thanks to Luca Sabbatini for this coe
  const int maxRingIdxOfRelevance = 2;
  const int maxChamberIdxOfRelevance = 35;
  const int totalNumberOfChambersOfRelevance = 540;

  Int_t thisRingIdx = 0;
  Int_t numChambersInRing = 0;

  float nchamber1Avg = 0.;
  float thisRingNchamber1Avg = 0;
  float thisRingNchamber1Max = 0;
  // Note that thisRingNchamber1Min is only the minumum *non-zero*
  // occupancy among all the cmabers for the given ring
  float thisRingNchamber1Min = 0;
  Int_t thisRingNchamber1Base = 0;

  // rFloatRGB, gFloatRGB, bFloatRGB, are from the RGB color space;
  // all three range between [0,1]
  // hFloatHLS, lFloatHLS, sFloatHLS, are from the HLS color space;
  // lFloatHLS and sFloatHLS range between [0,1], while hFloatHLS
  // ranges between [0,360]
  Float_t rFloatRGB = 0.;
  Float_t gFloatRGB = 0.;
  Float_t bFloatRGB = 0.;
  Float_t hFloatHLS = 0.;
  Float_t lFloatHLS = 0.5;
  Float_t sFloatHLS = 1.0;
  TColor tempColor;

  // compute average chamber occupancy over all CSCs
  for (int i = 0; i < maxRingIdxOfRelevance + 1; i++) {
    for (int j = 0; j < maxChamberIdxOfRelevance + 1; j++) {
      nchamber1Avg += nchamber1[i][j];
    }
  }
  nchamber1Avg = nchamber1Avg / totalNumberOfChambersOfRelevance;

  Float_t myFavoriteLineWidth = 2.0;
  gStyle->SetLineWidth(myFavoriteLineWidth);
  float pi = 3.14159;
  TVector3 x(0, 0, 1);
  int linecolor = 1;

  // emptyLineColor is the color of the outline of empty chambers.
  // if emptyLineColor is negative, no outline is drawn at all for
  // that chamber.
  Int_t emptyLineColor = kYellow;

  // note that these can be greater than 360 because I have included
  // a line down below that performs a modulus operation: %360, which
  // ensures that the resulting hue will indeed be with [0,360]
  //
  // Hue:    0   30      60      120     180     240 270 300     360
  //         |   |       |       |       |       |   |   |       |
  // Color: Red Orange Yellow   Green   Cyan  Blue Vilet Magenta Red
  Float_t lineHueMin = 240;
  Float_t lineHueMax = 360;

  TLine *line1;
  TLine *line2;
  TLine *line3;
  TLine *line4;

  if (station == 1) {
    // station 1, ring 1 (inner-most ring)
    thisRingIdx = 0;
    numChambersInRing = 36;
    TVector3 p1(101, 9.361, 0);
    TVector3 p2(101, -9.361, 0);
    TVector3 p3(260, -22.353, 0);
    TVector3 p4(260, 22.353, 0);

    // compute thisRingNchamber1Min, thisRingNchamber1Max,
    //  thisRingNchamber1Avg and thisRingNchamber1Base
    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      // set the line color
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS =
              int((nchamber1[0][i] - thisRingNchamber1Base) / (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                      (lineHueMax - lineHueMin) +
                  lineHueMin) %
              360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      // draw the chamber outline using the line color (so long
      // as the the chamber isn't empty *and* the user did not
      // set a negative line color for empty chambers)
      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
        line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
        line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
        line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
      }

      // Rotate coordinate by 1 chamber
      p1.Rotate(2 * pi / numChambersInRing, x);
      p2.Rotate(2 * pi / numChambersInRing, x);
      p3.Rotate(2 * pi / numChambersInRing, x);
      p4.Rotate(2 * pi / numChambersInRing, x);
    }

    // station 1, ring 2 (middle-ring)
    thisRingIdx = 1;
    numChambersInRing = 36;
    TVector3 q1(281.49, 25.5, 0);
    TVector3 q2(281.49, -25.5, 0);
    TVector3 q3(455.99, -41.87, 0);
    TVector3 q4(455.99, 41.87, 0);

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS = int((nchamber1[thisRingIdx][i] - thisRingNchamber1Base) /
                              (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                              (lineHueMax - lineHueMin) +
                          lineHueMin) %
                      360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
        line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
        line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
        line4 = new TLine(q4(0), q4(1), q1(0), q1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();

        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(q3(0) * q3(0) + q3(1) * q3(1));
        float theta = atan2(q3(1), q3(0));
        r = r + 25;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      q1.Rotate(2 * pi / numChambersInRing, x);
      q2.Rotate(2 * pi / numChambersInRing, x);
      q3.Rotate(2 * pi / numChambersInRing, x);
      q4.Rotate(2 * pi / numChambersInRing, x);
    }

    // station 1, ring 3 (outer-most ring)
    thisRingIdx = 2;
    numChambersInRing = 36;
    TVector3 r1(511.99, 31.7, 0);
    TVector3 r2(511.99, -31.7, 0);
    TVector3 r3(676.15, -46.05, 0);
    TVector3 r4(676.15, 46.05, 0);

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS = int((nchamber1[thisRingIdx][i] - thisRingNchamber1Base) /
                              (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                              (lineHueMax - lineHueMin) +
                          lineHueMin) %
                      360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(r1(0), r1(1), r2(0), r2(1));
        line2 = new TLine(r2(0), r2(1), r3(0), r3(1));
        line3 = new TLine(r3(0), r3(1), r4(0), r4(1));
        line4 = new TLine(r4(0), r4(1), r1(0), r1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(r3(0) * r3(0) + r3(1) * r3(1));
        float theta = atan2(r3(1), r3(0));
        r = r + 25;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      r1.Rotate(2 * pi / numChambersInRing, x);
      r2.Rotate(2 * pi / numChambersInRing, x);
      r3.Rotate(2 * pi / numChambersInRing, x);
      r4.Rotate(2 * pi / numChambersInRing, x);
    }
  }

  if (station == 2) {
    // station 2, ring 1 (inner ring)
    thisRingIdx = 0;
    numChambersInRing = 18;
    TVector3 p1(146.9, 27.0, 0);
    TVector3 p2(146.9, -27.0, 0);
    TVector3 p3(336.56, -62.855, 0);
    TVector3 p4(336.56, 62.855, 0);

    // must "pre-rotate" by one-fourth of a chamber
    p1.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p2.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p3.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p4.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS =
              int((nchamber1[0][i] - thisRingNchamber1Base) / (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                      (lineHueMax - lineHueMin) +
                  lineHueMin) %
              360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
        line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
        line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
        line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(p3(0) * p3(0) + p3(1) * p3(1));
        float theta = atan2(p3(1), p3(0));
        r = r + 10;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      p1.Rotate(2 * pi / numChambersInRing, x);
      p2.Rotate(2 * pi / numChambersInRing, x);
      p3.Rotate(2 * pi / numChambersInRing, x);
      p4.Rotate(2 * pi / numChambersInRing, x);
    }

    // station 2, ring 2 (outer ring)
    thisRingIdx = 1;
    numChambersInRing = 36;
    TVector3 q1(364.02, 33.23, 0);
    TVector3 q2(364.02, -33.23, 0);
    TVector3 q3(687.08, -63.575, 0);
    TVector3 q4(687.08, 63.575, 0);

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS = int((nchamber1[thisRingIdx][i] - thisRingNchamber1Base) /
                              (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                              (lineHueMax - lineHueMin) +
                          lineHueMin) %
                      360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
        line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
        line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
        line4 = new TLine(q4(0), q4(1), q1(0), q1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(q3(0) * q3(0) + q3(1) * q3(1));
        float theta = atan2(q3(1), q3(0));
        r = r + 25;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      q1.Rotate(2 * pi / numChambersInRing, x);
      q2.Rotate(2 * pi / numChambersInRing, x);
      q3.Rotate(2 * pi / numChambersInRing, x);
      q4.Rotate(2 * pi / numChambersInRing, x);
    }
  }

  if (station == 3) {
    // station 3, ring 1 (inner ring)
    thisRingIdx = 0;
    numChambersInRing = 18;
    TVector3 p1(166.89, 30.7, 0);
    TVector3 p2(166.89, -30.7, 0);
    TVector3 p3(336.59, -62.855, 0);
    TVector3 p4(336.59, 62.855, 0);

    // must "pre-rotate" by one-fourth of a chamber
    p1.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p2.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p3.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p4.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);

    TLine *line1;
    TLine *line2;
    TLine *line3;
    TLine *line4;

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS =
              int((nchamber1[0][i] - thisRingNchamber1Base) / (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                      (lineHueMax - lineHueMin) +
                  lineHueMin) %
              360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
        line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
        line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
        line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(p3(0) * p3(0) + p3(1) * p3(1));
        float theta = atan2(p3(1), p3(0));
        r = r + 10;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      p1.Rotate(2 * pi / numChambersInRing, x);
      p2.Rotate(2 * pi / numChambersInRing, x);
      p3.Rotate(2 * pi / numChambersInRing, x);
      p4.Rotate(2 * pi / numChambersInRing, x);
    }

    // station 3, ring 2 (outer ring)
    thisRingIdx = 1;
    numChambersInRing = 36;
    TVector3 q1(364.02, 33.23, 0);
    TVector3 q2(364.02, -33.23, 0);
    TVector3 q3(687.08, -63.575, 0);
    TVector3 q4(687.08, 63.575, 0);

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS = int((nchamber1[thisRingIdx][i] - thisRingNchamber1Base) /
                              (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                              (lineHueMax - lineHueMin) +
                          lineHueMin) %
                      360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
        line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
        line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
        line4 = new TLine(q4(0), q4(1), q1(0), q1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(q3(0) * q3(0) + q3(1) * q3(1));
        float theta = atan2(q3(1), q3(0));
        r = r + 25;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      q1.Rotate(2 * pi / numChambersInRing, x);
      q2.Rotate(2 * pi / numChambersInRing, x);
      q3.Rotate(2 * pi / numChambersInRing, x);
      q4.Rotate(2 * pi / numChambersInRing, x);
    }
  }

  if (station == 4) {
    // station 4, ring 1 (the only ring on station 4... so far)
    thisRingIdx = 0;
    numChambersInRing = 18;
    TVector3 p1(186.99, 34.505, 0);
    TVector3 p2(186.99, -34.505, 0);
    TVector3 p3(336.41, -62.825, 0);
    TVector3 p4(336.41, 62.825, 0);

    // must "pre-rotate" by one-fourth of a chamber
    p1.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p2.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p3.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);
    p4.Rotate((1 / 4.) * 2 * pi / numChambersInRing, x);

    TLine *line1;
    TLine *line2;
    TLine *line3;
    TLine *line4;

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS =
              int((nchamber1[0][i] - thisRingNchamber1Base) / (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                      (lineHueMax - lineHueMin) +
                  lineHueMin) %
              360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(p1(0), p1(1), p2(0), p2(1));
        line2 = new TLine(p2(0), p2(1), p3(0), p3(1));
        line3 = new TLine(p3(0), p3(1), p4(0), p4(1));
        line4 = new TLine(p4(0), p4(1), p1(0), p1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(p3(0) * p3(0) + p3(1) * p3(1));
        float theta = atan2(p3(1), p3(0));
        r = r + 10;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      p1.Rotate(2 * pi / numChambersInRing, x);
      p2.Rotate(2 * pi / numChambersInRing, x);
      p3.Rotate(2 * pi / numChambersInRing, x);
      p4.Rotate(2 * pi / numChambersInRing, x);
    }

    // station 4, ring 2 (outer ring) (This is just a copy of the ME3/2 code)
    thisRingIdx = 1;
    numChambersInRing = 36;
    TVector3 q1(364.02, 33.23, 0);
    TVector3 q2(364.02, -33.23, 0);
    TVector3 q3(687.08, -63.575, 0);
    TVector3 q4(687.08, 63.575, 0);

    thisRingNchamber1Avg = 0.;
    thisRingNchamber1Max = 0.;
    for (int j = 0; j < numChambersInRing; j++) {
      thisRingNchamber1Avg += nchamber1[thisRingIdx][j];
      if (thisRingNchamber1Max < nchamber1[thisRingIdx][j]) {
        thisRingNchamber1Max = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Min = thisRingNchamber1Max;
    for (int j = 0; j < numChambersInRing; j++) {
      if (thisRingNchamber1Min > nchamber1[thisRingIdx][j] &&  //
          nchamber1[thisRingIdx][j] != 0) {
        thisRingNchamber1Min = nchamber1[thisRingIdx][j];
      }
    }
    thisRingNchamber1Avg = thisRingNchamber1Avg / numChambersInRing;
    thisRingNchamber1Base = thisRingNchamber1Min;

    for (int i = 0; i < numChambersInRing; i++) {
      if (nchamber1[thisRingIdx][i] != 0) {
        if ((thisRingNchamber1Max - thisRingNchamber1Base) != 0)
          hFloatHLS = int((nchamber1[thisRingIdx][i] - thisRingNchamber1Base) /
                              (thisRingNchamber1Max - thisRingNchamber1Base) *  //
                              (lineHueMax - lineHueMin) +
                          lineHueMin) %
                      360;
        tempColor.HLS2RGB(hFloatHLS, lFloatHLS, sFloatHLS, rFloatRGB, gFloatRGB, bFloatRGB);
        linecolor = tempColor.GetColor(rFloatRGB, gFloatRGB, bFloatRGB);
      } else if (emptyLineColor >= 0) {
        linecolor = emptyLineColor;
      }

      if ((nchamber1[thisRingIdx][i] != 0) ||  //
          (nchamber1[thisRingIdx][i] == 0 && emptyLineColor >= 0)) {
        line1 = new TLine(q1(0), q1(1), q2(0), q2(1));
        line2 = new TLine(q2(0), q2(1), q3(0), q3(1));
        line3 = new TLine(q3(0), q3(1), q4(0), q4(1));
        line4 = new TLine(q4(0), q4(1), q1(0), q1(1));
        line1->SetLineColor(linecolor);
        line2->SetLineColor(linecolor);
        line3->SetLineColor(linecolor);
        line4->SetLineColor(linecolor);

        line1->Draw();
        line2->Draw();
        line3->Draw();
        line4->Draw();
        int chamberNum = i + 1;
        char chStr[3];
        sprintf(chStr, "%d", chamberNum);
        TLatex l;
        l.SetTextAlign(22);
        l.SetTextSize(0.02);
        float r = sqrt(q3(0) * q3(0) + q3(1) * q3(1));
        float theta = atan2(q3(1), q3(0));
        r = r + 25;
        theta = pi / numChambersInRing + theta;
        float xpos = r * cos(theta);
        float ypos = r * sin(theta);
        l.DrawLatex(xpos, ypos, chStr);
      }
      q1.Rotate(2 * pi / numChambersInRing, x);
      q2.Rotate(2 * pi / numChambersInRing, x);
      q3.Rotate(2 * pi / numChambersInRing, x);
      q4.Rotate(2 * pi / numChambersInRing, x);
    }
  }
}