Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:57 

0001 //////////////////////////////////////////////////////////////////////////////
0002 // Usage:
0003 // .L PlotMPV.C+g
0004 //             For carrying out exponential fit to the MPV's as a function
0005 //             of accumulated luminosity and plotting them on canvas
0006 //     plotMPV(infile, eta, phi, depth, first, drawleg, iyear, lumis, ener,
0007 //             save, npar, debug);
0008 //
0009 //             For combining the fit results for 2 iphi values (63, 65) by
0010 //             weighted mean, carrying out the fit and plotting them on canvas
0011 //     plot2MPV(infile, eta, depth, first, drawleg, iyear, lumis, ener, save,
0012 //              npar, debug)
0013 //
0014 //             Get truncated mean for a range of eta, depth, ...
0015 //     getTruncateMean(infile, frac, type, nvx, save);
0016 //     getTruncatedMeanX(infile, type, nvx, year, phi, save);
0017 //
0018 //             Draw slopes as a function of luminosity
0019 //     drawSlope(infile, type, phi, depth, drawleg, iyear, lumis, ener,
0020 //               save, debug)
0021 //
0022 //             Draw the 2D plot of light loss as a function of eta, depth
0023 //      plotLightLoss(infile, lumis, iyear, ener, save, debug)
0024 //
0025 //  where
0026 //   infile (const char*) Nme of the input file
0027 //   eta    (int)         ieta value of the tower
0028 //   phi    (int)         iphi value of the tower (0 if iphi's of RBX combined
0029 //                        for 2017 and of all iphi's for 2018 and beyond)
0030 //   depth  (int)         depth value of the tower (-1 if all depths combined)
0031 //   first  (uint32_t)    the starting lumi point to be used (0)
0032 //   drawleg (bool)       legend to be drawn or not (true)
0033 //   iyear  (int)         Year of data taking, if > 1000 to be shown (17)
0034 //   lumis  (double)      Total integrated luminosity (49.0)
0035 //   ener   (int)         C.M. energy (13)
0036 //   npar   (int)         Number of parameters in the fit (1)
0037 //   save   (bool)        if the plot to be saved as a pdf file (false)
0038 //   debug  (bool)        For controlling debug printouts (false)
0039 //   frac   (double)      Fraction of events to be used
0040 //   type   (string)      Specifies the range (e.g. Data2018D7)
0041 //   nvx    (int)         Vertex bin to be considered (e.g. 2)
0042 //   year   (int)         Cuts made for 2017 or 2018
0043 //
0044 ///////////////////////////////////////////////////////////////////////////////
0045 
0046 #include <TArrow.h>
0047 #include <TASImage.h>
0048 #include <TCanvas.h>
0049 #include <TChain.h>
0050 #include <TFile.h>
0051 #include <TF1.h>
0052 #include <TFitResult.h>
0053 #include <TFitResultPtr.h>
0054 #include <TGraph.h>
0055 #include <TGraph2D.h>
0056 #include <TGraphErrors.h>
0057 #include <TGraphAsymmErrors.h>
0058 #include <TH1D.h>
0059 #include <TH2D.h>
0060 #include <TLatex.h>
0061 #include <TLegend.h>
0062 #include <TPaveStats.h>
0063 #include <TPaveText.h>
0064 #include <TProfile.h>
0065 #include <TROOT.h>
0066 #include <TStyle.h>
0067 
0068 #include <cstdlib>
0069 #include <iomanip>
0070 #include <iostream>
0071 #include <fstream>
0072 #include <map>
0073 #include <string>
0074 #include <vector>
0075 
0076 const bool debugMode = false;
0077 
0078 std::vector<std::string> splitString(const std::string& fLine) {
0079   std::vector<std::string> result;
0080   int start = 0;
0081   bool empty = true;
0082   for (unsigned i = 0; i <= fLine.size(); i++) {
0083     if (fLine[i] == ' ' || i == fLine.size()) {
0084       if (!empty) {
0085         std::string item(fLine, start, i - start);
0086         result.push_back(item);
0087         empty = true;
0088       }
0089       start = i + 1;
0090     } else {
0091       if (empty)
0092         empty = false;
0093     }
0094   }
0095   return result;
0096 }
0097 
0098 void meanMax(const char* infile, bool debug) {
0099   std::map<std::string, std::pair<double, double> > means;
0100   std::ifstream fInput(infile);
0101   if (!fInput.good()) {
0102     std::cout << "Cannot open file " << infile << std::endl;
0103   } else {
0104     char buffer[1024];
0105     unsigned int all(0), good(0);
0106     while (fInput.getline(buffer, 1024)) {
0107       ++all;
0108       if (buffer[0] == '-')
0109         continue;  //ignore comment
0110       std::vector<std::string> items = splitString(std::string(buffer));
0111       if (items.size() < 5) {
0112         if (debug)
0113           std::cout << "Ignore " << items.size() << " in line: " << buffer << std::endl;
0114       } else {
0115         ++good;
0116         double mean = std::atof(items[4].c_str());
0117         std::map<std::string, std::pair<double, double> >::iterator itr = means.find(items[2]);
0118         if (itr == means.end()) {
0119           means[items[2]] = std::pair<double, double>(mean, mean);
0120         } else {
0121           double mmin = std::min(mean, itr->second.first);
0122           double mmax = std::max(mean, itr->second.second);
0123           means[items[2]] = std::pair<double, double>(mmin, mmax);
0124         }
0125       }
0126     }
0127     fInput.close();
0128     std::cout << "Reads total of " << all << " and " << good << " good "
0129               << " records from " << infile << std::endl;
0130     for (std::map<std::string, std::pair<double, double> >::iterator itr = means.begin(); itr != means.end(); ++itr)
0131       std::cout << itr->first << " Mean " << itr->second.first << ":" << itr->second.second << std::endl;
0132   }
0133 }
0134 
0135 void readMPVs(const char* infile,
0136               const int eta,
0137               const int phi,
0138               const int dep,
0139               const unsigned int maxperiod,
0140               std::map<int, std::pair<double, double> >& mpvs,
0141               bool debug) {
0142   mpvs.clear();
0143   std::ifstream fInput(infile);
0144   if (!fInput.good()) {
0145     std::cout << "Cannot open file " << infile << std::endl;
0146   } else {
0147     char buffer[1024];
0148     unsigned int all(0), good(0), select(0);
0149     while (fInput.getline(buffer, 1024)) {
0150       ++all;
0151       if (buffer[0] == '#')
0152         continue;  //ignore comment
0153       std::vector<std::string> items = splitString(std::string(buffer));
0154       if (items.size() < 7) {
0155         if (debug)
0156           std::cout << "Ignore  line: " << buffer << std::endl;
0157       } else {
0158         ++good;
0159         int period = std::atoi(items[0].c_str());
0160         int ieta = std::atoi(items[1].c_str());
0161         int iphi = std::atoi(items[2].c_str());
0162         int depth = std::atoi(items[3].c_str());
0163         double mpv = std::atof(items[5].c_str());
0164         double dmpv = std::atof(items[6].c_str());
0165         if ((ieta == eta) && (iphi == phi) && (depth == dep) && (period > 0) && (period <= (int)(maxperiod))) {
0166           ++select;
0167           mpvs[period] = std::pair<double, double>(mpv, dmpv);
0168         }
0169       }
0170     }
0171     fInput.close();
0172     std::cout << "Reads total of " << all << " and " << good << " good and " << select << " selected records from "
0173               << infile << std::endl;
0174   }
0175   if (debug) {
0176     unsigned int k(0);
0177     for (std::map<int, std::pair<double, double> >::const_iterator itr = mpvs.begin(); itr != mpvs.end(); ++itr, ++k) {
0178       std::cout << "[" << k << "] Period " << itr->first << " MPV " << (itr->second).first << " +- "
0179                 << (itr->second).second << std::endl;
0180     }
0181   }
0182 }
0183 
0184 struct rType {
0185   TCanvas* pad;
0186   int bin;
0187   double mean, error, chisq;
0188   rType() {
0189     pad = nullptr;
0190     bin = 0;
0191     mean = error = chisq = 0;
0192   }
0193 };
0194 
0195 rType drawMPV(const char* name,
0196               int eta,
0197               int phi,
0198               int depth,
0199               const unsigned int nmax,
0200               const std::map<int, std::pair<double, double> >& mpvs,
0201               unsigned int first = 0,
0202               bool drawleg = true,
0203               int iyear = 17,
0204               double lumis = 0,
0205               int ener = 13,
0206               int normalize = 0,
0207               bool save = false,
0208               int npar = 1,
0209               bool debug = false) {
0210   rType results;
0211   gStyle->SetCanvasBorderMode(0);
0212   gStyle->SetCanvasColor(kWhite);
0213   gStyle->SetPadColor(kWhite);
0214   gStyle->SetFillColor(kWhite);
0215   gStyle->SetOptTitle(0);
0216   gStyle->SetOptStat(0);
0217   gStyle->SetOptFit(0);
0218   results.pad = new TCanvas(name, name);
0219   results.pad->SetRightMargin(0.10);
0220   results.pad->SetTopMargin(0.10);
0221   TLegend* legend = new TLegend(0.25, 0.84, 0.59, 0.89);
0222   legend->SetFillColor(kWhite);
0223   int iy = iyear % 100;
0224   std::vector<float> lumi, dlumi;
0225   double xmin(0), xmax(0);
0226   if (iy == 17) {
0227     xmax = 50;
0228     float lum[13] = {
0229         1.361, 4.540, 8.315, 12.319, 16.381, 20.828, 25.192, 29.001, 32.206, 35.634, 39.739, 43.782, 47.393};
0230     for (unsigned int k = 0; k < nmax; ++k) {
0231       lumi.push_back(lum[k]);
0232       dlumi.push_back(0);
0233     }
0234   } else {
0235     xmax = 70;
0236     float lum[20] = {1.531,  4.401,  7.656,  11.217, 14.133, 17.382, 21.332, 25.347, 28.982, 32.102,
0237                      35.191, 38.326, 41.484, 44.645, 47.740, 50.790, 53.867, 56.918, 59.939, 64.453};
0238     for (unsigned int k = 0; k < nmax; ++k) {
0239       lumi.push_back(lum[k]);
0240       dlumi.push_back(0);
0241     }
0242   }
0243   std::cout << iyear << ":" << iy << " N " << nmax << " X " << xmin << ":" << xmax << std::endl;
0244   const int p1[6] = {2, 5, 6, 11, 15, 19};
0245   if (mpvs.size() > 0) {
0246     float lums[nmax], dlum[nmax], mpv[nmax], dmpv[nmax];
0247     unsigned int np = mpvs.size();
0248     if (np > nmax)
0249       np = nmax;
0250     float ymin(0), ymax(0);
0251     unsigned int k(0), k1(0);
0252     double startvalues[2];
0253     startvalues[0] = 0;
0254     for (std::map<int, std::pair<double, double> >::const_iterator itr = mpvs.begin(); itr != mpvs.end(); ++itr, ++k1) {
0255       if (k1 == 0)
0256         startvalues[1] = (itr->second).first;
0257       if (k1 >= first && k < np) {
0258         mpv[k] = (itr->second).first;
0259         dmpv[k] = (itr->second).second;
0260         if (mpvs.size() == 6) {
0261           lums[k] = lumi[p1[k1]];
0262         } else {
0263           lums[k] = lumi[k1];
0264         }
0265         dlum[k] = 0;
0266         if (k == 0) {
0267           ymin = mpv[k] - dmpv[k];
0268           ymax = mpv[k] + dmpv[k];
0269         } else {
0270           if (mpv[k] - dmpv[k] < ymin)
0271             ymin = mpv[k] - dmpv[k];
0272           if (mpv[k] + dmpv[k] > ymax)
0273             ymax = mpv[k] + dmpv[k];
0274         }
0275         if (debug)
0276           std::cout << "[" << k << "] = " << mpv[k] << " +- " << dmpv[k] << "\n";
0277         ++k;
0278       }
0279     }
0280     np = k;
0281     if (debug)
0282       std::cout << "YRange (Initlal) : " << ymin << ":" << ymax << ":" << startvalues[1];
0283     if (normalize % 10 > 0) {
0284       for (unsigned int k = 0; k < np; ++k) {
0285         mpv[k] /= startvalues[1];
0286         dmpv[k] /= startvalues[1];
0287         if (debug)
0288           std::cout << "[" << k << "] = " << mpv[k] << " +- " << dmpv[k] << "\n";
0289       }
0290       ymin /= startvalues[1];
0291       ymax /= startvalues[1];
0292       startvalues[1] = (mpvs.size() == 6) ? 5.0 : 1.0;
0293       if (debug)
0294         std::cout << "YRange (Initlal) : " << ymin << ":" << ymax << ":" << startvalues[1];
0295     }
0296     int imin = (int)(0.01 * ymin) - 2;
0297     int imax = (int)(0.01 * ymax) + 3;
0298     if (normalize % 10 > 0) {
0299       ymax = 1.2;
0300       ymin = 0.8;
0301     } else {
0302       ymin = 100 * imin;
0303       ymax = 100 * imax;
0304     }
0305     if (debug)
0306       std::cout << " (Final) : " << ymin << ":" << ymax << std::endl;
0307 
0308     TGraphErrors* g = new TGraphErrors(np, lums, mpv, dlum, dmpv);
0309     TF1* f = (npar == 2) ? (new TF1("f", "[1]*TMath::Exp(-x*[0])")) : (new TF1("f", "TMath::Exp(-x*[0])"));
0310     f->SetParameters(startvalues);
0311     f->SetLineColor(kRed);
0312 
0313     double value1(0), error1(0), factor(0);
0314     for (int j = 0; j < 2; ++j) {
0315       TFitResultPtr Fit = g->Fit(f, "QS");
0316       value1 = Fit->Value(0);
0317       error1 = Fit->FitResult::Error(0);
0318       factor = sqrt(f->GetChisquare() / np);
0319       if (factor > 1.2 && ((normalize / 10) % 10 > 0)) {
0320         for (unsigned int k = 0; k < np; ++k)
0321           dmpv[k] *= factor;
0322         delete g;
0323         g = new TGraphErrors(np, lums, mpv, dlum, dmpv);
0324       } else {
0325         break;
0326       }
0327     }
0328     results.mean = value1;
0329     results.error = error1;
0330     results.chisq = factor * factor;
0331     g->SetMarkerStyle(8);
0332     g->SetMarkerColor(kBlue);
0333 
0334     std::cout << "ieta " << eta << " iphi " << phi << " depth " << depth << " mu " << value1 << " +- " << error1
0335               << " Chisq " << f->GetChisquare() << "/" << np << std::endl;
0336     std::ofstream log1("slopes_newCode.txt", std::ios_base::app | std::ios_base::out);
0337     log1 << "ieta\t" << eta << "\tdepth\t" << depth << "\tmu\t" << value1 << "\terror\t " << error1 << "\tChisq\t"
0338          << f->GetChisquare() << std::endl;
0339     log1.close();
0340 
0341     g->GetXaxis()->SetRangeUser(xmin, xmax);
0342     g->GetYaxis()->SetRangeUser(ymin, ymax);
0343     g->SetMarkerSize(1.5);
0344     if (normalize % 10 > 0) {
0345       g->GetYaxis()->SetTitle("MPV_{Charge} (L) /MPV_{Charge} (0)");
0346     } else {
0347       g->GetYaxis()->SetTitle("MPV_{Charge} (fC)");
0348     }
0349     g->GetXaxis()->SetTitle("Delivered Luminosity (fb^{-1})");
0350     g->GetYaxis()->SetTitleSize(0.04);
0351     g->GetXaxis()->SetTitleSize(0.04);
0352     g->GetXaxis()->SetNdivisions(20);
0353     g->GetXaxis()->SetLabelSize(0.04);
0354     g->GetYaxis()->SetLabelSize(0.04);
0355     g->GetXaxis()->SetTitleOffset(1.0);
0356     g->GetYaxis()->SetTitleOffset(1.2);
0357     g->SetTitle("");
0358     g->Draw("AP");
0359 
0360     char namel[100];
0361     if (phi > 0) {
0362       if (depth > 0)
0363         sprintf(namel, "i#eta %d, i#phi %d, depth %d", eta, phi, depth);
0364       else
0365         sprintf(namel, "i#eta %d, i#phi %d (combined depths)", eta, phi);
0366     } else {
0367       if (depth > 0)
0368         sprintf(namel, "i#eta %d, depth %d", eta, depth);
0369       else
0370         sprintf(namel, "i#eta %d (combined depths)", eta);
0371     }
0372     legend->AddEntry(g, namel, "lp");
0373     if (drawleg) {
0374       legend->Draw("same");
0375       results.pad->Update();
0376 
0377       sprintf(namel, "Slope = %6.4f #pm %6.4f", value1, error1);
0378       TPaveText* txt0 = new TPaveText(0.60, 0.84, 0.89, 0.89, "blNDC");
0379       txt0->SetFillColor(0);
0380       txt0->AddText(namel);
0381       txt0->Draw("same");
0382     }
0383     if (iyear > 1000) {
0384       char txt[100];
0385       TPaveText* txt1 = new TPaveText(0.60, 0.91, 0.90, 0.96, "blNDC");
0386       txt1->SetFillColor(0);
0387       sprintf(txt, "%d, %d TeV %5.1f fb^{-1}", iyear, ener, lumis);
0388       txt1->AddText(txt);
0389       txt1->Draw("same");
0390       TPaveText* txt2 = new TPaveText(0.10, 0.91, 0.18, 0.96, "blNDC");
0391       txt2->SetFillColor(0);
0392       sprintf(txt, "CMS");
0393       txt2->AddText(txt);
0394       txt2->Draw("same");
0395     }
0396     results.pad->Modified();
0397     results.pad->Update();
0398     if (save) {
0399       sprintf(namel, "%s_comb.pdf", results.pad->GetName());
0400       results.pad->Print(namel);
0401     }
0402   }
0403   return results;
0404 }
0405 
0406 rType plotMPV(const char* infile,
0407               int eta,
0408               int phi,
0409               int depth,
0410               unsigned int first = 0,
0411               bool drawleg = true,
0412               int iyear = 17,
0413               double lumis = 0.0,
0414               int ener = 13,
0415               int normalize = 0,
0416               bool save = false,
0417               int npar = 1,
0418               bool debug = false) {
0419   char name[100];
0420   int iy = iyear % 100;
0421   const unsigned int nmax = (iy == 17) ? 13 : 20;
0422   double lumi = ((lumis > 0) ? lumis : ((iy == 17) ? 49.0 : 68.5));
0423   sprintf(name, "mpvE%dF%dD%d", eta, phi, depth);
0424 
0425   std::map<int, std::pair<double, double> > mpvs;
0426   readMPVs(infile, eta, phi, depth, nmax, mpvs, debug);
0427   rType results =
0428       drawMPV(name, eta, phi, depth, nmax, mpvs, first, drawleg, iyear, lumi, ener, normalize, save, npar, debug);
0429   return results;
0430 }
0431 
0432 rType plot2MPV(const char* infile,
0433                int eta,
0434                int depth,
0435                unsigned int first = 0,
0436                bool drawleg = true,
0437                int iyear = 17,
0438                double lumis = 0.0,
0439                int ener = 13,
0440                int normalize = 0,
0441                bool save = false,
0442                int npar = 1,
0443                bool debug = false) {
0444   char name[100];
0445   sprintf(name, "mpvE%dD%d", eta, depth);
0446 
0447   int iy = iyear % 100;
0448   const unsigned int nmax = (iy == 17) ? 13 : 20;
0449   double lumi = ((lumis > 0) ? lumis : ((iy == 17) ? 49.0 : 68.5));
0450   std::map<int, std::pair<double, double> > mpvs1, mpvs2, mpvs;
0451   readMPVs(infile, eta, 63, depth, nmax, mpvs1, debug);
0452   readMPVs(infile, eta, 65, depth, nmax, mpvs2, debug);
0453   if (mpvs1.size() == mpvs2.size()) {
0454     unsigned int k(0);
0455     for (std::map<int, std::pair<double, double> >::const_iterator itr = mpvs1.begin(); itr != mpvs1.end();
0456          ++itr, ++k) {
0457       int period = (itr->first);
0458       double mpv1 = (itr->second).first;
0459       double empv1 = (itr->second).second;
0460       double wt1 = 1.0 / (empv1 * empv1);
0461       double mpv2 = mpvs2[period].first;
0462       double empv2 = mpvs2[period].second;
0463       double wt2 = 1.0 / (empv2 * empv2);
0464       double mpv = (mpv1 * wt1 + mpv2 * wt2) / (wt1 + wt2);
0465       double empv = 1.0 / sqrt(wt1 + wt2);
0466       if (debug)
0467         std::cout << "Period " << period << " MPV1 " << mpv1 << " +- " << empv1 << " MPV2 " << mpv2 << " +- " << empv2
0468                   << " --> MPV = " << mpv << " +- " << empv << std::endl;
0469       mpvs[period] = std::pair<double, double>(mpv, empv);
0470     }
0471   }
0472   rType results =
0473       drawMPV(name, eta, 0, depth, nmax, mpvs, first, drawleg, iyear, lumi, ener, normalize, save, npar, debug);
0474   return results;
0475 }
0476 
0477 void plotMPVs(const char* infile,
0478               int phi,
0479               int nvx = 3,
0480               int normalize = 0,
0481               unsigned int first = 0,
0482               int iyear = 17,
0483               double lumis = 0.0,
0484               int ener = 13,
0485               bool save = false,
0486               int npar = 1,
0487               bool debug = false) {
0488   const unsigned int nEta = 22;
0489   int ieta[nEta] = {-26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26};
0490   int ndep[nEta] = {7, 6, 6, 6, 6, 6, 6, 6, 5, 3, 1, 1, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7};
0491   int idep[nEta] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 4, 4, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1};
0492 
0493   char name[100];
0494   int iy = iyear % 100;
0495   sprintf(name, "Data%d%d_trunc.txt", phi, iy);
0496   std::ofstream log(name, std::ios_base::app | std::ios_base::out);
0497   int k0 = (iy == 17) ? 10 : 0;
0498   const unsigned int nmax = (iy == 17) ? 13 : 20;
0499   double lumi = ((lumis > 0) ? lumis : ((iy == 17) ? 49.0 : 68.5));
0500 
0501   for (unsigned int k = k0; k < nEta; ++k) {
0502     int eta = ieta[k];
0503     for (int depth = idep[k]; depth <= ndep[k]; ++depth) {
0504       sprintf(name, "mpvE%dF%dD%d", eta, phi, depth);
0505       std::map<int, std::pair<double, double> > mpvs;
0506       readMPVs(infile, eta, phi, depth, nmax, mpvs, debug);
0507       rType rt =
0508           drawMPV(name, eta, phi, depth, nmax, mpvs, first, true, iyear, lumi, ener, normalize, save, npar, debug);
0509       if (rt.pad != nullptr) {
0510         char line[100];
0511         sprintf(line, "%3d %2d %d %d   %8.4f  %8.4f  %8.4f", eta, phi, depth, nvx, rt.mean, rt.error, rt.chisq);
0512         std::cout << line << std::endl;
0513         log << eta << '\t' << phi << '\t' << depth << '\t' << nvx << '\t' << rt.mean << '\t' << rt.error << '\t'
0514             << rt.chisq << std::endl;
0515       }
0516     }
0517   }
0518 }
0519 
0520 rType plotHist(
0521     const char* infile, double frac, int eta = 25, int depth = 1, int phi = 0, int nvx = 0, bool debug = true) {
0522   TFile* file = new TFile(infile);
0523   rType rt;
0524   gStyle->SetCanvasBorderMode(0);
0525   gStyle->SetCanvasColor(kWhite);
0526   gStyle->SetPadColor(kWhite);
0527   gStyle->SetFillColor(kWhite);
0528   gStyle->SetOptTitle(0);
0529   gStyle->SetOptStat(111110);
0530   if (file != nullptr) {
0531     char name[100];
0532     int fi = (phi == 0) ? 1 : phi;
0533     sprintf(name, "ChrgE%dF%dD%dV%dP0", eta, fi, depth - 1, nvx);
0534     TH1D* hist = (TH1D*)file->FindObjectAny(name);
0535     if (hist != nullptr) {
0536       sprintf(name, "c_E%dF%dD%dV%d", eta, fi, depth - 1, nvx);
0537       rt.pad = new TCanvas(name, name, 500, 500);
0538       rt.pad->SetRightMargin(0.10);
0539       rt.pad->SetTopMargin(0.10);
0540       hist->GetXaxis()->SetTitleSize(0.04);
0541       hist->GetXaxis()->SetTitle("Charge (fC)");
0542       hist->GetYaxis()->SetTitle("Tracks");
0543       hist->GetYaxis()->SetLabelOffset(0.005);
0544       hist->GetYaxis()->SetTitleSize(0.04);
0545       hist->GetYaxis()->SetLabelSize(0.035);
0546       hist->GetYaxis()->SetTitleOffset(1.10);
0547       hist->SetMarkerStyle(20);
0548       hist->SetMarkerColor(2);
0549       hist->SetLineColor(2);
0550       hist->Draw();
0551       rt.pad->Update();
0552       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0553       if (st1 != nullptr) {
0554         st1->SetY1NDC(0.75);
0555         st1->SetY2NDC(0.90);
0556         st1->SetX1NDC(0.65);
0557         st1->SetX2NDC(0.90);
0558       }
0559       TPaveText* txt1 = new TPaveText(0.60, 0.91, 0.90, 0.96, "blNDC");
0560       txt1->SetFillColor(0);
0561       sprintf(name, "#eta = %d, #phi = %d, depth = %d nvx = %d", eta, phi, depth, nvx);
0562       txt1->AddText(name);
0563       txt1->Draw("same");
0564       rt.pad->Modified();
0565       rt.pad->Update();
0566       int nbin = hist->GetNbinsX();
0567       int bins(0);
0568       double total(0), mom1(0), mom2(0);
0569       for (int k = 1; k < nbin; ++k) {
0570         double xx = hist->GetBinLowEdge(k) + 0.5 * hist->GetBinWidth(k);
0571         double en = hist->GetBinContent(k);
0572         total += en;
0573         mom1 += xx * en;
0574         mom2 += (xx * xx * en);
0575       }
0576       mom1 /= total;
0577       mom2 /= total;
0578       double err1 = sqrt((mom2 - mom1 * mom1) / total);
0579       double total0(0), mom10(0), mom20(0);
0580       for (int k = 1; k < nbin; ++k) {
0581         double xx = hist->GetBinLowEdge(k) + 0.5 * hist->GetBinWidth(k);
0582         double en = hist->GetBinContent(k);
0583         if (total0 < frac * total) {
0584           bins = k;
0585           total0 += en;
0586           mom10 += xx * en;
0587           mom20 += (xx * xx * en);
0588         }
0589       }
0590       mom10 /= total0;
0591       mom20 /= total0;
0592       rt.bin = bins;
0593       rt.mean = mom10;
0594       rt.error = sqrt((mom20 - mom10 * mom10) / total0);
0595       if (debug)
0596         std::cout << "Eta " << eta << " phi " << phi << " depth " << depth << " nvx " << nvx << " Mean = " << mom1
0597                   << " +- " << err1 << " " << frac * 100 << "% Truncated Mean = " << rt.mean << " +- " << rt.error
0598                   << "\n";
0599       sprintf(name, "%4.2f truncated mean = %8.2f +- %6.2f", frac, rt.mean, rt.error);
0600       TPaveText* txt2 = new TPaveText(0.10, 0.91, 0.55, 0.96, "blNDC");
0601       txt2->SetFillColor(0);
0602       txt2->AddText(name);
0603       txt2->Draw("same");
0604       TArrow* arrow = new TArrow(rt.bin, 0.032, rt.bin, 6.282, 0.05, "<");
0605       arrow->SetFillColor(kBlack);
0606       arrow->SetFillStyle(1001);
0607       arrow->SetLineColor(kBlack);
0608       arrow->SetLineStyle(2);
0609       arrow->SetLineWidth(2);
0610       arrow->Draw();
0611       rt.pad->Modified();
0612       rt.pad->Update();
0613     }
0614   }
0615   return rt;
0616 }
0617 
0618 void getTruncateMean(const char* infile, double frac, std::string type, int nvx = 2, bool save = false) {
0619   int ieta[22] = {-26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26};
0620   int ndep[22] = {7, 6, 6, 6, 6, 6, 6, 6, 5, 3, 1, 1, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7};
0621   int idep[22] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 4, 4, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1};
0622 
0623   for (int k = 0; k < 22; ++k) {
0624     for (int depth = idep[k]; depth <= ndep[k]; ++depth) {
0625       rType rt = plotHist(infile, frac, ieta[k], depth, 0, nvx, false);
0626       if (rt.pad != nullptr) {
0627         char line[100];
0628         sprintf(
0629             line, "%s %2d   0 %d %d   %8.4f  %8.4f  %d", type.c_str(), ieta[k], depth, nvx, rt.mean, rt.error, rt.bin);
0630         std::cout << line << std::endl;
0631         if (save) {
0632           char name[100];
0633           sprintf(name, "%s_comb.pdf", rt.pad->GetName());
0634           rt.pad->Print(name);
0635         }
0636       }
0637     }
0638   }
0639 }
0640 
0641 void getTruncatedMeanX(
0642     const char* infile, std::string type, int nvx = 2, int year = 18, int phi = 0, bool save = false) {
0643   const unsigned int nEta = 22;
0644   int ieta[nEta] = {-26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26};
0645   int ndep[nEta] = {7, 6, 6, 6, 6, 6, 6, 6, 6, 3, 4, 4, 3, 6, 6, 6, 6, 6, 6, 6, 6, 7};
0646   int idep[nEta] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 4, 4, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1};
0647   int ioff[nEta] = {0, 7, 13, 19, 25, 31, 37, 43, 49, 54, 56, 57, 58, 60, 65, 71, 77, 83, 89, 95, 101, 107};
0648   const unsigned int maxIndx = 114;
0649   double frac18[maxIndx] = {
0650       0.40, 0.55, 0.50, 0.55, 0.60, 0.60, 0.60, 0.40, 0.55, 0.60, 0.60, 0.60, 0.60, 0.40, 0.55, 0.60, 0.60, 0.60, 0.60,
0651       0.45, 0.55, 0.60, 0.60, 0.60, 0.60, 0.55, 0.60, 0.60, 0.60, 0.60, 0.60, 0.55, 0.60, 0.60, 0.60, 0.60, 0.65, 0.60,
0652       0.60, 0.60, 0.60, 0.60, 0.65, 0.60, 0.60, 0.60, 0.60, 0.60, 0.65, 0.60, 0.60, 0.60, 0.60, 0.60, 0.65, 0.60, 0.80,
0653       0.80, 0.65, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.65, 0.60, 0.60, 0.60, 0.60, 0.60,
0654       0.65, 0.60, 0.60, 0.60, 0.60, 0.60, 0.65, 0.55, 0.60, 0.60, 0.60, 0.60, 0.65, 0.55, 0.60, 0.60, 0.60, 0.60, 0.60,
0655       0.45, 0.55, 0.60, 0.60, 0.60, 0.60, 0.40, 0.55, 0.60, 0.60, 0.60, 0.60, 0.40, 0.55, 0.50, 0.55, 0.60, 0.60, 0.60};
0656   double frac17[maxIndx] = {
0657       0.70, 0.65, 0.55, 0.60, 0.60, 0.65, 0.65, 0.67, 0.55, 0.55, 0.57, 0.55, 0.62, 0.65, 0.57, 0.60, 0.55, 0.70, 0.65,
0658       0.65, 0.60, 0.55, 0.65, 0.75, 0.60, 0.60, 0.60, 0.65, 0.62, 0.75, 0.65, 0.70, 0.67, 0.60, 0.65, 0.55, 0.65, 0.60,
0659       0.55, 0.75, 0.65, 0.75, 0.70, 0.60, 0.55, 0.75, 0.67, 0.65, 0.60, 0.60, 0.45, 0.50, 0.52, 0.70, 0.65, 0.60, 0.80,
0660       0.80, 0.60, 0.65, 0.60, 0.55, 0.75, 0.67, 0.65, 0.60, 0.55, 0.75, 0.65, 0.75, 0.70, 0.70, 0.67, 0.60, 0.65, 0.55,
0661       0.65, 0.60, 0.60, 0.65, 0.62, 0.75, 0.65, 0.65, 0.60, 0.55, 0.65, 0.75, 0.60, 0.65, 0.60, 0.55, 0.65, 0.75, 0.60,
0662       0.65, 0.57, 0.60, 0.55, 0.70, 0.65, 0.67, 0.55, 0.55, 0.57, 0.55, 0.62, 0.70, 0.65, 0.55, 0.60, 0.60, 0.65, 0.65};
0663   double frac63[maxIndx] = {
0664       0.63, 0.60, 0.57, 0.56, 0.60, 0.61, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.61, 0.58, 0.60, 0.60, 0.60, 0.62, 0.63,
0665       0.57, 0.60, 0.60, 0.60, 0.65, 0.61, 0.57, 0.60, 0.60, 0.60, 0.62, 0.64, 0.60, 0.60, 0.60, 0.60, 0.63, 0.65, 0.65,
0666       0.60, 0.63, 0.60, 0.64, 0.65, 0.60, 0.60, 0.63, 0.62, 0.62, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.80,
0667       0.80, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.60, 0.63, 0.62, 0.62, 0.65, 0.65, 0.60, 0.63, 0.60, 0.64,
0668       0.65, 0.60, 0.60, 0.60, 0.60, 0.63, 0.65, 0.57, 0.60, 0.60, 0.60, 0.62, 0.64, 0.57, 0.60, 0.60, 0.60, 0.65, 0.61,
0669       0.58, 0.60, 0.60, 0.60, 0.62, 0.63, 0.60, 0.60, 0.60, 0.60, 0.60, 0.61, 0.63, 0.60, 0.57, 0.56, 0.60, 0.61, 0.60};
0670 
0671   char name[100];
0672   sprintf(name, "Data%d_trunc.txt", year);
0673   std::ofstream log(name, std::ios_base::app | std::ios_base::out);
0674   int k0 = (year == 17) ? (nEta / 2) : 0;
0675   for (unsigned int k = k0; k < nEta; ++k) {
0676     for (int depth = idep[k]; depth <= ndep[k]; ++depth) {
0677       int indx = ioff[k] + depth - idep[k];
0678       double frac = ((year == 17) ? ((phi == 0) ? frac17[indx] : frac63[indx]) : frac18[indx]);
0679       if (debugMode)
0680         std::cout << k << " E " << ieta[k] << " D " << depth << " I " << indx << " R " << frac << std::endl;
0681       rType rt = plotHist(infile, frac, ieta[k], depth, phi, nvx, false);
0682       if (rt.pad != nullptr) {
0683         char line[100];
0684         sprintf(line,
0685                 "%s %2d   %d %d %d   %8.4f  %8.4f  %d",
0686                 type.c_str(),
0687                 ieta[k],
0688                 phi,
0689                 depth,
0690                 nvx,
0691                 rt.mean,
0692                 rt.error,
0693                 rt.bin);
0694         std::cout << line << std::endl;
0695         log << type << '\t' << ieta[k] << '\t' << phi << '\t' << depth << '\t' << nvx << '\t' << rt.mean << '\t'
0696             << rt.error << std::endl;
0697         if (save) {
0698           sprintf(name, "%s_comb.pdf", rt.pad->GetName());
0699           rt.pad->Print(name);
0700         }
0701       }
0702     }
0703   }
0704 }
0705 
0706 std::map<int, std::pair<double, double> > readSlopes(
0707     const char* infile, const int typ, const int phi, const int dep, bool debug) {
0708   std::map<int, std::pair<double, double> > slopes;
0709   std::ifstream fInput(infile);
0710   if (!fInput.good()) {
0711     std::cout << "Cannot open file " << infile << std::endl;
0712   } else {
0713     char buffer[1024];
0714     unsigned int all(0), good(0), select(0);
0715     while (fInput.getline(buffer, 1024)) {
0716       ++all;
0717       if (buffer[0] == '#')
0718         continue;  //ignore comment
0719       std::vector<std::string> items = splitString(std::string(buffer));
0720       if (items.size() < 6) {
0721         if (debug)
0722           std::cout << "Ignore  line: " << buffer << std::endl;
0723       } else {
0724         ++good;
0725         int type = std::atoi(items[0].c_str());
0726         int ieta = std::atoi(items[1].c_str());
0727         int iphi = std::atoi(items[2].c_str());
0728         int depth = std::atoi(items[3].c_str());
0729         double val = std::atof(items[4].c_str());
0730         double err = std::atof(items[5].c_str());
0731         if ((type == typ) && (iphi == phi) && (depth == dep)) {
0732           ++select;
0733           slopes[ieta] = std::pair<double, double>(val, err);
0734         }
0735       }
0736     }
0737     fInput.close();
0738     std::cout << "Reads total of " << all << " and " << good << " good and " << select << " selected records from "
0739               << infile << std::endl;
0740   }
0741   if (debug) {
0742     unsigned int k(0);
0743     for (std::map<int, std::pair<double, double> >::const_iterator itr = slopes.begin(); itr != slopes.end();
0744          ++itr, ++k) {
0745       std::cout << "[" << k << "] ieta " << itr->first << " Slope " << (itr->second).first << " +- "
0746                 << (itr->second).second << std::endl;
0747     }
0748   }
0749   return slopes;
0750 }
0751 
0752 void drawSlope(const char* infile,
0753                int type,
0754                int phi,
0755                int depth,
0756                bool drawleg = true,
0757                int iyear = 2018,
0758                double lumis = 0,
0759                int ener = 13,
0760                bool save = false,
0761                bool debug = false) {
0762   std::map<int, std::pair<double, double> > slopes = readSlopes(infile, type, phi, depth, debug);
0763   if (slopes.size() > 0) {
0764     gStyle->SetCanvasBorderMode(0);
0765     gStyle->SetCanvasColor(kWhite);
0766     gStyle->SetPadColor(kWhite);
0767     gStyle->SetFillColor(kWhite);
0768     gStyle->SetOptTitle(0);
0769     gStyle->SetOptStat(0);
0770     gStyle->SetOptFit(0);
0771     char name[100], cname[100];
0772     std::string method[2] = {"Fit", "TM"};
0773     std::string methods[2] = {"Landau Fit", "Truncated Mean"};
0774     int iy = iyear % 100;
0775     sprintf(cname, "%s%d", method[type].c_str(), iy);
0776     TCanvas* pad = new TCanvas(cname, cname);
0777     pad->SetLeftMargin(0.12);
0778     pad->SetRightMargin(0.10);
0779     pad->SetTopMargin(0.10);
0780     TLegend* legend = new TLegend(0.2206304, 0.8259023, 0.760745, 0.8768577, NULL, "brNDC");
0781     legend->SetFillColor(kWhite);
0782     legend->SetBorderSize(1);
0783     unsigned int nmax = slopes.size();
0784     int np(0);
0785     double xv[nmax], dxv[nmax], yv[nmax], dyv[nmax];
0786     for (std::map<int, std::pair<double, double> >::const_iterator itr = slopes.begin(); itr != slopes.end(); ++itr) {
0787       xv[np] = itr->first;
0788       dxv[np] = 0;
0789       yv[np] = (itr->second).first;
0790       dyv[np] = (itr->second).second;
0791       ++np;
0792     }
0793     TGraphErrors* g = new TGraphErrors(np, xv, yv, dxv, dyv);
0794     g->SetMarkerStyle(8);
0795     g->SetMarkerColor(kBlue);
0796     g->GetXaxis()->SetRangeUser(-30.0, 30.0);
0797     g->GetYaxis()->SetRangeUser(-0.0005, 0.0075);
0798     g->SetMarkerSize(1.5);
0799     g->GetYaxis()->SetTitle("Slope");
0800     g->GetXaxis()->SetTitle("i#eta");
0801     g->GetYaxis()->SetTitleSize(0.04);
0802     g->GetXaxis()->SetTitleSize(0.04);
0803     g->GetXaxis()->SetLabelSize(0.04);
0804     g->GetYaxis()->SetLabelSize(0.04);
0805     g->GetXaxis()->SetTitleOffset(1.0);
0806     g->GetYaxis()->SetTitleOffset(1.5);
0807     g->SetTitle("");
0808     g->Draw("AP");
0809 
0810     if (phi > 0) {
0811       sprintf(name, "%s Method (i#phi %d, depth %d)", methods[type].c_str(), phi, depth);
0812     } else {
0813       if (iy == 17)
0814         sprintf(name, "%s Method (HEP17 depth %d)", methods[type].c_str(), depth);
0815       else
0816         sprintf(name, "%s Method (All #phi's, depth %d)", methods[type].c_str(), depth);
0817     }
0818     legend->AddEntry(g, name, "lp");
0819     if (drawleg) {
0820       legend->Draw("same");
0821       pad->Update();
0822     }
0823     if (iyear > 1000) {
0824       char txt[100];
0825       TPaveText* txt1 = new TPaveText(0.60, 0.91, 0.90, 0.96, "blNDC");
0826       txt1->SetFillColor(0);
0827       sprintf(txt, "%d, %d TeV %5.1f fb^{-1}", iyear, ener, lumis);
0828       txt1->AddText(txt);
0829       txt1->Draw("same");
0830       TPaveText* txt2 = new TPaveText(0.10, 0.91, 0.18, 0.96, "blNDC");
0831       txt2->SetFillColor(0);
0832       sprintf(txt, "CMS");
0833       txt2->AddText(txt);
0834       txt2->Draw("same");
0835     }
0836     pad->Modified();
0837     pad->Update();
0838     if (save) {
0839       sprintf(name, "%s_depth%d.pdf", pad->GetName(), depth);
0840       pad->Print(name);
0841     }
0842   }
0843 }
0844 
0845 TCanvas* plotLightLoss(std::string infile = "mu_HE_insitu_2018.txt",
0846                        double lumis = 66.5,
0847                        int iyear = 2018,
0848                        int ener = 13,
0849                        bool save = false,
0850                        bool debug = true) {
0851   std::map<std::pair<int, int>, double> losses;
0852   std::map<std::pair<int, int>, double>::const_iterator itr;
0853   std::ifstream fInput(infile.c_str());
0854   if (!fInput.good()) {
0855     std::cout << "Cannot open file " << infile << std::endl;
0856   } else {
0857     char buffer[1024];
0858     unsigned int all(0), good(0);
0859     while (fInput.getline(buffer, 1024)) {
0860       ++all;
0861       if (buffer[0] == '#')
0862         continue;  //ignore comment
0863       std::vector<std::string> items = splitString(std::string(buffer));
0864       if (items.size() < 5) {
0865         if (debug)
0866           std::cout << "Ignore  line: " << buffer << std::endl;
0867       } else {
0868         ++good;
0869         int ieta = std::atoi(items[2].c_str());
0870         int depth = std::atoi(items[1].c_str());
0871         double val = std::atof(items[3].c_str());
0872         double loss = exp(-val * lumis);
0873         losses[std::pair<int, int>(ieta, depth)] = loss;
0874       }
0875     }
0876     fInput.close();
0877     std::cout << "Reads total of " << all << " and " << good << " good records "
0878               << "from " << infile << std::endl;
0879     if (debug) {
0880       unsigned int k(0);
0881       for (itr = losses.begin(); itr != losses.end(); ++itr, ++k)
0882         std::cout << "[" << k << "] eta|depth " << (itr->first).first << "|" << (itr->first).second << " Loss "
0883                   << (itr->second) << std::endl;
0884     }
0885   }
0886   TCanvas* pad;
0887   if (losses.size() > 0) {
0888     gStyle->SetCanvasBorderMode(0);
0889     gStyle->SetCanvasColor(kWhite);
0890     gStyle->SetPadColor(kWhite);
0891     gStyle->SetFillColor(kWhite);
0892     gStyle->SetOptTitle(0);
0893     gStyle->SetOptStat(0);
0894     gStyle->SetOptFit(0);
0895     gStyle->SetPalette(1);
0896     char text[100], cname[100];
0897     int iy = iyear % 100;
0898     sprintf(cname, "light%d", iy);
0899     pad = new TCanvas(cname, cname, 600, 700);
0900     pad->SetLeftMargin(0.10);
0901     pad->SetRightMargin(0.14);
0902     pad->SetTopMargin(0.10);
0903     const int neta = 15;
0904     TH2D* h = new TH2D("cname", "cname", 8, 0, 8, neta, 0, neta);
0905     TGraph2D* dt = new TGraph2D();
0906     dt->SetNpx(8);
0907     dt->SetNpy(15);
0908     h->GetXaxis()->SetTitle("Detector depth number");
0909     h->GetYaxis()->SetTitle("Tower");
0910     h->GetZaxis()->SetRangeUser(0.75, 1.25);
0911     std::map<std::pair<int, int>, double>::const_iterator itr;
0912 
0913     int ieta, depth;
0914     for (itr = losses.begin(); itr != losses.end(); ++itr) {
0915       ieta = (itr->first).first;
0916       depth = (itr->first).second;
0917       h->Fill((depth + 0.5), (30.5 - ieta), (itr->second));
0918     }
0919     h->Draw("colz");
0920 
0921     gPad->Update();
0922     TAxis* a = h->GetYaxis();
0923     a->SetNdivisions(20);
0924     std::string val[neta] = {"30", "29", "28", "27", "26", "25", "24", "23", "22", "21", "20", "19", "18", "17", "16"};
0925     for (int i = 0; i < neta; i++)
0926       a->ChangeLabel((i + 1), -1, -1, -1, -1, -1, val[i]);
0927     a->CenterLabels(kTRUE);
0928     gPad->Modified();
0929     gPad->Update();
0930     TPaveText* txt1 = new TPaveText(0.60, 0.91, 0.90, 0.97, "blNDC");
0931     txt1->SetFillColor(0);
0932     sprintf(text, "%d, %d TeV %5.1f fb^{-1}", iyear, ener, lumis);
0933     txt1->AddText(text);
0934     txt1->Draw("same");
0935     TPaveText* txt2 = new TPaveText(0.11, 0.91, 0.18, 0.97, "blNDC");
0936     txt2->SetFillColor(0);
0937     sprintf(text, "CMS");
0938     txt2->AddText(text);
0939     txt2->Draw("same");
0940     pad->Modified();
0941     pad->Update();
0942     if (save) {
0943       sprintf(cname, "%s.pdf", pad->GetName());
0944       pad->Print(cname);
0945     }
0946   }
0947   return pad;
0948 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!