File indexing completed on 2023-11-18 03:07:10
0001 #include "Alignment/OfflineValidation/interface/PlotAlignmentValidation.h"
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 PlotAlignmentValidation::PlotAlignmentValidation(bool bigtext) : bigtext_(bigtext) {
0015 setOutputDir(".");
0016 setTreeBaseDir();
0017 sourcelist = nullptr;
0018
0019 moreThanOneSource = false;
0020 useFit_ = false;
0021
0022
0023 TGaxis::SetMaxDigits(4);
0024
0025
0026
0027
0028
0029 TH1::StatOverflows(kTRUE);
0030
0031
0032 legendOptions(TkAlStyle::legendoptions);
0033 }
0034
0035
0036
0037
0038
0039 PlotAlignmentValidation::PlotAlignmentValidation(
0040 const char* inputFile, std::string legendName, int lineColor, int lineStyle, bool bigtext)
0041 : PlotAlignmentValidation(bigtext) {
0042 loadFileList(inputFile, legendName, lineColor, lineStyle);
0043 }
0044
0045
0046
0047
0048
0049 PlotAlignmentValidation::~PlotAlignmentValidation() {
0050 for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0051 delete (*it);
0052 }
0053
0054 delete sourcelist;
0055 }
0056
0057
0058
0059
0060
0061
0062
0063 void PlotAlignmentValidation::openSummaryFile() {
0064 if (!openedsummaryfile) {
0065 openedsummaryfile = true;
0066 summaryfile.open(outputDir + "/" + summaryfilename + ".txt");
0067
0068 rootsummaryfile = new TFile(outputDir + "/" + summaryfilename + ".root", "RECREATE");
0069
0070 for (auto vars : sourceList) {
0071 summaryfile << "\t" << vars->getName();
0072 }
0073 summaryfile << "\tformat={}\n";
0074 } else {
0075
0076 if (!rootsummaryfile->IsOpen())
0077 rootsummaryfile->Open(outputDir + "/" + summaryfilename + ".root", "UPDATE");
0078 }
0079 }
0080
0081
0082
0083
0084
0085 void PlotAlignmentValidation::loadFileList(const char* inputFile, std::string legendName, int lineColor, int lineStyle) {
0086 if (openedsummaryfile) {
0087 std::cout << "Can't load a root file after opening the summary file!" << std::endl;
0088 assert(0);
0089 }
0090 sourceList.push_back(new TkOfflineVariables(inputFile, treeBaseDir, legendName, lineColor, lineStyle));
0091 }
0092
0093
0094
0095
0096
0097 void PlotAlignmentValidation::useFitForDMRplots(bool usefit) { useFit_ = usefit; }
0098
0099
0100
0101
0102
0103
0104 int PlotAlignmentValidation::numberOfLayers(int phase, int subdetector) {
0105 switch (phase) {
0106 case 0:
0107 switch (subdetector) {
0108 case 1:
0109 return 3;
0110 case 2:
0111 return 2;
0112 case 3:
0113 return 4;
0114 case 4:
0115 return 3;
0116 case 5:
0117 return 6;
0118 case 6:
0119 return 9;
0120 default:
0121 assert(false);
0122 }
0123 case 1:
0124 switch (subdetector) {
0125 case 1:
0126 return 4;
0127 case 2:
0128 return 3;
0129 case 3:
0130 return 4;
0131 case 4:
0132 return 3;
0133 case 5:
0134 return 6;
0135 case 6:
0136 return 9;
0137 default:
0138 assert(false);
0139 }
0140 default:
0141 assert(false);
0142 }
0143 return 0;
0144 }
0145
0146
0147
0148
0149
0150 int PlotAlignmentValidation::maxNumberOfLayers(int subdetector) {
0151 int result = 0;
0152 for (auto it = sourceList.begin(); it != sourceList.end(); ++it) {
0153 result = std::max(result, numberOfLayers((*it)->getPhase(), subdetector));
0154 }
0155 return result;
0156 }
0157
0158
0159
0160
0161
0162 void PlotAlignmentValidation::legendOptions(TString options) {
0163 showMean_ = false;
0164 showRMS_ = false;
0165 showMeanError_ = false;
0166 showRMSError_ = false;
0167 showModules_ = false;
0168 showUnderOverFlow_ = false;
0169 options.ReplaceAll(" ", "").ToLower();
0170 if (options.Contains("mean") || options.Contains("all"))
0171 showMean_ = true;
0172 if (options.Contains("meanerror") || options.Contains("all"))
0173 showMeanError_ = true;
0174 if (options.Contains("rms") || options.Contains("all"))
0175 showRMS_ = true;
0176 if (options.Contains("rmserror") || options.Contains("all"))
0177 showRMSError_ = true;
0178 if (options.Contains("modules") || options.Contains("all"))
0179 showModules_ = true;
0180 if (options.Contains("under") || options.Contains("over") || options.Contains("outside") || options.Contains("all"))
0181 showUnderOverFlow_ = true;
0182
0183 twolines_ = (showUnderOverFlow_ && (showMean_ + showMeanError_ + showRMS_ + showRMSError_ >= 1) && bigtext_);
0184 }
0185
0186
0187
0188
0189
0190 void PlotAlignmentValidation::setOutputDir(std::string dir) {
0191 if (openedsummaryfile) {
0192 std::cout << "Can't set the output dir after opening the summary file!" << std::endl;
0193 assert(0);
0194 }
0195 outputDir = dir;
0196 gSystem->mkdir(outputDir.data(), true);
0197 }
0198
0199
0200
0201
0202
0203 void PlotAlignmentValidation::plotSubDetResiduals(bool plotNormHisto, unsigned int subDetId) {
0204 gStyle->SetOptStat(11111);
0205 gStyle->SetOptFit(0000);
0206
0207 TCanvas* c = new TCanvas("c", "c");
0208 c->SetTopMargin(0.15);
0209 TString histoName = "";
0210 if (plotNormHisto) {
0211 histoName = "h_NormXprime";
0212 } else
0213 histoName = "h_Xprime_";
0214 switch (subDetId) {
0215 case 1:
0216 histoName += "TPBBarrel_0";
0217 break;
0218 case 2:
0219 histoName += "TPEendcap_1";
0220 break;
0221 case 3:
0222 histoName += "TPEendcap_2";
0223 break;
0224 case 4:
0225 histoName += "TIBBarrel_0";
0226 break;
0227 case 5:
0228 histoName += "TIDEndcap_1";
0229 break;
0230 case 6:
0231 histoName += "TIDEndcap_2";
0232 break;
0233 case 7:
0234 histoName += "TOBBarrel_3";
0235 break;
0236 case 8:
0237 histoName += "TECEndcap_4";
0238 break;
0239 case 9:
0240 histoName += "TECEndcap_5";
0241 break;
0242 }
0243 int tmpcounter = 0;
0244 TH1* sumHisto = nullptr;
0245 for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0246 if (tmpcounter == 0) {
0247 TFile* f = (*it)->getFile();
0248 sumHisto = (TH1*)f->FindKeyAny(histoName)->ReadObj();
0249 sumHisto->SetLineColor(tmpcounter + 1);
0250 sumHisto->SetLineStyle(tmpcounter + 1);
0251 sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
0252 sumHisto->Draw();
0253
0254
0255
0256
0257
0258
0259
0260 tmpcounter++;
0261 } else {
0262 sumHisto = (TH1*)(*it)->getFile()->FindObjectAny(histoName);
0263 sumHisto->SetLineColor(tmpcounter + 1);
0264 sumHisto->SetLineStyle(tmpcounter + 1);
0265 sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
0266
0267
0268 c->Update();
0269 tmpcounter++;
0270 }
0271 TObject* statObj = sumHisto->GetListOfFunctions()->FindObject("stats");
0272 if (statObj && statObj->InheritsFrom(TPaveStats::Class())) {
0273 TPaveStats* stats = static_cast<TPaveStats*>(statObj);
0274 stats->SetLineColor(tmpcounter + 1);
0275 stats->SetTextColor(tmpcounter + 1);
0276 stats->SetFillColor(10);
0277 stats->SetX1NDC(0.91 - tmpcounter * 0.1);
0278 stats->SetX2NDC(0.15);
0279 stats->SetY1NDC(1);
0280 stats->SetY2NDC(0.10);
0281 sumHisto->Draw("sames");
0282 }
0283 }
0284
0285 char PlotName[1000];
0286 sprintf(PlotName, "%s/%s.png", outputDir.c_str(), histoName.Data());
0287 c->Print(PlotName);
0288 sprintf(PlotName, "%s/%s.eps", outputDir.c_str(), histoName.Data());
0289 c->Print(PlotName);
0290 sprintf(PlotName, "%s/%s.pdf", outputDir.c_str(), histoName.Data());
0291 c->Print(PlotName);
0292 sprintf(PlotName, "%s/%s.root", outputDir.c_str(), histoName.Data());
0293 c->Print(PlotName);
0294
0295
0296 }
0297
0298
0299 void PlotAlignmentValidation::plotHitMaps() {
0300
0301
0302 TCanvas* c = new TCanvas("c", "c", 1200, 400);
0303 c->Divide(3, 1);
0304
0305
0306
0307
0308
0309 std::string histName_ = "Entriesprofile";
0310 c->cd(1);
0311 TTree* tree = (*sourceList.begin())->getTree();
0312 tree->Draw("entries:posR:posZ", "", "COLZ2Prof");
0313 c->cd(2);
0314 tree->Draw("entries:posY:posX", "", "COLZ2Prof");
0315 c->cd(3);
0316 tree->Draw("entries:posR:posPhi", "", "COLZ2Prof");
0317
0318 char PlotName[1000];
0319 sprintf(PlotName, "%s/%s.png", outputDir.c_str(), histName_.c_str());
0320 c->Print(PlotName);
0321 sprintf(PlotName, "%s/%s.eps", outputDir.c_str(), histName_.c_str());
0322 c->Print(PlotName);
0323 sprintf(PlotName, "%s/%s.pdf", outputDir.c_str(), histName_.c_str());
0324 c->Print(PlotName);
0325 sprintf(PlotName, "%s/%s.root", outputDir.c_str(), histName_.c_str());
0326 c->Print(PlotName);
0327
0328 c->Close();
0329
0330 }
0331
0332
0333 void PlotAlignmentValidation::plotOutlierModules(const char* outputFileName,
0334 std::string plotVariable,
0335 float plotVariable_cut,
0336 int unsigned minHits) {
0337 Int_t counter = 0;
0338
0339 gStyle->SetOptStat(111111);
0340 gStyle->SetStatY(0.9);
0341
0342
0343 TCanvas* c1 = new TCanvas("canv", "canv", 800, 500);
0344 outputFile = outputDir + '/' + outputFileName;
0345 c1->Print((outputFile + '[').Data());
0346
0347 c1->Divide(2, 1);
0348
0349 TTree* tree = (*sourceList.begin())->getTree();
0350 TkOffTreeVariables* treeMem = nullptr;
0351 tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
0352
0353 Long64_t nentries = tree->GetEntriesFast();
0354
0355 for (Long64_t i = 0; i < nentries; i++) {
0356 tree->GetEntry(i);
0357 float var = 0;
0358 if (plotVariable == "chi2PerDofX")
0359 var = treeMem->chi2PerDofX;
0360 else if (plotVariable == "chi2PerDofY")
0361 var = treeMem->chi2PerDofY;
0362 else if (plotVariable == "fitMeanX")
0363 var = treeMem->fitMeanX;
0364 else if (plotVariable == "fitMeanY")
0365 var = treeMem->fitMeanY;
0366 else if (plotVariable == "fitSigmaX")
0367 var = treeMem->fitSigmaX;
0368 else if (plotVariable == "fitSigmaY")
0369 var = treeMem->fitSigmaY;
0370 else {
0371 std::cout << "There is no variable " << plotVariable << " included in the tree." << std::endl;
0372 break;
0373 }
0374
0375
0376
0377
0378 if (var > plotVariable_cut && treeMem->entries > minHits) {
0379 TFile* f = (*sourceList.begin())->getFile();
0380
0381 if (f->FindKeyAny(treeMem->histNameX.c_str()) != nullptr) {
0382 TH1* h =
0383 (TH1*)f->FindKeyAny(treeMem->histNameX.c_str())->ReadObj();
0384 gStyle->SetOptFit(0111);
0385 std::cout << "hist name " << h->GetName() << std::endl;
0386
0387 TString path = (char*)strstr(gDirectory->GetPath(), "TrackerOfflineValidation");
0388
0389
0390 if (h)
0391 std::cout << h->GetEntries() << std::endl;
0392
0393
0394 c1->cd(0);
0395 TPaveText* text = new TPaveText(0, 0.95, 0.99, 0.99);
0396 text->AddText(path);
0397 text->SetFillColor(0);
0398 text->SetShadowColor(0);
0399 text->SetBorderSize(0);
0400 text->Draw();
0401
0402
0403 c1->cd(1);
0404 TPad* subpad = (TPad*)c1->GetPad(1);
0405 subpad->SetPad(0, 0, 0.5, 0.94);
0406 h->Draw();
0407
0408
0409 h = (TH1*)f->FindObjectAny(treeMem->histNameNormX.c_str());
0410 if (h)
0411 std::cout << h->GetEntries() << std::endl;
0412 c1->cd(2);
0413 TPad* subpad2 = (TPad*)c1->GetPad(2);
0414 subpad2->SetPad(0.5, 0, 0.99, 0.94);
0415 h->Draw();
0416
0417 c1->Print(outputFile);
0418 counter++;
0419 } else {
0420 std::cout << "There are no residual histograms on module level stored!" << std::endl;
0421 std::cout << "Please make sure that moduleLevelHistsTransient = cms.bool(False) in the validation job!"
0422 << std::endl;
0423 break;
0424 }
0425 }
0426 }
0427 c1->Print((outputFile + "]").Data());
0428 if (counter == 0)
0429 std::cout << "no bad modules found" << std::endl;
0430
0431
0432
0433
0434
0435
0436 }
0437
0438
0439
0440
0441
0442 TList* PlotAlignmentValidation::getTreeList() {
0443 TList* treeList = new TList();
0444 TFile* first_source = (TFile*)sourcelist->First();
0445 std::cout << first_source->GetName() << std::endl;
0446 TDirectoryFile* d = (TDirectoryFile*)first_source->Get(treeBaseDir.c_str());
0447 treeList->Add((TTree*)(*d).Get("TkOffVal"));
0448
0449 if (moreThanOneSource == true) {
0450 TFile* nextsource = (TFile*)sourcelist->After(first_source);
0451 while (nextsource) {
0452 std::cout << nextsource->GetName() << std::endl;
0453 d = (TDirectoryFile*)nextsource->Get("TrackerOfflineValidation");
0454
0455 treeList->Add((TTree*)(*d).Get("TkOffVal"));
0456
0457 nextsource = (TFile*)sourcelist->After(nextsource);
0458 }
0459 }
0460 return treeList;
0461 }
0462
0463
0464 void PlotAlignmentValidation::setTreeBaseDir(std::string dir) { treeBaseDir = dir; }
0465
0466
0467 void PlotAlignmentValidation::plotSurfaceShapes(const std::string& options, const std::string& residType) {
0468 std::cout << "-------- plotSurfaceShapes called with " << options << std::endl;
0469 if (options == "none")
0470 return;
0471 else if (options == "coarse") {
0472 plotSS("subdet=1");
0473 plotSS("subdet=2");
0474 plotSS("subdet=3");
0475 plotSS("subdet=4");
0476 plotSS("subdet=5");
0477 plotSS("subdet=6");
0478 }
0479
0480 else
0481 plotSS(options, residType);
0482
0483 return;
0484 }
0485
0486
0487 void PlotAlignmentValidation::plotSS(const std::string& options, const std::string& residType) {
0488 if (residType.empty()) {
0489 plotSS(options, "ResXvsXProfile");
0490 plotSS(options, "ResXvsYProfile");
0491 return;
0492 }
0493
0494 int bkperrorx = gStyle->GetErrorX();
0495 gStyle->SetErrorX(1);
0496
0497 int plotLayerN = 0;
0498
0499
0500 bool plotLayers = false;
0501
0502 int plotSubDetN = 0;
0503
0504 TRegexp layer_re("layer=[0-9]+");
0505 Ssiz_t index, len;
0506 if (options.find("layers") != std::string::npos) {
0507 plotLayers = true;
0508 }
0509 if ((index = layer_re.Index(options, &len)) != -1) {
0510 if (plotLayers) {
0511 std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
0512 } else {
0513 std::string substr = options.substr(index + 6, len - 6);
0514 plotLayerN = atoi(substr.c_str());
0515 }
0516 }
0517
0518 TRegexp subdet_re("subdet=[1-6]+");
0519 if ((index = subdet_re.Index(options, &len)) != -1) {
0520 std::string substr = options.substr(index + 7, len - 7);
0521 plotSubDetN = atoi(substr.c_str());
0522 }
0523
0524 gStyle->SetOptStat(0);
0525
0526 TCanvas c("canv", "canv");
0527
0528
0529
0530
0531 for (int iSubDet = 1; iSubDet <= 6; ++iSubDet) {
0532
0533 bool isTEC = (iSubDet == 6);
0534
0535
0536 if (plotSubDetN != 0 && iSubDet != plotSubDetN)
0537 continue;
0538
0539
0540
0541
0542 if (plotLayerN > maxNumberOfLayers(iSubDet)) {
0543 continue;
0544 }
0545
0546 int minlayer = plotLayers ? 1 : plotLayerN;
0547 int maxlayer = plotLayers ? maxNumberOfLayers(iSubDet) : plotLayerN;
0548
0549 int maxlayerphase0 = plotLayers ? numberOfLayers(0, iSubDet) : plotLayerN;
0550
0551 for (int layer = minlayer; layer <= maxlayer; layer++) {
0552
0553 for (int iTEC = 0; iTEC < 2; iTEC++) {
0554 if (!isTEC && iTEC == 0)
0555 continue;
0556
0557 char selection[1000];
0558 if (!isTEC) {
0559 if (layer == 0)
0560 sprintf(selection, "subDetId==%d", iSubDet);
0561 else
0562 sprintf(selection, "subDetId==%d && layer == %d", iSubDet, layer);
0563 } else {
0564 if (iTEC == 0)
0565 sprintf(selection, "subDetId==%d && ring <= 4", iSubDet);
0566 else
0567 sprintf(selection, "subDetId==%d && ring > 4", iSubDet);
0568 }
0569
0570
0571
0572 TString subDetName;
0573 switch (iSubDet) {
0574 case 1:
0575 subDetName = "BPIX";
0576 break;
0577 case 2:
0578 subDetName = "FPIX";
0579 break;
0580 case 3:
0581 subDetName = "TIB";
0582 break;
0583 case 4:
0584 subDetName = "TID";
0585 break;
0586 case 5:
0587 subDetName = "TOB";
0588 break;
0589 case 6:
0590 subDetName = "TEC";
0591 break;
0592 }
0593
0594 TString secondline = "";
0595 if (layer != 0) {
0596
0597 if (iSubDet == 4 || iSubDet == 6)
0598 secondline = "disc ";
0599 else {
0600 secondline = "layer ";
0601 }
0602 secondline += Form("%d", layer);
0603 secondline += " ";
0604 }
0605 if (isTEC && iTEC == 0)
0606 secondline += TString("R1-4");
0607 if (isTEC && iTEC > 0)
0608 secondline += TString("R5-7");
0609
0610
0611 TLegend* legend = nullptr;
0612
0613 THStack* hs = addHists(selection, residType, &legend, false, layer <= maxlayerphase0);
0614 if (!hs || hs->GetHists() == nullptr || hs->GetHists()->GetSize() == 0) {
0615 std::cout << "No histogram for " << subDetName << ", perhaps not enough data? Creating default histogram."
0616 << std::endl;
0617 if (hs == nullptr)
0618 hs = new THStack("hstack", "");
0619
0620 TProfile* defhist = new TProfile("defhist", "Empty default histogram", 100, -1, 1, -1, 1);
0621 hs->Add(defhist);
0622 hs->Draw();
0623 } else {
0624 hs->Draw("nostack PE");
0625 modifySSHistAndLegend(hs, legend);
0626 legend->Draw();
0627 setTitleStyle(*hs, "", "", iSubDet, true, secondline);
0628
0629
0630 TH1* firstHisto = (TH1*)hs->GetHists()->First();
0631 TString xName = firstHisto->GetXaxis()->GetTitle();
0632 TString yName = firstHisto->GetYaxis()->GetTitle();
0633 hs->GetHistogram()->GetXaxis()->SetTitleColor(kBlack);
0634 hs->GetHistogram()->GetXaxis()->SetTitle(xName);
0635 hs->GetHistogram()->GetYaxis()->SetTitleColor(kBlack);
0636
0637 yName.ReplaceAll("cm", "#mum");
0638 hs->GetHistogram()->GetYaxis()->SetTitle(yName);
0639 }
0640
0641
0642 std::ostringstream plotName;
0643 plotName << outputDir << "/SurfaceShape_" << subDetName << "_";
0644 plotName << residType;
0645 if (layer != 0) {
0646 plotName << "_";
0647
0648 if (iSubDet == 4 || iSubDet == 6)
0649 plotName << "disc";
0650 else {
0651 plotName << "layer";
0652 }
0653 plotName << layer;
0654 }
0655 if (isTEC && iTEC == 0)
0656 plotName << "_"
0657 << "R1-4";
0658 if (isTEC && iTEC > 0)
0659 plotName << "_"
0660 << "R5-7";
0661
0662
0663 c.Update();
0664 c.Print((plotName.str() + ".png").c_str());
0665 c.Print((plotName.str() + ".eps").c_str());
0666 c.Print((plotName.str() + ".pdf").c_str());
0667
0668
0669 TFile f((plotName.str() + ".root").c_str(), "recreate");
0670 c.Write();
0671 f.Close();
0672
0673 delete legend;
0674 delete hs;
0675 }
0676 }
0677 }
0678 gStyle->SetErrorX(bkperrorx);
0679
0680 return;
0681 }
0682
0683
0684
0685
0686
0687 void PlotAlignmentValidation::plotDMR(const std::string& variable,
0688 Int_t minHits,
0689 const std::string& options,
0690 const std::string& filterName) {
0691
0692
0693
0694
0695 std::size_t findres = variable.find(',');
0696 if (findres != std::string::npos) {
0697 std::string substring1 = variable.substr(0, findres);
0698 std::string substring2 = variable.substr(findres + 1, std::string::npos);
0699 plotDMR(substring1, minHits, options, filterName);
0700 plotDMR(substring2, minHits, options, filterName);
0701 return;
0702 }
0703
0704
0705
0706 if (variable == "mean" || variable == "median" || variable == "meanNorm" || variable == "rms" ||
0707 variable == "rmsNorm") {
0708 plotDMR(variable + "X", minHits, options, filterName);
0709 plotDMR(variable + "Y", minHits, options, filterName);
0710 return;
0711 }
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721 TRegexp layer_re("layer=[0-9]+");
0722 bool plotPlain = false, plotSplits = false, plotLayers = false;
0723 int plotLayerN = 0;
0724 Ssiz_t index, len;
0725 if (options.find("plain") != std::string::npos) {
0726 plotPlain = true;
0727 }
0728 if (options.find("split") != std::string::npos) {
0729 plotSplits = true;
0730 }
0731 if (options.find("layers") != std::string::npos) {
0732 plotLayers = true;
0733 }
0734 if ((index = layer_re.Index(options, &len)) != -1) {
0735 if (plotLayers) {
0736 std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
0737 } else {
0738 std::string substr = options.substr(index + 6, len - 6);
0739 plotLayerN = atoi(substr.c_str());
0740 }
0741 }
0742
0743
0744
0745 if (!plotPlain && !plotSplits) {
0746 plotPlain = true;
0747 }
0748
0749
0750
0751 static bool plotSplitsFor[6] = {true, true, true, false, true, false};
0752
0753 DMRPlotInfo plotinfo;
0754
0755 gStyle->SetOptStat(0);
0756
0757 TCanvas c("canv", "canv");
0758
0759 plotinfo.variable = variable;
0760 plotinfo.minHits = minHits;
0761 plotinfo.plotPlain = plotPlain;
0762 plotinfo.plotLayers = plotLayers;
0763 plotinfo.filterName = filterName;
0764
0765
0766
0767
0768 if (variable == "meanX") {
0769 plotinfo.nbins = 50;
0770 plotinfo.min = -0.001;
0771 plotinfo.max = 0.001;
0772 } else if (variable == "meanY") {
0773 plotinfo.nbins = 50;
0774 plotinfo.min = -0.005;
0775 plotinfo.max = 0.005;
0776 } else if (variable == "medianX")
0777 if (plotSplits) {
0778 plotinfo.nbins = 50;
0779 plotinfo.min = -0.0005;
0780 plotinfo.max = 0.0005;
0781 } else {
0782 plotinfo.nbins = 100;
0783 plotinfo.min = -0.001;
0784 plotinfo.max = 0.001;
0785 }
0786 else if (variable == "medianY")
0787 if (plotSplits) {
0788 plotinfo.nbins = 50;
0789 plotinfo.min = -0.0005;
0790 plotinfo.max = 0.0005;
0791 } else {
0792 plotinfo.nbins = 100;
0793 plotinfo.min = -0.001;
0794 plotinfo.max = 0.001;
0795 }
0796 else if (variable == "meanNormX") {
0797 plotinfo.nbins = 100;
0798 plotinfo.min = -2.0;
0799 plotinfo.max = 2.0;
0800 } else if (variable == "meanNormY") {
0801 plotinfo.nbins = 100;
0802 plotinfo.min = -2.0;
0803 plotinfo.max = 2.0;
0804 } else if (variable == "rmsX") {
0805 plotinfo.nbins = 100;
0806 plotinfo.min = 0.0;
0807 plotinfo.max = 0.1;
0808 } else if (variable == "rmsY") {
0809 plotinfo.nbins = 100;
0810 plotinfo.min = 0.0;
0811 plotinfo.max = 0.1;
0812 } else if (variable == "rmsNormX") {
0813 plotinfo.nbins = 100;
0814 plotinfo.min = 0.3;
0815 plotinfo.max = 1.8;
0816 } else if (variable == "rmsNormY") {
0817 plotinfo.nbins = 100;
0818 plotinfo.min = 0.3;
0819 plotinfo.max = 1.8;
0820 } else {
0821 std::cerr << "Unknown variable " << variable << std::endl;
0822 plotinfo.nbins = 100;
0823 plotinfo.min = -0.1;
0824 plotinfo.max = 0.1;
0825 }
0826
0827 for (int i = 1; i <= 6; ++i) {
0828
0829 if (i != 1 && i != 2 && variable.length() > 0 && variable[variable.length() - 1] == 'Y') {
0830 continue;
0831 }
0832
0833
0834 if (plotLayerN > maxNumberOfLayers(i)) {
0835 continue;
0836 }
0837
0838 plotinfo.plotSplits = plotSplits && plotSplitsFor[i - 1];
0839 if (!plotinfo.plotPlain && !plotinfo.plotSplits) {
0840 continue;
0841 }
0842
0843
0844
0845 bool hasheader = (TkAlStyle::legendheader != "");
0846
0847 int nPlots = 1;
0848 if (plotinfo.plotSplits) {
0849 nPlots = 3;
0850 }
0851
0852
0853 if (plotinfo.plotLayers) {
0854 nPlots *= maxNumberOfLayers(i);
0855 }
0856 nPlots *= sourceList.size();
0857 if (twolines_) {
0858 nPlots *= 2;
0859 }
0860 nPlots += hasheader;
0861
0862 double legendY = 0.80;
0863 if (nPlots > 3) {
0864 legendY -= 0.01 * (nPlots - 3);
0865 }
0866 if (bigtext_) {
0867 legendY -= 0.05;
0868 }
0869 if (legendY < 0.6) {
0870 std::cerr << "Warning: Huge legend!" << std::endl;
0871 legendY = 0.6;
0872 }
0873
0874 THStack hstack("hstack", "hstack");
0875 plotinfo.maxY = 0;
0876 plotinfo.subDetId = i;
0877 plotinfo.legend = new TLegend(0.17, legendY, 0.85, 0.88);
0878 plotinfo.legend->SetNColumns(2);
0879 if (hasheader)
0880 plotinfo.legend->SetHeader(TkAlStyle::legendheader);
0881 if (bigtext_)
0882 plotinfo.legend->SetTextSize(TkAlStyle::textSize);
0883 plotinfo.legend->SetFillStyle(0);
0884 plotinfo.hstack = &hstack;
0885 plotinfo.h = plotinfo.h1 = plotinfo.h2 = nullptr;
0886 plotinfo.firsthisto = true;
0887
0888 openSummaryFile();
0889 vmean.clear();
0890 vrms.clear();
0891 vdeltamean.clear();
0892 vmeanerror.clear();
0893 vPValueEqualSplitMeans.clear(), vAlignmentUncertainty.clear();
0894 vPValueMeanEqualIdeal.clear();
0895 vPValueRMSEqualIdeal.clear();
0896
0897 std::string stringsubdet;
0898 switch (i) {
0899 case 1:
0900 stringsubdet = "BPIX";
0901 break;
0902 case 2:
0903 stringsubdet = "FPIX";
0904 break;
0905 case 3:
0906 stringsubdet = "TIB";
0907 break;
0908 case 4:
0909 stringsubdet = "TID";
0910 break;
0911 case 5:
0912 stringsubdet = "TOB";
0913 break;
0914 case 6:
0915 stringsubdet = "TEC";
0916 break;
0917 }
0918
0919 for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0920 plotinfo.vars = *it;
0921 plotinfo.h1 = plotinfo.h2 = plotinfo.h = nullptr;
0922
0923 int minlayer = plotLayers ? 1 : plotLayerN;
0924
0925 if (plotinfo.plotPlain)
0926 minlayer = 0;
0927 int maxlayer = plotLayers ? numberOfLayers(plotinfo.vars->getPhase(), plotinfo.subDetId) : plotLayerN;
0928
0929 for (int layer = minlayer; layer <= maxlayer; layer++) {
0930 if (plotinfo.plotPlain) {
0931 plotDMRHistogram(plotinfo, 0, layer, stringsubdet);
0932 }
0933
0934 if (plotinfo.plotSplits) {
0935 plotDMRHistogram(plotinfo, -1, layer, stringsubdet);
0936 plotDMRHistogram(plotinfo, 1, layer, stringsubdet);
0937 }
0938
0939 if (plotinfo.plotPlain) {
0940 if (plotinfo.h) {
0941 setDMRHistStyleAndLegend(plotinfo.h, plotinfo, 0, layer);
0942 } else {
0943 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") &&
0944 layer == 0) {
0945 vmean.push_back(nan(""));
0946 vrms.push_back(nan(""));
0947 vmeanerror.push_back(nan(""));
0948 vAlignmentUncertainty.push_back(nan(""));
0949 vPValueMeanEqualIdeal.push_back(nan(""));
0950 vPValueRMSEqualIdeal.push_back(nan(""));
0951 if (plotinfo.plotSplits && plotinfo.plotPlain) {
0952 vdeltamean.push_back(nan(""));
0953 vPValueEqualSplitMeans.push_back(nan(""));
0954 }
0955 }
0956 }
0957 }
0958
0959 if (plotinfo.plotSplits) {
0960
0961 if (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && !plotinfo.plotPlain) {
0962 std::ostringstream legend;
0963 std::string unit = " #mum";
0964 legend.precision(3);
0965 legend << std::fixed;
0966 float factor = 10000.0f;
0967 if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" ||
0968 plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY") {
0969 factor = 1.0f;
0970 unit = "";
0971 }
0972 float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
0973 legend << plotinfo.vars->getName();
0974 if (layer > 0) {
0975
0976 if (i == 4 || i == 6)
0977 legend << ", disc ";
0978 else
0979 legend << ", layer ";
0980 legend << layer;
0981 }
0982 plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
0983 legend.str("");
0984 legend << "#Delta#mu = " << deltamu << unit;
0985 plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
0986
0987 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && !plotLayers && layer == 0) {
0988 vdeltamean.push_back(deltamu);
0989 }
0990 }
0991 if (plotinfo.h1) {
0992 setDMRHistStyleAndLegend(plotinfo.h1, plotinfo, -1, layer);
0993 }
0994 if (plotinfo.h2) {
0995 setDMRHistStyleAndLegend(plotinfo.h2, plotinfo, 1, layer);
0996 }
0997 }
0998 }
0999 }
1000
1001 if (hstack.GetHists() != nullptr && hstack.GetHists()->GetSize() != 0) {
1002 hstack.Draw("nostack");
1003 hstack.SetMaximum(plotinfo.maxY * 1.3);
1004 setTitleStyle(hstack, variable.c_str(), "#modules", plotinfo.subDetId);
1005 setHistStyle(*hstack.GetHistogram(), variable.c_str(), "#modules", 1);
1006
1007 plotinfo.legend->Draw();
1008 } else {
1009
1010 plotinfo.h = new TH1F("defhist", "Empty default histogram", plotinfo.nbins, plotinfo.min, plotinfo.max);
1011 plotinfo.h->SetMaximum(10);
1012 if (plotinfo.variable.find("Norm") == std::string::npos)
1013 scaleXaxis(plotinfo.h, 10000);
1014 setTitleStyle(*plotinfo.h, variable.c_str(), "#modules", plotinfo.subDetId);
1015 setHistStyle(*plotinfo.h, variable.c_str(), "#modules", 1);
1016 plotinfo.h->Draw();
1017 }
1018
1019 std::ostringstream plotName;
1020 plotName << outputDir << "/D";
1021
1022 if (variable == "medianX")
1023 plotName << "medianR_";
1024 else if (variable == "medianY")
1025 plotName << "medianYR_";
1026 else if (variable == "meanX")
1027 plotName << "meanR_";
1028 else if (variable == "meanY")
1029 plotName << "meanYR_";
1030 else if (variable == "meanNormX")
1031 plotName << "meanNR_";
1032 else if (variable == "meanNormY")
1033 plotName << "meanNYR_";
1034 else if (variable == "rmsX")
1035 plotName << "rmsR_";
1036 else if (variable == "rmsY")
1037 plotName << "rmsYR_";
1038 else if (variable == "rmsNormX")
1039 plotName << "rmsNR_";
1040 else if (variable == "rmsNormY")
1041 plotName << "rmsNYR_";
1042
1043 TString subdet;
1044 switch (i) {
1045 case 1:
1046 subdet = "BPIX";
1047 break;
1048 case 2:
1049 subdet = "FPIX";
1050 break;
1051 case 3:
1052 subdet = "TIB";
1053 break;
1054 case 4:
1055 subdet = "TID";
1056 break;
1057 case 5:
1058 subdet = "TOB";
1059 break;
1060 case 6:
1061 subdet = "TEC";
1062 break;
1063 }
1064
1065 plotName << subdet;
1066
1067 if (plotPlain && !plotSplits) {
1068 plotName << "_plain";
1069 } else if (!plotPlain && plotSplits) {
1070 plotName << "_split";
1071 }
1072 if (plotLayers) {
1073
1074 if (i == 4 || i == 6)
1075 plotName << "_discs";
1076 else
1077 plotName << "_layers";
1078 }
1079 if (plotLayerN > 0) {
1080
1081 if (i == 4 || i == 6)
1082 plotName << "_disc";
1083 else
1084 plotName << "_layer";
1085 plotName << plotLayerN;
1086 }
1087
1088
1089 c.Update();
1090 c.Print((plotName.str() + ".png").c_str());
1091 c.Print((plotName.str() + ".eps").c_str());
1092 c.Print((plotName.str() + ".pdf").c_str());
1093
1094
1095 TFile f((plotName.str() + ".root").c_str(), "recreate");
1096 c.Write();
1097 f.Close();
1098
1099
1100 delete plotinfo.h;
1101 delete plotinfo.h1;
1102 delete plotinfo.h2;
1103
1104 if (!vmean.empty()) {
1105 summaryfile << " mu_" << subdet;
1106 if (plotinfo.variable == "medianY")
1107 summaryfile << "_y";
1108 summaryfile << " (um)\t"
1109 << "latexname=$\\mu_\\text{" << subdet << "}";
1110 if (plotinfo.variable == "medianY")
1111 summaryfile << "^{y}";
1112 summaryfile << "$ ($\\mu$m)\t"
1113 << "format={:.3g}\t"
1114 << "latexformat=${:.3g}$";
1115 for (auto mu : vmean)
1116 summaryfile << "\t" << mu;
1117 summaryfile << "\n";
1118 }
1119 if (!vrms.empty()) {
1120 summaryfile << "sigma_" << subdet;
1121 if (plotinfo.variable == "medianY")
1122 summaryfile << "_y";
1123 summaryfile << " (um)\t"
1124 << "latexname=$\\sigma_\\text{" << subdet << "}";
1125 if (plotinfo.variable == "medianY")
1126 summaryfile << "^{y}";
1127 summaryfile << "$ ($\\mu$m)\t"
1128 << "format={:.3g}\t"
1129 << "latexformat=${:.3g}$";
1130 for (auto sigma : vrms)
1131 summaryfile << "\t" << sigma;
1132 summaryfile << "\n";
1133 }
1134 if (!vdeltamean.empty()) {
1135 summaryfile << " dmu_" << subdet;
1136 if (plotinfo.variable == "medianY")
1137 summaryfile << "_y";
1138 summaryfile << " (um)\t"
1139 << "latexname=$\\Delta\\mu_\\text{" << subdet << "}";
1140 if (plotinfo.variable == "medianY")
1141 summaryfile << "^{y}";
1142 summaryfile << "$ ($\\mu$m)\t"
1143 << "format={:.3g}\t"
1144 << "latexformat=${:.3g}$";
1145 for (auto dmu : vdeltamean)
1146 summaryfile << "\t" << dmu;
1147 summaryfile << "\n";
1148 }
1149 if (!vmeanerror.empty()) {
1150 summaryfile << " sigma_mu_" << subdet;
1151 if (plotinfo.variable == "medianY")
1152 summaryfile << "_y";
1153 summaryfile << " (um)\t"
1154 << "latexname=$\\sigma\\mu_\\text{" << subdet << "}";
1155 if (plotinfo.variable == "medianY")
1156 summaryfile << "^{y}";
1157 summaryfile << "$ ($\\mu$m)\t"
1158 << "format={:.3g}\t"
1159 << "latexformat=${:.3g}$";
1160 for (auto dmu : vmeanerror)
1161 summaryfile << "\t" << dmu;
1162 summaryfile << "\n";
1163 }
1164 if (!vPValueEqualSplitMeans.empty()) {
1165 summaryfile << " p_delta_mu_equal_zero_" << subdet;
1166 if (plotinfo.variable == "medianY")
1167 summaryfile << "_y";
1168 summaryfile << "\t"
1169 << "latexname=$P(\\delta\\mu_\\text{" << subdet << "}=0)";
1170 if (plotinfo.variable == "medianY")
1171 summaryfile << "^{y}";
1172 summaryfile << "$\t"
1173 << "format={:.3g}\t"
1174 << "latexformat=${:.3g}$";
1175 for (auto dmu : vPValueEqualSplitMeans)
1176 summaryfile << "\t" << dmu;
1177 summaryfile << "\n";
1178 }
1179 if (!vAlignmentUncertainty.empty()) {
1180 summaryfile << " alignment_uncertainty_" << subdet;
1181 if (plotinfo.variable == "medianY")
1182 summaryfile << "_y";
1183 summaryfile << " (um)\t"
1184 << "latexname=$\\sigma_\\text{align}_\\text{" << subdet << "}";
1185 if (plotinfo.variable == "medianY")
1186 summaryfile << "^{y}";
1187 summaryfile << "$ ($\\mu$m)\t"
1188 << "format={:.3g}\t"
1189 << "latexformat=${:.3g}$";
1190 for (auto dmu : vAlignmentUncertainty)
1191 summaryfile << "\t" << dmu;
1192 summaryfile << "\n";
1193 }
1194 if (!vPValueMeanEqualIdeal.empty()) {
1195 summaryfile << " p_mean_equals_ideal_" << subdet;
1196 if (plotinfo.variable == "medianY")
1197 summaryfile << "_y";
1198 summaryfile << "\t"
1199 << "latexname=$P(\\mu_\\text{" << subdet << "}=\\mu_\\text{ideal})";
1200 if (plotinfo.variable == "medianY")
1201 summaryfile << "^{y}";
1202 summaryfile << "$\t"
1203 << "format={:.3g}\t"
1204 << "latexformat=${:.3g}$";
1205 for (auto dmu : vPValueMeanEqualIdeal)
1206 summaryfile << "\t" << dmu;
1207 summaryfile << "\n";
1208 }
1209 if (!vPValueRMSEqualIdeal.empty()) {
1210 summaryfile << " p_RMS_equals_ideal_" << subdet;
1211 if (plotinfo.variable == "medianY")
1212 summaryfile << "_y";
1213 summaryfile << "\t"
1214 << "latexname=$P(\\sigma_\\text{" << subdet << "}=\\sigma_\\text{ideal})";
1215 if (plotinfo.variable == "medianY")
1216 summaryfile << "^{y}";
1217 summaryfile << "$\t"
1218 << "format={:.3g}\t"
1219 << "latexformat=${:.3g}$";
1220 for (auto dmu : vPValueRMSEqualIdeal)
1221 summaryfile << "\t" << dmu;
1222 summaryfile << "\n";
1223 }
1224 }
1225 }
1226
1227
1228 void PlotAlignmentValidation::plotChi2(const char* inputFile) {
1229
1230
1231
1232 Bool_t errorflag = kFALSE;
1233 TFile* fi1 = TFile::Open(inputFile, "read");
1234 if (fi1 != nullptr) {
1235 if (fi1->GetDirectory("TrackerOfflineValidation/GlobalTrackVariables") == nullptr) {
1236 errorflag = kTRUE;
1237 }
1238 } else {
1239 errorflag = kTRUE;
1240 }
1241 if (errorflag) {
1242 std::cout << "PlotAlignmentValidation::plotChi2: Can't find GlobalTrackVariables given file,"
1243 << " no chi^2-plots produced" << std::endl;
1244 return;
1245 }
1246
1247 auto normchi = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_normchi2");
1248
1249 normchi->GetFrame()->ResetBit(kCanDelete);
1250
1251 auto chiprob = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_chi2Prob");
1252
1253 chiprob->GetFrame()->ResetBit(kCanDelete);
1254
1255 if (normchi == nullptr || chiprob == nullptr) {
1256 errorflag = kTRUE;
1257 }
1258 if (errorflag) {
1259 std::cout << "PlotAlignmentValidation::plotChi2: Can't find data from given file,"
1260 << " no chi^2-plots produced" << std::endl;
1261 return;
1262 }
1263
1264 TLegend* legend = nullptr;
1265 for (auto primitive : *normchi->GetListOfPrimitives()) {
1266 legend = dynamic_cast<TLegend*>(primitive);
1267 if (legend)
1268 break;
1269 }
1270 if (legend) {
1271 openSummaryFile();
1272 summaryfile << "ntracks";
1273 for (auto alignment : sourceList) {
1274 summaryfile << "\t";
1275 TString title = alignment->getName();
1276 int color = alignment->getLineColor();
1277 int style = alignment->getLineStyle();
1278 bool foundit = false;
1279 for (auto entry : *legend->GetListOfPrimitives()) {
1280 TLegendEntry* legendentry = dynamic_cast<TLegendEntry*>(entry);
1281 assert(legendentry);
1282 TH1* h = dynamic_cast<TH1*>(legendentry->GetObject());
1283 if (!h)
1284 continue;
1285 if (legendentry->GetLabel() == title && h->GetLineColor() == color && h->GetLineStyle() == style) {
1286 foundit = true;
1287 summaryfile << h->GetEntries();
1288 break;
1289 }
1290 }
1291 if (!foundit) {
1292 summaryfile << 0;
1293 }
1294 }
1295 summaryfile << "\n";
1296 }
1297
1298 chiprob->Draw();
1299 normchi->Draw();
1300
1301
1302 normchi->Print((outputDir + "/h_normchi2.png").c_str());
1303 chiprob->Print((outputDir + "/h_chi2Prob.png").c_str());
1304 normchi->Print((outputDir + "/h_normchi2.eps").c_str());
1305 chiprob->Print((outputDir + "/h_chi2Prob.eps").c_str());
1306 normchi->Print((outputDir + "/h_normchi2.pdf").c_str());
1307 chiprob->Print((outputDir + "/h_chi2Prob.pdf").c_str());
1308
1309
1310 TFile fi2((outputDir + "/h_normchi2.root").c_str(), "recreate");
1311 normchi->Write();
1312 fi2.Close();
1313
1314 TFile fi3((outputDir + "/h_chi2Prob.root").c_str(), "recreate");
1315 chiprob->Write();
1316 fi3.Close();
1317
1318 fi1->Close();
1319 delete fi1;
1320 }
1321
1322
1323 THStack* PlotAlignmentValidation::addHists(
1324 const TString& selection, const TString& residType, TLegend** myLegend, bool printModuleIds, bool validforphase0) {
1325 enum ResidType {
1326 xPrimeRes,
1327 yPrimeRes,
1328 xPrimeNormRes,
1329 yPrimeNormRes,
1330 xRes,
1331 yRes,
1332 xNormRes,
1333 ResXvsXProfile,
1334 ResXvsYProfile,
1335 ResYvsXProfile,
1336 ResYvsYProfile
1337 };
1338 ResidType rType = xPrimeRes;
1339 if (residType == "xPrime")
1340 rType = xPrimeRes;
1341 else if (residType == "yPrime")
1342 rType = yPrimeRes;
1343 else if (residType == "xPrimeNorm")
1344 rType = xPrimeNormRes;
1345 else if (residType == "yPrimeNorm")
1346 rType = yPrimeNormRes;
1347 else if (residType == "x")
1348 rType = xRes;
1349 else if (residType == "y")
1350 rType = yRes;
1351 else if (residType == "xNorm")
1352 rType = xNormRes;
1353
1354 else if (residType == "ResXvsXProfile")
1355 rType = ResXvsXProfile;
1356 else if (residType == "ResYvsXProfile")
1357 rType = ResYvsXProfile;
1358 else if (residType == "ResXvsYProfile")
1359 rType = ResXvsYProfile;
1360 else if (residType == "ResYvsYProfile")
1361 rType = ResYvsYProfile;
1362 else {
1363 std::cout << "PlotAlignmentValidation::addHists: Unknown residual type " << residType << std::endl;
1364 return nullptr;
1365 }
1366
1367 std::cout << "PlotAlignmentValidation::addHists: using selection " << selection << std::endl;
1368 THStack* retHistoStack = new THStack("hstack", "");
1369 if (myLegend != nullptr)
1370 if (*myLegend == nullptr) {
1371 *myLegend = new TLegend(0.17, 0.80, 0.85, 0.88);
1372 }
1373
1374 for (std::vector<TkOfflineVariables*>::iterator itSourceFile = sourceList.begin(); itSourceFile != sourceList.end();
1375 ++itSourceFile) {
1376 std::vector<TString> histnames;
1377
1378 TFile* f = (*itSourceFile)->getFile();
1379 TTree* tree = (*itSourceFile)->getTree();
1380 int myLineColor = (*itSourceFile)->getLineColor();
1381 int myLineStyle = (*itSourceFile)->getLineStyle();
1382 TString myLegendName = (*itSourceFile)->getName();
1383 TH1* h = nullptr;
1384 UInt_t nEmpty = 0;
1385 Long64_t nentries = tree->GetEntriesFast();
1386 if (!f || !tree) {
1387 std::cout << "PlotAlignmentValidation::addHists: no tree or no file" << std::endl;
1388 return nullptr;
1389 }
1390
1391 bool histnamesfilled = false;
1392 int phase = (bool)(f->Get("TrackerOfflineValidation/Pixel/P1PXBBarrel_1"));
1393 if (residType.Contains("Res") && residType.Contains("Profile")) {
1394 TString basename = TString(residType)
1395 .ReplaceAll("Res", "p_res")
1396 .ReplaceAll("vs", "")
1397 .ReplaceAll("Profile", "_");
1398 if (selection == "subDetId==1") {
1399 if (phase == 1)
1400 histnames.push_back(TString(basename) += "P1PXBBarrel_1");
1401 else
1402 histnames.push_back(TString(basename) += "TPBBarrel_1");
1403 histnamesfilled = true;
1404 } else if (selection == "subDetId==2") {
1405 if (phase == 1) {
1406 histnames.push_back(TString(basename) += "P1PXECEndcap_2");
1407 histnames.push_back(TString(basename) += "P1PXECEndcap_3");
1408 } else {
1409 histnames.push_back(TString(basename) += "TPEEndcap_2");
1410 histnames.push_back(TString(basename) += "TPEEndcap_3");
1411 }
1412 histnamesfilled = true;
1413 } else if (selection == "subDetId==3") {
1414 histnames.push_back(TString(basename) += "TIBBarrel_1");
1415 histnamesfilled = true;
1416 } else if (selection == "subDetId==4") {
1417 histnames.push_back(TString(basename) += "TIDEndcap_2");
1418 histnames.push_back(TString(basename) += "TIDEndcap_3");
1419 histnamesfilled = true;
1420 } else if (selection == "subDetId==5") {
1421 histnames.push_back(TString(basename) += "TOBBarrel_4");
1422 histnamesfilled = true;
1423 } else if (selection == "subDetId==6") {
1424 histnames.push_back(TString(basename) += "TECEndcap_5");
1425 histnames.push_back(TString(basename) += "TECEndcap_6");
1426 histnamesfilled = true;
1427 } else if (selection == "subDetId==6 && ring <= 4") {
1428
1429 for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1430 for (int iDisk = 1; iDisk <= 9; iDisk++)
1431 for (int iSide = 1; iSide <= 2; iSide++)
1432 for (int iPetal = 1; iPetal <= 8; iPetal++)
1433 for (int iRing = 1; iRing <= 4 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9); iRing++)
1434
1435
1436 {
1437 std::stringstream s;
1438 s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
1439 << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1440 histnames.push_back(TString(s.str()));
1441 }
1442 histnamesfilled = true;
1443 } else if (selection == "subDetId==6 && ring > 4") {
1444
1445 for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1446 for (int iDisk = 1; iDisk <= 9; iDisk++)
1447 for (int iSide = 1; iSide <= 2; iSide++)
1448 for (int iPetal = 1; iPetal <= 8; iPetal++)
1449 for (int iRing = 5 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1450 iRing <= 7 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1451 iRing++)
1452
1453
1454 {
1455 std::stringstream s;
1456 s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
1457 << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1458 histnames.push_back(TString(s.str()));
1459 }
1460 histnamesfilled = true;
1461 }
1462 }
1463
1464 Long64_t nSel = 0;
1465 if (histnamesfilled && !histnames.empty()) {
1466 nSel = (Long64_t)histnames.size();
1467 }
1468 if (!histnamesfilled) {
1469
1470
1471 nSel = tree->Draw("Entry$", selection, "goff");
1472 if (nSel == -1)
1473 return nullptr;
1474 if (nSel == 0) {
1475 std::cout << "PlotAlignmentValidation::addHists: no selected module." << std::endl;
1476 return nullptr;
1477 }
1478
1479 const std::vector<double> selected(tree->GetV1(), tree->GetV1() + nSel);
1480
1481 std::vector<double>::const_iterator iterEnt = selected.begin();
1482
1483
1484
1485 TkOffTreeVariables* treeMem = nullptr;
1486 tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
1487 for (Long64_t i = 0; i < nentries; i++) {
1488 if (i < *iterEnt - 0.1
1489 || iterEnt == selected.end()) {
1490 continue;
1491 } else if (TMath::Abs(i - *iterEnt) < 0.11) {
1492 ++iterEnt;
1493 } else
1494 std::cout << "Must not happen: " << i << " " << *iterEnt << std::endl;
1495
1496 tree->GetEntry(i);
1497 if (printModuleIds) {
1498 std::cout << treeMem->moduleId << ": " << treeMem->entries << " entries" << std::endl;
1499 }
1500 if (treeMem->entries <= 0) {
1501 ++nEmpty;
1502 continue;
1503 }
1504 TString hName;
1505 switch (rType) {
1506 case xPrimeRes:
1507 hName = treeMem->histNameX.c_str();
1508 break;
1509 case yPrimeRes:
1510 hName = treeMem->histNameY.c_str();
1511 break;
1512 case xPrimeNormRes:
1513 hName = treeMem->histNameNormX.c_str();
1514 break;
1515 case yPrimeNormRes:
1516 hName = treeMem->histNameNormY.c_str();
1517 break;
1518 case xRes:
1519 hName = treeMem->histNameLocalX.c_str();
1520 break;
1521 case yRes:
1522 hName = treeMem->histNameLocalY.c_str();
1523 break;
1524 case xNormRes:
1525 hName = treeMem->histNameNormLocalX.c_str();
1526 break;
1527
1528 case ResXvsXProfile:
1529 hName = treeMem->profileNameResXvsX.c_str();
1530 break;
1531 case ResXvsYProfile:
1532 hName = treeMem->profileNameResXvsY.c_str();
1533 break;
1534 case ResYvsXProfile:
1535 hName = treeMem->profileNameResYvsX.c_str();
1536 break;
1537 case ResYvsYProfile:
1538 hName = treeMem->profileNameResYvsY.c_str();
1539 break;
1540 }
1541 histnames.push_back(hName);
1542 }
1543 }
1544
1545 for (std::vector<TString>::iterator ithistname = histnames.begin(); ithistname != histnames.end(); ++ithistname) {
1546 if (phase == 0 && !validforphase0)
1547 break;
1548 TH1* newHist;
1549 if (ithistname->Contains("/")) {
1550 newHist = (TH1*)f->Get(*ithistname);
1551 } else {
1552 TKey* histKey = f->FindKeyAny(*ithistname);
1553 newHist = (histKey ? static_cast<TH1*>(histKey->ReadObj()) : nullptr);
1554 }
1555 if (!newHist) {
1556 std::cout << "Hist " << *ithistname << " not found in file, break loop." << std::endl;
1557 break;
1558 }
1559 if (newHist->GetEntries() == 0) {
1560 nEmpty++;
1561 continue;
1562 }
1563 newHist->SetLineColor(myLineColor);
1564 newHist->SetLineStyle(myLineStyle);
1565 if (!h) {
1566 TString name(newHist->GetName());
1567 Ssiz_t pos_ = 0;
1568 for (UInt_t i2 = 0; i2 < 3; ++i2)
1569 pos_ = name.Index("_", pos_ + 1);
1570 name = name(0, pos_);
1571 h = static_cast<TH1*>(newHist->Clone("summed_" + name));
1572
1573
1574 } else {
1575 h->Add(newHist);
1576 }
1577 delete newHist;
1578 }
1579
1580 std::cout << "PlotAlignmentValidation::addHists"
1581 << "Result is merged from " << nSel - nEmpty << " hists, " << nEmpty << " hists were empty." << std::endl;
1582
1583 if (nSel - nEmpty == 0)
1584 continue;
1585
1586 if (myLegend != nullptr)
1587 (*myLegend)->AddEntry(h, myLegendName, "L");
1588
1589 retHistoStack->Add(h);
1590 }
1591
1592 return retHistoStack;
1593 }
1594
1595
1596
1597
1598
1599 TF1* PlotAlignmentValidation::fitGauss(TH1* hist, int color) {
1600
1601
1602
1603 if (!hist || hist->GetEntries() < 20)
1604 return nullptr;
1605
1606 float mean = hist->GetMean();
1607 float sigma = hist->GetRMS();
1608 std::string functionname = "gaussian_";
1609 functionname += hist->GetName();
1610 TF1* func = new TF1(functionname.c_str(), "gaus", mean - 2. * sigma, mean + 2. * sigma);
1611
1612 func->SetLineColor(color);
1613 func->SetLineStyle(2);
1614 if (0 == hist->Fit(func, "QNR")) {
1615 mean = func->GetParameter(1);
1616 sigma = func->GetParameter(2);
1617
1618 func->SetRange(mean - 3. * sigma, mean + 3. * sigma);
1619
1620
1621 if (0 == hist->Fit(func, "Q0ILR")) {
1622 if (hist->GetFunction(func->GetName())) {
1623
1624 }
1625 }
1626 }
1627 return func;
1628 }
1629
1630
1631
1632
1633
1634 void PlotAlignmentValidation::storeHistogramInRootfile(TH1* hist) {
1635
1636 rootsummaryfile->cd();
1637 hist->Write();
1638 }
1639
1640
1641 void PlotAlignmentValidation::scaleXaxis(TH1* hist, Int_t scale) {
1642 Double_t xmin = hist->GetXaxis()->GetXmin();
1643 Double_t xmax = hist->GetXaxis()->GetXmax();
1644 hist->GetXaxis()->SetLimits(xmin * scale, xmax * scale);
1645 }
1646
1647
1648 TObject* PlotAlignmentValidation::findObjectFromCanvas(TCanvas* canv, const char* className, Int_t n) {
1649
1650 TIter next(canv->GetListOfPrimitives());
1651 TObject* obj = nullptr;
1652 Int_t found = 0;
1653 while ((obj = next())) {
1654 if (strncmp(obj->ClassName(), className, 10) == 0) {
1655 if (++found == n)
1656 return obj;
1657 }
1658 }
1659
1660 return nullptr;
1661 }
1662
1663
1664 void PlotAlignmentValidation::setTitleStyle(
1665 TNamed& hist, const char* titleX, const char* titleY, int subDetId, bool isSurfaceDeformation, TString secondline) {
1666 std::stringstream title_Xaxis;
1667 std::stringstream title_Yaxis;
1668 TString titleXAxis = titleX;
1669 TString titleYAxis = titleY;
1670 if (titleXAxis != "" && titleYAxis != "")
1671 std::cout << "plot " << titleXAxis << " vs " << titleYAxis << std::endl;
1672
1673 hist.SetTitle("");
1674 TkAlStyle::drawStandardTitle();
1675
1676
1677 TString subD;
1678 switch (subDetId) {
1679 case 1:
1680 subD = "BPIX";
1681 break;
1682 case 2:
1683 subD = "FPIX";
1684 break;
1685 case 3:
1686 subD = "TIB";
1687 break;
1688 case 4:
1689 subD = "TID";
1690 break;
1691 case 5:
1692 subD = "TOB";
1693 break;
1694 case 6:
1695 subD = "TEC";
1696 break;
1697 }
1698
1699 TPaveText* text2;
1700 if (!isSurfaceDeformation) {
1701 text2 = new TPaveText(0.7, 0.3, 0.9, 0.6, "brNDC");
1702 } else {
1703 std::cout << "Surface Deformation" << std::endl;
1704 text2 = new TPaveText(0.8, 0.75, 0.9, 0.9, "brNDC");
1705 }
1706 text2->SetTextSize(0.06);
1707 text2->SetTextFont(42);
1708 text2->SetFillStyle(0);
1709 text2->SetBorderSize(0);
1710 text2->SetMargin(0.01);
1711 text2->SetTextAlign(12);
1712 text2->AddText(0.01, 0.75, subD);
1713 if (secondline != "") {
1714 text2->AddText(0.01, 0.25, secondline);
1715 }
1716 text2->Draw();
1717 }
1718
1719
1720
1721
1722
1723 void PlotAlignmentValidation::setHistStyle(TH1& hist, const char* titleX, const char* titleY, int color) {
1724 std::stringstream title_Xaxis;
1725 std::stringstream title_Yaxis;
1726 TString titleXAxis = titleX;
1727 TString titleYAxis = titleY;
1728
1729 if (titleXAxis.Contains("Phi"))
1730 title_Xaxis << titleX << "[rad]";
1731 else if (titleXAxis.Contains("meanX"))
1732 title_Xaxis << "#LTx'_{pred}-x'_{hit}#GT[#mum]";
1733 else if (titleXAxis.Contains("meanY"))
1734 title_Xaxis << "#LTy'_{pred}-y'_{hit}#GT[#mum]";
1735 else if (titleXAxis.Contains("rmsX"))
1736 title_Xaxis << "RMS(x'_{pred}-x'_{hit})[#mum]";
1737 else if (titleXAxis.Contains("rmsY"))
1738 title_Xaxis << "RMS(y'_{pred}-y'_{hit})[#mum]";
1739 else if (titleXAxis.Contains("meanNormX"))
1740 title_Xaxis << "#LTx'_{pred}-x'_{hit}/#sigma#GT";
1741 else if (titleXAxis.Contains("meanNormY"))
1742 title_Xaxis << "#LTy'_{pred}-y'_{hit}/#sigma#GT";
1743 else if (titleXAxis.Contains("rmsNormX"))
1744 title_Xaxis << "RMS(x'_{pred}-x'_{hit}/#sigma)";
1745 else if (titleXAxis.Contains("rmsNormY"))
1746 title_Xaxis << "RMS(y'_{pred}-y'_{hit}/#sigma)";
1747 else if (titleXAxis.Contains("meanLocalX"))
1748 title_Xaxis << "#LTx_{pred}-x_{hit}#GT[#mum]";
1749 else if (titleXAxis.Contains("rmsLocalX"))
1750 title_Xaxis << "RMS(x_{pred}-x_{hit})[#mum]";
1751 else if (titleXAxis.Contains("meanNormLocalX"))
1752 title_Xaxis << "#LTx_{pred}-x_{hit}/#sigma#GT[#mum]";
1753 else if (titleXAxis.Contains("rmsNormLocalX"))
1754 title_Xaxis << "RMS(x_{pred}-x_{hit}/#sigma)[#mum]";
1755 else if (titleXAxis.Contains("medianX"))
1756 title_Xaxis << "median(x'_{pred}-x'_{hit})[#mum]";
1757 else if (titleXAxis.Contains("medianY"))
1758 title_Xaxis << "median(y'_{pred}-y'_{hit})[#mum]";
1759 else
1760 title_Xaxis << titleX << "[cm]";
1761
1762 if (hist.IsA()->InheritsFrom(TH1F::Class()))
1763 hist.SetLineColor(color);
1764 if (hist.IsA()->InheritsFrom(TProfile::Class())) {
1765 hist.SetMarkerStyle(20);
1766 hist.SetMarkerSize(0.8);
1767 hist.SetMarkerColor(color);
1768 }
1769
1770 hist.GetXaxis()->SetTitle((title_Xaxis.str()).c_str());
1771
1772 double binning = (hist.GetXaxis()->GetXmax() - hist.GetXaxis()->GetXmin()) / hist.GetNbinsX();
1773 title_Yaxis.precision(2);
1774
1775 if (((titleYAxis.Contains("layer") || titleYAxis.Contains("ring")) && titleYAxis.Contains("subDetId")) ||
1776 titleYAxis.Contains("#modules")) {
1777 title_Yaxis << "number of modules";
1778 if (TString(title_Xaxis.str()).Contains("[#mum]"))
1779 title_Yaxis << " / " << binning << " #mum";
1780 else if (TString(title_Xaxis.str()).Contains("[cm]"))
1781 title_Yaxis << " / " << binning << " cm";
1782 else
1783 title_Yaxis << " / " << binning;
1784 } else
1785 title_Yaxis << titleY << "[cm]";
1786
1787 hist.GetYaxis()->SetTitle((title_Yaxis.str()).c_str());
1788
1789 hist.GetXaxis()->SetTitleFont(42);
1790 hist.GetYaxis()->SetTitleFont(42);
1791 }
1792
1793
1794 std::string PlotAlignmentValidation::getSelectionForDMRPlot(int minHits, int subDetId, int direction, int layer) {
1795 std::ostringstream builder;
1796 builder << "entries >= " << minHits;
1797 builder << " && subDetId == " << subDetId;
1798 if (direction != 0) {
1799 if (subDetId == 2) {
1800 builder << " && zDirection == " << direction;
1801 } else {
1802 builder << " && rDirection == " << direction;
1803 }
1804 }
1805 if (layer > 0) {
1806 builder << " && layer == " << layer;
1807 }
1808 return builder.str();
1809 }
1810
1811 std::string PlotAlignmentValidation::getVariableForDMRPlot(
1812 const std::string& histoname, const std::string& variable, int nbins, double min, double max) {
1813 std::ostringstream builder;
1814 builder << variable << ">>" << histoname << "(" << nbins << "," << min << "," << max << ")";
1815 return builder.str();
1816 }
1817
1818
1819 void PlotAlignmentValidation::setDMRHistStyleAndLegend(TH1F* h,
1820 PlotAlignmentValidation::DMRPlotInfo& plotinfo,
1821 int direction,
1822 int layer) {
1823 TF1* fitResults = nullptr;
1824
1825 h->SetDirectory(nullptr);
1826
1827
1828
1829 h->SetLineWidth((direction == 0 || (plotinfo.plotSplits && !plotinfo.plotPlain)) ? 2 : 1);
1830
1831
1832
1833
1834
1835 int linestyle = plotinfo.vars->getLineStyle() - 1, linestyleplus = 0;
1836 if (direction == 1) {
1837 linestyleplus = 1;
1838 }
1839 if (direction == -1) {
1840 linestyleplus = 2;
1841 }
1842 if (direction != 0 && plotinfo.plotSplits && !plotinfo.plotPlain) {
1843 linestyleplus--;
1844 }
1845 linestyle = (linestyle + linestyleplus) % 4 + 1;
1846
1847 int linecolor = plotinfo.vars->getLineColor();
1848 if (plotinfo.plotLayers && layer > 0) {
1849 linecolor += layer - 1;
1850 }
1851
1852 if (plotinfo.firsthisto) {
1853 setHistStyle(*h, plotinfo.variable.c_str(), "#modules", 1);
1854 plotinfo.firsthisto = false;
1855 }
1856
1857 h->SetLineColor(linecolor);
1858 h->SetLineStyle(linestyle);
1859
1860 if (plotinfo.maxY < h->GetMaximum()) {
1861 plotinfo.maxY = h->GetMaximum();
1862 }
1863
1864
1865 if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1866 plotinfo.variable == "meanY") {
1867 fitResults = fitGauss(h, linecolor);
1868 }
1869
1870 plotinfo.hstack->Add(h);
1871
1872 std::ostringstream legend;
1873 legend.precision(3);
1874 legend << std::fixed;
1875
1876
1877 if (direction == -1 && plotinfo.subDetId != 2) {
1878 legend << "rDirection < 0";
1879 } else if (direction == 1 && plotinfo.subDetId != 2) {
1880 legend << "rDirection > 0";
1881 } else if (direction == -1 && plotinfo.subDetId == 2) {
1882 legend << "zDirection < 0";
1883 } else if (direction == 1 && plotinfo.subDetId == 2) {
1884 legend << "zDirection > 0";
1885 } else {
1886 legend << plotinfo.vars->getName();
1887 if (layer > 0) {
1888
1889 if (plotinfo.subDetId == 4 || plotinfo.subDetId == 6)
1890 legend << ", disc ";
1891 else
1892 legend << ", layer ";
1893 legend << layer << "";
1894 }
1895 }
1896
1897 plotinfo.legend->AddEntry(h, legend.str().c_str(), "l");
1898 legend.str("");
1899
1900
1901 double mean = 0.0, meanerror = 0.0, rms = 0.0, rmserror = 0.0;
1902 TString rmsname, units;
1903 bool showdeltamu =
1904 (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && plotinfo.plotSplits && plotinfo.plotPlain && direction == 0);
1905 if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1906 plotinfo.variable == "meanY" || plotinfo.variable == "rmsX" || plotinfo.variable == "rmsY") {
1907 if (useFit_ && fitResults) {
1908 mean = fitResults->GetParameter(1) * 10000;
1909 meanerror = fitResults->GetParError(1) * 10000;
1910 rms = fitResults->GetParameter(2) * 10000;
1911 rmserror = fitResults->GetParError(2) * 10000;
1912 rmsname = "#sigma";
1913 delete fitResults;
1914 } else {
1915 mean = h->GetMean(1) * 10000;
1916 meanerror = h->GetMeanError(1) * 10000;
1917 rms = h->GetRMS(1) * 10000;
1918 rmserror = h->GetRMSError(1) * 10000;
1919 rmsname = "rms";
1920 }
1921 units = " #mum";
1922 } else if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
1923 plotinfo.variable == "rmsNormY") {
1924 mean = h->GetMean(1);
1925 meanerror = h->GetMeanError(1);
1926 rms = h->GetRMS(1);
1927 rmserror = h->GetRMSError(1);
1928 rmsname = "rms";
1929 units = "";
1930 }
1931 if (showMean_) {
1932 legend << " #mu = " << mean;
1933 if (showMeanError_)
1934 legend << " #pm " << meanerror;
1935 legend << units;
1936 if (showRMS_ || showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1937 legend << ", ";
1938 }
1939 if (showRMS_) {
1940 legend << " " << rmsname << " = " << rms;
1941 if (showRMSError_)
1942 legend << " #pm " << rmserror;
1943 legend << units;
1944 if (showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1945 legend << ", ";
1946 }
1947
1948 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && layer == 0 &&
1949 direction == 0) {
1950 vmean.push_back(mean);
1951 vrms.push_back(rms);
1952 vmeanerror.push_back(meanerror);
1953 TH1F* ideal = (TH1F*)plotinfo.hstack->GetHists()->At(0);
1954 TH1F* h = plotinfo.h;
1955 if (h->GetRMS() >= ideal->GetRMS()) {
1956 vAlignmentUncertainty.push_back(sqrt(pow(h->GetRMS(), 2) - pow(ideal->GetRMS(), 2)));
1957 } else {
1958 vAlignmentUncertainty.push_back(nan(""));
1959 }
1960 float p = (float)resampleTestOfEqualMeans(ideal, h, 10000);
1961 vPValueMeanEqualIdeal.push_back(p);
1962 p = resampleTestOfEqualRMS(ideal, h, 10000);
1963 vPValueRMSEqualIdeal.push_back(p);
1964 }
1965
1966
1967 if (showdeltamu) {
1968 float factor = 10000.0f;
1969 if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
1970 plotinfo.variable == "rmsNormY") {
1971 factor = 1.0f;
1972 }
1973 float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
1974 legend << "#Delta#mu = " << deltamu << units;
1975 if ((showModules_ || showUnderOverFlow_) && !twolines_)
1976 legend << ", ";
1977
1978 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && layer == 0 &&
1979 direction == 0) {
1980 vdeltamean.push_back(deltamu);
1981 if (plotinfo.h1->GetEntries() && plotinfo.h2->GetEntries()) {
1982 float p = (float)resampleTestOfEqualMeans(plotinfo.h1, plotinfo.h2, 10000);
1983 vPValueEqualSplitMeans.push_back(p);
1984 }
1985 }
1986 }
1987
1988 if (twolines_) {
1989 plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
1990 plotinfo.legend->AddEntry((TObject*)nullptr, "", "");
1991 legend.str("");
1992 }
1993
1994 if (!showUnderOverFlow_ && showModules_) {
1995 legend << (int)h->GetEntries() << " modules";
1996 }
1997 if (showUnderOverFlow_) {
1998 if (showModules_) {
1999 legend << (int)h->GetEntries() << " modules ("
2000 << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " outside range)";
2001 } else {
2002 legend << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " modules outside range";
2003 }
2004 }
2005 plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
2006
2007
2008 if (plotinfo.variable.find("Norm") == std::string::npos)
2009 scaleXaxis(h, 10000);
2010 }
2011
2012
2013
2014
2015
2016
2017 void PlotAlignmentValidation::plotDMRHistogram(PlotAlignmentValidation::DMRPlotInfo& plotinfo,
2018 int direction,
2019 int layer,
2020 std::string subdet) {
2021 TH1F* h = nullptr;
2022
2023
2024 TString histoname = "";
2025 if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY")
2026 histoname = "median";
2027 else if (plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY")
2028 histoname = "DrmsNR";
2029 histoname += "_";
2030 histoname += plotinfo.vars->getName();
2031 histoname.ReplaceAll(" ", "_");
2032 histoname += "_";
2033 histoname += subdet.c_str();
2034 if (plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormY")
2035 histoname += "_y";
2036 if (layer != 0) {
2037 if (subdet == "TID" || subdet == "TEC")
2038 histoname += "_disc";
2039 else
2040 histoname += "_layer";
2041 histoname += std::to_string(layer);
2042 }
2043 if (direction == -1) {
2044 histoname += "_minus";
2045 } else if (direction == 1) {
2046 histoname += "_plus";
2047 } else {
2048 histoname += "";
2049 }
2050
2051
2052 std::string plotVariable =
2053 getVariableForDMRPlot(histoname.Data(), plotinfo.variable, plotinfo.nbins, plotinfo.min, plotinfo.max);
2054 std::string selection = "";
2055 if (plotinfo.filterName.empty()) {
2056
2057 selection = getSelectionForDMRPlot(plotinfo.minHits, plotinfo.subDetId, direction, layer);
2058 plotinfo.vars->getTree()->Draw(plotVariable.c_str(), selection.c_str(), "goff");
2059 if (gDirectory)
2060 gDirectory->GetObject(histoname.Data(), h);
2061 if (h && h->GetEntries() > 0) {
2062 if (direction == -1) {
2063 plotinfo.h1 = h;
2064 } else if (direction == 1) {
2065 plotinfo.h2 = h;
2066 } else {
2067 plotinfo.h = h;
2068 }
2069 }
2070 } else {
2071 TTreeReader reader(plotinfo.vars->getTree());
2072 TTreeReaderValue<Float_t> varToPlot(reader, plotinfo.variable.c_str());
2073 TTreeReaderValue<unsigned int> _entries(reader, "entries");
2074 TTreeReaderValue<unsigned int> _subDetId(reader, "subDetId");
2075 TTreeReaderValue<unsigned int> _moduleId(reader, "moduleId");
2076 TTreeReaderValue<Float_t> _zDirection(reader, "zDirection");
2077 TTreeReaderValue<Float_t> _rDirection(reader, "rDirection");
2078 TTreeReaderValue<unsigned int> _layer(reader, "layer");
2079 std::string badModulesFile_ = plotinfo.filterName;
2080 TFile* fBadModules = new TFile(badModulesFile_.c_str(), "READ");
2081 TTree* tBadModules = (TTree*)fBadModules->Get("alignTree");
2082 TTreeReader readerBad(tBadModules);
2083 TTreeReaderValue<int> _valid(readerBad, "valid");
2084 TTreeReaderValue<int> _bad_id(readerBad, "id");
2085 TTreeReaderValue<double> _bad_lumi(readerBad, "lumi");
2086
2087
2088 std::ofstream fUsedModules;
2089 fUsedModules.open("usedModules.txt", std::ios::out | std::ios::app);
2090
2091
2092 for (uint i = 0; i < plotinfo.vars->getTree()->GetEntries(); i++) {
2093 reader.SetEntry(i);
2094 if (*_entries < uint(plotinfo.minHits))
2095 continue;
2096 if (*_subDetId != uint(plotinfo.subDetId))
2097 continue;
2098 if (direction != 0) {
2099 if (plotinfo.subDetId == 2) {
2100 if (*_zDirection != direction)
2101 continue;
2102 } else {
2103 if (*_rDirection != direction)
2104 continue;
2105 }
2106 }
2107 if (layer > 0) {
2108 if (*_layer != uint(layer))
2109 continue;
2110 }
2111 bool isBadModule = false;
2112 for (int ibad = 0; ibad < tBadModules->GetEntries(); ibad++) {
2113 readerBad.SetEntry(ibad);
2114
2115 if (subdet == "BPIX" || subdet == "FPIX") {
2116 if (*_bad_lumi <= 2.0)
2117 continue;
2118 } else {
2119 if (*_bad_lumi <= 7.0)
2120 continue;
2121 }
2122
2123 if (*_moduleId == uint(*_bad_id))
2124 isBadModule = true;
2125 }
2126
2127 if (isBadModule)
2128 continue;
2129 fUsedModules << *_moduleId << "\n";
2130 if (h) {
2131 h->Fill(*varToPlot);
2132 }
2133 }
2134
2135
2136 fUsedModules.close();
2137 fBadModules->Close();
2138 if (h) {
2139 h->SetName(histoname.Data());
2140 }
2141 }
2142
2143
2144 if (h && h->GetEntries() > 0) {
2145 if (direction == -1) {
2146 plotinfo.h1 = h;
2147 } else if (direction == 1) {
2148 plotinfo.h2 = h;
2149 } else {
2150 plotinfo.h = h;
2151 }
2152 }
2153 if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormX" ||
2154 plotinfo.variable == "rmsNormY") {
2155 storeHistogramInRootfile(h);
2156 }
2157 }
2158
2159 void PlotAlignmentValidation::modifySSHistAndLegend(THStack* hs, TLegend* legend) {
2160
2161
2162 Double_t legendY = 0.80;
2163 bool hasheader = (TkAlStyle::legendheader != "");
2164 if (hasheader)
2165 legend->SetHeader(TkAlStyle::legendheader);
2166 legend->SetFillStyle(0);
2167 int legendsize = hs->GetHists()->GetSize() + hasheader;
2168
2169 if (legendsize > 3)
2170 legendY -= 0.01 * (legendsize - 3);
2171 if (bigtext_) {
2172 legendY -= 0.05;
2173 }
2174 if (legendY < 0.6) {
2175 std::cerr << "Warning: Huge legend!" << std::endl;
2176 legendY = 0.6;
2177 }
2178 legend->SetY1(legendY);
2179 if (bigtext_)
2180 legend->SetTextSize(TkAlStyle::textSize);
2181
2182
2183 TProfile* prof = nullptr;
2184 TIter next(hs->GetHists());
2185 Int_t index = hasheader;
2186 while ((prof = (TProfile*)next())) {
2187
2188 Double_t scale = 10000;
2189 prof->Scale(scale);
2190
2191 Double_t stats[6] = {0};
2192 prof->GetStats(stats);
2193
2194 std::ostringstream legendtext;
2195 legendtext.precision(3);
2196 legendtext << std::fixed;
2197 legendtext << ": y mean = " << stats[4] / stats[0] * scale << " #mum";
2198
2199 TLegendEntry* entry = (TLegendEntry*)legend->GetListOfPrimitives()->At(index);
2200 if (entry == nullptr)
2201 std::cout << "PlotAlignmentValidation::PlotAlignmentValidation::modifySSLegend: Bad legend!" << std::endl;
2202 else
2203 entry->SetLabel((entry->GetLabel() + legendtext.str()).c_str());
2204 index++;
2205 }
2206
2207
2208 hs->SetMaximum(hs->GetMaximum("nostack PE") * 1.3);
2209 }
2210
2211
2212
2213
2214
2215
2216
2217 double PlotAlignmentValidation::resampleTestOfEqualRMS(TH1F* h1, TH1F* h2, int numSamples) {
2218
2219 std::vector<double> diff;
2220 diff.clear();
2221
2222 double rmsdiff = abs(h1->GetRMS() - h2->GetRMS());
2223
2224 double m1 = h1->GetMean();
2225 double m2 = h2->GetMean();
2226
2227 double d1 = 0;
2228 double d2 = 0;
2229
2230 double test_mean = 0;
2231 for (int i = 0; i < numSamples; i++) {
2232 d1 = 0;
2233 d2 = 0;
2234 for (int i = 0; i < h1->GetEntries(); i++) {
2235 d1 += h1->GetRandom() - m1;
2236 }
2237 for (int i = 0; i < h2->GetEntries(); i++) {
2238 d2 += h2->GetRandom() + m2;
2239 }
2240 d1 /= h1->GetEntries();
2241 d2 /= h2->GetEntries();
2242 diff.push_back(abs(d1 - d2 - rmsdiff));
2243 test_mean += abs(d1 - d2 - rmsdiff);
2244 }
2245 test_mean /= numSamples;
2246 edm::LogPrint("") << "test mean:" << test_mean;
2247
2248 double p = 0;
2249 for (double d : diff) {
2250 if (d > rmsdiff) {
2251 p += 1;
2252 }
2253 }
2254
2255 p /= numSamples;
2256 return p;
2257 }
2258
2259
2260
2261
2262
2263
2264
2265 double PlotAlignmentValidation::resampleTestOfEqualMeans(TH1F* h1, TH1F* h2, int numSamples) {
2266
2267 std::vector<double> diff;
2268 diff.clear();
2269
2270 double meandiff = abs(h1->GetMean() - h2->GetMean());
2271
2272 double d1 = 0;
2273 double d2 = 0;
2274
2275 double test_mean = 0;
2276 for (int i = 0; i < numSamples; i++) {
2277 d1 = 0;
2278 d2 = 0;
2279 for (int i = 0; i < h1->GetEntries(); i++) {
2280 d1 += h1->GetRandom();
2281 }
2282 for (int i = 0; i < h2->GetEntries(); i++) {
2283 d2 += h2->GetRandom();
2284 }
2285 d1 /= h1->GetEntries();
2286 d2 /= h2->GetEntries();
2287 diff.push_back(abs(d1 - d2 - meandiff));
2288 test_mean += abs(d1 - d2 - meandiff);
2289 }
2290 test_mean /= numSamples;
2291 edm::LogPrint("") << "test mean:" << test_mean;
2292
2293 double p = 0;
2294 for (double d : diff) {
2295 if (d > meandiff) {
2296 p += 1;
2297 }
2298 }
2299
2300 p /= numSamples;
2301 return p;
2302 }
2303
2304 float PlotAlignmentValidation::twotailedStudentTTestEqualMean(float t, float v) {
2305 return 2 * (1 - ROOT::Math::tdistribution_cdf(abs(t), v));
2306 }
2307
2308 const TString PlotAlignmentValidation::summaryfilename = "OfflineValidationSummary";
2309
2310 std::vector<TH1*> PlotAlignmentValidation::findmodule(TFile* f, unsigned int moduleid) {
2311
2312 TString histnamex;
2313 TString histnamey;
2314
2315 auto t = (TTree*)f->Get("TrackerOfflineValidation/TkOffVal");
2316
2317 TkOffTreeVariables* variables = nullptr;
2318 t->SetBranchAddress("TkOffTreeVariables", &variables);
2319 unsigned int number_of_entries = t->GetEntries();
2320 for (unsigned int i = 0; i < number_of_entries; i++) {
2321 t->GetEntry(i);
2322 if (variables->moduleId == moduleid) {
2323 histnamex = variables->histNameX;
2324 histnamey = variables->histNameY;
2325 break;
2326 }
2327 }
2328
2329 std::vector<TH1*> h;
2330
2331 auto h1 = (TH1*)f->FindObjectAny(histnamex);
2332 auto h2 = (TH1*)f->FindObjectAny(histnamey);
2333
2334 h1->SetDirectory(nullptr);
2335 h2->SetDirectory(nullptr);
2336
2337 h.push_back(h1);
2338 h.push_back(h2);
2339
2340 return h;
2341 }
2342
2343 void PlotAlignmentValidation::residual_by_moduleID(unsigned int moduleid) {
2344 TCanvas* cx = new TCanvas("x_residual");
2345 TCanvas* cy = new TCanvas("y_residual");
2346 TLegend* legendx = new TLegend(0.55, 0.7, 1, 0.9);
2347 TLegend* legendy = new TLegend(0.55, 0.7, 1, 0.9);
2348
2349 legendx->SetTextSize(0.016);
2350 legendx->SetTextAlign(12);
2351 legendy->SetTextSize(0.016);
2352 legendy->SetTextAlign(12);
2353
2354 for (auto it : sourceList) {
2355 TFile* file = it->getFile();
2356 int color = it->getLineColor();
2357 int linestyle = it->getLineStyle();
2358 TString legendname = it->getName();
2359 std::vector<TH1*> hist = findmodule(file, moduleid);
2360
2361 TString histnamex = legendname + " NEntries: " + std::to_string(int(hist[0]->GetEntries()));
2362 hist[0]->SetTitle(histnamex);
2363 hist[0]->SetStats(false);
2364 hist[0]->Rebin(50);
2365 hist[0]->SetBit(TH1::kNoTitle);
2366 hist[0]->SetLineColor(color);
2367 hist[0]->SetLineStyle(linestyle);
2368 cx->cd();
2369 hist[0]->Draw("Same");
2370 legendx->AddEntry(hist[0], histnamex, "l");
2371
2372 TString histnamey = legendname + " NEntries: " + std::to_string(int(hist[1]->GetEntries()));
2373 hist[1]->SetTitle(histnamey);
2374 hist[1]->SetStats(false);
2375 hist[1]->Rebin(50);
2376 hist[1]->SetBit(TH1::kNoTitle);
2377 hist[1]->SetLineColor(color);
2378 hist[1]->SetLineStyle(linestyle);
2379 cy->cd();
2380 hist[1]->Draw("Same");
2381 legendy->AddEntry(hist[1], histnamey, "l");
2382 }
2383
2384 TString filenamex = "x_residual_" + std::to_string(moduleid);
2385 TString filenamey = "y_residual_" + std::to_string(moduleid);
2386 cx->cd();
2387 legendx->Draw();
2388 cx->SaveAs(outputDir + "/" + filenamex + ".root");
2389 cx->SaveAs(outputDir + "/" + filenamex + ".pdf");
2390 cx->SaveAs(outputDir + "/" + filenamex + ".png");
2391 cx->SaveAs(outputDir + "/" + filenamex + ".eps");
2392
2393 cy->cd();
2394 legendy->Draw();
2395 cy->SaveAs(outputDir + "/" + filenamey + ".root");
2396 cy->SaveAs(outputDir + "/" + filenamey + ".pdf");
2397 cy->SaveAs(outputDir + "/" + filenamey + ".png");
2398 cy->SaveAs(outputDir + "/" + filenamey + ".eps");
2399 }