SiStripTools/bin/OccupancyPlotMacros.cc

0001 #include "TText.h"
0002 #include "TLatex.h"
0003 #include "TLine.h"
0004 #include "TGaxis.h"
0005 #include "TFile.h"
0006 #include "TDirectory.h"
0007 #include "TH1F.h"
0008 #include "TProfile.h"
0009 #include "TH1D.h"
0010 #include "TList.h"
0011 #include "TBox.h"
0012 #include "TFrame.h"
0013 #include "TStyle.h"
0014 #include "TCanvas.h"
0015 #include "TColor.h"
0016 #include "TROOT.h"
0017 #include <cstring>
0018 #include <iostream>
0019 #include <cmath>
0020 #include <algorithm>
0021 #include "TROOT.h"
0022 #include "DPGAnalysis/SiStripTools/interface/CommonAnalyzer.h"
0023 #include "OccupancyPlotMacros.h"
0024 //#include <vector>
0025
0026 float linear(float x) { return x; }
0027 float logarithm(float x) { return log(x); }
0028
0029 std::pair<float, float> presentbin(int i) {
0030   float dz = -1;
0031   float dr = -1;
0032
0033   if (i > 100 && i < 200) {
0034     dz = 3.33;
0035     dr = 0.4;
0036   }  // BPIX
0037
0038   if (i > 200 && i < 1000 && (i % 10 == 1 || i % 10 == 7)) {
0039     dz = 0.8;
0040     dr = 0.4;
0041   }  // FPIX
0042   if (i > 200 && i < 1000 && !(i % 10 == 1 || i % 10 == 7)) {
0043     dz = 0.8;
0044     dr = 0.8;
0045   }
0046
0047   if (i > 1000 && i < 2000) {
0048     dz = 5.948;
0049     dr = 0.4;
0050   }  // TIB
0051
0052   if (i > 3000 && i < 4000) {
0053     dz = 9.440;
0054     dr = 0.4;
0055   }  // TOB
0056
0057   if (i > 2000 && i < 3000 && (i % 1000) / 100 == 1) {
0058     dz = 0.8;
0059     dr = 5.647;
0060   }  // TID
0061   if (i > 2000 && i < 3000 && (i % 1000) / 100 == 2) {
0062     dz = 0.8;
0063     dr = 4.512;
0064   }
0065   if (i > 2000 && i < 3000 && (i % 1000) / 100 == 3) {
0066     dz = 0.8;
0067     dr = 5.637;
0068   }
0069
0070   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 1) {
0071     dz = 0.8;
0072     dr = 4.362;
0073   }  // TEC
0074   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 2) {
0075     dz = 0.8;
0076     dr = 4.512;
0077   }
0078   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 3) {
0079     dz = 0.8;
0080     dr = 5.637;
0081   }
0082   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 4) {
0083     dz = 0.8;
0084     dr = 5.862;
0085   }
0086   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 5) {
0087     dz = 0.8;
0088     dr = 7.501;
0089   }
0090   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 6) {
0091     dz = 0.8;
0092     dr = 9.336;
0093   }
0094   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 7) {
0095     dz = 0.8;
0096     dr = 10.373;
0097   }
0098
0099   std::pair<float, float> res(dz, dr);
0100   return res;
0101 }
0102
0103 std::pair<float, float> phase2bin(int i) {
0104   float dz = -1;
0105   float dr = -1;
0106
0107   if (i > 1000 && i < 1040) {
0108     dz = 3.33;
0109     dr = 0.4;
0110   }
0111   if (i > 1040 && i < 1130) {
0112     dr = 3.33;
0113     dz = 0.4;
0114   }
0115
0116   if (i > 2000 && i < 2550) {
0117     dz = 2.5;
0118     dr = 0.1;
0119   }
0120   if (i > 2550 && i < 3000) {
0121     dz = 5.0;
0122     dr = 0.1;
0123   }
0124
0125   if (i > 3000 && i < 5000) {
0126     dz = 0.2;
0127     dr = 5.0;
0128   }
0129
0130   if (i > 3100 && i < 3119) {
0131     dz = 0.2;
0132     dr = 2.5;
0133   }
0134   if (i > 3140 && i < 3159) {
0135     dz = 0.2;
0136     dr = 2.5;
0137   }
0138   if (i > 3180 && i < 3199) {
0139     dz = 0.2;
0140     dr = 2.5;
0141   }
0142   if (i > 3220 && i < 3239) {
0143     dz = 0.2;
0144     dr = 2.5;
0145   }
0146   if (i > 3260 && i < 3279) {
0147     dz = 0.2;
0148     dr = 2.5;
0149   }
0150
0151   if (i > 4100 && i < 4119) {
0152     dz = 0.2;
0153     dr = 2.5;
0154   }
0155   if (i > 4140 && i < 4159) {
0156     dz = 0.2;
0157     dr = 2.5;
0158   }
0159   if (i > 4180 && i < 4199) {
0160     dz = 0.2;
0161     dr = 2.5;
0162   }
0163   if (i > 4220 && i < 4239) {
0164     dz = 0.2;
0165     dr = 2.5;
0166   }
0167   if (i > 4260 && i < 4279) {
0168     dz = 0.2;
0169     dr = 2.5;
0170   }
0171
0172   std::pair<float, float> res(dz, dr);
0173   return res;
0174 }
0175
0176 std::pair<float, float> phase1bin(int i) {
0177   float dz = -1;
0178   float dr = -1;
0179
0180   if (i > 1000 && i < 1040) {
0181     dz = 3.33;
0182     dr = 0.4;
0183   }  // BPIX
0184   if (i > 1040 && i < 1080) {
0185     dr = 3.33;
0186     dz = 0.4;
0187   }  // FPIX
0188
0189   if (i > 1100 && i < 2000) {
0190     dz = 5.948;
0191     dr = 0.4;
0192   }  // TIB
0193
0194   if (i > 3000 && i < 4000) {
0195     dz = 9.440;
0196     dr = 0.4;
0197   }  // TOB
0198
0199   if (i > 2000 && i < 3000 && (i % 1000) / 100 == 1) {
0200     dz = 0.8;
0201     dr = 5.647;
0202   }  // TID
0203   if (i > 2000 && i < 3000 && (i % 1000) / 100 == 2) {
0204     dz = 0.8;
0205     dr = 4.512;
0206   }
0207   if (i > 2000 && i < 3000 && (i % 1000) / 100 == 3) {
0208     dz = 0.8;
0209     dr = 5.637;
0210   }
0211
0212   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 1) {
0213     dz = 0.8;
0214     dr = 4.362;
0215   }  // TEC
0216   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 2) {
0217     dz = 0.8;
0218     dr = 4.512;
0219   }
0220   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 3) {
0221     dz = 0.8;
0222     dr = 5.637;
0223   }
0224   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 4) {
0225     dz = 0.8;
0226     dr = 5.862;
0227   }
0228   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 5) {
0229     dz = 0.8;
0230     dr = 7.501;
0231   }
0232   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 6) {
0233     dz = 0.8;
0234     dr = 9.336;
0235   }
0236   if (i > 4000 && i < 6000 && (i % 1000) / 100 == 7) {
0237     dz = 0.8;
0238     dr = 10.373;
0239   }
0240
0241   std::pair<float, float> res(dz, dr);
0242   return res;
0243 }
0244
0245 void printFrame(TCanvas* c, TH1D* h, const char* label, const int frame, const int min, const int max, const bool same) {
0246   c->cd(frame);
0247   h->SetAxisRange(min, max);
0248   if (same) {
0249     h->DrawCopy("same");
0250   } else {
0251     h->DrawCopy();
0252     TText* t = new TText((max + min) / 2, h->GetMaximum(), label);
0253     t->SetTextAlign(22);
0254     t->DrawClone();
0255   }
0256 }
0257
0258 void PlotOccupancyMapGeneric(TFile* ff,
0259                              const char* module,
0260                              const float min,
0261                              const float max,
0262                              const float mmin,
0263                              const float mmax,
0264                              const int color,
0265                              std::pair<float, float> (*size)(int),
0266                              std::vector<SubDetParams>& vsub) {
0267   gROOT->SetStyle("Plain");
0268
0269   if (ff->cd(module)) {
0270     TProfile* aveoccu = (TProfile*)gDirectory->Get("aveoccu");
0271     TProfile* avemult = (TProfile*)gDirectory->Get("avemult");
0272     TH1F* nchannels = (TH1F*)gDirectory->Get("nchannels_real");
0273
0274     TProfile* averadius = (TProfile*)gDirectory->Get("averadius");
0275     TProfile* avez = (TProfile*)gDirectory->Get("avez");
0276
0277     std::cout << "pointers " << aveoccu << " " << avemult << " " << nchannels << " " << averadius << " " << avez
0278               << std::endl;
0279
0280     if (aveoccu && avemult && nchannels && averadius && avez) {
0281       nchannels->Sumw2();
0282       for (int i = 1; i < nchannels->GetNbinsX() + 1; ++i) {
0283         nchannels->SetBinError(i, 0.);
0284       }
0285
0286       TH1D* haveoccu = aveoccu->ProjectionX("haveoccu");
0287       haveoccu->SetDirectory(nullptr);
0288       haveoccu->Divide(nchannels);
0289
0290       TH1D* havemult = avemult->ProjectionX("havemult");
0291       havemult->SetDirectory(nullptr);
0292       havemult->Divide(nchannels);
0293
0294       TH1D* havewidth = (TH1D*)haveoccu->Clone("havewidth");
0295       havewidth->SetDirectory(nullptr);
0296       havewidth->SetTitle("Average Cluster Size");
0297       havewidth->Divide(havemult);
0298
0299       new TCanvas("occupancy", "occupancy", 1200, 500);
0300       gPad->SetLogy(1);
0301       haveoccu->SetStats(false);
0302       haveoccu->SetLineColor(kRed);
0303       haveoccu->SetMarkerColor(kRed);
0304       haveoccu->DrawCopy();
0305
0306       new TCanvas("multiplicity", "multiplicity", 1200, 500);
0307       gPad->SetLogy(1);
0308       havemult->SetStats(false);
0309       havemult->SetLineColor(kRed);
0310       havemult->SetMarkerColor(kRed);
0311       havemult->DrawCopy();
0312
0313       new TCanvas("width", "width", 1200, 500);
0314       havewidth->SetStats(false);
0315       havewidth->SetLineColor(kRed);
0316       havewidth->SetMarkerColor(kRed);
0317       havewidth->DrawCopy();
0318
0319       bool same = false;
0320       TCanvas* o2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("occupancy2");
0321       if (o2) {
0322         same = true;
0323         haveoccu->SetLineColor(kBlue);
0324         haveoccu->SetMarkerColor(kBlue);
0325       } else {
0326         o2 = new TCanvas("occupancy2", "occupancy2", 1200, 800);
0327         o2->Divide(3, 2);
0328       }
0329       for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
0330         printFrame(o2, haveoccu, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
0331       }
0332
0333       same = false;
0334       TCanvas* m2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("multiplicity2");
0335       if (m2) {
0336         same = true;
0337         havemult->SetLineColor(kBlue);
0338         havemult->SetMarkerColor(kBlue);
0339       } else {
0340         m2 = new TCanvas("multiplicity2", "multiplicity2", 1200, 800);
0341         m2->Divide(3, 2);
0342       }
0343       for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
0344         printFrame(m2, havemult, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
0345       }
0346
0347       same = false;
0348       TCanvas* w2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("width2");
0349       if (w2) {
0350         same = true;
0351         havewidth->SetLineColor(kBlue);
0352         havewidth->SetMarkerColor(kBlue);
0353       } else {
0354         w2 = new TCanvas("width2", "width2", 1200, 800);
0355         w2->Divide(3, 2);
0356       }
0357       for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
0358         printFrame(w2, havewidth, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
0359       }
0360
0361       float (*scale)(float);
0362       scale = &logarithm;
0363
0364       drawMap("multmap", havemult, averadius, avez, mmin, mmax, size, scale, color);
0365       drawMap("occumap", haveoccu, averadius, avez, min, max, size, scale, color, "channel occupancy");
0366     }
0367   }
0368 }
0369
0370 float combinedOccupancy(TFile* ff, const char* module, const int lowerbin, const int upperbin) {
0371   float cumoccu = -2.;
0372   double cumerr = -2;
0373
0374   if (ff->cd(module)) {
0375     TProfile* aveoccu = (TProfile*)gDirectory->Get("aveoccu");
0376     //    TProfile* avemult= (TProfile*)gDirectory->Get("avemult");
0377     TH1F* nchannels = (TH1F*)gDirectory->Get("nchannels_real");
0378
0379     float sumoccu = 0.;
0380     float sumnchannels = 0;
0381     double sumerrsq = 0;
0382
0383     for (int i = lowerbin; i < upperbin + 1; ++i) {
0384       std::cout << "processing bin " << i << " " << aveoccu->GetBinContent(i) << "+/-" << aveoccu->GetBinError(i)
0385                 << std::endl;
0386       sumoccu += aveoccu->GetBinContent(i);
0387       sumnchannels += nchannels->GetBinContent(i);
0388       sumerrsq += aveoccu->GetBinError(i) * aveoccu->GetBinError(i);
0389     }
0390     cumoccu = sumnchannels != 0 ? sumoccu / sumnchannels : -1;
0391     cumerr = sumnchannels != 0 ? sqrt(sumerrsq) / sumnchannels : -1;
0392     std::cout << "Cumulative occupancy: " << sumoccu << " " << sumnchannels << " " << cumoccu << "+/-" << cumerr;
0393   }
0394
0395   return cumoccu;
0396 }
0397
0398 void PlotOccupancyMap(TFile* ff,
0399                       const char* module,
0400                       const float min,
0401                       const float max,
0402                       const float mmin,
0403                       const float mmax,
0404                       const int color) {
0405   std::pair<float, float> (*size)(int);
0406   size = &presentbin;
0407
0408   std::vector<SubDetParams> vsub;
0409   SubDetParams ppix = {"BPIX+FPIX", 100, 270};
0410   vsub.push_back(ppix);
0411   SubDetParams ptib = {"TIB", 1050, 1450};
0412   vsub.push_back(ptib);
0413   SubDetParams ptid = {"TID", 2070, 2400};
0414   vsub.push_back(ptid);
0415   SubDetParams ptob = {"TOB", 3000, 3700};
0416   vsub.push_back(ptob);
0417   SubDetParams ptecm = {"TEC-", 4000, 4850};
0418   vsub.push_back(ptecm);
0419   SubDetParams ptecp = {"TEC+", 5000, 5850};
0420   vsub.push_back(ptecp);
0421
0422   PlotOccupancyMapGeneric(ff, module, min, max, mmin, mmax, color, size, vsub);
0423 }
0424
0425 void PlotOccupancyMapPhase1(TFile* ff,
0426                             const char* module,
0427                             const float min,
0428                             const float max,
0429                             const float mmin,
0430                             const float mmax,
0431                             const int color) {
0432   std::pair<float, float> (*size)(int);
0433   size = &phase1bin;
0434
0435   std::vector<SubDetParams> vsub;
0436   SubDetParams ppix = {"BPIX+FPIX", 1000, 1080};
0437   vsub.push_back(ppix);
0438   SubDetParams ptib = {"TIB", 1090, 1450};
0439   vsub.push_back(ptib);
0440   SubDetParams ptid = {"TID", 2070, 2400};
0441   vsub.push_back(ptid);
0442   SubDetParams ptob = {"TOB", 3000, 3700};
0443   vsub.push_back(ptob);
0444   SubDetParams ptecm = {"TEC-", 4000, 4850};
0445   vsub.push_back(ptecm);
0446   SubDetParams ptecp = {"TEC+", 5000, 5850};
0447   vsub.push_back(ptecp);
0448
0449   PlotOccupancyMapGeneric(ff, module, min, max, mmin, mmax, color, size, vsub);
0450 }
0451
0452 void PlotOccupancyMapPhase2(TFile* ff,
0453                             const char* module,
0454                             const float min,
0455                             const float max,
0456                             const float mmin,
0457                             const float mmax,
0458                             const int color) {
0459   std::pair<float, float> (*size)(int);
0460   size = &phase2bin;
0461
0462   std::vector<SubDetParams> vsub;
0463   SubDetParams ppix = {"BPIX+FPIX", 1000, 1090};
0464   vsub.push_back(ppix);
0465   SubDetParams ptob = {"TOB", 2000, 2900};
0466   vsub.push_back(ptob);
0467   SubDetParams ptecm = {"TEC-", 3100, 3300};
0468   vsub.push_back(ptecm);
0469   SubDetParams ptecp = {"TEC+", 4100, 4300};
0470   vsub.push_back(ptecp);
0471
0472   PlotOccupancyMapGeneric(ff, module, min, max, mmin, mmax, color, size, vsub);
0473 }
0474
0475 void PlotOnTrackOccupancy(
0476     TFile* ff, const char* module, const char* ontrkmod, const float mmin, const float mmax, const int color) {
0477   std::pair<float, float> (*size)(int);
0478   size = &presentbin;
0479
0480   std::vector<SubDetParams> vsub;
0481   SubDetParams ppix = {"BPIX+FPIX", 100, 270};
0482   vsub.push_back(ppix);
0483   SubDetParams ptib = {"TIB", 1050, 1450};
0484   vsub.push_back(ptib);
0485   SubDetParams ptid = {"TID", 2070, 2400};
0486   vsub.push_back(ptid);
0487   SubDetParams ptob = {"TOB", 3000, 3700};
0488   vsub.push_back(ptob);
0489   SubDetParams ptecm = {"TEC-", 4000, 4850};
0490   vsub.push_back(ptecm);
0491   SubDetParams ptecp = {"TEC+", 5000, 5850};
0492   vsub.push_back(ptecp);
0493
0494   PlotOnTrackOccupancyGeneric(ff, module, ontrkmod, mmin, mmax, color, size, vsub);
0495 }
0496
0497 void PlotOnTrackOccupancyPhase1(
0498     TFile* ff, const char* module, const char* ontrkmod, const float mmin, const float mmax, const int color) {
0499   std::pair<float, float> (*size)(int);
0500   size = &phase1bin;
0501
0502   std::vector<SubDetParams> vsub;
0503   SubDetParams ppix = {"BPIX+FPIX", 1000, 1080};
0504   vsub.push_back(ppix);
0505   SubDetParams ptib = {"TIB", 1090, 1450};
0506   vsub.push_back(ptib);
0507   SubDetParams ptid = {"TID", 2070, 2400};
0508   vsub.push_back(ptid);
0509   SubDetParams ptob = {"TOB", 3000, 3700};
0510   vsub.push_back(ptob);
0511   SubDetParams ptecm = {"TEC-", 4000, 4850};
0512   vsub.push_back(ptecm);
0513   SubDetParams ptecp = {"TEC+", 5000, 5850};
0514   vsub.push_back(ptecp);
0515
0516   PlotOnTrackOccupancyGeneric(ff, module, ontrkmod, mmin, mmax, color, size, vsub);
0517 }
0518
0519 void PlotOnTrackOccupancyPhase2(
0520     TFile* ff, const char* module, const char* ontrkmod, const float mmin, const float mmax, const int color) {
0521   std::pair<float, float> (*size)(int);
0522   size = &phase2bin;
0523
0524   std::vector<SubDetParams> vsub;
0525   SubDetParams ppix = {"BPIX+FPIX", 1000, 1090};
0526   vsub.push_back(ppix);
0527   SubDetParams ptob = {"TOB", 2000, 2900};
0528   vsub.push_back(ptob);
0529   SubDetParams ptecm = {"TEC-", 3100, 3300};
0530   vsub.push_back(ptecm);
0531   SubDetParams ptecp = {"TEC+", 4100, 4300};
0532   vsub.push_back(ptecp);
0533
0534   PlotOnTrackOccupancyGeneric(ff, module, ontrkmod, mmin, mmax, color, size, vsub);
0535 }
0536
0537 void PlotOnTrackOccupancyGeneric(TFile* ff,
0538                                  const char* module,
0539                                  const char* ontrkmod,
0540                                  const float mmin,
0541                                  const float mmax,
0542                                  const int color,
0543                                  std::pair<float, float> (*size)(int),
0544                                  const std::vector<SubDetParams>& vsub) {
0545   gROOT->SetStyle("Plain");
0546
0547   TProfile* avemult = nullptr;
0548   TProfile* aveontrkmult = nullptr;
0549   TProfile* averadius = nullptr;
0550   TProfile* avez = nullptr;
0551
0552   if (ff->cd(module)) {
0553     avemult = (TProfile*)gDirectory->Get("avemult");
0554     averadius = (TProfile*)gDirectory->Get("averadius");
0555     avez = (TProfile*)gDirectory->Get("avez");
0556   }
0557   if (ff->cd(ontrkmod))
0558     aveontrkmult = (TProfile*)gDirectory->Get("avemult");
0559
0560   std::cout << "pointers " << avemult << " " << aveontrkmult << " " << averadius << " " << avez << std::endl;
0561
0562   if (averadius && avez && avemult && aveontrkmult) {
0563     TH1D* havemult = avemult->ProjectionX("havemult");
0564     TH1D* haveontrkmult = aveontrkmult->ProjectionX("haveontrkmult");
0565     havemult->SetDirectory(nullptr);
0566     haveontrkmult->SetDirectory(nullptr);
0567     haveontrkmult->Divide(havemult);
0568
0569     new TCanvas("ontrkmult", "ontrkmult", 1200, 500);
0570     gPad->SetLogy(1);
0571     haveontrkmult->SetStats(false);
0572     haveontrkmult->SetLineColor(kRed);
0573     haveontrkmult->SetMarkerColor(kRed);
0574     haveontrkmult->SetMarkerSize(.5);
0575     haveontrkmult->SetMarkerStyle(20);
0576     haveontrkmult->DrawCopy();
0577
0578     bool same = false;
0579     TCanvas* o2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("ontrkmult2");
0580     if (o2) {
0581       same = true;
0582       haveontrkmult->SetLineColor(kBlue);
0583       haveontrkmult->SetMarkerColor(kBlue);
0584     } else {
0585       o2 = new TCanvas("ontrkmult2", "ontrkmult2", 1200, 800);
0586       o2->Divide(3, 2);
0587     }
0588     for (unsigned int isub = 0; isub < vsub.size(); ++isub) {
0589       printFrame(o2, haveontrkmult, vsub[isub].label.c_str(), isub + 1, vsub[isub].min, vsub[isub].max, same);
0590     }
0591
0592     float (*scale)(float);
0593     scale = &linear;
0594
0595     drawMap("ontrkmultmap", haveontrkmult, averadius, avez, mmin, mmax, size, scale, color);
0596   }
0597 }
0598
0599 TCanvas* drawMap(const char* cname,
0600                  const TH1* hval,
0601                  const TProfile* averadius,
0602                  const TProfile* avez,
0603                  const float mmin,
0604                  const float mmax,
0605                  std::pair<float, float> (*size)(int),
0606                  float (*scale)(float),
0607                  const int color,
0608                  const char* ptitle) {
0609   if (color == 1) {
0610     // A not-so-great color version
0611     const Int_t NRGBs = 5;
0612     const Int_t NCont = 255;
0613     Double_t stops[NRGBs] = {0.00, 0.25, 0.50, 0.75, 1.00};
0614     Double_t red[NRGBs] = {0.00, 0.00, 0.40, 1.00, 1.00};
0615     Double_t green[NRGBs] = {0.00, 0.40, 0.70, 0.60, 1.00};
0616     Double_t blue[NRGBs] = {0.30, 0.60, 0.00, 0.00, 0.20};
0617     TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
0618     gStyle->SetNumberContours(NCont);
0619   } else if (color == 2) {
0620     // Gray scale
0621     const Int_t NRGBs = 3;
0622     const Int_t NCont = 255;
0623     Double_t stops[NRGBs] = {0.00, 0.50, 1.00};
0624     Double_t red[NRGBs] = {0.90, 0.50, 0.00};
0625     Double_t green[NRGBs] = {0.90, 0.50, 0.00};
0626     Double_t blue[NRGBs] = {0.90, 0.50, 0.00};
0627     TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
0628     gStyle->SetNumberContours(NCont);
0629   } else if (color == 3) {
0630     // used by Kevin in the TRK-11-001 paper
0631     const Int_t NRGBs = 7;
0632     const Int_t NCont = 255;
0633     Double_t stops[NRGBs] = {0.00, 0.15, 0.30, 0.45, 0.65, 0.85, 1.00};
0634     Double_t red[NRGBs] = {0.60, 0.30, 0.00, 0.00, 0.60, 0.40, 0.00};
0635     Double_t green[NRGBs] = {1.00, 0.90, 0.80, 0.75, 0.20, 0.00, 0.00};
0636     Double_t blue[NRGBs] = {1.00, 1.00, 1.00, 0.30, 0.00, 0.00, 0.00};
0637     TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
0638     gStyle->SetNumberContours(NCont);
0639   }
0640
0641   int ncol = gStyle->GetNumberOfColors();
0642   std::cout << "Number of colors " << ncol << std::endl;
0643   // Loop on bins and creation of boxes
0644
0645   TList modulesmult;
0646
0647   for (int i = 1; i < hval->GetNbinsX(); ++i) {
0648     if ((averadius->GetBinEntries(i) * avez->GetBinEntries(i)) != 0) {
0649       double dz = -1.;
0650       double dr = -1.;
0651       // determine module size
0652
0653       dz = (*size)(i).first;
0654       dr = (*size)(i).second;
0655
0656       if (dz < 0 && dr < 0)
0657         continue;
0658
0659       {
0660         TBox* modmult = new TBox(avez->GetBinContent(i) - dz,
0661                                  averadius->GetBinContent(i) - dr,
0662                                  avez->GetBinContent(i) + dz,
0663                                  averadius->GetBinContent(i) + dr);
0664         modmult->SetFillStyle(1001);
0665         if (color < 0) {
0666           modmult->SetFillColor(kBlack);
0667         } else {
0668           int icol = int(ncol * (scale(hval->GetBinContent(i)) - scale(mmin)) / (scale(mmax) - scale(mmin)));
0669           if (icol < 0)
0670             icol = 0;
0671           if (icol > (ncol - 1))
0672             icol = (ncol - 1);
0673           std::cout << i << " " << icol << " " << hval->GetBinContent(i) << std::endl;
0674           modmult->SetFillColor(gStyle->GetColorPalette(icol));
0675         }
0676         modulesmult.Add(modmult);
0677       }
0678     }
0679   }
0680   // eta boundaries lines
0681   double etavalext[] = {4., 3.5, 3., 2.8, 2.6, 2.4, 2.2, 2.0, 1.8, 1.6};
0682   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};
0683   TList etalines;
0684   TList etalabels;
0685   TList paperlabels;
0686   for (int i = 0; i < 10; ++i) {
0687     //  double eta = 3.0-i*0.2;
0688     double eta = etavalext[i];
0689     TLine* lin = new TLine(295, 2 * 295 / (exp(eta) - exp(-eta)), 305, 2 * 305 / (exp(eta) - exp(-eta)));
0690     etalines.Add(lin);
0691     char lab[100];
0692     sprintf(lab, "%3.1f", eta);
0693     TText* label = new TText(285, 2 * 285 / (exp(eta) - exp(-eta)), lab);
0694     label->SetTextSize(.03);
0695     label->SetTextAlign(22);
0696     etalabels.Add(label);
0697   }
0698   for (int i = 0; i < 10; ++i) {
0699     //  double eta = -3.0+i*0.2;
0700     double eta = -1 * etavalext[i];
0701     TLine* lin = new TLine(-295, -2 * 295 / (exp(eta) - exp(-eta)), -305, -2 * 305 / (exp(eta) - exp(-eta)));
0702     etalines.Add(lin);
0703     char lab[100];
0704     sprintf(lab, "%3.1f", eta);
0705     TText* label = new TText(-285, -2 * 285 / (exp(eta) - exp(-eta)), lab);
0706     label->SetTextSize(.03);
0707     label->SetTextAlign(22);
0708     etalabels.Add(label);
0709   }
0710   for (int i = 0; i < 15; ++i) {
0711     //  double eta = -1.4+i*0.2;
0712     double eta = etavalint[i];
0713     TLine* lin = new TLine(130. * (exp(eta) - exp(-eta)) / 2., 130, 138. * (exp(eta) - exp(-eta)) / 2., 138);
0714     etalines.Add(lin);
0715     char lab[100];
0716     sprintf(lab, "%3.1f", eta);
0717     TText* label = new TText(125. * (exp(eta) - exp(-eta)) / 2., 125, lab);
0718     label->SetTextSize(.03);
0719     label->SetTextAlign(22);
0720     etalabels.Add(label);
0721   }
0722   TLatex* etalab = new TLatex(0, 115, "#eta");
0723   etalab->SetTextSize(.03);
0724   etalab->SetTextAlign(22);
0725   etalabels.Add(etalab);
0726
0727   // CMS label
0728   TLatex* cmslab = new TLatex(0.15, 0.965, "CMS");
0729   cmslab->SetNDC();
0730   cmslab->SetTextSize(0.04);
0731   cmslab->SetTextAlign(31);
0732   paperlabels.Add(cmslab);
0733   TLatex* enelab = new TLatex(0.92, 0.965, "#sqrt{s} = 7 TeV");
0734   enelab->SetNDC();
0735   enelab->SetTextSize(0.04);
0736   enelab->SetTextAlign(31);
0737   paperlabels.Add(enelab);
0738   /*
0739     TLatex *lumilab = new TLatex(0.6,0.965,Form("L = %.1f  fb^{-1}",19.7));
0740     lumilab->SetNDC();
0741     lumilab->SetTextSize(0.04);
0742     lumilab->SetTextAlign(31);
0743     paperlabels.Add(lumilab);
0744   */
0745
0746   TGaxis* raxis = new TGaxis(-310, 0, -310, 140, 0, 140, 10, "S");
0747   TGaxis* zaxis = new TGaxis(-310, 0, 310, 0, -310, 310, 10, "S");
0748   raxis->SetTickSize(.01);
0749   zaxis->SetTickSize(.01);
0750   raxis->SetTitle("r (cm)");
0751   zaxis->SetTitle("z (cm)");
0752
0753   TList mpalette;
0754
0755   for (int i = 0; i < ncol; ++i) {
0756     TBox* box = new TBox(315, 0 + 140. / ncol * i, 330, 0 + 140. / ncol * (i + 1));
0757     box->SetFillStyle(1001);
0758     box->SetFillColor(gStyle->GetColorPalette(i));
0759     mpalette.Add(box);
0760   }
0761
0762   TGaxis* mpaxis = nullptr;
0763   if (scale(1) != 1) {
0764     mpaxis = new TGaxis(330, 0, 330, 140, mmin, mmax, 510, "SLG+");
0765   } else {
0766     mpaxis = new TGaxis(330, 0, 330, 140, mmin, mmax, 510, "SL+");
0767   }
0768   mpaxis->SetTickSize(.02);
0769   mpaxis->SetLabelOffset(mpaxis->GetLabelOffset() * 0.5);
0770   mpaxis->SetTitle(ptitle);
0771   mpalette.Add(mpaxis);
0772
0773   TCanvas* cc2 = new TCanvas(cname, cname, 1000, 500);
0774   cc2->Range(-370., -20., 390., 150.);
0775   TFrame* fr2 = new TFrame(-310, 0, 310, 140);
0776   fr2->UseCurrentStyle();
0777   fr2->Draw();
0778   raxis->Draw();
0779   zaxis->Draw();
0780   std::cout << modulesmult.GetSize() << std::endl;
0781   etalines.Draw();
0782   etalabels.Draw();
0783   if (color >= 0)
0784     mpalette.Draw();
0785   modulesmult.Draw();
0786
0787   return cc2;
0788 }
0789
0790 void PlotDebugFPIX_XYMap(TFile* ff, const char* module, const unsigned int ioffset, const char* name) {
0791   gROOT->SetStyle("Plain");
0792
0793   TCanvas* cc = new TCanvas(name, name, 750, 750);
0794   cc->Range(-25, -25, 25, 25);
0795   TFrame* fr1 = new TFrame(-20, -20, 20, 20);
0796   fr1->UseCurrentStyle();
0797   fr1->Draw();
0798   ff->cd(module);
0799   gDirectory->ls();
0800   TProfile* avex = (TProfile*)gDirectory->Get("avex");
0801   TProfile* avey = (TProfile*)gDirectory->Get("avey");
0802   TProfile* avez = (TProfile*)gDirectory->Get("avez");
0803
0804   if (avex && avey && avez) {
0805     TText* tittext = new TText(0, 0, name);
0806     tittext->SetTextSize(.04);
0807     tittext->SetTextAlign(22);
0808     tittext->Draw();
0809     for (unsigned int mod = ioffset + 1; mod < ioffset + 57; ++mod) {
0810       double x = avex->GetBinContent(mod);
0811       double y = avey->GetBinContent(mod);
0812       //      TBox* modbox = new TBox(x-1,y-1,x+1,y+1);
0813       char modstring[30];
0814       sprintf(modstring, "%d", mod % 100);
0815       TText* modtext = new TText(x, y, modstring);
0816       modtext->SetTextAngle(atan(y / x) * 180 / 3.14159);
0817       modtext->SetTextSize(.02);
0818       modtext->SetTextAlign(22);
0819       modtext->SetTextColor(kRed);
0820       std::cout << mod << " " << x << " " << y << std::endl;
0821       //      modbox->Draw();
0822       modtext->Draw();
0823     }
0824     for (unsigned int mod = ioffset + 101; mod < ioffset + 157; ++mod) {
0825       double x = avex->GetBinContent(mod);
0826       double y = avey->GetBinContent(mod);
0827       //      TBox* modbox = new TBox(x-1,y-1,x+1,y+1);
0828       char modstring[30];
0829       sprintf(modstring, "%d", mod % 100);
0830       TText* modtext = new TText(x, y, modstring);
0831       modtext->SetTextAngle(atan(y / x) * 180 / 3.14159);
0832       modtext->SetTextSize(.02);
0833       modtext->SetTextAlign(22);
0834       modtext->SetTextColor(kBlue);
0835       std::cout << mod << " " << x << " " << y << " " << atan(y / x) << std::endl;
0836       //      modbox->Draw();
0837       modtext->Draw();
0838     }
0839   }
0840 }
0841
0842 void PlotTrackerXsect(TFile* ff, const char* module) {
0843   gROOT->SetStyle("Plain");
0844
0845   if (ff->cd(module)) {
0846     TProfile* averadius = (TProfile*)gDirectory->Get("averadius");
0847     TProfile* avez = (TProfile*)gDirectory->Get("avez");
0848
0849     std::cout << "pointers " << averadius << " " << avez << std::endl;
0850
0851     if (averadius && avez) {
0852       std::pair<float, float> (*size)(int);
0853       size = &presentbin;
0854       float (*scale)(float);
0855       scale = &linear;
0856
0857       drawMap("trackermap", averadius, averadius, avez, 0, 0, size, scale, -1);
0858     }
0859   }
0860 }
0861
0862 TH1D* TrendPlotSingleBin(TFile* ff, const char* module, const char* hname, const int bin) {
0863   CommonAnalyzer caoccu(ff, "", module);
0864
0865   TH1D* occutrend = new TH1D("occutrend", "Average number of clusters vs run", 10, 0., 10.);
0866   occutrend->SetCanExtend(TH1::kXaxis);
0867   occutrend->Sumw2();
0868
0869   std::vector<unsigned int> runs = caoccu.getRunList();
0870   std::sort(runs.begin(), runs.end());
0871
0872   {
0873     for (unsigned int i = 0; i < runs.size(); ++i) {
0874       char runlabel[100];
0875       sprintf(runlabel, "%d", runs[i]);
0876       char runpath[100];
0877       sprintf(runpath, "run_%d", runs[i]);
0878       caoccu.setPath(runpath);
0879
0880       TProfile* occu = nullptr;
0881       if (occu == nullptr)
0882         occu = (TProfile*)caoccu.getObject(hname);
0883       if (occu) {
0884         const int ibin = occu->FindBin(bin);
0885         std::cout << runlabel << " "
0886                   << " " << ibin << " " << occu->GetBinContent(ibin) << " " << occu->GetBinError(ibin) << std::endl;
0887         const int jbin = occutrend->Fill(runlabel, occu->GetBinContent(ibin));
0888         occutrend->SetBinError(jbin, occu->GetBinError(ibin));
0889       }
0890     }
0891   }
0892   return occutrend;
0893 }