File indexing completed on 2025-01-04 00:29:07
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 Float_t maxBadLumiPixel,
0692 Float_t maxBadLumiStrip) {
0693
0694
0695
0696
0697 std::size_t findres = variable.find(',');
0698 if (findres != std::string::npos) {
0699 std::string substring1 = variable.substr(0, findres);
0700 std::string substring2 = variable.substr(findres + 1, std::string::npos);
0701 plotDMR(substring1, minHits, options, filterName, maxBadLumiPixel, maxBadLumiStrip);
0702 plotDMR(substring2, minHits, options, filterName, maxBadLumiPixel, maxBadLumiStrip);
0703 return;
0704 }
0705
0706
0707
0708 if (variable == "mean" || variable == "median" || variable == "meanNorm" || variable == "rms" ||
0709 variable == "rmsNorm") {
0710 plotDMR(variable + "X", minHits, options, filterName, maxBadLumiPixel, maxBadLumiStrip);
0711 plotDMR(variable + "Y", minHits, options, filterName, maxBadLumiPixel, maxBadLumiStrip);
0712 return;
0713 }
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723 TRegexp layer_re("layer=[0-9]+");
0724 bool plotPlain = false, plotSplits = false, plotLayers = false;
0725 int plotLayerN = 0;
0726 Ssiz_t index, len;
0727 if (options.find("plain") != std::string::npos) {
0728 plotPlain = true;
0729 }
0730 if (options.find("split") != std::string::npos) {
0731 plotSplits = true;
0732 }
0733 if (options.find("layers") != std::string::npos) {
0734 plotLayers = true;
0735 }
0736 if ((index = layer_re.Index(options, &len)) != -1) {
0737 if (plotLayers) {
0738 std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
0739 } else {
0740 std::string substr = options.substr(index + 6, len - 6);
0741 plotLayerN = atoi(substr.c_str());
0742 }
0743 }
0744
0745
0746
0747 if (!plotPlain && !plotSplits) {
0748 plotPlain = true;
0749 }
0750
0751
0752
0753 static bool plotSplitsFor[6] = {true, true, true, false, true, false};
0754
0755 DMRPlotInfo plotinfo;
0756
0757 gStyle->SetOptStat(0);
0758
0759 TCanvas c("canv", "canv");
0760
0761 plotinfo.variable = variable;
0762 plotinfo.minHits = minHits;
0763 plotinfo.plotPlain = plotPlain;
0764 plotinfo.plotLayers = plotLayers;
0765 plotinfo.filterName = filterName;
0766 plotinfo.maxBadLumiPixel = maxBadLumiPixel;
0767 plotinfo.maxBadLumiStrip = maxBadLumiStrip;
0768
0769
0770
0771
0772 if (variable == "meanX") {
0773 plotinfo.nbins = 50;
0774 plotinfo.min = -0.001;
0775 plotinfo.max = 0.001;
0776 } else if (variable == "meanY") {
0777 plotinfo.nbins = 50;
0778 plotinfo.min = -0.005;
0779 plotinfo.max = 0.005;
0780 } else if (variable == "medianX")
0781 if (plotSplits) {
0782 plotinfo.nbins = 50;
0783 plotinfo.min = -0.0005;
0784 plotinfo.max = 0.0005;
0785 } else {
0786 plotinfo.nbins = 100;
0787 plotinfo.min = -0.001;
0788 plotinfo.max = 0.001;
0789 }
0790 else if (variable == "medianY")
0791 if (plotSplits) {
0792 plotinfo.nbins = 50;
0793 plotinfo.min = -0.0005;
0794 plotinfo.max = 0.0005;
0795 } else {
0796 plotinfo.nbins = 100;
0797 plotinfo.min = -0.001;
0798 plotinfo.max = 0.001;
0799 }
0800 else if (variable == "meanNormX") {
0801 plotinfo.nbins = 100;
0802 plotinfo.min = -2.0;
0803 plotinfo.max = 2.0;
0804 } else if (variable == "meanNormY") {
0805 plotinfo.nbins = 100;
0806 plotinfo.min = -2.0;
0807 plotinfo.max = 2.0;
0808 } else if (variable == "rmsX") {
0809 plotinfo.nbins = 100;
0810 plotinfo.min = 0.0;
0811 plotinfo.max = 0.1;
0812 } else if (variable == "rmsY") {
0813 plotinfo.nbins = 100;
0814 plotinfo.min = 0.0;
0815 plotinfo.max = 0.1;
0816 } else if (variable == "rmsNormX") {
0817 plotinfo.nbins = 100;
0818 plotinfo.min = 0.3;
0819 plotinfo.max = 1.8;
0820 } else if (variable == "rmsNormY") {
0821 plotinfo.nbins = 100;
0822 plotinfo.min = 0.3;
0823 plotinfo.max = 1.8;
0824 } else {
0825 std::cerr << "Unknown variable " << variable << std::endl;
0826 plotinfo.nbins = 100;
0827 plotinfo.min = -0.1;
0828 plotinfo.max = 0.1;
0829 }
0830
0831 for (int i = 1; i <= 6; ++i) {
0832
0833 if (!plotinfo.filterName.empty()) {
0834 if (i == 1 || i == 2) {
0835 if (variable == "medianX") {
0836 if (plotSplits) {
0837 plotinfo.nbins = 50;
0838 } else {
0839 plotinfo.nbins = 50;
0840 }
0841 } else if (variable == "medianY") {
0842 if (plotSplits) {
0843 plotinfo.nbins = 50;
0844 } else {
0845 plotinfo.nbins = 25;
0846 }
0847 }
0848 } else if (i == 3 || i == 4 || i == 5 || i == 6) {
0849 if (variable == "medianX" || variable == "medianY") {
0850 if (plotSplits) {
0851 plotinfo.nbins = 50;
0852 } else {
0853 plotinfo.nbins = 25;
0854 }
0855 }
0856 }
0857 }
0858
0859
0860 if (i != 1 && i != 2 && !variable.empty() && variable[variable.length() - 1] == 'Y') {
0861 continue;
0862 }
0863
0864
0865 if (plotLayerN > maxNumberOfLayers(i)) {
0866 continue;
0867 }
0868
0869 plotinfo.plotSplits = plotSplits && plotSplitsFor[i - 1];
0870 if (!plotinfo.plotPlain && !plotinfo.plotSplits) {
0871 continue;
0872 }
0873
0874
0875
0876 bool hasheader = (TkAlStyle::legendheader != "");
0877
0878 int nPlots = 1;
0879 if (plotinfo.plotSplits) {
0880 nPlots = 3;
0881 }
0882
0883
0884 if (plotinfo.plotLayers) {
0885 nPlots *= maxNumberOfLayers(i);
0886 }
0887 nPlots *= sourceList.size();
0888 if (twolines_) {
0889 nPlots *= 2;
0890 }
0891 nPlots += hasheader;
0892
0893 double legendY = 0.80;
0894 if (nPlots > 3) {
0895 legendY -= 0.01 * (nPlots - 3);
0896 }
0897 if (bigtext_) {
0898 legendY -= 0.05;
0899 }
0900 if (legendY < 0.6) {
0901 std::cerr << "Warning: Huge legend!" << std::endl;
0902 legendY = 0.6;
0903 }
0904
0905 THStack hstack("hstack", "hstack");
0906 plotinfo.maxY = 0;
0907 plotinfo.subDetId = i;
0908 plotinfo.legend = new TLegend(0.17, legendY, 0.85, 0.88);
0909 plotinfo.legend->SetNColumns(2);
0910 if (hasheader)
0911 plotinfo.legend->SetHeader(TkAlStyle::legendheader);
0912 if (bigtext_)
0913 plotinfo.legend->SetTextSize(TkAlStyle::textSize);
0914 plotinfo.legend->SetFillStyle(0);
0915 plotinfo.hstack = &hstack;
0916 plotinfo.h = plotinfo.h1 = plotinfo.h2 = nullptr;
0917 plotinfo.firsthisto = true;
0918
0919 openSummaryFile();
0920 vmean.clear();
0921 vrms.clear();
0922 vdeltamean.clear();
0923 vmeanerror.clear();
0924 vPValueEqualSplitMeans.clear(), vAlignmentUncertainty.clear();
0925 vPValueMeanEqualIdeal.clear();
0926 vPValueRMSEqualIdeal.clear();
0927
0928 std::string stringsubdet;
0929 switch (i) {
0930 case 1:
0931 stringsubdet = "BPIX";
0932 break;
0933 case 2:
0934 stringsubdet = "FPIX";
0935 break;
0936 case 3:
0937 stringsubdet = "TIB";
0938 break;
0939 case 4:
0940 stringsubdet = "TID";
0941 break;
0942 case 5:
0943 stringsubdet = "TOB";
0944 break;
0945 case 6:
0946 stringsubdet = "TEC";
0947 break;
0948 }
0949
0950 for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0951 plotinfo.vars = *it;
0952 plotinfo.h1 = plotinfo.h2 = plotinfo.h = nullptr;
0953
0954 int minlayer = plotLayers ? 1 : plotLayerN;
0955
0956 if (plotinfo.plotPlain)
0957 minlayer = 0;
0958 int maxlayer = plotLayers ? numberOfLayers(plotinfo.vars->getPhase(), plotinfo.subDetId) : plotLayerN;
0959
0960 for (int layer = minlayer; layer <= maxlayer; layer++) {
0961 if (plotinfo.plotPlain) {
0962 plotDMRHistogram(plotinfo, 0, layer, stringsubdet);
0963 }
0964
0965 if (plotinfo.plotSplits) {
0966 plotDMRHistogram(plotinfo, -1, layer, stringsubdet);
0967 plotDMRHistogram(plotinfo, 1, layer, stringsubdet);
0968 }
0969
0970 if (plotinfo.plotPlain) {
0971 if (plotinfo.h) {
0972 setDMRHistStyleAndLegend(plotinfo.h, plotinfo, 0, layer);
0973 } else {
0974 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") &&
0975 layer == 0) {
0976 vmean.push_back(nan(""));
0977 vrms.push_back(nan(""));
0978 vmeanerror.push_back(nan(""));
0979 vAlignmentUncertainty.push_back(nan(""));
0980 vPValueMeanEqualIdeal.push_back(nan(""));
0981 vPValueRMSEqualIdeal.push_back(nan(""));
0982 if (plotinfo.plotSplits && plotinfo.plotPlain) {
0983 vdeltamean.push_back(nan(""));
0984 vPValueEqualSplitMeans.push_back(nan(""));
0985 }
0986 }
0987 }
0988 }
0989
0990 if (plotinfo.plotSplits) {
0991
0992 if (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && !plotinfo.plotPlain) {
0993 std::ostringstream legend;
0994 std::string unit = " #mum";
0995 legend.precision(3);
0996 legend << std::fixed;
0997 float factor = 10000.0f;
0998 if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" ||
0999 plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY") {
1000 factor = 1.0f;
1001 unit = "";
1002 }
1003 float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
1004 legend << plotinfo.vars->getName();
1005 if (layer > 0) {
1006
1007 if (i == 4 || i == 6)
1008 legend << ", disc ";
1009 else
1010 legend << ", layer ";
1011 legend << layer;
1012 }
1013 plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
1014 legend.str("");
1015 legend << "#Delta#mu = " << deltamu << unit;
1016 plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
1017
1018 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && !plotLayers && layer == 0) {
1019 vdeltamean.push_back(deltamu);
1020 }
1021 }
1022 if (plotinfo.h1) {
1023 setDMRHistStyleAndLegend(plotinfo.h1, plotinfo, -1, layer);
1024 }
1025 if (plotinfo.h2) {
1026 setDMRHistStyleAndLegend(plotinfo.h2, plotinfo, 1, layer);
1027 }
1028 }
1029 }
1030 }
1031
1032 if (hstack.GetHists() != nullptr && hstack.GetHists()->GetSize() != 0) {
1033 hstack.Draw("nostack");
1034 hstack.SetMaximum(plotinfo.maxY * 1.3);
1035 setTitleStyle(hstack, variable.c_str(), "#modules", plotinfo.subDetId);
1036 setHistStyle(*hstack.GetHistogram(), variable.c_str(), "#modules", 1);
1037
1038 plotinfo.legend->Draw();
1039 } else {
1040
1041 plotinfo.h = new TH1F("defhist", "Empty default histogram", plotinfo.nbins, plotinfo.min, plotinfo.max);
1042 plotinfo.h->SetMaximum(10);
1043 if (plotinfo.variable.find("Norm") == std::string::npos)
1044 scaleXaxis(plotinfo.h, 10000);
1045 setTitleStyle(*plotinfo.h, variable.c_str(), "#modules", plotinfo.subDetId);
1046 setHistStyle(*plotinfo.h, variable.c_str(), "#modules", 1);
1047 plotinfo.h->Draw();
1048 }
1049
1050 std::ostringstream plotName;
1051 plotName << outputDir << "/D";
1052
1053 if (variable == "medianX")
1054 plotName << "medianR_";
1055 else if (variable == "medianY")
1056 plotName << "medianYR_";
1057 else if (variable == "meanX")
1058 plotName << "meanR_";
1059 else if (variable == "meanY")
1060 plotName << "meanYR_";
1061 else if (variable == "meanNormX")
1062 plotName << "meanNR_";
1063 else if (variable == "meanNormY")
1064 plotName << "meanNYR_";
1065 else if (variable == "rmsX")
1066 plotName << "rmsR_";
1067 else if (variable == "rmsY")
1068 plotName << "rmsYR_";
1069 else if (variable == "rmsNormX")
1070 plotName << "rmsNR_";
1071 else if (variable == "rmsNormY")
1072 plotName << "rmsNYR_";
1073
1074 TString subdet;
1075 switch (i) {
1076 case 1:
1077 subdet = "BPIX";
1078 break;
1079 case 2:
1080 subdet = "FPIX";
1081 break;
1082 case 3:
1083 subdet = "TIB";
1084 break;
1085 case 4:
1086 subdet = "TID";
1087 break;
1088 case 5:
1089 subdet = "TOB";
1090 break;
1091 case 6:
1092 subdet = "TEC";
1093 break;
1094 }
1095
1096 plotName << subdet;
1097
1098 if (plotPlain && !plotSplits) {
1099 plotName << "_plain";
1100 } else if (!plotPlain && plotSplits) {
1101 plotName << "_split";
1102 }
1103 if (plotLayers) {
1104
1105 if (i == 4 || i == 6)
1106 plotName << "_discs";
1107 else
1108 plotName << "_layers";
1109 }
1110 if (plotLayerN > 0) {
1111
1112 if (i == 4 || i == 6)
1113 plotName << "_disc";
1114 else
1115 plotName << "_layer";
1116 plotName << plotLayerN;
1117 }
1118
1119
1120 c.Update();
1121 c.Print((plotName.str() + ".png").c_str());
1122 c.Print((plotName.str() + ".eps").c_str());
1123 c.Print((plotName.str() + ".pdf").c_str());
1124
1125
1126 TFile f((plotName.str() + ".root").c_str(), "recreate");
1127 c.Write();
1128 f.Close();
1129
1130
1131 delete plotinfo.h;
1132 delete plotinfo.h1;
1133 delete plotinfo.h2;
1134
1135 if (!vmean.empty()) {
1136 summaryfile << " mu_" << subdet;
1137 if (plotinfo.variable == "medianY")
1138 summaryfile << "_y";
1139 summaryfile << " (um)\t"
1140 << "latexname=$\\mu_\\text{" << subdet << "}";
1141 if (plotinfo.variable == "medianY")
1142 summaryfile << "^{y}";
1143 summaryfile << "$ ($\\mu$m)\t"
1144 << "format={:.3g}\t"
1145 << "latexformat=${:.3g}$";
1146 for (auto mu : vmean)
1147 summaryfile << "\t" << mu;
1148 summaryfile << "\n";
1149 }
1150 if (!vrms.empty()) {
1151 summaryfile << "sigma_" << subdet;
1152 if (plotinfo.variable == "medianY")
1153 summaryfile << "_y";
1154 summaryfile << " (um)\t"
1155 << "latexname=$\\sigma_\\text{" << subdet << "}";
1156 if (plotinfo.variable == "medianY")
1157 summaryfile << "^{y}";
1158 summaryfile << "$ ($\\mu$m)\t"
1159 << "format={:.3g}\t"
1160 << "latexformat=${:.3g}$";
1161 for (auto sigma : vrms)
1162 summaryfile << "\t" << sigma;
1163 summaryfile << "\n";
1164 }
1165 if (!vdeltamean.empty()) {
1166 summaryfile << " dmu_" << subdet;
1167 if (plotinfo.variable == "medianY")
1168 summaryfile << "_y";
1169 summaryfile << " (um)\t"
1170 << "latexname=$\\Delta\\mu_\\text{" << subdet << "}";
1171 if (plotinfo.variable == "medianY")
1172 summaryfile << "^{y}";
1173 summaryfile << "$ ($\\mu$m)\t"
1174 << "format={:.3g}\t"
1175 << "latexformat=${:.3g}$";
1176 for (auto dmu : vdeltamean)
1177 summaryfile << "\t" << dmu;
1178 summaryfile << "\n";
1179 }
1180 if (!vmeanerror.empty()) {
1181 summaryfile << " sigma_mu_" << subdet;
1182 if (plotinfo.variable == "medianY")
1183 summaryfile << "_y";
1184 summaryfile << " (um)\t"
1185 << "latexname=$\\sigma\\mu_\\text{" << subdet << "}";
1186 if (plotinfo.variable == "medianY")
1187 summaryfile << "^{y}";
1188 summaryfile << "$ ($\\mu$m)\t"
1189 << "format={:.3g}\t"
1190 << "latexformat=${:.3g}$";
1191 for (auto dmu : vmeanerror)
1192 summaryfile << "\t" << dmu;
1193 summaryfile << "\n";
1194 }
1195 if (!vPValueEqualSplitMeans.empty()) {
1196 summaryfile << " p_delta_mu_equal_zero_" << subdet;
1197 if (plotinfo.variable == "medianY")
1198 summaryfile << "_y";
1199 summaryfile << "\t"
1200 << "latexname=$P(\\delta\\mu_\\text{" << subdet << "}=0)";
1201 if (plotinfo.variable == "medianY")
1202 summaryfile << "^{y}";
1203 summaryfile << "$\t"
1204 << "format={:.3g}\t"
1205 << "latexformat=${:.3g}$";
1206 for (auto dmu : vPValueEqualSplitMeans)
1207 summaryfile << "\t" << dmu;
1208 summaryfile << "\n";
1209 }
1210 if (!vAlignmentUncertainty.empty()) {
1211 summaryfile << " alignment_uncertainty_" << subdet;
1212 if (plotinfo.variable == "medianY")
1213 summaryfile << "_y";
1214 summaryfile << " (um)\t"
1215 << "latexname=$\\sigma_\\text{align}_\\text{" << subdet << "}";
1216 if (plotinfo.variable == "medianY")
1217 summaryfile << "^{y}";
1218 summaryfile << "$ ($\\mu$m)\t"
1219 << "format={:.3g}\t"
1220 << "latexformat=${:.3g}$";
1221 for (auto dmu : vAlignmentUncertainty)
1222 summaryfile << "\t" << dmu;
1223 summaryfile << "\n";
1224 }
1225 if (!vPValueMeanEqualIdeal.empty()) {
1226 summaryfile << " p_mean_equals_ideal_" << subdet;
1227 if (plotinfo.variable == "medianY")
1228 summaryfile << "_y";
1229 summaryfile << "\t"
1230 << "latexname=$P(\\mu_\\text{" << subdet << "}=\\mu_\\text{ideal})";
1231 if (plotinfo.variable == "medianY")
1232 summaryfile << "^{y}";
1233 summaryfile << "$\t"
1234 << "format={:.3g}\t"
1235 << "latexformat=${:.3g}$";
1236 for (auto dmu : vPValueMeanEqualIdeal)
1237 summaryfile << "\t" << dmu;
1238 summaryfile << "\n";
1239 }
1240 if (!vPValueRMSEqualIdeal.empty()) {
1241 summaryfile << " p_RMS_equals_ideal_" << subdet;
1242 if (plotinfo.variable == "medianY")
1243 summaryfile << "_y";
1244 summaryfile << "\t"
1245 << "latexname=$P(\\sigma_\\text{" << subdet << "}=\\sigma_\\text{ideal})";
1246 if (plotinfo.variable == "medianY")
1247 summaryfile << "^{y}";
1248 summaryfile << "$\t"
1249 << "format={:.3g}\t"
1250 << "latexformat=${:.3g}$";
1251 for (auto dmu : vPValueRMSEqualIdeal)
1252 summaryfile << "\t" << dmu;
1253 summaryfile << "\n";
1254 }
1255 }
1256 }
1257
1258
1259 void PlotAlignmentValidation::plotChi2(const char* inputFile) {
1260
1261
1262
1263 Bool_t errorflag = kFALSE;
1264 TFile* fi1 = TFile::Open(inputFile, "read");
1265 if (fi1 != nullptr) {
1266 if (fi1->GetDirectory("TrackerOfflineValidation/GlobalTrackVariables") == nullptr) {
1267 errorflag = kTRUE;
1268 }
1269 } else {
1270 errorflag = kTRUE;
1271 }
1272 if (errorflag) {
1273 std::cout << "PlotAlignmentValidation::plotChi2: Can't find GlobalTrackVariables given file,"
1274 << " no chi^2-plots produced" << std::endl;
1275 return;
1276 }
1277
1278 auto normchi = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_normchi2");
1279
1280 normchi->GetFrame()->ResetBit(kCanDelete);
1281
1282 auto chiprob = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_chi2Prob");
1283
1284 chiprob->GetFrame()->ResetBit(kCanDelete);
1285
1286 if (normchi == nullptr || chiprob == nullptr) {
1287 errorflag = kTRUE;
1288 }
1289 if (errorflag) {
1290 std::cout << "PlotAlignmentValidation::plotChi2: Can't find data from given file,"
1291 << " no chi^2-plots produced" << std::endl;
1292 return;
1293 }
1294
1295 TLegend* legend = nullptr;
1296 for (auto primitive : *normchi->GetListOfPrimitives()) {
1297 legend = dynamic_cast<TLegend*>(primitive);
1298 if (legend)
1299 break;
1300 }
1301 if (legend) {
1302 openSummaryFile();
1303 summaryfile << "ntracks";
1304 for (auto alignment : sourceList) {
1305 summaryfile << "\t";
1306 TString title = alignment->getName();
1307 int color = alignment->getLineColor();
1308 int style = alignment->getLineStyle();
1309 bool foundit = false;
1310 for (auto entry : *legend->GetListOfPrimitives()) {
1311 TLegendEntry* legendentry = dynamic_cast<TLegendEntry*>(entry);
1312 assert(legendentry);
1313 TH1* h = dynamic_cast<TH1*>(legendentry->GetObject());
1314 if (!h)
1315 continue;
1316 if (legendentry->GetLabel() == title && h->GetLineColor() == color && h->GetLineStyle() == style) {
1317 foundit = true;
1318 summaryfile << h->GetEntries();
1319 break;
1320 }
1321 }
1322 if (!foundit) {
1323 summaryfile << 0;
1324 }
1325 }
1326 summaryfile << "\n";
1327 }
1328
1329 chiprob->Draw();
1330 normchi->Draw();
1331
1332
1333 normchi->Print((outputDir + "/h_normchi2.png").c_str());
1334 chiprob->Print((outputDir + "/h_chi2Prob.png").c_str());
1335 normchi->Print((outputDir + "/h_normchi2.eps").c_str());
1336 chiprob->Print((outputDir + "/h_chi2Prob.eps").c_str());
1337 normchi->Print((outputDir + "/h_normchi2.pdf").c_str());
1338 chiprob->Print((outputDir + "/h_chi2Prob.pdf").c_str());
1339
1340
1341 TFile fi2((outputDir + "/h_normchi2.root").c_str(), "recreate");
1342 normchi->Write();
1343 fi2.Close();
1344
1345 TFile fi3((outputDir + "/h_chi2Prob.root").c_str(), "recreate");
1346 chiprob->Write();
1347 fi3.Close();
1348
1349 fi1->Close();
1350 delete fi1;
1351 }
1352
1353
1354 THStack* PlotAlignmentValidation::addHists(
1355 const TString& selection, const TString& residType, TLegend** myLegend, bool printModuleIds, bool validforphase0) {
1356 enum ResidType {
1357 xPrimeRes,
1358 yPrimeRes,
1359 xPrimeNormRes,
1360 yPrimeNormRes,
1361 xRes,
1362 yRes,
1363 xNormRes,
1364 ResXvsXProfile,
1365 ResXvsYProfile,
1366 ResYvsXProfile,
1367 ResYvsYProfile
1368 };
1369 ResidType rType = xPrimeRes;
1370 if (residType == "xPrime")
1371 rType = xPrimeRes;
1372 else if (residType == "yPrime")
1373 rType = yPrimeRes;
1374 else if (residType == "xPrimeNorm")
1375 rType = xPrimeNormRes;
1376 else if (residType == "yPrimeNorm")
1377 rType = yPrimeNormRes;
1378 else if (residType == "x")
1379 rType = xRes;
1380 else if (residType == "y")
1381 rType = yRes;
1382 else if (residType == "xNorm")
1383 rType = xNormRes;
1384
1385 else if (residType == "ResXvsXProfile")
1386 rType = ResXvsXProfile;
1387 else if (residType == "ResYvsXProfile")
1388 rType = ResYvsXProfile;
1389 else if (residType == "ResXvsYProfile")
1390 rType = ResXvsYProfile;
1391 else if (residType == "ResYvsYProfile")
1392 rType = ResYvsYProfile;
1393 else {
1394 std::cout << "PlotAlignmentValidation::addHists: Unknown residual type " << residType << std::endl;
1395 return nullptr;
1396 }
1397
1398 std::cout << "PlotAlignmentValidation::addHists: using selection " << selection << std::endl;
1399 THStack* retHistoStack = new THStack("hstack", "");
1400 if (myLegend != nullptr)
1401 if (*myLegend == nullptr) {
1402 *myLegend = new TLegend(0.17, 0.80, 0.85, 0.88);
1403 }
1404
1405 for (std::vector<TkOfflineVariables*>::iterator itSourceFile = sourceList.begin(); itSourceFile != sourceList.end();
1406 ++itSourceFile) {
1407 std::vector<TString> histnames;
1408
1409 TFile* f = (*itSourceFile)->getFile();
1410 TTree* tree = (*itSourceFile)->getTree();
1411 int myLineColor = (*itSourceFile)->getLineColor();
1412 int myLineStyle = (*itSourceFile)->getLineStyle();
1413 TString myLegendName = (*itSourceFile)->getName();
1414 TH1* h = nullptr;
1415 UInt_t nEmpty = 0;
1416 Long64_t nentries = tree->GetEntriesFast();
1417 if (!f || !tree) {
1418 std::cout << "PlotAlignmentValidation::addHists: no tree or no file" << std::endl;
1419 return nullptr;
1420 }
1421
1422 bool histnamesfilled = false;
1423 int phase = (bool)(f->Get("TrackerOfflineValidation/Pixel/P1PXBBarrel_1"));
1424 if (residType.Contains("Res") && residType.Contains("Profile")) {
1425 TString basename = TString(residType)
1426 .ReplaceAll("Res", "p_res")
1427 .ReplaceAll("vs", "")
1428 .ReplaceAll("Profile", "_");
1429 if (selection == "subDetId==1") {
1430 if (phase == 1)
1431 histnames.push_back(TString(basename) += "P1PXBBarrel_1");
1432 else
1433 histnames.push_back(TString(basename) += "TPBBarrel_1");
1434 histnamesfilled = true;
1435 } else if (selection == "subDetId==2") {
1436 if (phase == 1) {
1437 histnames.push_back(TString(basename) += "P1PXECEndcap_2");
1438 histnames.push_back(TString(basename) += "P1PXECEndcap_3");
1439 } else {
1440 histnames.push_back(TString(basename) += "TPEEndcap_2");
1441 histnames.push_back(TString(basename) += "TPEEndcap_3");
1442 }
1443 histnamesfilled = true;
1444 } else if (selection == "subDetId==3") {
1445 histnames.push_back(TString(basename) += "TIBBarrel_1");
1446 histnamesfilled = true;
1447 } else if (selection == "subDetId==4") {
1448 histnames.push_back(TString(basename) += "TIDEndcap_2");
1449 histnames.push_back(TString(basename) += "TIDEndcap_3");
1450 histnamesfilled = true;
1451 } else if (selection == "subDetId==5") {
1452 histnames.push_back(TString(basename) += "TOBBarrel_4");
1453 histnamesfilled = true;
1454 } else if (selection == "subDetId==6") {
1455 histnames.push_back(TString(basename) += "TECEndcap_5");
1456 histnames.push_back(TString(basename) += "TECEndcap_6");
1457 histnamesfilled = true;
1458 } else if (selection == "subDetId==6 && ring <= 4") {
1459
1460 for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1461 for (int iDisk = 1; iDisk <= 9; iDisk++)
1462 for (int iSide = 1; iSide <= 2; iSide++)
1463 for (int iPetal = 1; iPetal <= 8; iPetal++)
1464 for (int iRing = 1; iRing <= 4 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9); iRing++)
1465
1466
1467 {
1468 std::stringstream s;
1469 s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
1470 << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1471 histnames.push_back(TString(s.str()));
1472 }
1473 histnamesfilled = true;
1474 } else if (selection == "subDetId==6 && ring > 4") {
1475
1476 for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1477 for (int iDisk = 1; iDisk <= 9; iDisk++)
1478 for (int iSide = 1; iSide <= 2; iSide++)
1479 for (int iPetal = 1; iPetal <= 8; iPetal++)
1480 for (int iRing = 5 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1481 iRing <= 7 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1482 iRing++)
1483
1484
1485 {
1486 std::stringstream s;
1487 s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
1488 << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1489 histnames.push_back(TString(s.str()));
1490 }
1491 histnamesfilled = true;
1492 }
1493 }
1494
1495 Long64_t nSel = 0;
1496 if (histnamesfilled && !histnames.empty()) {
1497 nSel = (Long64_t)histnames.size();
1498 }
1499 if (!histnamesfilled) {
1500
1501
1502 nSel = tree->Draw("Entry$", selection, "goff");
1503 if (nSel == -1)
1504 return nullptr;
1505 if (nSel == 0) {
1506 std::cout << "PlotAlignmentValidation::addHists: no selected module." << std::endl;
1507 return nullptr;
1508 }
1509
1510 const std::vector<double> selected(tree->GetV1(), tree->GetV1() + nSel);
1511
1512 std::vector<double>::const_iterator iterEnt = selected.begin();
1513
1514
1515
1516 TkOffTreeVariables* treeMem = nullptr;
1517 tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
1518 for (Long64_t i = 0; i < nentries; i++) {
1519 if (i < *iterEnt - 0.1
1520 || iterEnt == selected.end()) {
1521 continue;
1522 } else if (TMath::Abs(i - *iterEnt) < 0.11) {
1523 ++iterEnt;
1524 } else
1525 std::cout << "Must not happen: " << i << " " << *iterEnt << std::endl;
1526
1527 tree->GetEntry(i);
1528 if (printModuleIds) {
1529 std::cout << treeMem->moduleId << ": " << treeMem->entries << " entries" << std::endl;
1530 }
1531 if (treeMem->entries <= 0) {
1532 ++nEmpty;
1533 continue;
1534 }
1535 TString hName;
1536 switch (rType) {
1537 case xPrimeRes:
1538 hName = treeMem->histNameX.c_str();
1539 break;
1540 case yPrimeRes:
1541 hName = treeMem->histNameY.c_str();
1542 break;
1543 case xPrimeNormRes:
1544 hName = treeMem->histNameNormX.c_str();
1545 break;
1546 case yPrimeNormRes:
1547 hName = treeMem->histNameNormY.c_str();
1548 break;
1549 case xRes:
1550 hName = treeMem->histNameLocalX.c_str();
1551 break;
1552 case yRes:
1553 hName = treeMem->histNameLocalY.c_str();
1554 break;
1555 case xNormRes:
1556 hName = treeMem->histNameNormLocalX.c_str();
1557 break;
1558
1559 case ResXvsXProfile:
1560 hName = treeMem->profileNameResXvsX.c_str();
1561 break;
1562 case ResXvsYProfile:
1563 hName = treeMem->profileNameResXvsY.c_str();
1564 break;
1565 case ResYvsXProfile:
1566 hName = treeMem->profileNameResYvsX.c_str();
1567 break;
1568 case ResYvsYProfile:
1569 hName = treeMem->profileNameResYvsY.c_str();
1570 break;
1571 }
1572 histnames.push_back(hName);
1573 }
1574 }
1575
1576 for (std::vector<TString>::iterator ithistname = histnames.begin(); ithistname != histnames.end(); ++ithistname) {
1577 if (phase == 0 && !validforphase0)
1578 break;
1579 TH1* newHist;
1580 if (ithistname->Contains("/")) {
1581 newHist = (TH1*)f->Get(*ithistname);
1582 } else {
1583 TKey* histKey = f->FindKeyAny(*ithistname);
1584 newHist = (histKey ? static_cast<TH1*>(histKey->ReadObj()) : nullptr);
1585 }
1586 if (!newHist) {
1587 std::cout << "Hist " << *ithistname << " not found in file, break loop." << std::endl;
1588 break;
1589 }
1590 if (newHist->GetEntries() == 0) {
1591 nEmpty++;
1592 continue;
1593 }
1594 newHist->SetLineColor(myLineColor);
1595 newHist->SetLineStyle(myLineStyle);
1596 if (!h) {
1597 TString name(newHist->GetName());
1598 Ssiz_t pos_ = 0;
1599 for (UInt_t i2 = 0; i2 < 3; ++i2)
1600 pos_ = name.Index("_", pos_ + 1);
1601 name = name(0, pos_);
1602 h = static_cast<TH1*>(newHist->Clone("summed_" + name));
1603
1604
1605 } else {
1606 h->Add(newHist);
1607 }
1608 delete newHist;
1609 }
1610
1611 std::cout << "PlotAlignmentValidation::addHists"
1612 << "Result is merged from " << nSel - nEmpty << " hists, " << nEmpty << " hists were empty." << std::endl;
1613
1614 if (nSel - nEmpty == 0)
1615 continue;
1616
1617 if (myLegend != nullptr)
1618 (*myLegend)->AddEntry(h, myLegendName, "L");
1619
1620 retHistoStack->Add(h);
1621 }
1622
1623 return retHistoStack;
1624 }
1625
1626
1627
1628
1629
1630 TF1* PlotAlignmentValidation::fitGauss(TH1* hist, int color) {
1631
1632
1633
1634 if (!hist || hist->GetEntries() < 20)
1635 return nullptr;
1636
1637 float mean = hist->GetMean();
1638 float sigma = hist->GetRMS();
1639 std::string functionname = "gaussian_";
1640 functionname += hist->GetName();
1641 TF1* func = new TF1(functionname.c_str(), "gaus", mean - 2. * sigma, mean + 2. * sigma);
1642
1643 func->SetLineColor(color);
1644 func->SetLineStyle(2);
1645 if (0 == hist->Fit(func, "QNR")) {
1646 mean = func->GetParameter(1);
1647 sigma = func->GetParameter(2);
1648
1649 func->SetRange(mean - 3. * sigma, mean + 3. * sigma);
1650
1651
1652 if (0 == hist->Fit(func, "Q0ILR")) {
1653 if (hist->GetFunction(func->GetName())) {
1654
1655 }
1656 }
1657 }
1658 return func;
1659 }
1660
1661
1662
1663
1664
1665 void PlotAlignmentValidation::storeHistogramInRootfile(TH1* hist) {
1666
1667 rootsummaryfile->cd();
1668 hist->Write();
1669 }
1670
1671
1672 void PlotAlignmentValidation::scaleXaxis(TH1* hist, Int_t scale) {
1673 Double_t xmin = hist->GetXaxis()->GetXmin();
1674 Double_t xmax = hist->GetXaxis()->GetXmax();
1675 hist->GetXaxis()->SetLimits(xmin * scale, xmax * scale);
1676 }
1677
1678
1679 TObject* PlotAlignmentValidation::findObjectFromCanvas(TCanvas* canv, const char* className, Int_t n) {
1680
1681 TIter next(canv->GetListOfPrimitives());
1682 TObject* obj = nullptr;
1683 Int_t found = 0;
1684 while ((obj = next())) {
1685 if (strncmp(obj->ClassName(), className, 10) == 0) {
1686 if (++found == n)
1687 return obj;
1688 }
1689 }
1690
1691 return nullptr;
1692 }
1693
1694
1695 void PlotAlignmentValidation::setTitleStyle(
1696 TNamed& hist, const char* titleX, const char* titleY, int subDetId, bool isSurfaceDeformation, TString secondline) {
1697 std::stringstream title_Xaxis;
1698 std::stringstream title_Yaxis;
1699 TString titleXAxis = titleX;
1700 TString titleYAxis = titleY;
1701 if (titleXAxis != "" && titleYAxis != "")
1702 std::cout << "plot " << titleXAxis << " vs " << titleYAxis << std::endl;
1703
1704 hist.SetTitle("");
1705 TkAlStyle::drawStandardTitle();
1706
1707
1708 TString subD;
1709 switch (subDetId) {
1710 case 1:
1711 subD = "BPIX";
1712 break;
1713 case 2:
1714 subD = "FPIX";
1715 break;
1716 case 3:
1717 subD = "TIB";
1718 break;
1719 case 4:
1720 subD = "TID";
1721 break;
1722 case 5:
1723 subD = "TOB";
1724 break;
1725 case 6:
1726 subD = "TEC";
1727 break;
1728 }
1729
1730 TPaveText* text2;
1731 if (!isSurfaceDeformation) {
1732 text2 = new TPaveText(0.7, 0.3, 0.9, 0.6, "brNDC");
1733 } else {
1734 std::cout << "Surface Deformation" << std::endl;
1735 text2 = new TPaveText(0.8, 0.75, 0.9, 0.9, "brNDC");
1736 }
1737 text2->SetTextSize(0.06);
1738 text2->SetTextFont(42);
1739 text2->SetFillStyle(0);
1740 text2->SetBorderSize(0);
1741 text2->SetMargin(0.01);
1742 text2->SetTextAlign(12);
1743 text2->AddText(0.01, 0.75, subD);
1744 if (secondline != "") {
1745 text2->AddText(0.01, 0.25, secondline);
1746 }
1747 text2->Draw();
1748 }
1749
1750
1751
1752
1753
1754 void PlotAlignmentValidation::setHistStyle(TH1& hist, const char* titleX, const char* titleY, int color) {
1755 std::stringstream title_Xaxis;
1756 std::stringstream title_Yaxis;
1757 TString titleXAxis = titleX;
1758 TString titleYAxis = titleY;
1759
1760 if (titleXAxis.Contains("Phi"))
1761 title_Xaxis << titleX << "[rad]";
1762 else if (titleXAxis.Contains("meanX"))
1763 title_Xaxis << "#LTx'_{pred}-x'_{hit}#GT[#mum]";
1764 else if (titleXAxis.Contains("meanY"))
1765 title_Xaxis << "#LTy'_{pred}-y'_{hit}#GT[#mum]";
1766 else if (titleXAxis.Contains("rmsX"))
1767 title_Xaxis << "RMS(x'_{pred}-x'_{hit})[#mum]";
1768 else if (titleXAxis.Contains("rmsY"))
1769 title_Xaxis << "RMS(y'_{pred}-y'_{hit})[#mum]";
1770 else if (titleXAxis.Contains("meanNormX"))
1771 title_Xaxis << "#LTx'_{pred}-x'_{hit}/#sigma#GT";
1772 else if (titleXAxis.Contains("meanNormY"))
1773 title_Xaxis << "#LTy'_{pred}-y'_{hit}/#sigma#GT";
1774 else if (titleXAxis.Contains("rmsNormX"))
1775 title_Xaxis << "RMS(x'_{pred}-x'_{hit}/#sigma)";
1776 else if (titleXAxis.Contains("rmsNormY"))
1777 title_Xaxis << "RMS(y'_{pred}-y'_{hit}/#sigma)";
1778 else if (titleXAxis.Contains("meanLocalX"))
1779 title_Xaxis << "#LTx_{pred}-x_{hit}#GT[#mum]";
1780 else if (titleXAxis.Contains("rmsLocalX"))
1781 title_Xaxis << "RMS(x_{pred}-x_{hit})[#mum]";
1782 else if (titleXAxis.Contains("meanNormLocalX"))
1783 title_Xaxis << "#LTx_{pred}-x_{hit}/#sigma#GT[#mum]";
1784 else if (titleXAxis.Contains("rmsNormLocalX"))
1785 title_Xaxis << "RMS(x_{pred}-x_{hit}/#sigma)[#mum]";
1786 else if (titleXAxis.Contains("medianX"))
1787 title_Xaxis << "median(x'_{pred}-x'_{hit})[#mum]";
1788 else if (titleXAxis.Contains("medianY"))
1789 title_Xaxis << "median(y'_{pred}-y'_{hit})[#mum]";
1790 else
1791 title_Xaxis << titleX << "[cm]";
1792
1793 if (hist.IsA()->InheritsFrom(TH1F::Class()))
1794 hist.SetLineColor(color);
1795 if (hist.IsA()->InheritsFrom(TProfile::Class())) {
1796 hist.SetMarkerStyle(20);
1797 hist.SetMarkerSize(0.8);
1798 hist.SetMarkerColor(color);
1799 }
1800
1801 hist.GetXaxis()->SetTitle((title_Xaxis.str()).c_str());
1802
1803 double binning = (hist.GetXaxis()->GetXmax() - hist.GetXaxis()->GetXmin()) / hist.GetNbinsX();
1804 title_Yaxis.precision(2);
1805
1806 if (((titleYAxis.Contains("layer") || titleYAxis.Contains("ring")) && titleYAxis.Contains("subDetId")) ||
1807 titleYAxis.Contains("#modules")) {
1808 title_Yaxis << "number of modules";
1809 if (TString(title_Xaxis.str()).Contains("[#mum]"))
1810 title_Yaxis << " / " << binning << " #mum";
1811 else if (TString(title_Xaxis.str()).Contains("[cm]"))
1812 title_Yaxis << " / " << binning << " cm";
1813 else
1814 title_Yaxis << " / " << binning;
1815 } else
1816 title_Yaxis << titleY << "[cm]";
1817
1818 hist.GetYaxis()->SetTitle((title_Yaxis.str()).c_str());
1819
1820 hist.GetXaxis()->SetTitleFont(42);
1821 hist.GetYaxis()->SetTitleFont(42);
1822 }
1823
1824
1825 std::string PlotAlignmentValidation::getSelectionForDMRPlot(int minHits, int subDetId, int direction, int layer) {
1826 std::ostringstream builder;
1827 builder << "entries >= " << minHits;
1828 builder << " && subDetId == " << subDetId;
1829 if (direction != 0) {
1830 if (subDetId == 2) {
1831 builder << " && zDirection == " << direction;
1832 } else {
1833 builder << " && rDirection == " << direction;
1834 }
1835 }
1836 if (layer > 0) {
1837 builder << " && layer == " << layer;
1838 }
1839 return builder.str();
1840 }
1841
1842 std::string PlotAlignmentValidation::getVariableForDMRPlot(
1843 const std::string& histoname, const std::string& variable, int nbins, double min, double max) {
1844 std::ostringstream builder;
1845 builder << variable << ">>" << histoname << "(" << nbins << "," << min << "," << max << ")";
1846 return builder.str();
1847 }
1848
1849
1850 void PlotAlignmentValidation::setDMRHistStyleAndLegend(TH1F* h,
1851 PlotAlignmentValidation::DMRPlotInfo& plotinfo,
1852 int direction,
1853 int layer) {
1854 TF1* fitResults = nullptr;
1855
1856 h->SetDirectory(nullptr);
1857
1858
1859
1860 h->SetLineWidth((direction == 0 || (plotinfo.plotSplits && !plotinfo.plotPlain)) ? 2 : 1);
1861
1862
1863
1864
1865
1866 int linestyle = plotinfo.vars->getLineStyle() - 1, linestyleplus = 0;
1867 if (direction == 1) {
1868 linestyleplus = 1;
1869 }
1870 if (direction == -1) {
1871 linestyleplus = 2;
1872 }
1873 if (direction != 0 && plotinfo.plotSplits && !plotinfo.plotPlain) {
1874 linestyleplus--;
1875 }
1876 linestyle = (linestyle + linestyleplus) % 4 + 1;
1877
1878 int linecolor = plotinfo.vars->getLineColor();
1879 if (plotinfo.plotLayers && layer > 0) {
1880 linecolor += layer - 1;
1881 }
1882
1883 if (plotinfo.firsthisto) {
1884 setHistStyle(*h, plotinfo.variable.c_str(), "#modules", 1);
1885 plotinfo.firsthisto = false;
1886 }
1887
1888 h->SetLineColor(linecolor);
1889 h->SetLineStyle(linestyle);
1890
1891 if (plotinfo.maxY < h->GetMaximum()) {
1892 plotinfo.maxY = h->GetMaximum();
1893 }
1894
1895
1896 if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1897 plotinfo.variable == "meanY") {
1898 fitResults = fitGauss(h, linecolor);
1899 }
1900
1901 plotinfo.hstack->Add(h);
1902
1903 std::ostringstream legend;
1904 legend.precision(3);
1905 legend << std::fixed;
1906
1907
1908 if (direction == -1 && plotinfo.subDetId != 2) {
1909 legend << "rDirection < 0";
1910 } else if (direction == 1 && plotinfo.subDetId != 2) {
1911 legend << "rDirection > 0";
1912 } else if (direction == -1 && plotinfo.subDetId == 2) {
1913 legend << "zDirection < 0";
1914 } else if (direction == 1 && plotinfo.subDetId == 2) {
1915 legend << "zDirection > 0";
1916 } else {
1917 legend << plotinfo.vars->getName();
1918 if (layer > 0) {
1919
1920 if (plotinfo.subDetId == 4 || plotinfo.subDetId == 6)
1921 legend << ", disc ";
1922 else
1923 legend << ", layer ";
1924 legend << layer << "";
1925 }
1926 }
1927
1928 plotinfo.legend->AddEntry(h, legend.str().c_str(), "l");
1929 legend.str("");
1930
1931
1932 double mean = 0.0, meanerror = 0.0, rms = 0.0, rmserror = 0.0;
1933 TString rmsname, units;
1934 bool showdeltamu =
1935 (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && plotinfo.plotSplits && plotinfo.plotPlain && direction == 0);
1936 if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1937 plotinfo.variable == "meanY" || plotinfo.variable == "rmsX" || plotinfo.variable == "rmsY") {
1938 if (useFit_ && fitResults) {
1939 mean = fitResults->GetParameter(1) * 10000;
1940 meanerror = fitResults->GetParError(1) * 10000;
1941 rms = fitResults->GetParameter(2) * 10000;
1942 rmserror = fitResults->GetParError(2) * 10000;
1943 rmsname = "#sigma";
1944 delete fitResults;
1945 } else {
1946 mean = h->GetMean(1) * 10000;
1947 meanerror = h->GetMeanError(1) * 10000;
1948 rms = h->GetRMS(1) * 10000;
1949 rmserror = h->GetRMSError(1) * 10000;
1950 rmsname = "rms";
1951 }
1952 units = " #mum";
1953 } else if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
1954 plotinfo.variable == "rmsNormY") {
1955 mean = h->GetMean(1);
1956 meanerror = h->GetMeanError(1);
1957 rms = h->GetRMS(1);
1958 rmserror = h->GetRMSError(1);
1959 rmsname = "rms";
1960 units = "";
1961 }
1962 if (showMean_) {
1963 legend << " #mu = " << mean;
1964 if (showMeanError_)
1965 legend << " #pm " << meanerror;
1966 legend << units;
1967 if (showRMS_ || showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1968 legend << ", ";
1969 }
1970 if (showRMS_) {
1971 legend << " " << rmsname << " = " << rms;
1972 if (showRMSError_)
1973 legend << " #pm " << rmserror;
1974 legend << units;
1975 if (showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1976 legend << ", ";
1977 }
1978
1979 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && layer == 0 &&
1980 direction == 0) {
1981 vmean.push_back(mean);
1982 vrms.push_back(rms);
1983 vmeanerror.push_back(meanerror);
1984 TH1F* ideal = (TH1F*)plotinfo.hstack->GetHists()->At(0);
1985 TH1F* h = plotinfo.h;
1986 if (h->GetRMS() >= ideal->GetRMS()) {
1987 vAlignmentUncertainty.push_back(sqrt(pow(h->GetRMS(), 2) - pow(ideal->GetRMS(), 2)));
1988 } else {
1989 vAlignmentUncertainty.push_back(nan(""));
1990 }
1991 float p = (float)resampleTestOfEqualMeans(ideal, h, 10000);
1992 vPValueMeanEqualIdeal.push_back(p);
1993 p = resampleTestOfEqualRMS(ideal, h, 10000);
1994 vPValueRMSEqualIdeal.push_back(p);
1995 }
1996
1997
1998 if (showdeltamu) {
1999 float factor = 10000.0f;
2000 if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
2001 plotinfo.variable == "rmsNormY") {
2002 factor = 1.0f;
2003 }
2004 float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
2005 legend << "#Delta#mu = " << deltamu << units;
2006 if ((showModules_ || showUnderOverFlow_) && !twolines_)
2007 legend << ", ";
2008
2009 if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && layer == 0 &&
2010 direction == 0) {
2011 vdeltamean.push_back(deltamu);
2012 if (plotinfo.h1->GetEntries() && plotinfo.h2->GetEntries()) {
2013 float p = (float)resampleTestOfEqualMeans(plotinfo.h1, plotinfo.h2, 10000);
2014 vPValueEqualSplitMeans.push_back(p);
2015 }
2016 }
2017 }
2018
2019 if (twolines_) {
2020 plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
2021 plotinfo.legend->AddEntry((TObject*)nullptr, "", "");
2022 legend.str("");
2023 }
2024
2025 if (!showUnderOverFlow_ && showModules_) {
2026 legend << (int)h->GetEntries() << " modules";
2027 }
2028 if (showUnderOverFlow_) {
2029 if (showModules_) {
2030 legend << (int)h->GetEntries() << " modules ("
2031 << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " outside range)";
2032 } else {
2033 legend << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " modules outside range";
2034 }
2035 }
2036 plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
2037
2038
2039 if (plotinfo.variable.find("Norm") == std::string::npos)
2040 scaleXaxis(h, 10000);
2041 }
2042
2043
2044
2045
2046
2047
2048 void PlotAlignmentValidation::plotDMRHistogram(PlotAlignmentValidation::DMRPlotInfo& plotinfo,
2049 int direction,
2050 int layer,
2051 std::string subdet) {
2052 TH1F* h = nullptr;
2053
2054
2055 TString histoname = "";
2056 if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY")
2057 histoname = "median";
2058 else if (plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY")
2059 histoname = "DrmsNR";
2060 histoname += "_";
2061 histoname += plotinfo.vars->getName();
2062 histoname.ReplaceAll(" ", "_");
2063 histoname += "_";
2064 histoname += subdet.c_str();
2065 if (plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormY")
2066 histoname += "_y";
2067 if (layer != 0) {
2068 if (subdet == "TID" || subdet == "TEC")
2069 histoname += "_disc";
2070 else
2071 histoname += "_layer";
2072 histoname += std::to_string(layer);
2073 }
2074 if (direction == -1) {
2075 histoname += "_minus";
2076 } else if (direction == 1) {
2077 histoname += "_plus";
2078 } else {
2079 histoname += "";
2080 }
2081
2082
2083 std::string plotVariable =
2084 getVariableForDMRPlot(histoname.Data(), plotinfo.variable, plotinfo.nbins, plotinfo.min, plotinfo.max);
2085 std::string selection = "";
2086 if (plotinfo.filterName.empty()) {
2087
2088 selection = getSelectionForDMRPlot(plotinfo.minHits, plotinfo.subDetId, direction, layer);
2089 plotinfo.vars->getTree()->Draw(plotVariable.c_str(), selection.c_str(), "goff");
2090 if (gDirectory)
2091 gDirectory->GetObject(histoname.Data(), h);
2092 if (h && h->GetEntries() > 0) {
2093 if (direction == -1) {
2094 plotinfo.h1 = h;
2095 } else if (direction == 1) {
2096 plotinfo.h2 = h;
2097 } else {
2098 plotinfo.h = h;
2099 }
2100 }
2101 } else {
2102 plotinfo.vars->getTree()->Draw(
2103 plotVariable.c_str(), selection.c_str(), "goff");
2104 if (gDirectory)
2105 gDirectory->GetObject(histoname.Data(), h);
2106 h->Reset();
2107 std::cout << "Module filter enabled." << std::endl;
2108 TTreeReader reader(plotinfo.vars->getTree());
2109 TTreeReaderValue<Float_t> varToPlot(reader, plotinfo.variable.c_str());
2110 TTreeReaderValue<unsigned int> _entries(reader, "entries");
2111 TTreeReaderValue<unsigned int> _subDetId(reader, "subDetId");
2112 TTreeReaderValue<unsigned int> _moduleId(reader, "moduleId");
2113 TTreeReaderValue<Float_t> _zDirection(reader, "zDirection");
2114 TTreeReaderValue<Float_t> _rDirection(reader, "rDirection");
2115 TTreeReaderValue<unsigned int> _layer(reader, "layer");
2116 std::string badModulesFile_ = plotinfo.filterName;
2117 TFile* fBadModules = new TFile(badModulesFile_.c_str(), "READ");
2118 TTree* tBadModules = (TTree*)fBadModules->Get("alignTree");
2119 TTreeReader readerBad(tBadModules);
2120 TTreeReaderValue<int> _valid(readerBad, "valid");
2121 TTreeReaderValue<int> _bad_id(readerBad, "id");
2122 TTreeReaderValue<double> _bad_lumi(readerBad, "lumi");
2123
2124
2125 std::ofstream fUsedModules;
2126 std::ofstream fNotUsedModules;
2127 fUsedModules.open("usedModules.txt", std::ios::out | std::ios::app);
2128 fNotUsedModules.open("notUsedModules.txt", std::ios::out | std::ios::app);
2129
2130
2131 for (uint i = 0; i < plotinfo.vars->getTree()->GetEntries(); i++) {
2132 reader.SetEntry(i);
2133 if (*_entries < uint(plotinfo.minHits))
2134 continue;
2135 if (*_subDetId != uint(plotinfo.subDetId))
2136 continue;
2137 if (direction != 0) {
2138 if (plotinfo.subDetId == 2) {
2139 if (*_zDirection != direction)
2140 continue;
2141 } else {
2142 if (*_rDirection != direction)
2143 continue;
2144 }
2145 }
2146 if (layer > 0) {
2147 if (*_layer != uint(layer))
2148 continue;
2149 }
2150 bool isBadModule = false;
2151 for (int ibad = 0; ibad < tBadModules->GetEntries(); ibad++) {
2152 readerBad.SetEntry(ibad);
2153
2154 if (subdet == "BPIX" || subdet == "FPIX") {
2155 if (*_bad_lumi <= plotinfo.maxBadLumiPixel)
2156 continue;
2157 } else {
2158 if (*_bad_lumi <= plotinfo.maxBadLumiStrip)
2159 continue;
2160 }
2161
2162 if (*_moduleId == uint(*_bad_id))
2163 isBadModule = true;
2164 }
2165
2166 if (isBadModule) {
2167 fNotUsedModules << *_moduleId << "\n";
2168 continue;
2169 }
2170 fUsedModules << *_moduleId << "\n";
2171 if (h) {
2172 h->Fill(*varToPlot);
2173 }
2174 }
2175
2176
2177 fUsedModules.close();
2178 fNotUsedModules.close();
2179 fBadModules->Close();
2180 if (h) {
2181 h->SetName(histoname.Data());
2182 }
2183 }
2184
2185
2186 if (h && h->GetEntries() > 0) {
2187 if (direction == -1) {
2188 plotinfo.h1 = h;
2189 } else if (direction == 1) {
2190 plotinfo.h2 = h;
2191 } else {
2192 plotinfo.h = h;
2193 }
2194 }
2195 if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormX" ||
2196 plotinfo.variable == "rmsNormY") {
2197 storeHistogramInRootfile(h);
2198 }
2199 }
2200
2201 void PlotAlignmentValidation::modifySSHistAndLegend(THStack* hs, TLegend* legend) {
2202
2203
2204 Double_t legendY = 0.80;
2205 bool hasheader = (TkAlStyle::legendheader != "");
2206 if (hasheader)
2207 legend->SetHeader(TkAlStyle::legendheader);
2208 legend->SetFillStyle(0);
2209 int legendsize = hs->GetHists()->GetSize() + hasheader;
2210
2211 if (legendsize > 3)
2212 legendY -= 0.01 * (legendsize - 3);
2213 if (bigtext_) {
2214 legendY -= 0.05;
2215 }
2216 if (legendY < 0.6) {
2217 std::cerr << "Warning: Huge legend!" << std::endl;
2218 legendY = 0.6;
2219 }
2220 legend->SetY1(legendY);
2221 if (bigtext_)
2222 legend->SetTextSize(TkAlStyle::textSize);
2223
2224
2225 TProfile* prof = nullptr;
2226 TIter next(hs->GetHists());
2227 Int_t index = hasheader;
2228 while ((prof = (TProfile*)next())) {
2229
2230 Double_t scale = 10000;
2231 prof->Scale(scale);
2232
2233 Double_t stats[6] = {0};
2234 prof->GetStats(stats);
2235
2236 std::ostringstream legendtext;
2237 legendtext.precision(3);
2238 legendtext << std::fixed;
2239 legendtext << ": y mean = " << stats[4] / stats[0] * scale << " #mum";
2240
2241 TLegendEntry* entry = (TLegendEntry*)legend->GetListOfPrimitives()->At(index);
2242 if (entry == nullptr)
2243 std::cout << "PlotAlignmentValidation::PlotAlignmentValidation::modifySSLegend: Bad legend!" << std::endl;
2244 else
2245 entry->SetLabel((entry->GetLabel() + legendtext.str()).c_str());
2246 index++;
2247 }
2248
2249
2250 hs->SetMaximum(hs->GetMaximum("nostack PE") * 1.3);
2251 }
2252
2253
2254
2255
2256
2257
2258
2259 double PlotAlignmentValidation::resampleTestOfEqualRMS(TH1F* h1, TH1F* h2, int numSamples) {
2260
2261 std::vector<double> diff;
2262 diff.clear();
2263
2264 double rmsdiff = std::abs(h1->GetRMS() - h2->GetRMS());
2265
2266 double m1 = h1->GetMean();
2267 double m2 = h2->GetMean();
2268
2269 double d1 = 0;
2270 double d2 = 0;
2271
2272 double test_mean = 0;
2273 for (int i = 0; i < numSamples; i++) {
2274 d1 = 0;
2275 d2 = 0;
2276 for (int i = 0; i < h1->GetEntries(); i++) {
2277 d1 += h1->GetRandom() - m1;
2278 }
2279 for (int i = 0; i < h2->GetEntries(); i++) {
2280 d2 += h2->GetRandom() + m2;
2281 }
2282 d1 /= h1->GetEntries();
2283 d2 /= h2->GetEntries();
2284 diff.push_back(std::abs(d1 - d2 - rmsdiff));
2285 test_mean += std::abs(d1 - d2 - rmsdiff);
2286 }
2287 test_mean /= numSamples;
2288 edm::LogPrint("") << "test mean:" << test_mean;
2289
2290 double p = 0;
2291 for (double d : diff) {
2292 if (d > rmsdiff) {
2293 p += 1;
2294 }
2295 }
2296
2297 p /= numSamples;
2298 return p;
2299 }
2300
2301
2302
2303
2304
2305
2306
2307 double PlotAlignmentValidation::resampleTestOfEqualMeans(TH1F* h1, TH1F* h2, int numSamples) {
2308
2309 std::vector<double> diff;
2310 diff.clear();
2311
2312 double meandiff = std::abs(h1->GetMean() - h2->GetMean());
2313
2314 double d1 = 0;
2315 double d2 = 0;
2316
2317 double test_mean = 0;
2318 for (int i = 0; i < numSamples; i++) {
2319 d1 = 0;
2320 d2 = 0;
2321 for (int i = 0; i < h1->GetEntries(); i++) {
2322 d1 += h1->GetRandom();
2323 }
2324 for (int i = 0; i < h2->GetEntries(); i++) {
2325 d2 += h2->GetRandom();
2326 }
2327 d1 /= h1->GetEntries();
2328 d2 /= h2->GetEntries();
2329 diff.push_back(std::abs(d1 - d2 - meandiff));
2330 test_mean += std::abs(d1 - d2 - meandiff);
2331 }
2332 test_mean /= numSamples;
2333 edm::LogPrint("") << "test mean:" << test_mean;
2334
2335 double p = 0;
2336 for (double d : diff) {
2337 if (d > meandiff) {
2338 p += 1;
2339 }
2340 }
2341
2342 p /= numSamples;
2343 return p;
2344 }
2345
2346 float PlotAlignmentValidation::twotailedStudentTTestEqualMean(float t, float v) {
2347 return 2 * (1 - ROOT::Math::tdistribution_cdf(std::abs(t), v));
2348 }
2349
2350 const TString PlotAlignmentValidation::summaryfilename = "OfflineValidationSummary";
2351
2352 std::vector<TH1*> PlotAlignmentValidation::findmodule(TFile* f, unsigned int moduleid) {
2353
2354 TString histnamex;
2355 TString histnamey;
2356
2357 auto t = (TTree*)f->Get("TrackerOfflineValidation/TkOffVal");
2358
2359 TkOffTreeVariables* variables = nullptr;
2360 t->SetBranchAddress("TkOffTreeVariables", &variables);
2361 unsigned int number_of_entries = t->GetEntries();
2362 for (unsigned int i = 0; i < number_of_entries; i++) {
2363 t->GetEntry(i);
2364 if (variables->moduleId == moduleid) {
2365 histnamex = variables->histNameX;
2366 histnamey = variables->histNameY;
2367 break;
2368 }
2369 }
2370
2371 std::vector<TH1*> h;
2372
2373 auto h1 = (TH1*)f->FindObjectAny(histnamex);
2374 auto h2 = (TH1*)f->FindObjectAny(histnamey);
2375
2376 h1->SetDirectory(nullptr);
2377 h2->SetDirectory(nullptr);
2378
2379 h.push_back(h1);
2380 h.push_back(h2);
2381
2382 return h;
2383 }
2384
2385 void PlotAlignmentValidation::residual_by_moduleID(unsigned int moduleid) {
2386 TCanvas* cx = new TCanvas("x_residual");
2387 TCanvas* cy = new TCanvas("y_residual");
2388 TLegend* legendx = new TLegend(0.55, 0.7, 1, 0.9);
2389 TLegend* legendy = new TLegend(0.55, 0.7, 1, 0.9);
2390
2391 legendx->SetTextSize(0.016);
2392 legendx->SetTextAlign(12);
2393 legendy->SetTextSize(0.016);
2394 legendy->SetTextAlign(12);
2395
2396 for (auto it : sourceList) {
2397 TFile* file = it->getFile();
2398 int color = it->getLineColor();
2399 int linestyle = it->getLineStyle();
2400 TString legendname = it->getName();
2401 std::vector<TH1*> hist = findmodule(file, moduleid);
2402
2403 TString histnamex = legendname + " NEntries: " + std::to_string(int(hist[0]->GetEntries()));
2404 hist[0]->SetTitle(histnamex);
2405 hist[0]->SetStats(false);
2406 hist[0]->Rebin(50);
2407 hist[0]->SetBit(TH1::kNoTitle);
2408 hist[0]->SetLineColor(color);
2409 hist[0]->SetLineStyle(linestyle);
2410 cx->cd();
2411 hist[0]->Draw("Same");
2412 legendx->AddEntry(hist[0], histnamex, "l");
2413
2414 TString histnamey = legendname + " NEntries: " + std::to_string(int(hist[1]->GetEntries()));
2415 hist[1]->SetTitle(histnamey);
2416 hist[1]->SetStats(false);
2417 hist[1]->Rebin(50);
2418 hist[1]->SetBit(TH1::kNoTitle);
2419 hist[1]->SetLineColor(color);
2420 hist[1]->SetLineStyle(linestyle);
2421 cy->cd();
2422 hist[1]->Draw("Same");
2423 legendy->AddEntry(hist[1], histnamey, "l");
2424 }
2425
2426 TString filenamex = "x_residual_" + std::to_string(moduleid);
2427 TString filenamey = "y_residual_" + std::to_string(moduleid);
2428 cx->cd();
2429 legendx->Draw();
2430 cx->SaveAs(outputDir + "/" + filenamex + ".root");
2431 cx->SaveAs(outputDir + "/" + filenamex + ".pdf");
2432 cx->SaveAs(outputDir + "/" + filenamex + ".png");
2433 cx->SaveAs(outputDir + "/" + filenamex + ".eps");
2434
2435 cy->cd();
2436 legendy->Draw();
2437 cy->SaveAs(outputDir + "/" + filenamey + ".root");
2438 cy->SaveAs(outputDir + "/" + filenamey + ".pdf");
2439 cy->SaveAs(outputDir + "/" + filenamey + ".png");
2440 cy->SaveAs(outputDir + "/" + filenamey + ".eps");
2441 }