Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:47

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