OccupancyPlotMacros.cc

Go to the documentation of this file.

#include "TText.h"
#include "TLatex.h"
#include "TLine.h"
#include "TGaxis.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TH1F.h"
#include "TProfile.h"
#include "TH1D.h"
#include "TList.h"
#include "TBox.h"
#include "TFrame.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TROOT.h"
#include <cstring>
#include <iostream>
#include <cmath>
#include <algorithm>
#include "TROOT.h"
#include "DPGAnalysis/SiStripTools/interface/CommonAnalyzer.h"
#include "OccupancyPlotMacros.h"
//#include <vector>
 
float linear(float x) { return x; }
float logarithm(float x) { return log(x); }
 
std::pair<float, float> presentbin(int i) {
  float dz = -1;
  float dr = -1;
 
  if (i > 100 && i < 200) {
    dz = 3.33;
    dr = 0.4;
  }  // BPIX
 
  if (i > 200 && i < 1000 && (i % 10 == 1 || i % 10 == 7)) {
    dz = 0.8;
    dr = 0.4;
  }  // FPIX
  if (i > 200 && i < 1000 && !(i % 10 == 1 || i % 10 == 7)) {
    dz = 0.8;
    dr = 0.8;
  }
 
  if (i > 1000 && i < 2000) {
    dz = 5.948;
    dr = 0.4;
  }  // TIB
 
  if (i > 3000 && i < 4000) {
    dz = 9.440;
    dr = 0.4;
  }  // TOB
 
  if (i > 2000 && i < 3000 && (i % 1000) / 100 == 1) {
    dz = 0.8;
    dr = 5.647;
  }  // TID
  if (i > 2000 && i < 3000 && (i % 1000) / 100 == 2) {
    dz = 0.8;
    dr = 4.512;
  }
  if (i > 2000 && i < 3000 && (i % 1000) / 100 == 3) {
    dz = 0.8;
    dr = 5.637;
  }
 
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 1) {
    dz = 0.8;
    dr = 4.362;
  }  // TEC
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 2) {
    dz = 0.8;
    dr = 4.512;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 3) {
    dz = 0.8;
    dr = 5.637;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 4) {
    dz = 0.8;
    dr = 5.862;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 5) {
    dz = 0.8;
    dr = 7.501;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 6) {
    dz = 0.8;
    dr = 9.336;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 7) {
    dz = 0.8;
    dr = 10.373;
  }
 
  std::pair<float, float> res(dz, dr);
  return res;
}
 
std::pair<float, float> phase2bin(int i) {
  float dz = -1;
  float dr = -1;
 
  if (i > 1000 && i < 1040) {
    dz = 3.33;
    dr = 0.4;
  }
  if (i > 1040 && i < 1130) {
    dr = 3.33;
    dz = 0.4;
  }
 
  if (i > 2000 && i < 2550) {
    dz = 2.5;
    dr = 0.1;
  }
  if (i > 2550 && i < 3000) {
    dz = 5.0;
    dr = 0.1;
  }
 
  if (i > 3000 && i < 5000) {
    dz = 0.2;
    dr = 5.0;
  }
 
  if (i > 3100 && i < 3119) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 3140 && i < 3159) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 3180 && i < 3199) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 3220 && i < 3239) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 3260 && i < 3279) {
    dz = 0.2;
    dr = 2.5;
  }
 
  if (i > 4100 && i < 4119) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 4140 && i < 4159) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 4180 && i < 4199) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 4220 && i < 4239) {
    dz = 0.2;
    dr = 2.5;
  }
  if (i > 4260 && i < 4279) {
    dz = 0.2;
    dr = 2.5;
  }
 
  std::pair<float, float> res(dz, dr);
  return res;
}
 
std::pair<float, float> phase1bin(int i) {
  float dz = -1;
  float dr = -1;
 
  if (i > 1000 && i < 1040) {
    dz = 3.33;
    dr = 0.4;
  }  // BPIX
  if (i > 1040 && i < 1080) {
    dr = 3.33;
    dz = 0.4;
  }  // FPIX
 
  if (i > 1100 && i < 2000) {
    dz = 5.948;
    dr = 0.4;
  }  // TIB
 
  if (i > 3000 && i < 4000) {
    dz = 9.440;
    dr = 0.4;
  }  // TOB
 
  if (i > 2000 && i < 3000 && (i % 1000) / 100 == 1) {
    dz = 0.8;
    dr = 5.647;
  }  // TID
  if (i > 2000 && i < 3000 && (i % 1000) / 100 == 2) {
    dz = 0.8;
    dr = 4.512;
  }
  if (i > 2000 && i < 3000 && (i % 1000) / 100 == 3) {
    dz = 0.8;
    dr = 5.637;
  }
 
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 1) {
    dz = 0.8;
    dr = 4.362;
  }  // TEC
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 2) {
    dz = 0.8;
    dr = 4.512;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 3) {
    dz = 0.8;
    dr = 5.637;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 4) {
    dz = 0.8;
    dr = 5.862;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 5) {
    dz = 0.8;
    dr = 7.501;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 6) {
    dz = 0.8;
    dr = 9.336;
  }
  if (i > 4000 && i < 6000 && (i % 1000) / 100 == 7) {
    dz = 0.8;
    dr = 10.373;
  }
 
  std::pair<float, float> res(dz, dr);
  return res;
}
 
void printFrame(TCanvas* c, TH1D* h, const char* label, const int frame, const int min, const int max, const bool same) {
  c->cd(frame);
  h->SetAxisRange(min, max);
  if (same) {
    h->DrawCopy("same");
  } else {
    h->DrawCopy();
    TText* t = new TText((max + min) / 2, h->GetMaximum(), label);
    t->SetTextAlign(22);
    t->DrawClone();
  }
}
 
void PlotOccupancyMapGeneric(TFile* ff,
                             const char* module,
                             const float min,
                             const float max,
                             const float mmin,
                             const float mmax,
                             const int color,
                             std::pair<float, float> (*size)(int),
                             std::vector<SubDetParams>& vsub) {
  gROOT->SetStyle("Plain");
 
  if (ff->cd(module)) {
    TProfile* aveoccu = (TProfile*)gDirectory->Get("aveoccu");
    TProfile* avemult = (TProfile*)gDirectory->Get("avemult");
    TH1F* nchannels = (TH1F*)gDirectory->Get("nchannels_real");
 
    TProfile* averadius = (TProfile*)gDirectory->Get("averadius");
    TProfile* avez = (TProfile*)gDirectory->Get("avez");
 
    std::cout << "pointers " << aveoccu << " " << avemult << " " << nchannels << " " << averadius << " " << avez
              << std::endl;
 
    if (aveoccu && avemult && nchannels && averadius && avez) {
      nchannels->Sumw2();
      for (int i = 1; i < nchannels->GetNbinsX() + 1; ++i) {
        nchannels->SetBinError(i, 0.);
      }
 
      TH1D* haveoccu = aveoccu->ProjectionX("haveoccu");
      haveoccu->SetDirectory(nullptr);
      haveoccu->Divide(nchannels);
 
      TH1D* havemult = avemult->ProjectionX("havemult");
      havemult->SetDirectory(nullptr);
      havemult->Divide(nchannels);
 
      TH1D* havewidth = (TH1D*)haveoccu->Clone("havewidth");
      havewidth->SetDirectory(nullptr);
      havewidth->SetTitle("Average Cluster Size");
      havewidth->Divide(havemult);
 
      new TCanvas("occupancy", "occupancy", 1200, 500);
      gPad->SetLogy(1);
      haveoccu->SetStats(false);
      haveoccu->SetLineColor(kRed);
      haveoccu->SetMarkerColor(kRed);
      haveoccu->DrawCopy();
 
      new TCanvas("multiplicity", "multiplicity", 1200, 500);
      gPad->SetLogy(1);
      havemult->SetStats(false);
      havemult->SetLineColor(kRed);
      havemult->SetMarkerColor(kRed);
      havemult->DrawCopy();
 
      new TCanvas("width", "width", 1200, 500);
      havewidth->SetStats(false);
      havewidth->SetLineColor(kRed);
      havewidth->SetMarkerColor(kRed);
      havewidth->DrawCopy();
 
      bool same = false;
      TCanvas* o2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("occupancy2");
      if (o2) {
        same = true;
        haveoccu->SetLineColor(kBlue);
        haveoccu->SetMarkerColor(kBlue);
      } else {
        o2 = new TCanvas("occupancy2", "occupancy2", 1200, 800);
        o2->Divide(3, 2);
      }
      for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
        printFrame(o2, haveoccu, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
      }
 
      same = false;
      TCanvas* m2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("multiplicity2");
      if (m2) {
        same = true;
        havemult->SetLineColor(kBlue);
        havemult->SetMarkerColor(kBlue);
      } else {
        m2 = new TCanvas("multiplicity2", "multiplicity2", 1200, 800);
        m2->Divide(3, 2);
      }
      for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
        printFrame(m2, havemult, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
      }
 
      same = false;
      TCanvas* w2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("width2");
      if (w2) {
        same = true;
        havewidth->SetLineColor(kBlue);
        havewidth->SetMarkerColor(kBlue);
      } else {
        w2 = new TCanvas("width2", "width2", 1200, 800);
        w2->Divide(3, 2);
      }
      for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
        printFrame(w2, havewidth, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
      }
 
      float (*scale)(float);
      scale = &logarithm;
 
      drawMap("multmap", havemult, averadius, avez, mmin, mmax, size, scale, color);
      drawMap("occumap", haveoccu, averadius, avez, min, max, size, scale, color, "channel occupancy");
    }
  }
}
 
float combinedOccupancy(TFile* ff, const char* module, const int lowerbin, const int upperbin) {
  float cumoccu = -2.;
  double cumerr = -2;
 
  if (ff->cd(module)) {
    TProfile* aveoccu = (TProfile*)gDirectory->Get("aveoccu");
    //    TProfile* avemult= (TProfile*)gDirectory->Get("avemult");
    TH1F* nchannels = (TH1F*)gDirectory->Get("nchannels_real");
 
    float sumoccu = 0.;
    float sumnchannels = 0;
    double sumerrsq = 0;
 
    for (int i = lowerbin; i < upperbin + 1; ++i) {
      std::cout << "processing bin " << i << " " << aveoccu->GetBinContent(i) << "+/-" << aveoccu->GetBinError(i)
                << std::endl;
      sumoccu += aveoccu->GetBinContent(i);
      sumnchannels += nchannels->GetBinContent(i);
      sumerrsq += aveoccu->GetBinError(i) * aveoccu->GetBinError(i);
    }
    cumoccu = sumnchannels != 0 ? sumoccu / sumnchannels : -1;
    cumerr = sumnchannels != 0 ? sqrt(sumerrsq) / sumnchannels : -1;
    std::cout << "Cumulative occupancy: " << sumoccu << " " << sumnchannels << " " << cumoccu << "+/-" << cumerr;
  }
 
  return cumoccu;
}
 
void PlotOccupancyMap(TFile* ff,
                      const char* module,
                      const float min,
                      const float max,
                      const float mmin,
                      const float mmax,
                      const int color) {
  std::pair<float, float> (*size)(int);
  size = &presentbin;
 
  std::vector<SubDetParams> vsub;
  SubDetParams ppix = {"BPIX+FPIX", 100, 270};
  vsub.push_back(ppix);
  SubDetParams ptib = {"TIB", 1050, 1450};
  vsub.push_back(ptib);
  SubDetParams ptid = {"TID", 2070, 2400};
  vsub.push_back(ptid);
  SubDetParams ptob = {"TOB", 3000, 3700};
  vsub.push_back(ptob);
  SubDetParams ptecm = {"TEC-", 4000, 4850};
  vsub.push_back(ptecm);
  SubDetParams ptecp = {"TEC+", 5000, 5850};
  vsub.push_back(ptecp);
 
  PlotOccupancyMapGeneric(ff, module, min, max, mmin, mmax, color, size, vsub);
}
 
void PlotOccupancyMapPhase1(TFile* ff,
                            const char* module,
                            const float min,
                            const float max,
                            const float mmin,
                            const float mmax,
                            const int color) {
  std::pair<float, float> (*size)(int);
  size = &phase1bin;
 
  std::vector<SubDetParams> vsub;
  SubDetParams ppix = {"BPIX+FPIX", 1000, 1080};
  vsub.push_back(ppix);
  SubDetParams ptib = {"TIB", 1090, 1450};
  vsub.push_back(ptib);
  SubDetParams ptid = {"TID", 2070, 2400};
  vsub.push_back(ptid);
  SubDetParams ptob = {"TOB", 3000, 3700};
  vsub.push_back(ptob);
  SubDetParams ptecm = {"TEC-", 4000, 4850};
  vsub.push_back(ptecm);
  SubDetParams ptecp = {"TEC+", 5000, 5850};
  vsub.push_back(ptecp);
 
  PlotOccupancyMapGeneric(ff, module, min, max, mmin, mmax, color, size, vsub);
}
 
void PlotOccupancyMapPhase2(TFile* ff,
                            const char* module,
                            const float min,
                            const float max,
                            const float mmin,
                            const float mmax,
                            const int color) {
  std::pair<float, float> (*size)(int);
  size = &phase2bin;
 
  std::vector<SubDetParams> vsub;
  SubDetParams ppix = {"BPIX+FPIX", 1000, 1090};
  vsub.push_back(ppix);
  SubDetParams ptob = {"TOB", 2000, 2900};
  vsub.push_back(ptob);
  SubDetParams ptecm = {"TEC-", 3100, 3300};
  vsub.push_back(ptecm);
  SubDetParams ptecp = {"TEC+", 4100, 4300};
  vsub.push_back(ptecp);
 
  PlotOccupancyMapGeneric(ff, module, min, max, mmin, mmax, color, size, vsub);
}
 
void PlotOnTrackOccupancy(
    TFile* ff, const char* module, const char* ontrkmod, const float mmin, const float mmax, const int color) {
  std::pair<float, float> (*size)(int);
  size = &presentbin;
 
  std::vector<SubDetParams> vsub;
  SubDetParams ppix = {"BPIX+FPIX", 100, 270};
  vsub.push_back(ppix);
  SubDetParams ptib = {"TIB", 1050, 1450};
  vsub.push_back(ptib);
  SubDetParams ptid = {"TID", 2070, 2400};
  vsub.push_back(ptid);
  SubDetParams ptob = {"TOB", 3000, 3700};
  vsub.push_back(ptob);
  SubDetParams ptecm = {"TEC-", 4000, 4850};
  vsub.push_back(ptecm);
  SubDetParams ptecp = {"TEC+", 5000, 5850};
  vsub.push_back(ptecp);
 
  PlotOnTrackOccupancyGeneric(ff, module, ontrkmod, mmin, mmax, color, size, vsub);
}
 
void PlotOnTrackOccupancyPhase1(
    TFile* ff, const char* module, const char* ontrkmod, const float mmin, const float mmax, const int color) {
  std::pair<float, float> (*size)(int);
  size = &phase1bin;
 
  std::vector<SubDetParams> vsub;
  SubDetParams ppix = {"BPIX+FPIX", 1000, 1080};
  vsub.push_back(ppix);
  SubDetParams ptib = {"TIB", 1090, 1450};
  vsub.push_back(ptib);
  SubDetParams ptid = {"TID", 2070, 2400};
  vsub.push_back(ptid);
  SubDetParams ptob = {"TOB", 3000, 3700};
  vsub.push_back(ptob);
  SubDetParams ptecm = {"TEC-", 4000, 4850};
  vsub.push_back(ptecm);
  SubDetParams ptecp = {"TEC+", 5000, 5850};
  vsub.push_back(ptecp);
 
  PlotOnTrackOccupancyGeneric(ff, module, ontrkmod, mmin, mmax, color, size, vsub);
}
 
void PlotOnTrackOccupancyPhase2(
    TFile* ff, const char* module, const char* ontrkmod, const float mmin, const float mmax, const int color) {
  std::pair<float, float> (*size)(int);
  size = &phase2bin;
 
  std::vector<SubDetParams> vsub;
  SubDetParams ppix = {"BPIX+FPIX", 1000, 1090};
  vsub.push_back(ppix);
  SubDetParams ptob = {"TOB", 2000, 2900};
  vsub.push_back(ptob);
  SubDetParams ptecm = {"TEC-", 3100, 3300};
  vsub.push_back(ptecm);
  SubDetParams ptecp = {"TEC+", 4100, 4300};
  vsub.push_back(ptecp);
 
  PlotOnTrackOccupancyGeneric(ff, module, ontrkmod, mmin, mmax, color, size, vsub);
}
 
void PlotOnTrackOccupancyGeneric(TFile* ff,
                                 const char* module,
                                 const char* ontrkmod,
                                 const float mmin,
                                 const float mmax,
                                 const int color,
                                 std::pair<float, float> (*size)(int),
                                 const std::vector<SubDetParams>& vsub) {
  gROOT->SetStyle("Plain");
 
  TProfile* avemult = nullptr;
  TProfile* aveontrkmult = nullptr;
  TProfile* averadius = nullptr;
  TProfile* avez = nullptr;
 
  if (ff->cd(module)) {
    avemult = (TProfile*)gDirectory->Get("avemult");
    averadius = (TProfile*)gDirectory->Get("averadius");
    avez = (TProfile*)gDirectory->Get("avez");
  }
  if (ff->cd(ontrkmod))
    aveontrkmult = (TProfile*)gDirectory->Get("avemult");
 
  std::cout << "pointers " << avemult << " " << aveontrkmult << " " << averadius << " " << avez << std::endl;
 
  if (averadius && avez && avemult && aveontrkmult) {
    TH1D* havemult = avemult->ProjectionX("havemult");
    TH1D* haveontrkmult = aveontrkmult->ProjectionX("haveontrkmult");
    havemult->SetDirectory(nullptr);
    haveontrkmult->SetDirectory(nullptr);
    haveontrkmult->Divide(havemult);
 
    new TCanvas("ontrkmult", "ontrkmult", 1200, 500);
    gPad->SetLogy(1);
    haveontrkmult->SetStats(false);
    haveontrkmult->SetLineColor(kRed);
    haveontrkmult->SetMarkerColor(kRed);
    haveontrkmult->SetMarkerSize(.5);
    haveontrkmult->SetMarkerStyle(20);
    haveontrkmult->DrawCopy();
 
    bool same = false;
    TCanvas* o2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("ontrkmult2");
    if (o2) {
      same = true;
      haveontrkmult->SetLineColor(kBlue);
      haveontrkmult->SetMarkerColor(kBlue);
    } else {
      o2 = new TCanvas("ontrkmult2", "ontrkmult2", 1200, 800);
      o2->Divide(3, 2);
    }
    for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
      printFrame(o2, haveontrkmult, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
    }
 
    float (*scale)(float);
    scale = &linear;
 
    drawMap("ontrkmultmap", haveontrkmult, averadius, avez, mmin, mmax, size, scale, color);
  }
}
 
TCanvas* drawMap(const char* cname,
                 const TH1* hval,
                 const TProfile* averadius,
                 const TProfile* avez,
                 const float mmin,
                 const float mmax,
                 std::pair<float, float> (*size)(int),
                 float (*scale)(float),
                 const int color,
                 const char* ptitle) {
  if (color == 1) {
    // A not-so-great color version
    const Int_t NRGBs = 5;
    const Int_t NCont = 255;
    Double_t stops[NRGBs] = {0.00, 0.25, 0.50, 0.75, 1.00};
    Double_t red[NRGBs] = {0.00, 0.00, 0.40, 1.00, 1.00};
    Double_t green[NRGBs] = {0.00, 0.40, 0.70, 0.60, 1.00};
    Double_t blue[NRGBs] = {0.30, 0.60, 0.00, 0.00, 0.20};
    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
    gStyle->SetNumberContours(NCont);
  } else if (color == 2) {
    // Gray scale
    const Int_t NRGBs = 3;
    const Int_t NCont = 255;
    Double_t stops[NRGBs] = {0.00, 0.50, 1.00};
    Double_t red[NRGBs] = {0.90, 0.50, 0.00};
    Double_t green[NRGBs] = {0.90, 0.50, 0.00};
    Double_t blue[NRGBs] = {0.90, 0.50, 0.00};
    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
    gStyle->SetNumberContours(NCont);
  } else if (color == 3) {
    // used by Kevin in the TRK-11-001 paper
    const Int_t NRGBs = 7;
    const Int_t NCont = 255;
    Double_t stops[NRGBs] = {0.00, 0.15, 0.30, 0.45, 0.65, 0.85, 1.00};
    Double_t red[NRGBs] = {0.60, 0.30, 0.00, 0.00, 0.60, 0.40, 0.00};
    Double_t green[NRGBs] = {1.00, 0.90, 0.80, 0.75, 0.20, 0.00, 0.00};
    Double_t blue[NRGBs] = {1.00, 1.00, 1.00, 0.30, 0.00, 0.00, 0.00};
    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
    gStyle->SetNumberContours(NCont);
  }
 
  int ncol = gStyle->GetNumberOfColors();
  std::cout << "Number of colors " << ncol << std::endl;
  // Loop on bins and creation of boxes
 
  TList modulesmult;
 
  for (int i = 1; i < hval->GetNbinsX(); ++i) {
    if ((averadius->GetBinEntries(i) * avez->GetBinEntries(i)) != 0) {
      double dz = -1.;
      double dr = -1.;
      // determine module size
 
      dz = (*size)(i).first;
      dr = (*size)(i).second;
 
      if (dz < 0 && dr < 0)
        continue;
 
      {
        TBox* modmult = new TBox(avez->GetBinContent(i) - dz,
                                 averadius->GetBinContent(i) - dr,
                                 avez->GetBinContent(i) + dz,
                                 averadius->GetBinContent(i) + dr);
        modmult->SetFillStyle(1001);
        if (color < 0) {
          modmult->SetFillColor(kBlack);
        } else {
          int icol = int(ncol * (scale(hval->GetBinContent(i)) - scale(mmin)) / (scale(mmax) - scale(mmin)));
          if (icol < 0)
            icol = 0;
          if (icol > (ncol - 1))
            icol = (ncol - 1);
          std::cout << i << " " << icol << " " << hval->GetBinContent(i) << std::endl;
          modmult->SetFillColor(gStyle->GetColorPalette(icol));
        }
        modulesmult.Add(modmult);
      }
    }
  }
  // eta boundaries lines
  double etavalext[] = {4., 3.5, 3., 2.8, 2.6, 2.4, 2.2, 2.0, 1.8, 1.6};
  double etavalint[] = {-1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4};
  TList etalines;
  TList etalabels;
  TList paperlabels;
  for (int i = 0; i < 10; ++i) {
    //  double eta = 3.0-i*0.2;
    double eta = etavalext[i];
    TLine* lin = new TLine(295, 2 * 295 / (exp(eta) - exp(-eta)), 305, 2 * 305 / (exp(eta) - exp(-eta)));
    etalines.Add(lin);
    char lab[100];
    sprintf(lab, "%3.1f", eta);
    TText* label = new TText(285, 2 * 285 / (exp(eta) - exp(-eta)), lab);
    label->SetTextSize(.03);
    label->SetTextAlign(22);
    etalabels.Add(label);
  }
  for (int i = 0; i < 10; ++i) {
    //  double eta = -3.0+i*0.2;
    double eta = -1 * etavalext[i];
    TLine* lin = new TLine(-295, -2 * 295 / (exp(eta) - exp(-eta)), -305, -2 * 305 / (exp(eta) - exp(-eta)));
    etalines.Add(lin);
    char lab[100];
    sprintf(lab, "%3.1f", eta);
    TText* label = new TText(-285, -2 * 285 / (exp(eta) - exp(-eta)), lab);
    label->SetTextSize(.03);
    label->SetTextAlign(22);
    etalabels.Add(label);
  }
  for (int i = 0; i < 15; ++i) {
    //  double eta = -1.4+i*0.2;
    double eta = etavalint[i];
    TLine* lin = new TLine(130. * (exp(eta) - exp(-eta)) / 2., 130, 138. * (exp(eta) - exp(-eta)) / 2., 138);
    etalines.Add(lin);
    char lab[100];
    sprintf(lab, "%3.1f", eta);
    TText* label = new TText(125. * (exp(eta) - exp(-eta)) / 2., 125, lab);
    label->SetTextSize(.03);
    label->SetTextAlign(22);
    etalabels.Add(label);
  }
  TLatex* etalab = new TLatex(0, 115, "#eta");
  etalab->SetTextSize(.03);
  etalab->SetTextAlign(22);
  etalabels.Add(etalab);
 
  // CMS label
  TLatex* cmslab = new TLatex(0.15, 0.965, "CMS");
  cmslab->SetNDC();
  cmslab->SetTextSize(0.04);
  cmslab->SetTextAlign(31);
  paperlabels.Add(cmslab);
  TLatex* enelab = new TLatex(0.92, 0.965, "#sqrt{s} = 7 TeV");
  enelab->SetNDC();
  enelab->SetTextSize(0.04);
  enelab->SetTextAlign(31);
  paperlabels.Add(enelab);
  /*
    TLatex *lumilab = new TLatex(0.6,0.965,Form("L = %.1f  fb^{-1}",19.7));
    lumilab->SetNDC();
    lumilab->SetTextSize(0.04);
    lumilab->SetTextAlign(31);
    paperlabels.Add(lumilab);
  */
 
  TGaxis* raxis = new TGaxis(-310, 0, -310, 140, 0, 140, 10, "S");
  TGaxis* zaxis = new TGaxis(-310, 0, 310, 0, -310, 310, 10, "S");
  raxis->SetTickSize(.01);
  zaxis->SetTickSize(.01);
  raxis->SetTitle("r (cm)");
  zaxis->SetTitle("z (cm)");
 
  TList mpalette;
 
  for (int i = 0; i < ncol; ++i) {
    TBox* box = new TBox(315, 0 + 140. / ncol * i, 330, 0 + 140. / ncol * (i + 1));
    box->SetFillStyle(1001);
    box->SetFillColor(gStyle->GetColorPalette(i));
    mpalette.Add(box);
  }
 
  TGaxis* mpaxis = nullptr;
  if (scale(1) != 1) {
    mpaxis = new TGaxis(330, 0, 330, 140, mmin, mmax, 510, "SLG+");
  } else {
    mpaxis = new TGaxis(330, 0, 330, 140, mmin, mmax, 510, "SL+");
  }
  mpaxis->SetTickSize(.02);
  mpaxis->SetLabelOffset(mpaxis->GetLabelOffset() * 0.5);
  mpaxis->SetTitle(ptitle);
  mpalette.Add(mpaxis);
 
  TCanvas* cc2 = new TCanvas(cname, cname, 1000, 500);
  cc2->Range(-370., -20., 390., 150.);
  TFrame* fr2 = new TFrame(-310, 0, 310, 140);
  fr2->UseCurrentStyle();
  fr2->Draw();
  raxis->Draw();
  zaxis->Draw();
  std::cout << modulesmult.GetSize() << std::endl;
  etalines.Draw();
  etalabels.Draw();
  if (color >= 0)
    mpalette.Draw();
  modulesmult.Draw();
 
  return cc2;
}
 
void PlotDebugFPIX_XYMap(TFile* ff, const char* module, const unsigned int ioffset, const char* name) {
  gROOT->SetStyle("Plain");
 
  TCanvas* cc = new TCanvas(name, name, 750, 750);
  cc->Range(-25, -25, 25, 25);
  TFrame* fr1 = new TFrame(-20, -20, 20, 20);
  fr1->UseCurrentStyle();
  fr1->Draw();
  ff->cd(module);
  gDirectory->ls();
  TProfile* avex = (TProfile*)gDirectory->Get("avex");
  TProfile* avey = (TProfile*)gDirectory->Get("avey");
  TProfile* avez = (TProfile*)gDirectory->Get("avez");
 
  if (avex && avey && avez) {
    TText* tittext = new TText(0, 0, name);
    tittext->SetTextSize(.04);
    tittext->SetTextAlign(22);
    tittext->Draw();
    for (unsigned int mod = ioffset + 1; mod < ioffset + 57; ++mod) {
      double x = avex->GetBinContent(mod);
      double y = avey->GetBinContent(mod);
      //      TBox* modbox = new TBox(x-1,y-1,x+1,y+1);
      char modstring[30];
      sprintf(modstring, "%d", mod % 100);
      TText* modtext = new TText(x, y, modstring);
      modtext->SetTextAngle(atan(y / x) * 180 / 3.14159);
      modtext->SetTextSize(.02);
      modtext->SetTextAlign(22);
      modtext->SetTextColor(kRed);
      std::cout << mod << " " << x << " " << y << std::endl;
      //      modbox->Draw();
      modtext->Draw();
    }
    for (unsigned int mod = ioffset + 101; mod < ioffset + 157; ++mod) {
      double x = avex->GetBinContent(mod);
      double y = avey->GetBinContent(mod);
      //      TBox* modbox = new TBox(x-1,y-1,x+1,y+1);
      char modstring[30];
      sprintf(modstring, "%d", mod % 100);
      TText* modtext = new TText(x, y, modstring);
      modtext->SetTextAngle(atan(y / x) * 180 / 3.14159);
      modtext->SetTextSize(.02);
      modtext->SetTextAlign(22);
      modtext->SetTextColor(kBlue);
      std::cout << mod << " " << x << " " << y << " " << atan(y / x) << std::endl;
      //      modbox->Draw();
      modtext->Draw();
    }
  }
}
 
void PlotTrackerXsect(TFile* ff, const char* module) {
  gROOT->SetStyle("Plain");
 
  if (ff->cd(module)) {
    TProfile* averadius = (TProfile*)gDirectory->Get("averadius");
    TProfile* avez = (TProfile*)gDirectory->Get("avez");
 
    std::cout << "pointers " << averadius << " " << avez << std::endl;
 
    if (averadius && avez) {
      std::pair<float, float> (*size)(int);
      size = &presentbin;
      float (*scale)(float);
      scale = &linear;
 
      drawMap("trackermap", averadius, averadius, avez, 0, 0, size, scale, -1);
    }
  }
}
 
TH1D* TrendPlotSingleBin(TFile* ff, const char* module, const char* hname, const int bin) {
  CommonAnalyzer caoccu(ff, "", module);
 
  TH1D* occutrend = new TH1D("occutrend", "Average number of clusters vs run", 10, 0., 10.);
  occutrend->SetCanExtend(TH1::kXaxis);
  occutrend->Sumw2();
 
  std::vector<unsigned int> runs = caoccu.getRunList();
  std::sort(runs.begin(), runs.end());
 
  {
    for (unsigned int i = 0; i < runs.size(); ++i) {
      char runlabel[100];
      sprintf(runlabel, "%d", runs[i]);
      char runpath[100];
      sprintf(runpath, "run_%d", runs[i]);
      caoccu.setPath(runpath);
 
      TProfile* occu = nullptr;
      if (occu == nullptr)
        occu = (TProfile*)caoccu.getObject(hname);
      if (occu) {
        const int ibin = occu->FindBin(bin);
        std::cout << runlabel << " "
                  << " " << ibin << " " << occu->GetBinContent(ibin) << " " << occu->GetBinError(ibin) << std::endl;
        const int jbin = occutrend->Fill(runlabel, occu->GetBinContent(ibin));
        occutrend->SetBinError(jbin, occu->GetBinError(ibin));
      }
    }
  }
  return occutrend;
}