OfflineValidation/macros/trackSplitPlot.C

0001 /***********************************
0002 Table Of Contents
0003 0. Track Split Plot
0004 1. Misalignment Dependence
0005 2. Make Plots
0006 3. Axis Label
0007 4. Axis Limits
0008 5. Place Legend
0009 ***********************************/
0010
0011 using namespace std;
0012 #include "trackSplitPlot.h"
0013 #include "Alignment/OfflineValidation/interface/TkAlStyle.h"
0014
0015 //===================
0016 //0. Track Split Plot
0017 //===================
0018
0019 TCanvas *trackSplitPlot(Int_t nFiles,
0020                         TString *files,
0021                         TString *names,
0022                         TString xvar,
0023                         TString yvar,
0024                         Bool_t relative,
0025                         Bool_t resolution,
0026                         Bool_t pull,
0027                         TString saveas,
0028                         ostream &summaryfile) {
0029   if (TkAlStyle::status() == NO_STATUS)
0030     TkAlStyle::set(INTERNAL);
0031   TString legendOptions = TkAlStyle::legendoptions;
0032   legendOptions.ReplaceAll("all", "meanerror,rmserror").ToLower();
0033   if (outliercut < 0)
0034     outliercut = -1;
0035   gStyle->SetMarkerSize(1.5);
0036   setupcolors();
0037   stufftodelete->SetOwner(true);
0038   cout << xvar << " " << yvar << endl;
0039   if (xvar == "" && yvar == "")
0040     return nullptr;
0041
0042   PlotType type;
0043   if (xvar == "")
0044     type = Histogram;
0045   else if (yvar == "")
0046     type = OrgHistogram;
0047   else if (resolution)
0048     type = Resolution;
0049   else if (nFiles < 1)
0050     type = ScatterPlot;
0051   else
0052     type = Profile;
0053   if (nFiles < 1)
0054     nFiles = 1;
0055
0056   const Int_t n = nFiles;
0057
0058   vector<TH1 *> p;
0059   Int_t lengths[n];
0060
0061   stringstream sx, sy, srel, ssigma1, ssigma2, ssigmaorg;
0062
0063   sx << xvar << "_org";
0064   TString xvariable = sx.str();
0065   TString xvariable2 = "";
0066   if (xvar == "runNumber")
0067     xvariable = "runNumber";
0068   if (xvar.BeginsWith("nHits")) {
0069     xvariable = xvar;
0070     xvariable2 = xvar;
0071     xvariable.Append("1_spl");
0072     xvariable2.Append("2_spl");
0073   }
0074
0075   sy << "Delta_" << yvar;
0076   TString yvariable = sy.str();
0077
0078   TString relvariable = "1";
0079   if (relative) {
0080     srel << yvar << "_org";
0081     relvariable = srel.str();
0082   }
0083
0084   TString sigma1variable = "", sigma2variable = "";
0085   if (pull) {
0086     ssigma1 << yvar << "1Err_spl";
0087     ssigma2 << yvar << "2Err_spl";
0088   }
0089   sigma1variable = ssigma1.str();
0090   sigma2variable = ssigma2.str();
0091
0092   TString sigmaorgvariable = "";
0093   if (pull && relative)
0094     ssigmaorg << yvar << "Err_org";
0095   sigmaorgvariable = ssigmaorg.str();
0096
0097   Double_t xmin = -1, xmax = 1, ymin = -1, ymax = 1, xbins = -1, ybins;
0098   if (type == Profile || type == ScatterPlot || type == OrgHistogram || type == Resolution)
0099     axislimits(nFiles, files, xvar, 'x', relative, pull, xmin, xmax, xbins);
0100   if (type == Profile || type == ScatterPlot || type == Histogram || type == Resolution)
0101     axislimits(nFiles, files, yvar, 'y', relative, pull, ymin, ymax, ybins);
0102
0103   std::vector<TString> meansrmss(n);
0104   std::vector<double> means(n);
0105   std::vector<double> rmss(n);
0106   //a file is not "used" if it's MC data and the x variable is run number, or if the filename is blank
0107   std::vector<bool> used(n);
0108
0109   for (Int_t i = 0; i < n; i++) {
0110     stringstream sid;
0111     sid << "p" << i;
0112     TString id = sid.str();
0113
0114     //for a profile or resolution, it fills a histogram, q[j], for each bin, then gets the mean and width from there.
0115     vector<TH1F *> q;
0116
0117     if (type == ScatterPlot)
0118       p.push_back(new TH2F(id, "", xbins, xmin, xmax, ybins, ymin, ymax));
0119     if (type == Histogram)
0120       p.push_back(new TH1F(id, "", ybins, ymin, ymax));
0121     if (type == OrgHistogram)
0122       p.push_back(new TH1F(id, "", xbins, xmin, xmax));
0123     if (type == Resolution || type == Profile) {
0124       p.push_back(new TH1F(id, "", xbins, xmin, xmax));
0125       for (Int_t j = 0; j < xbins; j++) {
0126         stringstream sid2;
0127         sid2 << "q" << i << j;
0128         TString id2 = sid2.str();
0129         q.push_back(new TH1F(id2, "", 1000, ymin * 10, ymax * 10));
0130       }
0131     }
0132
0133     p[i]->SetLineColor(colors[i]);
0134     if (type == Resolution || type == Profile) {
0135       p[i]->SetMarkerStyle(styles[i] / 100);
0136       p[i]->SetMarkerColor(colors[i]);
0137       p[i]->SetLineStyle(styles[i] % 100);
0138     } else {
0139       if (styles[i] >= 100) {
0140         p[i]->SetMarkerStyle(styles[i] / 100);
0141         p[i]->SetMarkerColor(colors[i]);
0142         p[i]->Sumw2();
0143       }
0144       p[i]->SetLineStyle(styles[i] % 100);
0145     }
0146
0147     stufftodelete->Add(p[i]);
0148     p[i]->SetBit(kCanDelete, true);
0149
0150     used[i] = true;
0151     //if it's MC data (run 1), the run number is meaningless
0152     if ((xvar == "runNumber" && findMax(files[i], "runNumber", 'x') < 2) || files[i] == "") {
0153       used[i] = false;
0154       p[i]->SetLineColor(kWhite);
0155       p[i]->SetMarkerColor(kWhite);
0156       for (unsigned int j = 0; j < q.size(); j++)
0157         delete q[j];
0158       continue;
0159     }
0160
0161     TFile *f = TFile::Open(files[i]);
0162     TTree *tree = (TTree *)f->Get("cosmicValidation/splitterTree");
0163     if (tree == nullptr)
0164       tree = (TTree *)f->Get("splitterTree");
0165
0166     lengths[i] = tree->GetEntries();
0167
0168     Double_t x = 0, y = 0, rel = 1, sigma1 = 1;
0169     Double_t sigma2 = 1;  //if !pull, we want to divide by sqrt(2) because we want the error from 1 track
0170     Double_t sigmaorg = 0;
0171
0172     Int_t xint = 0, xint2 = 0;
0173     Int_t runNumber = 0;
0174     double pt1 = 0, maxpt1 = 0;
0175
0176     if (!relative && !pull && (yvar == "dz" || yvar == "dxy"))
0177       rel = 1e-4;  //it's in cm but we want it in um, so divide by 1e-4
0178     if (!relative && !pull && (yvar == "phi" || yvar == "theta" || yvar == "qoverpt"))
0179       rel = 1e-3;  //make the axis labels manageable
0180
0181     tree->SetBranchAddress("runNumber", &runNumber);
0182     if (type == Profile || type == ScatterPlot || type == Resolution || type == OrgHistogram) {
0183       if (xvar == "runNumber")
0184         tree->SetBranchAddress(xvariable, &xint);
0185       else if (xvar.BeginsWith("nHits")) {
0186         tree->SetBranchAddress(xvariable, &xint);
0187         tree->SetBranchAddress(xvariable2, &xint2);
0188       } else
0189         tree->SetBranchAddress(xvariable, &x);
0190     }
0191     if (type == Profile || type == ScatterPlot || type == Resolution || type == Histogram) {
0192       int branchexists = tree->SetBranchAddress(yvariable, &y);
0193       if (branchexists == -5)  //i.e. it doesn't exist
0194       {
0195         yvariable.ReplaceAll("Delta_", "d");
0196         yvariable.Append("_spl");
0197         tree->SetBranchAddress(yvariable, &y);
0198       }
0199     }
0200     if (relative && xvar != yvar)  //if xvar == yvar, setting the branch here will undo setting it to x 2 lines earlier
0201       tree->SetBranchAddress(relvariable, &rel);  //setting the value of rel is then taken care of later: rel = x
0202     if (pull) {
0203       tree->SetBranchAddress(sigma1variable, &sigma1);
0204       tree->SetBranchAddress(sigma2variable, &sigma2);
0205     }
0206     if (relative && pull)
0207       tree->SetBranchAddress(sigmaorgvariable, &sigmaorg);
0208     if (xvar == "pt" || yvar == "pt" || xvar == "qoverpt" || yvar == "qoverpt") {
0209       tree->SetBranchAddress("pt1_spl", &pt1);
0210     } else {
0211       maxpt1 = 999;
0212     }
0213
0214     Int_t notincluded = 0;  //this counts the number that aren't in the right run range.
0215                             //it's subtracted from lengths[i] in order to normalize the histograms
0216
0217     for (Int_t j = 0; j < lengths[i]; j++) {
0218       tree->GetEntry(j);
0219       if (xvar == "runNumber" || xvar.BeginsWith("nHits"))
0220         x = xint;
0221       if (xvar == "runNumber")
0222         runNumber = x;
0223       if (yvar == "phi" && y >= pi)
0224         y -= 2 * pi;
0225       if (yvar == "phi" && y <= -pi)
0226         y += 2 * pi;
0227       if ((runNumber < minrun && runNumber > 1) ||
0228           (runNumber > maxrun && maxrun > 0))  //minrun and maxrun are global variables.
0229       {
0230         notincluded++;
0231         continue;
0232       }
0233       if (relative && xvar == yvar)
0234         rel = x;
0235       Double_t error = 0;
0236       if (relative && pull)
0237         error = sqrt((sigma1 / rel) * (sigma1 / rel) + (sigma2 / rel) * (sigma2 / rel) +
0238                      (sigmaorg * y / (rel * rel)) * (sigmaorg * x / (rel * rel)));
0239       else
0240         error = sqrt(sigma1 * sigma1 +
0241                      sigma2 * sigma2);  // = sqrt(2) if !pull; this divides by sqrt(2) to get the error in 1 track
0242       y /= (rel * error);
0243
0244       if (pt1 > maxpt1)
0245         maxpt1 = pt1;
0246
0247       if (ymin <= y && y < ymax && xmin <= x && x < xmax) {
0248         if (type == Histogram)
0249           p[i]->Fill(y);
0250         if (type == ScatterPlot)
0251           p[i]->Fill(x, y);
0252         if (type == Resolution || type == Profile) {
0253           int which = (p[i]->Fill(x, 0)) - 1;
0254           //get which q[j] by filling p[i] with nothing.  (TH1F::Fill returns the bin number)
0255           //p[i]'s actual contents are set later.
0256           if (which >= 0 && (unsigned)which < q.size())
0257             q[which]->Fill(y);
0258         }
0259         if (type == OrgHistogram)
0260           p[i]->Fill(x);
0261       }
0262
0263       if (xvar.BeginsWith("nHits")) {
0264         x = xint2;
0265         if (ymin <= y && y < ymax && xmin <= x && x < xmax) {
0266           if (type == Histogram)
0267             p[i]->Fill(y);
0268           if (type == ScatterPlot)
0269             p[i]->Fill(x, y);
0270           if (type == Resolution || type == Profile) {
0271             int which = (p[i]->Fill(x, 0)) - 1;
0272             if (which >= 0)
0273               q[which]->Fill(y);  //get which q[j] by filling p[i] (with nothing), which returns the bin number
0274           }
0275           if (type == OrgHistogram)
0276             p[i]->Fill(x);
0277         }
0278       }
0279
0280       if (lengths[i] < 10 ? true
0281                           : (((j + 1) / (int)(pow(10, (int)(log10(lengths[i])) - 1))) *
0282                                      (int)(pow(10, (int)(log10(lengths[i])) - 1)) ==
0283                                  j + 1 ||
0284                              j + 1 == lengths[i]))
0285       //print when j+1 is a multiple of 10^x, where 10^x has 1 less digit than lengths[i]
0286       // and when it's finished
0287       //For example, if lengths[i] = 123456, it will print this when j+1 = 10000, 20000, ..., 120000, 123456
0288       //So it will print between 10 and 100 times: 10 when lengths[i] = 10^x and 100 when lengths[i] = 10^x - 1
0289       {
0290         cout << j + 1 << "/" << lengths[i] << ": ";
0291         if (type == Profile || type == ScatterPlot || type == Resolution)
0292           cout << x << ", " << y << endl;
0293         if (type == OrgHistogram)
0294           cout << x << endl;
0295         if (type == Histogram)
0296           cout << y << endl;
0297       }
0298     }
0299     lengths[i] -= notincluded;
0300
0301     if (maxpt1 < 6) {  //0T
0302       used[i] = false;
0303       p[i]->SetLineColor(kWhite);
0304       p[i]->SetMarkerColor(kWhite);
0305       for (unsigned int j = 0; j < q.size(); j++)
0306         delete q[j];
0307       continue;
0308     }
0309
0310     meansrmss[i] = "";
0311     if (type == Histogram || type == OrgHistogram) {
0312       stringstream meanrms;
0313       meanrms.precision(3);
0314
0315       double average = -1e99;
0316       double rms = -1e99;
0317
0318       TString var = (type == Histogram ? yvar : xvar);
0319       char axis = (type == Histogram ? 'y' : 'x');
0320       TString varunits = "";
0321       if (!relative && !pull)
0322         varunits = units(var, axis);
0323       if (legendOptions.Contains("mean")) {
0324         if (outliercut < 0)
0325           average = p[i]->GetMean();
0326         else
0327           average = findAverage(files[i], var, axis, relative, pull);
0328         cout << "Average = " << average;
0329         meanrms << "#mu = " << average;
0330         means[i] = average;
0331         if (legendOptions.Contains("meanerror")) {
0332           if (outliercut < 0)
0333             rms = p[i]->GetRMS();
0334           else
0335             rms = findRMS(files[i], var, axis, relative, pull);
0336           meanrms << " #pm " << rms / TMath::Sqrt(lengths[i] * abs(outliercut));
0337           cout << " +/- " << rms / TMath::Sqrt(lengths[i] * abs(outliercut));
0338         }
0339         if (varunits != "") {
0340           meanrms << " " << varunits;
0341           cout << " " << varunits;
0342         }
0343         cout << endl;
0344         if (legendOptions.Contains("rms"))
0345           meanrms << ", ";
0346       }
0347       if (legendOptions.Contains("rms")) {
0348         if (rms < -1e98) {
0349           if (outliercut < 0)
0350             rms = p[i]->GetRMS();
0351           else
0352             rms = findRMS(files[i], var, axis, relative, pull);
0353         }
0354         cout << "RMS     = " << rms;
0355         meanrms << "rms = " << rms;
0356         rmss[i] = rms;
0357         if (legendOptions.Contains("rmserror")) {
0358           //https://root.cern.ch/root/html/src/TH1.cxx.html#7076
0359           meanrms << " #pm " << rms / TMath::Sqrt(2 * lengths[i] * abs(outliercut));
0360           cout << " +/- " << rms / TMath::Sqrt(2 * lengths[i] * abs(outliercut));
0361         }
0362         if (varunits != "") {
0363           meanrms << " " << varunits;
0364           cout << " " << varunits;
0365         }
0366         cout << endl;
0367       }
0368       meansrmss[i] = meanrms.str();
0369     }
0370
0371     if (type == Resolution) {
0372       for (Int_t j = 0; j < xbins; j++) {
0373         p[i]->SetBinContent(j + 1, q[j]->GetRMS());
0374         p[i]->SetBinError(j + 1, q[j]->GetRMSError());
0375         delete q[j];
0376       }
0377     }
0378
0379     if (type == Profile) {
0380       for (Int_t j = 0; j < xbins; j++) {
0381         p[i]->SetBinContent(j + 1, q[j]->GetMean());
0382         p[i]->SetBinError(j + 1, q[j]->GetMeanError());
0383         delete q[j];
0384       }
0385     }
0386
0387     setAxisLabels(p[i], type, xvar, yvar, relative, pull);
0388   }
0389
0390   if (type == Histogram && !pull && any_of(begin(used), end(used), ::identity<bool>)) {
0391     if (legendOptions.Contains("mean")) {
0392       summaryfile << "   mu_Delta" << yvar;
0393       if (relative)
0394         summaryfile << "/" << yvar;
0395       if (pull)
0396         summaryfile << "_pull";
0397       if (!pull && !relative && plainunits(yvar, 'y') != "")
0398         summaryfile << " (" << plainunits(yvar, 'y') << ")";
0399       summaryfile << "\t"
0400                   << "latexname=$\\mu_{" << latexlabel(yvar, 'y', relative, resolution, pull) << "}$";
0401       if (!pull && !relative && plainunits(yvar, 'y') != "")
0402         summaryfile << " (" << latexunits(yvar, 'y') << ")";
0403       summaryfile << "\t"
0404                   << "format={:.3g}\t"
0405                   << "latexformat=${:.3g}$";
0406       for (int i = 0; i < n; i++) {
0407         if (used[i]) {
0408           summaryfile << "\t" << means[i];
0409         } else {
0410           summaryfile << "\t" << nan("");
0411         }
0412       }
0413       summaryfile << "\n";
0414     }
0415     if (legendOptions.Contains("rms")) {
0416       summaryfile << "sigma_Delta" << yvar;
0417       if (relative)
0418         summaryfile << "/" << yvar;
0419       if (pull)
0420         summaryfile << "_pull";
0421       if (!pull && !relative && plainunits(yvar, 'y') != "")
0422         summaryfile << " (" << plainunits(yvar, 'y') << ")";
0423       summaryfile << "\t"
0424                   << "latexname=$\\sigma_{" << latexlabel(yvar, 'y', relative, resolution, pull) << "}$";
0425       if (!pull && !relative && latexunits(yvar, 'y') != "")
0426         summaryfile << " (" << latexunits(yvar, 'y') << ")";
0427       summaryfile << "\t"
0428                   << "format={:.3g}\t"
0429                   << "latexformat=${:.3g}$";
0430       for (int i = 0; i < n; i++) {
0431         if (used[i]) {
0432           summaryfile << "\t" << rmss[i];
0433         } else {
0434           summaryfile << "\t" << nan("");
0435         }
0436       }
0437       summaryfile << "\n";
0438     }
0439   }
0440
0441   TH1 *firstp = nullptr;
0442   for (int i = 0; i < n; i++) {
0443     if (used[i]) {
0444       firstp = p[i];
0445       break;
0446     }
0447   }
0448   if (firstp == nullptr) {
0449     stufftodelete->Clear();
0450     return nullptr;
0451   }
0452
0453   TCanvas *c1 = TCanvas::MakeDefCanvas();
0454
0455   TH1 *maxp = firstp;
0456   if (type == ScatterPlot)
0457     firstp->Draw("COLZ");
0458   else if (type == Resolution || type == Profile) {
0459     vector<TGraphErrors *> g;
0460     TMultiGraph *list = new TMultiGraph();
0461     for (Int_t i = 0, ii = 0; i < n; i++, ii++) {
0462       if (!used[i]) {
0463         ii--;
0464         continue;
0465       }
0466       g.push_back(new TGraphErrors(p[i]));
0467       for (Int_t j = 0; j < g[ii]->GetN(); j++) {
0468         if (g[ii]->GetY()[j] == 0 && g[ii]->GetEY()[j] == 0) {
0469           g[ii]->RemovePoint(j);
0470           j--;
0471         }
0472       }
0473       list->Add(g[ii]);
0474     }
0475     list->Draw("AP");
0476     Double_t yaxismax = list->GetYaxis()->GetXmax();
0477     Double_t yaxismin = list->GetYaxis()->GetXmin();
0478     delete list;  //automatically deletes g[i]
0479     if (yaxismin > 0) {
0480       yaxismax += yaxismin;
0481       yaxismin = 0;
0482     }
0483     firstp->GetYaxis()->SetRangeUser(yaxismin, yaxismax);
0484     if (xvar == "runNumber")
0485       firstp->GetXaxis()->SetNdivisions(505);
0486   } else if (type == Histogram || type == OrgHistogram) {
0487     Bool_t allthesame = true;
0488     for (Int_t i = 1; i < n && allthesame; i++) {
0489       if (lengths[i] != lengths[0])
0490         allthesame = false;
0491     }
0492     if (!allthesame && xvar != "runNumber")
0493       for (Int_t i = 0; i < n; i++) {
0494         //This does NOT include events that are out of the run number range (minrun and maxrun).
0495         //It DOES include events that are out of the histogram range.
0496         p[i]->Scale(1.0 / lengths[i]);
0497       }
0498     maxp = (TH1F *)firstp->Clone("maxp");
0499     stufftodelete->Add(maxp);
0500     maxp->SetBit(kCanDelete, true);
0501     maxp->SetLineColor(kWhite);
0502     for (Int_t i = 1; i <= maxp->GetNbinsX(); i++) {
0503       for (Int_t j = 0; j < n; j++) {
0504         if (!used[j])
0505           continue;
0506         maxp->SetBinContent(i, TMath::Max(maxp->GetBinContent(i), p[j]->GetBinContent(i)));
0507       }
0508     }
0509     maxp->SetMarkerStyle(0);
0510     maxp->SetMinimum(0);
0511     maxp->Draw("");
0512     if (xvar == "runNumber") {
0513       maxp->GetXaxis()->SetNdivisions(505);
0514       maxp->Draw("");
0515     }
0516   }
0517
0518   int nEntries = 0;
0519   for (int i = 0; i < n; i++)
0520     if (used[i])
0521       nEntries++;
0522   double width = 0.5;
0523   if (type == Histogram || type == OrgHistogram)
0524     width *= 2;
0525   TLegend *legend = TkAlStyle::legend(nEntries, width);
0526   legend->SetTextSize(0);
0527   if (type == Histogram || type == OrgHistogram)
0528     legend->SetNColumns(2);
0529   stufftodelete->Add(legend);
0530   legend->SetBit(kCanDelete, true);
0531
0532   for (Int_t i = 0; i < n; i++) {
0533     if (!used[i])
0534       continue;
0535     if (type == Resolution || type == Profile) {
0536       if (p[i] == firstp)
0537         p[i]->Draw("P");
0538       else
0539         p[i]->Draw("same P");
0540       legend->AddEntry(p[i], names[i], "pl");
0541     } else if (type == Histogram || type == OrgHistogram) {
0542       if (styles[i] >= 100) {
0543         p[i]->Draw("same P0E");
0544         legend->AddEntry(p[i], names[i], "pl");
0545       } else {
0546         p[i]->Draw("same hist");
0547         legend->AddEntry(p[i], names[i], "l");
0548       }
0549       legend->AddEntry((TObject *)nullptr, meansrmss[i], "");
0550     }
0551   }
0552   if (legend->GetListOfPrimitives()->At(0) == nullptr) {
0553     stufftodelete->Clear();
0554     deleteCanvas(c1);
0555     return nullptr;
0556   }
0557
0558   c1->Update();
0559   legend->Draw();
0560
0561   double legendfraction =
0562       legend->GetY2() -
0563       legend->GetY1();  //apparently GetY1 and GetY2 give NDC coordinates.  This is not a mistake on my part
0564   double padheight = gPad->GetUymax() - gPad->GetUymin();
0565   //legendfraction = legendheight / padheight = newlegendheight / newpadheight
0566   //newpadheight = padheight + x
0567   //newlegendheight = newpadheight - padheight = x so it doesn't cover anything
0568   //==>legendfraction = x/(padheight+x)
0569   /* ==> */ double x = padheight * legendfraction / (1 - legendfraction) * 1.5;  //1.5 to give extra room
0570   maxp->GetYaxis()->SetRangeUser(gPad->GetUymin(), gPad->GetUymax() + x);
0571
0572   TkAlStyle::drawStandardTitle();
0573
0574   c1->Update();
0575
0576   if (saveas != "")
0577     saveplot(c1, saveas);
0578
0579   return c1;
0580 }
0581
0582 //make a 1D histogram of Delta_yvar
0583
0584 TCanvas *trackSplitPlot(Int_t nFiles,
0585                         TString *files,
0586                         TString *names,
0587                         TString var,
0588                         Bool_t relative,
0589                         Bool_t pull,
0590                         TString saveas,
0591                         ostream &summaryfile) {
0592   return trackSplitPlot(nFiles, files, names, "", var, relative, false, pull, saveas, summaryfile);
0593 }
0594
0595 //For 1 file
0596
0597 TCanvas *trackSplitPlot(TString file,
0598                         TString xvar,
0599                         TString yvar,
0600                         Bool_t profile,
0601                         Bool_t relative,
0602                         Bool_t resolution,
0603                         Bool_t pull,
0604                         TString saveas,
0605                         ostream &summaryfile) {
0606   Int_t nFiles = 0;
0607   if (profile)  //it interprets nFiles < 1 as 1 file, make a scatterplot
0608     nFiles = 1;
0609   TString *files = &file;
0610   TString name = "";
0611   TString *names = &name;
0612   return trackSplitPlot(nFiles, files, names, xvar, yvar, relative, resolution, pull, saveas, summaryfile);
0613 }
0614
0615 //make a 1D histogram of Delta_yvar
0616
0617 TCanvas *trackSplitPlot(TString file, TString var, Bool_t relative, Bool_t pull, TString saveas, ostream &summaryfile) {
0618   Int_t nFiles = 1;
0619   TString *files = &file;
0620   TString name = "";
0621   TString *names = &name;
0622   return trackSplitPlot(nFiles, files, names, var, relative, pull, saveas, summaryfile);
0623 }
0624
0625 void saveplot(TCanvas *c1, TString saveas) {
0626   if (saveas == "")
0627     return;
0628   TString saveas2 = saveas, saveas3 = saveas;
0629   saveas2.ReplaceAll(".pngepsroot", "");
0630   saveas3.Remove(saveas3.Length() - 11);
0631   if (saveas2 == saveas3) {
0632     c1->SaveAs(saveas.ReplaceAll(".pngepsroot", ".png"));
0633     c1->SaveAs(saveas.ReplaceAll(".png", ".eps"));
0634     c1->SaveAs(saveas.ReplaceAll(".eps", ".root"));
0635     c1->SaveAs(saveas.ReplaceAll(".root", ".pdf"));
0636   } else {
0637     c1->SaveAs(saveas);
0638   }
0639 }
0640
0641 void deleteCanvas(TObject *canvas) {
0642   if (canvas == nullptr)
0643     return;
0644   if (!canvas->InheritsFrom("TCanvas")) {
0645     delete canvas;
0646     return;
0647   }
0648   TCanvas *c1 = (TCanvas *)canvas;
0649   delete c1;
0650 }
0651
0652 void setupcolors() {
0653   if (colorsset)
0654     return;
0655   colorsset = true;
0656   colors.clear();
0657   styles.clear();
0658   Color_t array[15] = {
0659       1, 2, 3, 4, 6, 7, 8, 9, kYellow + 3, kOrange + 10, kPink - 2, kTeal + 9, kAzure - 8, kViolet - 6, kSpring - 1};
0660   for (int i = 0; i < 15; i++) {
0661     colors.push_back(array[i]);
0662     styles.push_back(1);  //Set the default to 1
0663                           //This is to be consistent with the other validation
0664   }
0665 }
0666
0667 //This makes a plot, of Delta_yvar vs. runNumber, zoomed in to between firstrun and lastrun.
0668 //Each bin contains 1 run.
0669 //Before interpreting the results, make sure to look at the histogram of run number (using yvar = "")
0670 //There might be bins with very few events => big error bars,
0671 //or just 1 event => no error bar
0672
0673 void runNumberZoomed(Int_t nFiles,
0674                      TString *files,
0675                      TString *names,
0676                      TString yvar,
0677                      Bool_t relative,
0678                      Bool_t resolution,
0679                      Bool_t pull,
0680                      Int_t firstRun,
0681                      Int_t lastRun,
0682                      TString saveas) {
0683   Int_t tempminrun = minrun;
0684   Int_t tempmaxrun = maxrun;
0685   minrun = firstRun;
0686   maxrun = lastRun;
0687   trackSplitPlot(nFiles, files, names, "runNumber", yvar, relative, resolution, pull, saveas);
0688   minrun = tempminrun;
0689   maxrun = tempmaxrun;
0690 }
0691
0692 //==========================
0693 //1. Misalignment Dependence
0694 //==========================
0695
0696 //This can do three different things:
0697 // (1) if xvar == "", it will plot the mean (if !resolution) or width (if resolution) of Delta_yvar as a function
0698 //     of the misalignment values, as given in values.  misalignment (e.g. sagitta, elliptical) will be used as the
0699 //     x axis label.
0700 // (2) if xvar != "", it will fit the profile/resolution to a function.  If parameter > 0, it will plot the parameter given by parameter as
0701 //     a function of the misalignment.  parametername is used as the y axis label.  You can put a semicolon in parametername
0702 //     to separate the name from the units.  Functionname describes the funciton, and is put in brackets in the y axis label.
0703 //     For example, to fit to Delta_pt = [0]*(eta_org-[1]), you could use functionname = "#Deltap_{T} = A(#eta-B)",
0704 //     parameter = 0, and parametername = "A;GeV".
0705 // (3) if parameter < 0, it will draw the profile/resolution along with the fitted functions.
0706 //     The parameter of interest is still indicated by parameter, which is transformed to -parameter - 1.
0707 //     For example, -1 --> 0, -2 --> 1, -3 --> 2, ...
0708 //     This parameter's value and error will be in the legend.  You still need to enter parametername and functionname,
0709 //     because they will be used for labels.
0710
0711 //The best way to run misalignmentDependence is through makePlots.  If you want to run misalignmentDependence directly,
0712 //the LAST function, all the way at the bottom of this file, is probably the most practical to use (for all three of these).
0713
0714 // The first function takes a canvas as its argument.  This canvas needs to have been produced with trackSplitPlot using
0715 // the same values of xvar, yvar, relative, resolution, and pull or something strange could happen.
0716
0717 void misalignmentDependence(TCanvas *c1old,
0718                             Int_t nFiles,
0719                             TString *names,
0720                             TString misalignment,
0721                             Double_t *values,
0722                             Double_t *phases,
0723                             TString xvar,
0724                             TString yvar,
0725                             TF1 *function,
0726                             Int_t parameter,
0727                             TString parametername,
0728                             TString functionname,
0729                             Bool_t relative,
0730                             Bool_t resolution,
0731                             Bool_t pull,
0732                             TString saveas) {
0733   if (c1old == nullptr)
0734     return;
0735   c1old = (TCanvas *)c1old->Clone("c1old");
0736   if (misalignment == "" || yvar == "")
0737     return;
0738   Bool_t drawfits = (parameter < 0);
0739   if (parameter < 0)
0740     parameter = -parameter - 1;  //-1 --> 0, -2 --> 1, -3 --> 2, ...
0741   TString yaxislabel = nPart(1, parametername);
0742   TString parameterunits = nPart(2, parametername);
0743   if (parameterunits != "")
0744     yaxislabel.Append(" (").Append(parameterunits).Append(")");
0745   TList *list = c1old->GetListOfPrimitives();
0746   //const int n = list->GetEntries() - 2 - (xvar == "");
0747   const int n = nFiles;
0748
0749   gStyle->SetOptStat(0);
0750   gStyle->SetOptFit(0);
0751   gStyle->SetFitFormat("5.4g");
0752   gStyle->SetFuncColor(2);
0753   gStyle->SetFuncStyle(1);
0754   gStyle->SetFuncWidth(1);
0755
0756   TH1 **p = new TH1 *[n];
0757   TF1 **f = new TF1 *[n];
0758   bool used[n];
0759   for (Int_t i = 0; i < n; i++) {
0760     stringstream s0;
0761     s0 << "p" << i;
0762     TString pname = s0.str();
0763     p[i] = (TH1 *)list->/*At(i+1+(xvar == ""))*/ FindObject(pname);
0764     used[i] = (p[i] != nullptr);
0765     if (used[i])
0766       p[i]->SetDirectory(nullptr);
0767     if (xvar == "")
0768       continue;
0769     stringstream s;
0770     s << function->GetName() << i;
0771     TString newname = s.str();
0772     f[i] = (TF1 *)function->Clone(newname);
0773     stufftodelete->Add(f[i]);
0774   }
0775
0776   Double_t *result = new Double_t[nFiles];
0777   Double_t *error = new Double_t[nFiles];
0778   if (xvar == "") {
0779     yaxislabel = axislabel(yvar, 'y', relative, resolution, pull);
0780     for (Int_t i = 0; i < nFiles; i++) {
0781       if (!used[i])
0782         continue;
0783       if (!resolution) {
0784         result[i] = p[i]->GetMean();
0785         error[i] = p[i]->GetMeanError();
0786       } else {
0787         result[i] = p[i]->GetRMS();
0788         error[i] = p[i]->GetRMSError();
0789       }
0790       cout << result[i] << " +/- " << error[i] << endl;
0791     }
0792   } else {
0793     for (int i = 0; i < n; i++) {
0794       if (!used[i])
0795         continue;
0796       f[i]->SetLineColor(colors[i]);
0797       f[i]->SetLineStyle(styles[i]);
0798       f[i]->SetLineWidth(1);
0799       p[i]->SetMarkerColor(colors[i]);
0800       p[i]->SetMarkerStyle(20 + i);
0801       p[i]->SetLineColor(colors[i]);
0802       p[i]->SetLineStyle(styles[i]);
0803       p[i]->Fit(f[i], "IM");
0804       error[i] = f[i]->GetParError(parameter);
0805       //the fits sometimes don't work if the parameters are constrained.
0806       //take care of the constraining here.
0807       //for sine, make the amplitude positive and the phase between 0 and 2pi.
0808       //unless the amplitude is the only parameter (eg sagitta theta theta)
0809       if (function->GetName() == TString("sine") && function->GetNumberFreeParameters() >= 2) {
0810         if (f[i]->GetParameter(0) < 0) {
0811           f[i]->SetParameter(0, -f[i]->GetParameter(0));
0812           f[i]->SetParameter(2, f[i]->GetParameter(2) + pi);
0813         }
0814         while (f[i]->GetParameter(2) >= 2 * pi)
0815           f[i]->SetParameter(2, f[i]->GetParameter(2) - 2 * pi);
0816         while (f[i]->GetParameter(2) < 0)
0817           f[i]->SetParameter(2, f[i]->GetParameter(2) + 2 * pi);
0818       }
0819       result[i] = f[i]->GetParameter(parameter);
0820     }
0821   }
0822
0823   TCanvas *c1 = TCanvas::MakeDefCanvas();
0824
0825   if (drawfits && xvar != "" && yvar != "") {
0826     TString legendtitle = "[";
0827     legendtitle.Append(functionname);
0828     legendtitle.Append("]");
0829     TLegend *legend = new TLegend(.7, .7, .9, .9, legendtitle, "br");
0830     stufftodelete->Add(legend);
0831     TString drawoption = "";
0832     for (int i = 0; i < n; i++) {
0833       if (!used[i])
0834         continue;
0835       p[i]->Draw(drawoption);
0836       f[i]->Draw("same");
0837       drawoption = "same";
0838
0839       stringstream s;
0840       s.precision(3);
0841       s << nPart(1, parametername) << " = " << result[i] << " #pm " << error[i];
0842       if (parameterunits != "")
0843         s << " " << parameterunits;
0844       TString str = s.str();
0845       legend->AddEntry(p[i], names[i], "pl");
0846       legend->AddEntry(f[i], str, "l");
0847     }
0848     c1->Update();
0849     Double_t x1min = .98 * gPad->GetUxmin() + .02 * gPad->GetUxmax();
0850     Double_t x2max = .02 * gPad->GetUxmin() + .98 * gPad->GetUxmax();
0851     Double_t y1min = .98 * gPad->GetUymin() + .02 * gPad->GetUymax();
0852     Double_t y2max = .02 * gPad->GetUymin() + .98 * gPad->GetUymax();
0853     Double_t width = .4 * (x2max - x1min);
0854     Double_t height = (1. / 20) * legend->GetListOfPrimitives()->GetEntries() * (y2max - y1min);
0855     width *= 2;
0856     height /= 2;
0857     legend->SetNColumns(2);
0858
0859     Double_t newy2max = placeLegend(legend, width, height, x1min, y1min, x2max, y2max);
0860     p[0]->GetYaxis()->SetRangeUser(gPad->GetUymin(), (newy2max - .02 * gPad->GetUymin()) / .98);
0861
0862     legend->SetFillStyle(0);
0863     legend->Draw();
0864   } else {
0865     if (values == nullptr)
0866       return;
0867
0868     Bool_t phasesmatter = false;
0869     if (misalignment == "elliptical" || misalignment == "sagitta" || misalignment == "skew") {
0870       if (phases == nullptr) {
0871         cout << "This misalignment has a phase, but you didn't supply the phases!" << endl
0872              << "Can't produce plots depending on the misalignment value." << endl;
0873         return;
0874       }
0875       int firstnonzero = -1;
0876       for (Int_t i = 0; i < nFiles; i++) {
0877         if (values[i] == 0)
0878           continue;  //if the amplitude is 0 the phase is arbitrary
0879         if (firstnonzero == -1)
0880           firstnonzero = i;
0881         if (phases[i] != phases[firstnonzero])
0882           phasesmatter = true;
0883       }
0884     }
0885
0886     if (!phasesmatter) {
0887       TGraphErrors *g = new TGraphErrors(nFiles, values, result, (Double_t *)nullptr, error);
0888       g->SetName("");
0889       stufftodelete->Add(g);
0890
0891       TString xaxislabel = "#epsilon_{";
0892       xaxislabel.Append(misalignment);
0893       xaxislabel.Append("}");
0894       g->GetXaxis()->SetTitle(xaxislabel);
0895       if (xvar != "") {
0896         yaxislabel.Append("   [");
0897         yaxislabel.Append(functionname);
0898         yaxislabel.Append("]");
0899       }
0900       g->GetYaxis()->SetTitle(yaxislabel);
0901
0902       g->SetMarkerColor(colors[0]);
0903       g->SetMarkerStyle(20);
0904
0905       g->Draw("AP");
0906       Double_t yaxismax = g->GetYaxis()->GetXmax();
0907       Double_t yaxismin = g->GetYaxis()->GetXmin();
0908       if (yaxismin > 0) {
0909         yaxismax += yaxismin;
0910         yaxismin = 0;
0911       }
0912       g->GetYaxis()->SetRangeUser(yaxismin, yaxismax);
0913       g->Draw("AP");
0914     } else {
0915       double *xvalues = new double[nFiles];
0916       double *yvalues = new double[nFiles];  //these are not physically x and y (except in the case of skew)
0917       for (int i = 0; i < nFiles; i++) {
0918         xvalues[i] = values[i] * cos(phases[i]);
0919         yvalues[i] = values[i] * sin(phases[i]);
0920       }
0921       TGraph2DErrors *g =
0922           new TGraph2DErrors(nFiles, xvalues, yvalues, result, (Double_t *)nullptr, (Double_t *)nullptr, error);
0923       g->SetName("");
0924       stufftodelete->Add(g);
0925       delete[] xvalues;  //A TGraph2DErrors has its own copy of xvalues and yvalues, so it's ok to delete these copies.
0926       delete[] yvalues;
0927
0928       TString xaxislabel = "#epsilon_{";
0929       xaxislabel.Append(misalignment);
0930       xaxislabel.Append("}cos(#delta)");
0931       TString realyaxislabel = xaxislabel;
0932       realyaxislabel.ReplaceAll("cos(#delta)", "sin(#delta)");
0933       g->GetXaxis()->SetTitle(xaxislabel);
0934       g->GetYaxis()->SetTitle(realyaxislabel);
0935       TString zaxislabel = /*"fake"*/ yaxislabel;  //yaxislabel is defined earlier
0936       if (xvar != "") {
0937         zaxislabel.Append("   [");
0938         zaxislabel.Append(functionname);
0939         zaxislabel.Append("]");
0940       }
0941       g->GetZaxis()->SetTitle(zaxislabel);
0942       g->SetMarkerStyle(20);
0943       g->Draw("pcolerr");
0944     }
0945   }
0946
0947   if (saveas != "") {
0948     saveplot(c1, saveas);
0949     delete[] p;
0950     delete[] f;
0951     delete[] result;
0952     delete[] error;
0953     delete c1old;
0954   }
0955 }
0956
0957 //This version allows you to show multiple parameters.  It runs the previous version multiple times, once for each parameter.
0958 //saveas will be modified to indicate which parameter is being used each time.
0959
0960 void misalignmentDependence(TCanvas *c1old,
0961                             Int_t nFiles,
0962                             TString *names,
0963                             TString misalignment,
0964                             Double_t *values,
0965                             Double_t *phases,
0966                             TString xvar,
0967                             TString yvar,
0968                             TF1 *function,
0969                             Int_t nParameters,
0970                             Int_t *parameters,
0971                             TString *parameternames,
0972                             TString functionname,
0973                             Bool_t relative,
0974                             Bool_t resolution,
0975                             Bool_t pull,
0976                             TString saveas) {
0977   for (int i = 0; i < nParameters; i++) {
0978     TString saveasi = saveas;
0979     TString insert = nPart(1, parameternames[i]);
0980     insert.Prepend(".");
0981     saveasi.Insert(saveasi.Last('.'), insert);  //insert the parameter name before the file extension
0982     misalignmentDependence(c1old,
0983                            nFiles,
0984                            names,
0985                            misalignment,
0986                            values,
0987                            phases,
0988                            xvar,
0989                            yvar,
0990                            function,
0991                            parameters[i],
0992                            parameternames[i],
0993                            functionname,
0994                            relative,
0995                            resolution,
0996                            pull,
0997                            saveasi);
0998   }
0999 }
1000
1001 //This version does not take a canvas as its argument.  It runs trackSplitPlot to produce the canvas.
1002
1003 void misalignmentDependence(Int_t nFiles,
1004                             TString *files,
1005                             TString *names,
1006                             TString misalignment,
1007                             Double_t *values,
1008                             Double_t *phases,
1009                             TString xvar,
1010                             TString yvar,
1011                             TF1 *function,
1012                             Int_t parameter,
1013                             TString parametername,
1014                             TString functionname,
1015                             Bool_t relative,
1016                             Bool_t resolution,
1017                             Bool_t pull,
1018                             TString saveas) {
1019   misalignmentDependence(trackSplitPlot(nFiles, files, names, xvar, yvar, relative, resolution, pull, ""),
1020                          nFiles,
1021                          names,
1022                          misalignment,
1023                          values,
1024                          phases,
1025                          xvar,
1026                          yvar,
1027                          function,
1028                          parameter,
1029                          parametername,
1030                          functionname,
1031                          relative,
1032                          resolution,
1033                          pull,
1034                          saveas);
1035 }
1036
1037 void misalignmentDependence(Int_t nFiles,
1038                             TString *files,
1039                             TString *names,
1040                             TString misalignment,
1041                             Double_t *values,
1042                             Double_t *phases,
1043                             TString xvar,
1044                             TString yvar,
1045                             TF1 *function,
1046                             Int_t nParameters,
1047                             Int_t *parameters,
1048                             TString *parameternames,
1049                             TString functionname,
1050                             Bool_t relative,
1051                             Bool_t resolution,
1052                             Bool_t pull,
1053                             TString saveas) {
1054   for (int i = 0; i < nParameters; i++) {
1055     TString saveasi = saveas;
1056     TString insert = nPart(1, parameternames[i]);
1057     insert.Prepend(".");
1058     saveasi.Insert(saveasi.Last('.'), insert);  //insert the parameter name before the file extension
1059     misalignmentDependence(nFiles,
1060                            files,
1061                            names,
1062                            misalignment,
1063                            values,
1064                            phases,
1065                            xvar,
1066                            yvar,
1067                            function,
1068                            parameters[i],
1069                            parameternames[i],
1070                            functionname,
1071                            relative,
1072                            resolution,
1073                            pull,
1074                            saveasi);
1075   }
1076 }
1077
1078 // This version allows you to use a string for the function.  It creates a TF1 using this string and uses this TF1
1079
1080 void misalignmentDependence(TCanvas *c1old,
1081                             Int_t nFiles,
1082                             TString *names,
1083                             TString misalignment,
1084                             Double_t *values,
1085                             Double_t *phases,
1086                             TString xvar,
1087                             TString yvar,
1088                             TString function,
1089                             Int_t parameter,
1090                             TString parametername,
1091                             TString functionname,
1092                             Bool_t relative,
1093                             Bool_t resolution,
1094                             Bool_t pull,
1095                             TString saveas) {
1096   TF1 *f = new TF1("func", function);
1097   misalignmentDependence(c1old,
1098                          nFiles,
1099                          names,
1100                          misalignment,
1101                          values,
1102                          phases,
1103                          xvar,
1104                          yvar,
1105                          f,
1106                          parameter,
1107                          parametername,
1108                          functionname,
1109                          relative,
1110                          resolution,
1111                          pull,
1112                          saveas);
1113   delete f;
1114 }
1115
1116 void misalignmentDependence(TCanvas *c1old,
1117                             Int_t nFiles,
1118                             TString *names,
1119                             TString misalignment,
1120                             Double_t *values,
1121                             Double_t *phases,
1122                             TString xvar,
1123                             TString yvar,
1124                             TString function,
1125                             Int_t nParameters,
1126                             Int_t *parameters,
1127                             TString *parameternames,
1128                             TString functionname,
1129                             Bool_t relative,
1130                             Bool_t resolution,
1131                             Bool_t pull,
1132                             TString saveas) {
1133   for (int i = 0; i < nParameters; i++) {
1134     TString saveasi = saveas;
1135     TString insert = nPart(1, parameternames[i]);
1136     insert.Prepend(".");
1137     saveasi.Insert(saveasi.Last('.'), insert);  //insert the parameter name before the file extension
1138     misalignmentDependence(c1old,
1139                            nFiles,
1140                            names,
1141                            misalignment,
1142                            values,
1143                            phases,
1144                            xvar,
1145                            yvar,
1146                            function,
1147                            parameters[i],
1148                            parameternames[i],
1149                            functionname,
1150                            relative,
1151                            resolution,
1152                            pull,
1153                            saveasi);
1154   }
1155 }
1156
1157 void misalignmentDependence(Int_t nFiles,
1158                             TString *files,
1159                             TString *names,
1160                             TString misalignment,
1161                             Double_t *values,
1162                             Double_t *phases,
1163                             TString xvar,
1164                             TString yvar,
1165                             TString function,
1166                             Int_t parameter,
1167                             TString parametername,
1168                             TString functionname,
1169                             Bool_t relative,
1170                             Bool_t resolution,
1171                             Bool_t pull,
1172                             TString saveas) {
1173   TF1 *f = new TF1("func", function);
1174   misalignmentDependence(nFiles,
1175                          files,
1176                          names,
1177                          misalignment,
1178                          values,
1179                          phases,
1180                          xvar,
1181                          yvar,
1182                          f,
1183                          parameter,
1184                          parametername,
1185                          functionname,
1186                          relative,
1187                          resolution,
1188                          pull,
1189                          saveas);
1190   delete f;
1191 }
1192
1193 void misalignmentDependence(Int_t nFiles,
1194                             TString *files,
1195                             TString *names,
1196                             TString misalignment,
1197                             Double_t *values,
1198                             Double_t *phases,
1199                             TString xvar,
1200                             TString yvar,
1201                             TString function,
1202                             Int_t nParameters,
1203                             Int_t *parameters,
1204                             TString *parameternames,
1205                             TString functionname,
1206                             Bool_t relative,
1207                             Bool_t resolution,
1208                             Bool_t pull,
1209                             TString saveas) {
1210   for (int i = 0; i < nParameters; i++) {
1211     TString saveasi = saveas;
1212     TString insert = nPart(1, parameternames[i]);
1213     insert.Prepend(".");
1214     saveasi.Insert(saveasi.Last('.'), insert);  //insert the parameter name before the file extension
1215     misalignmentDependence(nFiles,
1216                            files,
1217                            names,
1218                            misalignment,
1219                            values,
1220                            phases,
1221                            xvar,
1222                            yvar,
1223                            function,
1224                            parameters[i],
1225                            parameternames[i],
1226                            functionname,
1227                            relative,
1228                            resolution,
1229                            pull,
1230                            saveasi);
1231   }
1232 }
1233
1234 //This version does not take a function as its argument.  It automatically determines what function, parameter,
1235 //functionname, and parametername to use based on misalignment, xvar, yvar, relative, resolution, and pull.
1236 //However, you have to manually put into the function which plots to fit to what shapes.
1237 //The 2012A data, using the prompt geometry, is a nice example if you want to see an elliptical misalignment.
1238 //If drawfits is true, it draws the fits; otherwise it plots the parameter as a function of misalignment as given by values.
1239
1240 //If the combination of misalignment, xvar, yvar, relative, resolution, pull has a default function to use, it returns true,
1241 // otherwise it returns false.
1242
1243 //This is the version called by makeThesePlots.C
1244
1245 Bool_t misalignmentDependence(TCanvas *c1old,
1246                               Int_t nFiles,
1247                               TString *names,
1248                               TString misalignment,
1249                               Double_t *values,
1250                               Double_t *phases,
1251                               TString xvar,
1252                               TString yvar,
1253                               Bool_t drawfits,
1254                               Bool_t relative,
1255                               Bool_t resolution,
1256                               Bool_t pull,
1257                               TString saveas) {
1258   if (xvar == "") {
1259     if (c1old == nullptr || misalignment == "" || values == nullptr)
1260       return false;
1261     misalignmentDependence(c1old,
1262                            nFiles,
1263                            names,
1264                            misalignment,
1265                            values,
1266                            phases,
1267                            xvar,
1268                            yvar,
1269                            (TF1 *)nullptr,
1270                            0,
1271                            "",
1272                            "",
1273                            relative,
1274                            resolution,
1275                            pull,
1276                            saveas);
1277     return true;
1278   }
1279   TF1 *f = nullptr;
1280   TString functionname = "";
1281
1282   //if only one parameter is of interest
1283   TString parametername = "";
1284   Int_t parameter = 9999;
1285
1286   //if multiple parameters are of interest
1287   Int_t nParameters = -1;
1288   TString *parameternames = nullptr;
1289   Int_t *parameters = nullptr;
1290
1291   if (misalignment == "sagitta") {
1292     if (xvar == "phi" && yvar == "phi" && !resolution && !pull) {
1293       f = new TF1("sine", "-[0]*cos([1]*x+[2])");
1294       f->FixParameter(1, 1);
1295       f->SetParameter(0, 6e-4);
1296       nParameters = 2;
1297       Int_t tempParameters[2] = {0, 2};
1298       TString tempParameterNames[2] = {"A;mrad", "B"};
1299       parameters = tempParameters;
1300       parameternames = tempParameterNames;
1301       functionname = "#Delta#phi=-Acos(#phi+B)";
1302     }
1303     if (xvar == "theta" && yvar == "theta" && !resolution && pull) {
1304       f = new TF1("line", "-[0]*(x+[1])");
1305       f->FixParameter(1, -pi / 2);
1306       parametername = "A";
1307       functionname = "#Delta#theta/#delta(#Delta#theta)=-A(#theta-#pi/2)";
1308       parameter = 0;
1309     }
1310     if (xvar == "theta" && yvar == "theta" && !resolution && !pull) {
1311       f = new TF1("sine", "[0]*sin([1]*x+[2])");
1312       f->FixParameter(1, 2);
1313       f->FixParameter(2, 0);
1314       parametername = "A;mrad";
1315       functionname = "#Delta#theta=-Asin(2#theta)";
1316       parameter = 0;
1317     }
1318   }
1319   if (misalignment == "elliptical") {
1320     if (xvar == "phi" && yvar == "dxy" && !resolution && !pull) {
1321       f = new TF1("sine", "[0]*sin([1]*x-[2])");
1322       //f = new TF1("sine","[0]*sin([1]*x-[2]) + [3]");
1323       f->FixParameter(1, -2);
1324       f->SetParameter(0, 5e-4);
1325
1326       nParameters = 2;
1327       Int_t tempParameters[2] = {0, 2};
1328       TString tempParameterNames[2] = {"A;#mum", "B"};
1329       //nParameters = 3;
1330       //Int_t tempParameters[3] = {0,2,3};
1331       //TString tempParameterNames[3] = {"A;#mum","B","C;#mum"};
1332
1333       parameters = tempParameters;
1334       parameternames = tempParameterNames;
1335       functionname = "#Deltad_{xy}=-Asin(2#phi+B)";
1336       //functionname = "#Deltad_{xy}=-Asin(2#phi+B)+C";
1337     }
1338     if (xvar == "phi" && yvar == "dxy" && !resolution && pull) {
1339       f = new TF1("sine", "[0]*sin([1]*x-[2])");
1340       //f = new TF1("sine","[0]*sin([1]*x-[2]) + [3]");
1341
1342       f->FixParameter(1, -2);
1343
1344       nParameters = 2;
1345       Int_t tempParameters[2] = {0, 2};
1346       TString tempParameterNames[2] = {"A", "B"};
1347       //nParameters = 3;
1348       //Int_t tempParameters[3] = {0,2,3};
1349       //TString tempParameterNames[3] = {"A","B","C"};
1350
1351       parameters = tempParameters;
1352       parameternames = tempParameterNames;
1353
1354       functionname = "#Deltad_{xy}/#delta(#Deltad_{xy})=-Asin(2#phi+B)";
1355       //functionname = "#Deltad_{xy}/#delta(#Deltad_{xy})=-Asin(2#phi+B)+C";
1356     }
1357
1358     if (xvar == "theta" && yvar == "dz" && !resolution && !pull) {
1359       f = new TF1("line", "-[0]*(x-[1])");
1360       f->FixParameter(1, pi / 2);
1361       parametername = "A;#mum";
1362       functionname = "#Deltad_{z}=-A(#theta-#pi/2)";
1363       parameter = 0;
1364     }
1365     /*
1366         This fit doesn't work
1367         if (xvar == "theta" && yvar == "dz" && !resolution && pull)
1368         {
1369             f = new TF1("sine","[0]*sin([1]*x+[2])");
1370             f->FixParameter(2,-pi/2);
1371             f->FixParameter(1,1);
1372             parametername = "A";
1373             functionname = "#Deltad_{z}/#delta(#Deltad_{z})=Acos(#theta)";
1374             parameter = 0;
1375         }
1376         */
1377     if (xvar == "dxy" && yvar == "phi" && !resolution && !pull) {
1378       f = new TF1("line", "-[0]*(x-[1])");
1379       f->FixParameter(1, 0);
1380       parametername = "A;mrad/cm";
1381       functionname = "#Delta#phi=-A(d_{xy})";
1382       parameter = 0;
1383     }
1384     if (xvar == "dxy" && yvar == "phi" && !resolution && pull) {
1385       f = new TF1("line", "-[0]*(x-[1])");
1386       f->FixParameter(1, 0);
1387       parametername = "A;cm^{-1}";
1388       functionname = "#Delta#phi/#delta(#Delta#phi)=-A(d_{xy})";
1389       parameter = 0;
1390     }
1391   }
1392   if (misalignment == "skew") {
1393     if (xvar == "phi" && yvar == "theta" && resolution && !pull) {
1394       f = new TF1("sine", "[0]*sin([1]*x+[2])+[3]");
1395       f->FixParameter(1, 2);
1396       nParameters = 3;
1397       Int_t tempParameters[3] = {0, 2, 3};
1398       TString tempParameterNames[3] = {"A;mrad", "B", "C;mrad"};
1399       parameters = tempParameters;
1400       parameternames = tempParameterNames;
1401       functionname = "#sigma(#Delta#theta)=Asin(2#phi+B)+C";
1402     }
1403     if (xvar == "phi" && yvar == "eta" && resolution && !pull) {
1404       f = new TF1("sine", "[0]*sin([1]*x+[2])+[3]");
1405       f->FixParameter(1, 2);
1406       nParameters = 3;
1407       Int_t tempParameters[3] = {0, 2, 3};
1408       TString tempParameterNames[3] = {"A;mrad", "B", "C;mrad"};
1409       parameters = tempParameters;
1410       parameternames = tempParameterNames;
1411       functionname = "#sigma(#Delta#eta)=Asin(2#phi+B)+C";
1412     }
1413     if (xvar == "phi" && yvar == "theta" && resolution && pull) {
1414       f = new TF1("sine", "[0]*sin([1]*x+[2])+[3]");
1415       f->FixParameter(1, 2);
1416       nParameters = 3;
1417       Int_t tempParameters[3] = {0, 2, 3};
1418       TString tempParameterNames[3] = {"A", "B", "C"};
1419       parameters = tempParameters;
1420       parameternames = tempParameterNames;
1421       functionname = "#sigma(#Delta#theta/#delta(#Delta#theta))=Asin(2#phi+B)+C";
1422     }
1423     if (xvar == "phi" && yvar == "eta" && resolution && pull) {
1424       f = new TF1("sine", "[0]*sin([1]*x+[2])+[3]");
1425       f->FixParameter(1, 2);
1426       nParameters = 3;
1427       Int_t tempParameters[3] = {0, 2, 3};
1428       TString tempParameterNames[3] = {"A", "B", "C"};
1429       parameters = tempParameters;
1430       parameternames = tempParameterNames;
1431       functionname = "#sigma(#Delta#eta/#delta(#Delta#eta))=Asin(2#phi+B)+C";
1432     }
1433     if (xvar == "phi" && yvar == "dz" && !resolution && !pull) {
1434       f = new TF1("tanh", "[0]*(tanh([1]*(x+[2]))   )");  // - tanh(([3]-[1])*x+[2]) + 1)");
1435       //f = new TF1("tanh","[0]*(tanh([1]*(x+[2])) + tanh([1]*([3]-[2]-x)) - 1)");
1436       f->SetParameter(0, 100);
1437       f->SetParLimits(1, -20, 20);
1438       f->SetParLimits(2, 0, pi);
1439       f->FixParameter(3, pi);
1440       nParameters = 3;
1441       Int_t tempParameters[3] = {0, 1, 2};
1442       TString tempParameterNames[3] = {"A;#mum", "B", "C"};
1443       parameters = tempParameters;
1444       parameternames = tempParameterNames;
1445       functionname = "#Deltad_{z}=Atanh(B(#phi+C))";
1446       //functionname = "#Deltad_{z}=A(tanh(B(#phi+C)) + tanh(B(#pi-#phi-C)) - 1";
1447     }
1448   }
1449   if (misalignment == "layerRot") {
1450     if (xvar == "qoverpt" && yvar == "qoverpt" && !relative && !resolution && !pull) {
1451       f = new TF1("sech", "[0]/cosh([1]*(x+[2]))+[3]");
1452       //f = new TF1("gauss","[0]/exp(([1]*(x+[2]))^2)+[3]");   //sech works better than a gaussian
1453       f->SetParameter(0, 1);
1454       f->SetParameter(1, 1);
1455       f->SetParLimits(1, 0, 10);
1456       f->FixParameter(2, 0);
1457       f->FixParameter(3, 0);
1458       nParameters = 2;
1459       Int_t tempParameters[2] = {0, 1};
1460       TString tempParameterNames[2] = {"A;10^{-3}e/GeV", "B;GeV/e"};
1461       parameters = tempParameters;
1462       parameternames = tempParameterNames;
1463       functionname = "#Delta(q/p_{T})=Asech(B(q/p_{T}))";
1464     }
1465   }
1466   if (misalignment == "telescope") {
1467     if (xvar == "theta" && yvar == "theta" && !relative && !resolution && !pull) {
1468       f = new TF1("gauss", "[0]/exp(([1]*(x+[2]))^2)+[3]");
1469       f->SetParameter(0, 1);
1470       f->SetParameter(1, 1);
1471       f->SetParLimits(1, 0, 10);
1472       f->FixParameter(2, -pi / 2);
1473       f->FixParameter(3, 0);
1474       nParameters = 2;
1475       Int_t tempParameters[2] = {0, 1};
1476       TString tempParameterNames[2] = {"A;mrad", "B"};
1477       parameters = tempParameters;
1478       parameternames = tempParameterNames;
1479       functionname = "#Delta#theta=Aexp(-(B(#theta-#pi/2))^{2})";
1480     }
1481   }
1482   if (functionname == "")
1483     return false;
1484   if (drawfits) {
1485     parameter = -parameter - 1;
1486     for (int i = 0; i < nParameters; i++)
1487       parameters[i] = -parameters[i] - 1;
1488   }
1489   if (nParameters > 0)
1490     misalignmentDependence(c1old,
1491                            nFiles,
1492                            names,
1493                            misalignment,
1494                            values,
1495                            phases,
1496                            xvar,
1497                            yvar,
1498                            f,
1499                            nParameters,
1500                            parameters,
1501                            parameternames,
1502                            functionname,
1503                            relative,
1504                            resolution,
1505                            pull,
1506                            saveas);
1507   else
1508     misalignmentDependence(c1old,
1509                            nFiles,
1510                            names,
1511                            misalignment,
1512                            values,
1513                            phases,
1514                            xvar,
1515                            yvar,
1516                            f,
1517                            parameter,
1518                            parametername,
1519                            functionname,
1520                            relative,
1521                            resolution,
1522                            pull,
1523                            saveas);
1524   delete f;
1525   return true;
1526 }
1527
1528 //This is the most practically useful version.  It does not take a canvas, but produces it automatically and then determines what
1529 //function to fit it to.
1530
1531 Bool_t misalignmentDependence(Int_t nFiles,
1532                               TString *files,
1533                               TString *names,
1534                               TString misalignment,
1535                               Double_t *values,
1536                               Double_t *phases,
1537                               TString xvar,
1538                               TString yvar,
1539                               Bool_t drawfits,
1540                               Bool_t relative,
1541                               Bool_t resolution,
1542                               Bool_t pull,
1543                               TString saveas) {
1544   return misalignmentDependence(trackSplitPlot(nFiles, files, names, xvar, yvar, relative, resolution, pull, ""),
1545                                 nFiles,
1546                                 names,
1547                                 misalignment,
1548                                 values,
1549                                 phases,
1550                                 xvar,
1551                                 yvar,
1552                                 drawfits,
1553                                 relative,
1554                                 resolution,
1555                                 pull,
1556                                 saveas);
1557 }
1558
1559 Bool_t hasFit(TString misalignment, TString xvar, TString yvar, Bool_t relative, Bool_t resolution, Bool_t pull) {
1560   return misalignmentDependence((TCanvas *)nullptr,
1561                                 0,
1562                                 (TString *)nullptr,
1563                                 misalignment,
1564                                 (Double_t *)nullptr,
1565                                 (Double_t *)nullptr,
1566                                 xvar,
1567                                 yvar,
1568                                 false,
1569                                 relative,
1570                                 resolution,
1571                                 pull,
1572                                 TString(""));
1573 }
1574
1575 //=============
1576 //2. Make Plots
1577 //=============
1578
1579 void makePlots(Int_t nFiles,
1580                TString *files,
1581                TString *names,
1582                TString misalignment,
1583                Double_t *values,
1584                Double_t *phases,
1585                TString directory,
1586                Bool_t matrix[xsize][ysize]) {
1587   stufftodelete->SetOwner(true);
1588
1589   for (Int_t i = 0, totaltime = 0; i < nFiles; i++) {
1590     TFile *f = nullptr;
1591     bool exists = false;
1592     if (files[i] == "")
1593       exists = true;
1594
1595     for (int j = 1; j <= 60 * 24 && !exists; j++, totaltime++)  //wait up to 1 day for the validation to be finished
1596     {
1597       f = TFile::Open(files[i]);
1598       if (f != nullptr)
1599         exists = f->IsOpen();
1600       delete f;
1601       if (exists)
1602         continue;
1603       gSystem->Sleep(60000);
1604       cout << "It's been ";
1605       if (j >= 60)
1606         cout << j / 60 << " hour";
1607       if (j >= 120)
1608         cout << "s";
1609       if (j % 60 != 0 && j >= 60)
1610         cout << " and ";
1611       if (j % 60 != 0)
1612         cout << j % 60 << " minute";
1613       if (j % 60 >= 2)
1614         cout << "s";
1615       cout << endl;
1616     }
1617     if (!exists)
1618       return;
1619     if (i == nFiles - 1 && totaltime > nFiles)
1620       gSystem->Sleep(60000);
1621   }
1622
1623   TString directorytomake = directory;
1624   gSystem->mkdir(directorytomake, true);
1625
1626   ofstream summaryfile(directorytomake + "/TrackSplittingValidationSummary.txt");
1627   for (int i = 0; i < nFiles; i++) {
1628     summaryfile << "\t" << TString(names[i]).ReplaceAll("#", "\\");
1629   }
1630   summaryfile << "\tformat={}\tlatexformat={}\n";
1631
1632   if (misalignment != "") {
1633     directorytomake.Append("/fits");
1634     gSystem->mkdir(directorytomake);
1635   }
1636
1637   for (Int_t x = 0; x < xsize; x++) {
1638     for (Int_t y = 0; y < ysize; y++) {
1639       for (Int_t pull = 0; pull == 0 || (pull == 1 && yvariables[y] != ""); pull++) {
1640         if (false)
1641           continue;  //this line is to make it easier to do e.g. all plots involving Delta eta
1642                      //(replace false with yvariables[y] != "eta")
1643
1644         if (!matrix[x][y])
1645           continue;
1646
1647         if (xvariables[x] == "" && yvariables[y] == "")
1648           continue;
1649
1650         Int_t nPlots =
1651             nFiles + 4;  //scatterplot for each (if you uncomment it), profile, resolution, and fits for each.
1652         vector<TString> s;
1653
1654         TString slashstring = "";
1655         if (directory.Last('/') != directory.Length() - 1)
1656           slashstring = "/";
1657
1658         vector<TString> plotnames;
1659         for (Int_t i = 0; i < nFiles; i++) {
1660           plotnames.push_back(names[i]);  //this is plotnames[i]
1661           plotnames[i].ReplaceAll(" ", "");
1662         }
1663
1664         plotnames.push_back("");  //this is plotnames[nFiles], but gets changed
1665         if (yvariables[y] == "")
1666           plotnames[nFiles] = "orghist";
1667         else if (xvariables[x] == "")
1668           plotnames[nFiles] = "hist";
1669         else
1670           plotnames[nFiles] = "profile";
1671
1672         plotnames.push_back("resolution");  //this is plotnames[nFiles+1]
1673
1674         plotnames.push_back("");  //this is plotnames[nFiles+2]
1675         plotnames.push_back("");  //this is plotnames[nFiles+3]
1676         if (plotnames[nFiles] == "profile") {
1677           plotnames[nFiles + 2] = ".profile";
1678           plotnames[nFiles + 2].Prepend(misalignment);
1679           plotnames[nFiles + 3] = ".resolution";
1680           plotnames[nFiles + 3].Prepend(misalignment);
1681           plotnames[nFiles + 2].Prepend("fits/");
1682           plotnames[nFiles + 3].Prepend("fits/");
1683         } else {
1684           plotnames[nFiles + 2] = "profile.";
1685           plotnames[nFiles + 2].Append(misalignment);
1686           plotnames[nFiles + 3] = "resolution.";
1687           plotnames[nFiles + 3].Append(misalignment);
1688         }
1689
1690         TString pullstring = "";
1691         if (pull)
1692           pullstring = "pull.";
1693
1694         TString xvarstring = xvariables[x];
1695         if (xvariables[x] != "runNumber" && !xvariables[x].BeginsWith("nHits") && xvariables[x] != "")
1696           xvarstring.Append("_org");
1697         if (xvariables[x] != "" && yvariables[y] != "")
1698           xvarstring.Append(".");
1699
1700         TString yvarstring = yvariables[y];
1701         if (yvariables[y] != "")
1702           yvarstring.Prepend("Delta_");
1703
1704         TString relativestring = "";
1705         if (relativearray[y])
1706           relativestring = ".relative";
1707
1708         for (Int_t i = 0; i < nPlots; i++) {
1709           stringstream ss;
1710           ss << directory << slashstring << plotnames[i] << "." << pullstring << xvarstring << yvarstring
1711              << relativestring << ".pngepsroot";
1712           s.push_back(ss.str());
1713           if (misalignment != "") {
1714             TString wrongway = misalignment;
1715             TString rightway = misalignment;
1716             wrongway.Append(".pull");
1717             rightway.Prepend("pull.");
1718             s[i].ReplaceAll(wrongway, rightway);
1719           }
1720         }
1721
1722         Int_t i;
1723         for (i = 0; i < nFiles; i++) {
1724           if (xvariables[x] == "" || yvariables[y] == "")
1725             continue;
1726           //uncomment this section to make scatterplots
1727           /*
1728                     trackSplitPlot(files[i],xvariables[x],yvariables[y],false,relativearray[y],false,(bool)pull,s[i]);
1729                     stufftodelete->Clear();
1730                     for ( ; gROOT->GetListOfCanvases()->GetEntries() > 0; )
1731                         deleteCanvas( gROOT->GetListOfCanvases()->Last());
1732                     for ( ; gROOT->GetListOfFiles()->GetEntries() > 0; )
1733                         delete (TFile*)gROOT->GetListOfFiles()->Last();
1734                     */
1735         }
1736
1737         if (xvariables[x] != "" && yvariables[y] != "") {
1738           //make profile
1739           TCanvas *c1 = trackSplitPlot(
1740               nFiles, files, names, xvariables[x], yvariables[y], relativearray[y], false, (bool)pull, s[i], summaryfile);
1741           if (misalignmentDependence(c1,
1742                                      nFiles,
1743                                      names,
1744                                      misalignment,
1745                                      values,
1746                                      phases,
1747                                      xvariables[x],
1748                                      yvariables[y],
1749                                      true,
1750                                      relativearray[y],
1751                                      false,
1752                                      (bool)pull,
1753                                      s[i + 2])) {
1754             s[i + 2].ReplaceAll(".png", ".parameter.png");
1755             misalignmentDependence(c1,
1756                                    nFiles,
1757                                    names,
1758                                    misalignment,
1759                                    values,
1760                                    phases,
1761                                    xvariables[x],
1762                                    yvariables[y],
1763                                    false,
1764                                    relativearray[y],
1765                                    false,
1766                                    (bool)pull,
1767                                    s[i + 2]);
1768           }
1769           stufftodelete->Clear();
1770           for (; gROOT->GetListOfCanvases()->GetEntries() > 0;)
1771             deleteCanvas(gROOT->GetListOfCanvases()->Last());
1772           for (; gROOT->GetListOfFiles()->GetEntries() > 0;)
1773             delete (TFile *)gROOT->GetListOfFiles()->Last();
1774
1775           //make resolution plot
1776           TCanvas *c2 = trackSplitPlot(nFiles,
1777                                        files,
1778                                        names,
1779                                        xvariables[x],
1780                                        yvariables[y],
1781                                        relativearray[y],
1782                                        true,
1783                                        (bool)pull,
1784                                        s[i + 1],
1785                                        summaryfile);
1786           if (misalignmentDependence(c2,
1787                                      nFiles,
1788                                      names,
1789                                      misalignment,
1790                                      values,
1791                                      phases,
1792                                      xvariables[x],
1793                                      yvariables[y],
1794                                      true,
1795                                      relativearray[y],
1796                                      true,
1797                                      (bool)pull,
1798                                      s[i + 3])) {
1799             s[i + 3].ReplaceAll(".png", ".parameter.png");
1800             misalignmentDependence(c2,
1801                                    nFiles,
1802                                    names,
1803                                    misalignment,
1804                                    values,
1805                                    phases,
1806                                    xvariables[x],
1807                                    yvariables[y],
1808                                    false,
1809                                    relativearray[y],
1810                                    true,
1811                                    (bool)pull,
1812                                    s[i + 3]);
1813           }
1814           stufftodelete->Clear();
1815           for (; gROOT->GetListOfCanvases()->GetEntries() > 0;)
1816             deleteCanvas(gROOT->GetListOfCanvases()->Last());
1817           for (; gROOT->GetListOfFiles()->GetEntries() > 0;)
1818             delete (TFile *)gROOT->GetListOfFiles()->Last();
1819         } else {
1820           //make histogram
1821           TCanvas *c1 = trackSplitPlot(
1822               nFiles, files, names, xvariables[x], yvariables[y], relativearray[y], false, (bool)pull, s[i], summaryfile);
1823           if (misalignmentDependence(c1,
1824                                      nFiles,
1825                                      names,
1826                                      misalignment,
1827                                      values,
1828                                      phases,
1829                                      xvariables[x],
1830                                      yvariables[y],
1831                                      true,
1832                                      relativearray[y],
1833                                      false,
1834                                      (bool)pull,
1835                                      s[i + 2])) {
1836             misalignmentDependence(c1,
1837                                    nFiles,
1838                                    names,
1839                                    misalignment,
1840                                    values,
1841                                    phases,
1842                                    xvariables[x],
1843                                    yvariables[y],
1844                                    true,
1845                                    relativearray[y],
1846                                    true,
1847                                    (bool)pull,
1848                                    s[i + 3]);
1849           }
1850           stufftodelete->Clear();
1851           for (; gROOT->GetListOfCanvases()->GetEntries() > 0;)
1852             deleteCanvas(gROOT->GetListOfCanvases()->Last());
1853           for (; gROOT->GetListOfFiles()->GetEntries() > 0;)
1854             delete (TFile *)gROOT->GetListOfFiles()->Last();
1855         }
1856       }
1857       cout << y + ysize * x + 1 << "/" << xsize * ysize << endl;
1858     }
1859   }
1860 }
1861
1862 void makePlots(Int_t nFiles, TString *files, TString *names, TString directory, Bool_t matrix[xsize][ysize]) {
1863   makePlots(nFiles, files, names, "", (Double_t *)nullptr, (Double_t *)nullptr, directory, matrix);
1864 }
1865
1866 void makePlots(TString file,
1867                TString misalignment,
1868                Double_t *values,
1869                Double_t *phases,
1870                TString directory,
1871                Bool_t matrix[xsize][ysize]) {
1872   setupcolors();
1873   file.Remove(TString::kTrailing, ',');
1874   unsigned int n = file.CountChar(',') + 1;
1875   TString *files = new TString[n];
1876   TString *names = new TString[n];
1877   vector<Color_t> tempcolors = colors;
1878   vector<Style_t> tempstyles = styles;
1879   for (unsigned int i = 0; i < n; i++) {
1880     TString thisfile = nPart(i + 1, file, ",");
1881     int numberofpipes = thisfile.CountChar('|');
1882     if (numberofpipes >= 0 && nPart(numberofpipes + 1, thisfile, "|").IsDigit()) {
1883       if (numberofpipes >= 1 && nPart(numberofpipes, thisfile, "|").IsDigit()) {
1884         colors[i] = nPart(numberofpipes, thisfile, "|").Atoi();
1885         styles[i] = nPart(numberofpipes + 1, thisfile, "|").Atoi();
1886         thisfile.Remove(thisfile.Length() - nPart(numberofpipes, thisfile, "|").Length() -
1887                         nPart(numberofpipes + 1, thisfile, "|").Length() - 2);
1888       } else {
1889         colors[i] = nPart(numberofpipes + 1, thisfile, "|").Atoi();
1890         thisfile.Remove(thisfile.Length() - nPart(numberofpipes + 1, thisfile, "|").Length() - 2);
1891       }
1892     }
1893     files[i] = nPart(1, thisfile, "=", true);
1894     names[i] = nPart(2, thisfile, "=", false);
1895   }
1896   if (n == 1 && names[0] == "")
1897     names[0] =
1898         "scatterplot";  //With 1 file there's no legend, so this is only used in the filename of the scatterplots, if made
1899   makePlots(n, files, names, misalignment, values, phases, directory, matrix);
1900   delete[] files;
1901   delete[] names;
1902   colors = tempcolors;
1903   styles = tempstyles;
1904 }
1905
1906 void makePlots(TString file, TString directory, Bool_t matrix[xsize][ysize]) {
1907   makePlots(file, "", (Double_t *)nullptr, (Double_t *)nullptr, directory, matrix);
1908 }
1909
1910 //***************************************************************************
1911 //functions to make plots for 1 row, column, or cell of the matrix
1912 //examples:
1913 //   xvar = "nHits", yvar = "ptrel" - makes plots of nHits vs Delta_pt/pt_org
1914 //   xvar = "all",   yvar = "pt"    - makes all plots involving Delta_pt
1915 //                                    (not Delta_pt/pt_org)
1916 //   xvar = "",      yvar = "all"   - makes all histograms of Delta_???
1917 //                                    (including Delta_pt/pt_org)
1918 //***************************************************************************
1919
1920 void makePlots(Int_t nFiles,
1921                TString *files,
1922                TString *names,
1923                TString misalignment,
1924                Double_t *values,
1925                Double_t *phases,
1926                TString directory,
1927                TString xvar,
1928                TString yvar) {
1929   Bool_t matrix[xsize][ysize];
1930   for (int x = 0; x < xsize; x++)
1931     for (int y = 0; y < ysize; y++) {
1932       bool xmatch = (xvar == "all" || xvar == xvariables[x]);
1933       bool ymatch = (yvar == "all" || yvar == yvariables[y]);
1934       if (yvar == "pt" && yvariables[y] == "pt" && relativearray[y] == true)
1935         ymatch = false;
1936       if (yvar == "ptrel" && yvariables[y] == "pt" && relativearray[y] == true)
1937         ymatch = true;
1938       matrix[x][y] = (xmatch && ymatch);
1939     }
1940   makePlots(nFiles, files, names, misalignment, values, phases, directory, matrix);
1941 }
1942
1943 void makePlots(Int_t nFiles, TString *files, TString *names, TString directory, TString xvar, TString yvar) {
1944   makePlots(nFiles, files, names, "", (Double_t *)nullptr, (Double_t *)nullptr, directory, xvar, yvar);
1945 }
1946
1947 void makePlots(TString file,
1948                TString misalignment,
1949                Double_t *values,
1950                Double_t *phases,
1951                TString directory,
1952                TString xvar,
1953                TString yvar) {
1954   setupcolors();
1955   file.Remove(TString::kTrailing, ',');
1956   unsigned int n = file.CountChar(',') + 1;
1957   TString *files = new TString[n];
1958   TString *names = new TString[n];
1959   vector<Color_t> tempcolors = colors;
1960   vector<Style_t> tempstyles = styles;
1961   for (unsigned int i = 0; i < n; i++) {
1962     TString thisfile = nPart(i + 1, file, ",");
1963     int numberofpipes = thisfile.CountChar('|');
1964     if (numberofpipes >= 0 && nPart(numberofpipes + 1, thisfile, "|").IsDigit()) {
1965       if (numberofpipes >= 1 && nPart(numberofpipes, thisfile, "|").IsDigit()) {
1966         colors[i] = nPart(numberofpipes, thisfile, "|").Atoi();
1967         styles[i] = nPart(numberofpipes + 1, thisfile, "|").Atoi();
1968         thisfile.Remove(thisfile.Length() - nPart(numberofpipes, thisfile, "|").Length() -
1969                         nPart(numberofpipes + 1, thisfile, "|").Length() - 2);
1970       } else {
1971         colors[i] = nPart(numberofpipes + 1, thisfile, "|").Atoi();
1972         thisfile.Remove(thisfile.Length() - nPart(numberofpipes + 1, thisfile, "|").Length() - 2);
1973       }
1974     }
1975     files[i] = nPart(1, thisfile, "=", true);
1976     names[i] = nPart(2, thisfile, "=", false);
1977   }
1978   if (n == 1 && names[0] == "")
1979     names[0] =
1980         "scatterplot";  //With 1 file there's no legend, so this is only used in the filename of the scatterplots, if made
1981   makePlots(n, files, names, misalignment, values, phases, directory, xvar, yvar);
1982   delete[] files;
1983   delete[] names;
1984   colors = tempcolors;
1985   styles = tempstyles;
1986 }
1987
1988 void makePlots(TString file, TString directory, TString xvar, TString yvar) {
1989   makePlots(file, "", (Double_t *)nullptr, (Double_t *)nullptr, directory, xvar, yvar);
1990 }
1991
1992 //***************************
1993 //functions to make all plots
1994 //***************************
1995
1996 void makePlots(Int_t nFiles,
1997                TString *files,
1998                TString *names,
1999                TString misalignment,
2000                Double_t *values,
2001                Double_t *phases,
2002                TString directory) {
2003   makePlots(nFiles, files, names, misalignment, values, phases, directory, "all", "all");
2004 }
2005
2006 void makePlots(Int_t nFiles, TString *files, TString *names, TString directory) {
2007   makePlots(nFiles, files, names, "", (Double_t *)nullptr, (Double_t *)nullptr, directory);
2008 }
2009
2010 void makePlots(TString file, TString misalignment, Double_t *values, Double_t *phases, TString directory) {
2011   setupcolors();
2012   file.Remove(TString::kTrailing, ',');
2013   unsigned int n = file.CountChar(',') + 1;
2014   TString *files = new TString[n];
2015   TString *names = new TString[n];
2016   vector<Color_t> tempcolors = colors;
2017   vector<Style_t> tempstyles = styles;
2018   for (unsigned int i = 0; i < n; i++) {
2019     TString thisfile = nPart(i + 1, file, ",");
2020     int numberofpipes = thisfile.CountChar('|');
2021     if (numberofpipes >= 0 && nPart(numberofpipes + 1, thisfile, "|").IsDigit()) {
2022       if (numberofpipes >= 1 && nPart(numberofpipes, thisfile, "|").IsDigit()) {
2023         colors[i] = nPart(numberofpipes, thisfile, "|").Atoi();
2024         styles[i] = nPart(numberofpipes + 1, thisfile, "|").Atoi();
2025         thisfile.Remove(thisfile.Length() - nPart(numberofpipes, thisfile, "|").Length() -
2026                         nPart(numberofpipes + 1, thisfile, "|").Length() - 2);
2027       } else {
2028         colors[i] = nPart(numberofpipes + 1, thisfile, "|").Atoi();
2029         thisfile.Remove(thisfile.Length() - nPart(numberofpipes + 1, thisfile, "|").Length() - 2);
2030       }
2031     }
2032     files[i] = nPart(1, thisfile, "=", true);
2033     names[i] = nPart(2, thisfile, "=", false);
2034   }
2035   if (n == 1 && names[0] == "")
2036     names[0] =
2037         "scatterplot";  //With 1 file there's no legend, so this is only used in the filename of the scatterplots, if made
2038   makePlots(n, files, names, misalignment, values, phases, directory);
2039   delete[] files;
2040   delete[] names;
2041   colors = tempcolors;
2042   styles = tempstyles;
2043 }
2044
2045 void makePlots(TString file, TString directory) {
2046   makePlots(file, "", (Double_t *)nullptr, (Double_t *)nullptr, directory);
2047 }
2048
2049 //=============
2050 //3. Axis Label
2051 //=============
2052
2053 TString fancyname(TString variable) {
2054   if (variable == "pt")
2055     return "p_{T}";
2056   else if (variable == "phi")
2057     return "#phi";
2058   else if (variable == "eta")
2059     return "#eta";
2060   else if (variable == "theta")
2061     return "#theta";
2062   else if (variable == "qoverpt")
2063     return "q/p_{T}";
2064   else if (variable == "runNumber")
2065     return "run number";
2066   else if (variable == "dxy" || variable == "dz")
2067     return variable.ReplaceAll("d", "d_{").Append("}");
2068   else
2069     return variable;
2070 }
2071
2072 //this gives the units, to be put in the axis label
2073 TString units(TString variable, Char_t axis) {
2074   if (variable == "pt")
2075     return "GeV";
2076   if (variable == "dxy" || variable == "dz") {
2077     if (axis == 'y')
2078       return "#mum";  //in the tree, it's listed in centimeters, but in trackSplitPlot the value is divided by 1e4
2079     if (axis == 'x')
2080       return "cm";
2081   }
2082   if (variable == "qoverpt") {
2083     if (axis == 'y')
2084       return "#times10^{-3}e/GeV";  //e/TeV is not particularly intuitive
2085     if (axis == 'x')
2086       return "e/GeV";
2087   }
2088   if (axis == 'y' && (variable == "phi" || variable == "theta"))
2089     return "mrad";
2090   return "";
2091 }
2092
2093 TString plainunits(TString variable, char axis) {
2094   TString result = units(variable, axis);
2095   result.ReplaceAll("#mu", "u");
2096   result.ReplaceAll("#times10^{-3}", "* 1e-3 ");
2097   return result;
2098 }
2099
2100 TString latexunits(TString variable, char axis) {
2101   TString result = units(variable, axis);
2102   result.ReplaceAll("#", "\\")
2103       .ReplaceAll("{", "{{")
2104       .ReplaceAll("}", "}}")
2105       .ReplaceAll("\\mum", "$\\mu$m")
2106       .ReplaceAll("\\times10^{{-3}}", "$\\times10^{{-3}}$");
2107   return result;
2108 }
2109
2110 //this gives the full axis label, including units.  It can handle any combination of relative, resolution, and pull.
2111 TString axislabel(TString variable, Char_t axis, Bool_t relative, Bool_t resolution, Bool_t pull) {
2112   if (axis == 'X' || axis == 'Y') {
2113     double min, max, bins;
2114     axislimits(0, nullptr, variable, tolower(axis), relative, pull, min, max, bins);
2115
2116     if (variable.BeginsWith("nHits"))
2117       return "fraction of tracks";
2118     if (variable == "runNumber")
2119       return "number of tracks";
2120
2121     stringstream s;
2122     s << "fraction of tracks / " << (max - min) / bins;
2123     if (!pull && !relative) {
2124       TString varunits = units(variable, tolower(axis));
2125       if (varunits != "")
2126         s << " " << varunits;
2127     }
2128     TString result = s.str();
2129     result.ReplaceAll(" #times", "#times");
2130     return result;
2131   }
2132
2133   stringstream s;
2134   if (resolution && axis == 'y')
2135     s << "#sigma(";
2136   if (axis == 'y')
2137     s << "#Delta";
2138   s << fancyname(variable);
2139   if (relative && axis == 'y') {
2140     s << " / ";
2141     if (!pull)
2142       s << "(";
2143     s << fancyname(variable);
2144   }
2145   if (axis == 'y') {
2146     if (pull) {
2147       s << " / #delta(#Delta" << fancyname(variable);
2148       if (relative)
2149         s << " / " << fancyname(variable);
2150       s << ")";
2151     } else {
2152       if (!relative)
2153         s << " / ";
2154       s << "#sqrt{2}";
2155       if (relative)
2156         s << ")";
2157     }
2158   }
2159   if (resolution && axis == 'y')
2160     s << ")";
2161   if (((!relative && !pull) || axis == 'x') && units(variable, axis) != "")
2162     s << " (" << units(variable, axis) << ")";
2163   TString result = s.str();
2164   result.ReplaceAll("#Deltaq/p_{T}", "#Delta(q/p_{T})");
2165   return result;
2166 }
2167
2168 TString latexlabel(TString variable, Char_t axis, Bool_t relative, Bool_t resolution, Bool_t pull) {
2169   TString result = axislabel(variable, axis, relative, resolution, pull);
2170   result.ReplaceAll(" (" + units(variable, axis) + ")", "");
2171   result.ReplaceAll("#", "\\").ReplaceAll("\\Delta", "\\Delta ");
2172   return result;
2173 }
2174
2175 void setAxisLabels(TH1 *p, PlotType type, TString xvar, TString yvar, Bool_t relative, Bool_t pull) {
2176   if (type == Histogram)
2177     p->SetXTitle(axislabel(yvar, 'y', relative, false, pull));
2178   if (type == ScatterPlot || type == Profile || type == Resolution || type == OrgHistogram)
2179     p->SetXTitle(axislabel(xvar, 'x'));
2180
2181   if (type == Histogram)
2182     p->SetYTitle(axislabel(yvar, 'Y', relative, false, pull));
2183   if (type == OrgHistogram)
2184     p->SetYTitle(axislabel(xvar, 'X', relative, false, pull));
2185   if (type == ScatterPlot || type == Profile)
2186     p->SetYTitle(axislabel(yvar, 'y', relative, false, pull));
2187   if (type == Resolution)
2188     p->SetYTitle(axislabel(yvar, 'y', relative, true, pull));
2189 }
2190
2191 void setAxisLabels(TMultiGraph *p, PlotType type, TString xvar, TString yvar, Bool_t relative, Bool_t pull) {
2192   if (type == Histogram)
2193     p->GetXaxis()->SetTitle(axislabel(yvar, 'y', relative, false, pull));
2194   if (type == ScatterPlot || type == Profile || type == Resolution || type == OrgHistogram)
2195     p->GetXaxis()->SetTitle(axislabel(xvar, 'x'));
2196
2197   if (type == Histogram)
2198     p->GetYaxis()->SetTitle(axislabel(yvar, 'Y', relative, false, pull));
2199   if (type == OrgHistogram)
2200     p->GetYaxis()->SetTitle(axislabel(xvar, 'X', relative, false, pull));
2201   if (type == ScatterPlot || type == Profile)
2202     p->GetYaxis()->SetTitle(axislabel(yvar, 'y', relative, false, pull));
2203   if (type == Resolution)
2204     p->GetYaxis()->SetTitle(axislabel(yvar, 'y', relative, true, pull));
2205 }
2206
2207 TString nPart(Int_t part, TString string, TString delimit, Bool_t removerest) {
2208   if (part <= 0)
2209     return "";
2210   for (int i = 1; i < part; i++)  //part-1 times
2211   {
2212     if (string.Index(delimit) < 0)
2213       return "";
2214     string.Replace(0, string.Index(delimit) + 1, "", 0);
2215   }
2216   if (string.Index(delimit) >= 0 && removerest)
2217     string.Remove(string.Index(delimit));
2218   return string;
2219 }
2220
2221 //==============
2222 //4. Axis Limits
2223 //==============
2224
2225 Double_t findStatistic(
2226     Statistic what, Int_t nFiles, TString *files, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2227   Double_t x = 0,  //if axis == 'x', var_org goes in x; if axis == 'y', Delta_var goes in x
2228       rel = 1,     //if relative, var_org goes in rel.  x is divided by rel, so you get Delta_var/var_org
2229       sigma1 = 1,  //if pull, the error for split track 1 goes in sigma1 and the error for split track 2 goes in sigma2.
2230       sigma2 = 1,  //x is divided by sqrt(sigma1^2+sigma2^2).  If !pull && axis == 'y', this divides by sqrt(2)
2231       sigmaorg = 0;  // because we want the error in one track.  sigmaorg is used when relative && pull
2232   Int_t xint = 0,
2233         xint2 = 0;  //xint is used for run number and nHits.  xint2 is used for nHits because each event has 2 values.
2234
2235   Int_t runNumber = 0;  //this is used to make sure the run number is between minrun and maxrun
2236
2237   if (axis == 'x') {
2238     sigma1 = 1 / sqrt(2);  //if axis == 'x' don't divide by sqrt(2)
2239     sigma2 = 1 / sqrt(2);
2240   }
2241
2242   Double_t totallength = 0;
2243   vector<double> xvect;
2244   Double_t result = 0;
2245   if (what == Minimum)
2246     result = 1e100;
2247   if (what == Maximum)
2248     result = -1e100;
2249
2250   stringstream sx, srel, ssigma1, ssigma2, ssigmaorg;
2251
2252   if (axis == 'y')
2253     sx << "Delta_";
2254   sx << var;
2255   if (axis == 'x' && var != "runNumber" && !var.BeginsWith("nHits"))
2256     sx << "_org";
2257   if (axis == 'x' && var.BeginsWith("nHits"))
2258     sx << "1_spl";
2259   TString variable = sx.str(), variable2 = variable;
2260   variable2.ReplaceAll("1_spl", "2_spl");
2261
2262   TString relvariable = "1";
2263   if (relative) {
2264     srel << var << "_org";
2265     relvariable = srel.str();
2266   }
2267
2268   if (pull) {
2269     ssigma1 << var << "1Err_spl";
2270     ssigma2 << var << "2Err_spl";
2271   }
2272   TString sigma1variable = ssigma1.str();
2273   TString sigma2variable = ssigma2.str();
2274
2275   if (pull && relative)
2276     ssigmaorg << var << "Err_org";
2277   TString sigmaorgvariable = ssigmaorg.str();
2278
2279   if (!relative && !pull && (variable == "Delta_dxy" || variable == "Delta_dz"))
2280     rel = 1e-4;  //it's in cm but we want um
2281   if (!relative && !pull && (variable == "Delta_phi" || variable == "Delta_theta" || variable == "Delta_qoverpt"))
2282     rel = 1e-3;  //make the axis labels manageable
2283
2284   for (Int_t j = 0; j < nFiles; j++) {
2285     if (((var == "runNumber" && what != Maximum) ? findMax(files[j], "runNumber", 'x') < 2 : false) ||
2286         files[j] == "")  //if it's MC data (run 1), the run number is meaningless
2287       continue;
2288     TFile *f = TFile::Open(files[j]);
2289     TTree *tree = (TTree *)f->Get("cosmicValidation/splitterTree");
2290     if (tree == nullptr)
2291       tree = (TTree *)f->Get("splitterTree");
2292     Int_t length = tree->GetEntries();
2293
2294     tree->SetBranchAddress("runNumber", &runNumber);
2295     if (var == "runNumber")
2296       tree->SetBranchAddress(variable, &xint);
2297     else if (var.BeginsWith("nHits")) {
2298       tree->SetBranchAddress(variable, &xint);
2299       tree->SetBranchAddress(variable2, &xint2);
2300     } else
2301       tree->SetBranchAddress(variable, &x);
2302
2303     if (relative)
2304       tree->SetBranchAddress(relvariable, &rel);
2305     if (pull) {
2306       tree->SetBranchAddress(sigma1variable, &sigma1);
2307       tree->SetBranchAddress(sigma2variable, &sigma2);
2308     }
2309     if (relative && pull)
2310       tree->SetBranchAddress(sigmaorgvariable, &sigmaorg);
2311
2312     for (Int_t i = 0; i < length; i++) {
2313       tree->GetEntry(i);
2314       if (var == "runNumber" || var.BeginsWith("nHits"))
2315         x = xint;
2316       if (var == "runNumber")
2317         runNumber = x;
2318       if (var == "phi" && x >= pi)
2319         x -= 2 * pi;
2320       if (var == "phi" && x <= -pi)
2321         x += 2 * pi;
2322       if ((runNumber < minrun && runNumber > 1) || (runNumber > maxrun && maxrun > 0))
2323         continue;
2324
2325       totallength++;
2326
2327       Double_t error;
2328       if (relative && pull)
2329         error = sqrt((sigma1 / rel) * (sigma1 / rel) + (sigma2 / rel) * (sigma2 / rel) +
2330                      (sigmaorg * x / (rel * rel)) * (sigmaorg * x / (rel * rel)));
2331       else
2332         error = sqrt(sigma1 * sigma1 + sigma2 * sigma2);  // = 1 if axis == 'x' && !pull
2333             // = sqrt(2) if axis == 'y' && !pull, so that you get the error in 1 track
2334             //       when you divide by it
2335       x /= (rel * error);
2336       if (!std::isfinite(x))  //e.g. in data with no pixels, the error occasionally comes out to be NaN
2337         continue;  //Filling a histogram with NaN is irrelevant, but here it would cause the whole result to be NaN
2338
2339       if (what == Minimum && x < result)
2340         result = x;
2341       if (what == Maximum && x > result)
2342         result = x;
2343       xvect.push_back(x);
2344       if (var.BeginsWith("nHits")) {
2345         x = xint2;
2346         if (what == Minimum && x < result)
2347           result = x;
2348         if (what == Maximum && x > result)
2349           result = x;
2350         xvect.push_back(x);
2351       }
2352     }
2353     delete f;  //automatically closes the file
2354   }
2355
2356   if (what == Minimum || what == Maximum)
2357     return result;
2358
2359   sort(xvect.begin(), xvect.end());
2360
2361   for (unsigned int i = (unsigned int)(xvect.size() * (1 - outliercut) / 2);
2362        i <= (unsigned int)(xvect.size() * (1 + outliercut) / 2 + .999);
2363        i++, totallength++)
2364     result += xvect[i];
2365
2366   result /= totallength;
2367
2368   if (what == RMS) {
2369     double average = result;
2370     result = 0;
2371     for (unsigned int i = (unsigned int)(xvect.size() * (1 - outliercut) / 2);
2372          i <= (unsigned int)(xvect.size() * (1 + outliercut) / 2 + .999);
2373          i++)
2374       result += (x - average) * (x - average);
2375     result = sqrt(result / (totallength - 1));
2376   }
2377   return result;
2378 }
2379
2380 Double_t findAverage(Int_t nFiles, TString *files, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2381   return findStatistic(Average, nFiles, files, var, axis, relative, pull);
2382 }
2383
2384 Double_t findMin(Int_t nFiles, TString *files, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2385   return findStatistic(Minimum, nFiles, files, var, axis, relative, pull);
2386 }
2387
2388 Double_t findMax(Int_t nFiles, TString *files, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2389   return findStatistic(Maximum, nFiles, files, var, axis, relative, pull);
2390 }
2391
2392 Double_t findRMS(Int_t nFiles, TString *files, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2393   return findStatistic(RMS, nFiles, files, var, axis, relative, pull);
2394 }
2395
2396 //These functions are for 1 file
2397
2398 Double_t findStatistic(Statistic what, TString file, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2399   return findStatistic(what, 1, &file, var, axis, relative, pull);
2400 }
2401
2402 Double_t findAverage(TString file, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2403   return findStatistic(Average, file, var, axis, relative, pull);
2404 }
2405
2406 Double_t findMin(TString file, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2407   return findStatistic(Minimum, file, var, axis, relative, pull);
2408 }
2409
2410 Double_t findMax(TString file, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2411   return findStatistic(Maximum, file, var, axis, relative, pull);
2412 }
2413
2414 Double_t findRMS(TString file, TString var, Char_t axis, Bool_t relative, Bool_t pull) {
2415   return findStatistic(RMS, file, var, axis, relative, pull);
2416 }
2417
2418 //This puts the axis limits that should be used for trackSplitPlot in min and max.
2419 //Default axis limits are defined for pt, qoverpt, dxy, dz, theta, eta, and phi.
2420 //For run number and nHits, the minimum and maximum are used.
2421 //For any other variable, average +/- 5*rms are used.
2422 //To use this instead of the default values, just comment out the part that says [else] if (var == "?") {min = ?; max = ?;}
2423
2424 void axislimits(Int_t nFiles,
2425                 TString *files,
2426                 TString var,
2427                 Char_t axis,
2428                 Bool_t relative,
2429                 Bool_t pull,
2430                 Double_t &min,
2431                 Double_t &max,
2432                 Double_t &bins) {
2433   bool pixel = subdetector.Contains("PIX");
2434   if (axis == 'x') {
2435     if (var == "pt") {
2436       min = 5;
2437       max = 100;
2438       bins = 38;
2439     } else if (var == "qoverpt") {
2440       min = -.35;
2441       max = .35;
2442       bins = 35;
2443     } else if (var == "dxy") {
2444       min = -100;
2445       max = 100;
2446       if (pixel) {
2447         min = -10;
2448         max = 10;
2449       }
2450       bins = 20;
2451     } else if (var == "dz") {
2452       min = -250;
2453       max = 250;
2454       if (pixel) {
2455         min = -25;
2456         max = 25;
2457       }
2458       bins = 25;
2459     } else if (var == "theta") {
2460       min = .5;
2461       max = 2.5;
2462       bins = 40;
2463     } else if (var == "eta") {
2464       min = -1.2;
2465       max = 1.2;
2466       bins = 40;
2467     } else if (var == "phi") {
2468       min = -3;
2469       max = 0;
2470       bins = 30;
2471     } else if (var == "runNumber" || var.BeginsWith("nHits")) {
2472       min = findMin(nFiles, files, var, 'x') - .5;
2473       max = findMax(nFiles, files, var, 'x') + .5;
2474       bins = max - min;
2475     } else {
2476       cout << "No x axis limits for " << var << ".  Using average +/- 5*rms" << endl;
2477       Double_t average = findAverage(nFiles, files, var, 'x');
2478       Double_t rms = findRMS(nFiles, files, var, 'x');
2479       max = TMath::Min(average + 5 * rms, findMax(nFiles, files, var, 'x'));
2480       min = TMath::Max(average - 5 * rms, findMin(nFiles, files, var, 'x'));
2481       bins = 50;
2482     }
2483   }
2484   if (axis == 'y') {
2485     if (pull) {
2486       min = -5;
2487       max = 5;
2488       bins = 40;
2489     } else if (var == "pt" && relative) {
2490       min = -.06;
2491       max = .06;
2492       bins = 30;
2493     } else if (var == "pt" && !relative) {
2494       min = -.8;
2495       max = .8;
2496       bins = 40;
2497     } else if (var == "qoverpt") {
2498       min = -2.5;
2499       max = 2.5;
2500       bins = 50;
2501     } else if (var == "dxy") {
2502       min = -1250;
2503       max = 1250;
2504       if (pixel) {
2505         min = -125;
2506         max = 125;
2507       }
2508       bins = 50;
2509     } else if (var == "dz") {
2510       min = -2000;
2511       max = 2000;
2512       if (pixel) {
2513         min = -200;
2514         max = 200;
2515       }
2516       bins = 40;
2517     } else if (var == "theta") {
2518       min = -10;
2519       max = 10;
2520       if (pixel) {
2521         min = -5;
2522         max = 5;
2523       }
2524       bins = 50;
2525     } else if (var == "eta") {
2526       min = -.007;
2527       max = .007;
2528       if (pixel) {
2529         min = -.003;
2530         max = .003;
2531       }
2532       bins = 30;
2533     } else if (var == "phi") {
2534       min = -2;
2535       max = 2;
2536       bins = 40;
2537     } else {
2538       cout << "No y axis limits for " << var << ".  Using average +/- 5 * rms." << endl;
2539       Double_t average = 0 /*findAverage(nFiles,files,var,'y',relative,pull)*/;
2540       Double_t rms = findRMS(nFiles, files, var, 'y', relative, pull);
2541       min = TMath::Max(TMath::Max(-TMath::Abs(average) - 5 * rms, findMin(nFiles, files, var, 'y', relative, pull)),
2542                        -findMax(nFiles, files, var, 'y', relative, pull));
2543       max = -min;
2544       bins = 50;
2545     }
2546   }
2547 }
2548
2549 //===============
2550 //5. Place Legend
2551 //===============
2552
2553 Double_t placeLegend(
2554     TLegend *l, Double_t width, Double_t height, Double_t x1min, Double_t y1min, Double_t x2max, Double_t y2max) {
2555   for (int i = legendGrid; i >= 0; i--) {
2556     for (int j = legendGrid; j >= 0; j--) {
2557       Double_t x1 = x1min * (1 - (double)i / legendGrid) + (x2max - width) * (double)i / legendGrid - margin * width;
2558       Double_t y1 = y1min * (1 - (double)j / legendGrid) + (y2max - height) * (double)j / legendGrid - margin * height;
2559       Double_t x2 = x1 + (1 + 2 * margin) * width;
2560       Double_t y2 = y1 + (1 + 2 * margin) * height;
2561       if (fitsHere(l, x1, y1, x2, y2)) {
2562         x1 += margin * width;
2563         y1 += margin * height;
2564         x2 -= margin * width;
2565         y2 -= margin * height;
2566         l->SetX1(x1);
2567         l->SetY1(y1);
2568         l->SetX2(x2);
2569         l->SetY2(y2);
2570         return y2max;
2571       }
2572     }
2573   }
2574   Double_t newy2max = y2max + increaseby * (y2max - y1min);
2575   Double_t newheight = height * (newy2max - y1min) / (y2max - y1min);
2576   return placeLegend(l, width, newheight, x1min, y1min, x2max, newy2max);
2577 }
2578
2579 Bool_t fitsHere(TLegend *l, Double_t x1, Double_t y1, Double_t x2, Double_t y2) {
2580   Bool_t fits = true;
2581   TList *list = l->GetListOfPrimitives();
2582   for (Int_t k = 0; list->At(k) != nullptr && fits; k++) {
2583     TObject *obj = ((TLegendEntry *)(list->At(k)))->GetObject();
2584     if (obj == nullptr)
2585       continue;
2586     TClass *cl = obj->IsA();
2587
2588     //Histogram, drawn as a histogram
2589     if (cl->InheritsFrom("TH1") && !cl->InheritsFrom("TH2") && !cl->InheritsFrom("TH3") && cl != TProfile::Class() &&
2590         ((TH1 *)obj)->GetMarkerColor() == kWhite) {
2591       Int_t where = 0;
2592       TH1 *h = (TH1 *)obj;
2593       for (Int_t i = 1; i <= h->GetNbinsX() && fits; i++) {
2594         if (h->GetBinLowEdge(i) + h->GetBinWidth(i) < x1)
2595           continue;  //to the left of the legend
2596         if (h->GetBinLowEdge(i) > x2)
2597           continue;  //to the right of the legend
2598         if (h->GetBinContent(i) > y1 && h->GetBinContent(i) < y2)
2599           fits = false;  //inside the legend
2600         if (h->GetBinContent(i) < y1) {
2601           if (where == 0)
2602             where = -1;  //below the legend
2603           if (where == 1)
2604             fits = false;  //a previous bin was above it so there's a vertical line through it
2605         }
2606         if (h->GetBinContent(i) > y2) {
2607           if (where == 0)
2608             where = 1;  //above the legend
2609           if (where == -1)
2610             fits = false;  //a previous bin was below it so there's a vertical line through it
2611         }
2612       }
2613       continue;
2614     }
2615     //Histogram, drawn with Draw("P")
2616     else if (cl->InheritsFrom("TH1") && !cl->InheritsFrom("TH2") && !cl->InheritsFrom("TH3") && cl != TProfile::Class())
2617     //Probably TProfile would be the same but I haven't tested it
2618     {
2619       TH1 *h = (TH1 *)obj;
2620       for (Int_t i = 1; i <= h->GetNbinsX() && fits; i++) {
2621         if (h->GetBinLowEdge(i) + h->GetBinWidth(i) / 2 < x1)
2622           continue;
2623         if (h->GetBinLowEdge(i) > x2)
2624           continue;
2625         if (h->GetBinContent(i) > y1 && h->GetBinContent(i) < y2)
2626           fits = false;
2627         if (h->GetBinContent(i) + h->GetBinError(i) > y2 && h->GetBinContent(i) - h->GetBinError(i) < y2)
2628           fits = false;
2629         if (h->GetBinContent(i) + h->GetBinError(i) > y1 && h->GetBinContent(i) - h->GetBinError(i) < y1)
2630           fits = false;
2631       }
2632     } else if (cl->InheritsFrom("TF1") && !cl->InheritsFrom("TF2")) {
2633       TF1 *f = (TF1 *)obj;
2634       Double_t max = f->GetMaximum(x1, x2);
2635       Double_t min = f->GetMinimum(x1, x2);
2636       if (min < y2 && max > y1)
2637         fits = false;
2638     }
2639     // else if (cl->InheritsFrom(...... add more objects here
2640     else {
2641       cout << "Don't know how to place the legend around objects of type " << obj->ClassName() << "." << endl
2642            << "Add this class into fitsHere() if you want it to work properly." << endl
2643            << "The legend will still be placed around any other objects." << endl;
2644     }
2645   }
2646   return fits;
2647 }