File indexing completed on 2021-02-14 12:45:47
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "trackSplitPlot.h"
0012 #include "Alignment/OfflineValidation/macros/TkAlStyle.cc"
0013
0014
0015
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);
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
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] == "")
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,
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;
0174 if (!relative && !pull && (yvar == "phi" || yvar == "theta" || yvar == "qoverpt"))
0175 rel = 1e-3;
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)
0194 {
0195 yvariable.ReplaceAll("Delta_","d");
0196 yvariable.Append("_spl");
0197 tree->SetBranchAddress(yvariable,&y);
0198 }
0199 }
0200 if (relative && xvar != yvar)
0201 tree->SetBranchAddress(relvariable,&rel);
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;
0216
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))
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);
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
0255
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);
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
0284
0285
0286
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) {
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
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;
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]);
0510
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();
0587 double padheight = gPad->GetUymax() - gPad->GetUymin();
0588
0589
0590
0591
0592 double x = padheight*legendfraction / (1-legendfraction) * 1.5;
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
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
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)
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
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);
0692
0693 }
0694 }
0695
0696
0697
0698
0699
0700
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
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
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;
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
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->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
0825
0826
0827
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;
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];
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;
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 = yaxislabel;
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
0990
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);
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
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);
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
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);
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);
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
1111
1112
1113
1114
1115
1116
1117
1118
1119
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
1137 TString parametername = "";
1138 Int_t parameter = 9999;
1139
1140
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
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
1190
1191
1192
1193 parameters = tempParameters;
1194 parameternames = tempParameterNames;
1195 functionname = "#Deltad_{xy}=-Asin(2#phi+B)";
1196
1197 }
1198 if (xvar == "phi" && yvar == "dxy" && !resolution && pull)
1199 {
1200 f = new TF1("sine","[0]*sin([1]*x-[2])");
1201
1202
1203 f->FixParameter(1,-2);
1204
1205 nParameters = 2;
1206 Int_t tempParameters[2] = {0,2};
1207 TString tempParameterNames[2] = {"A","B"};
1208
1209
1210
1211
1212 parameters = tempParameters;
1213 parameternames = tempParameterNames;
1214
1215 functionname = "#Deltad_{xy}/#delta(#Deltad_{xy})=-Asin(2#phi+B)";
1216
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
1229
1230
1231
1232
1233
1234
1235
1236
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])) )");
1305
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
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
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
1376
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
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++)
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;
1460
1461
1462 if (!matrix[x][y]) continue;
1463
1464 if (xvariables[x] == "" && yvariables[y] == "") continue;
1465
1466 Int_t nPlots = nFiles+4;
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]);
1476 plotnames[i].ReplaceAll(" ","");
1477 }
1478
1479 plotnames.push_back("");
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");
1488
1489 plotnames.push_back("");
1490 plotnames.push_back("");
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
1542
1543
1544
1545
1546
1547
1548
1549
1550 }
1551
1552 if (xvariables[x] != "" && yvariables[y] != "")
1553 {
1554
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
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
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";
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
1658
1659
1660
1661
1662
1663
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";
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
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";
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
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
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";
1827 if (axis == 'x')
1828 return "cm";
1829 }
1830 if (variable == "qoverpt")
1831 {
1832 if (axis == 'y')
1833 return "#times10^{-3}e/GeV";
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
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++)
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
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,
1986 rel = 1,
1987 sigma1 = 1,
1988 sigma2 = 1,
1989 sigmaorg = 0;
1990 Int_t xint = 0, xint2 = 0;
1991
1992 Int_t runNumber = 0;
1993
1994 if (axis == 'x')
1995 {
1996 sigma1 = 1/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;
2040 if (!relative && !pull && (variable == "Delta_phi" || variable == "Delta_theta" || variable == "Delta_qoverpt"))
2041 rel = 1e-3;
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] == "")
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);
2094
2095
2096 x /= (rel * error);
2097 if (!std::isfinite(x))
2098 continue;
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;
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
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
2191
2192
2193
2194
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 ;
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
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
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;
2411 if (h->GetBinLowEdge(i) > x2) continue;
2412 if (h->GetBinContent(i) > y1 && h->GetBinContent(i) < y2) fits = false;
2413 if (h->GetBinContent(i) < y1)
2414 {
2415 if (where == 0) where = -1;
2416 if (where == 1) fits = false;
2417 }
2418 if (h->GetBinContent(i) > y2)
2419 {
2420 if (where == 0) where = 1;
2421 if (where == -1) fits = false;
2422 }
2423 }
2424 continue;
2425 }
2426
2427 else if (cl->InheritsFrom("TH1") && !cl->InheritsFrom("TH2") && !cl->InheritsFrom("TH3")
2428 && cl != TProfile::Class())
2429
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
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 }