Project CMSSW displayed by LXR CMSSW_14_0_X_2023-12-08-2300/​Alignment/​OfflineValidation/​src/​PlotAlignmentValidation.cc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_14_0_X_2023-12-08-2300 CMSSW_14_0_X_2023-12-07-2300 CMSSW_14_0_X_2023-12-06-2300 CMSSW_14_0_X_2023-12-05-2300 CMSSW_14_0_X_2023-12-04-2300 CMSSW_14_0_X_2023-12-03-2300 Revert   Change

File indexing completed on 2023-11-18 03:07:10

 #include "Alignment/OfflineValidation/interface/PlotAlignmentValidation.h"
 
 /*! \class PlotAlignmentValidation
  *  \brief Class PlotAlignmentValidation
  *         Class used as the last step for Offline Track Validation tool.
  *         The main goal of this class is creating the plots regarding DMRs and Surface Deformations for modules and substructures.
  */
 
 //------------------------------------------------------------------------------
 /*! \fn PlotAlignmentValidation
  *  \brief Constructor for the class
  */
 
 PlotAlignmentValidation::PlotAlignmentValidation(bool bigtext) : bigtext_(bigtext) {
   setOutputDir(".");
   setTreeBaseDir();
   sourcelist = nullptr;
 
   moreThanOneSource = false;
   useFit_ = false;
 
   // Force ROOT to use scientific notation even with smaller datasets
   TGaxis::SetMaxDigits(4);
   // (This sets a static variable: correct in .eps images but must be set
   // again manually when viewing the .root files)
 
   // Make ROOT calculate histogram statistics using all data,
   // regardless of displayed range
   TH1::StatOverflows(kTRUE);
 
   //show all information in the legend by default
   legendOptions(TkAlStyle::legendoptions);
 }
 
 //------------------------------------------------------------------------------
 /*! \fn PlotAlignmentValidation
  *  \brief Constructor for the class. This function also retrieves the list of root files used to produce DMRs and Surface Deformations
  */
 PlotAlignmentValidation::PlotAlignmentValidation(
     const char* inputFile, std::string legendName, int lineColor, int lineStyle, bool bigtext)
     : PlotAlignmentValidation(bigtext) {
   loadFileList(inputFile, legendName, lineColor, lineStyle);
 }
 
 //------------------------------------------------------------------------------
 /*! \fn ~PlotAlignmentValidation
  *  \brief Default destructor
  */
 PlotAlignmentValidation::~PlotAlignmentValidation() {
   for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
     delete (*it);
   }
 
   delete sourcelist;
 }
 
 //------------------------------------------------------------------------------
 /*!
  * \fn openSummaryFile
  * \brief Create/open the root and txt summary files, where the DMR histograms and the associtated mean/sigma are stored respectively
  */
 
 void PlotAlignmentValidation::openSummaryFile() {
   if (!openedsummaryfile) {
     openedsummaryfile = true;
     summaryfile.open(outputDir + "/" + summaryfilename + ".txt");
     //Rootfile introduced to store the DMR histograms
     rootsummaryfile = new TFile(outputDir + "/" + summaryfilename + ".root", "RECREATE");
 
     for (auto vars : sourceList) {
       summaryfile << "\t" << vars->getName();
     }
     summaryfile << "\tformat={}\n";
   } else {
     //Check for the rootfile to be open, and open it in case it is not already.
     if (!rootsummaryfile->IsOpen())
       rootsummaryfile->Open(outputDir + "/" + summaryfilename + ".root", "UPDATE");
   }
 }
 
 //------------------------------------------------------------------------------
 /*! \fn loadFileList
  *  \brief Add to the list of sources the rootfile associated to a particular geometry
  */
 void PlotAlignmentValidation::loadFileList(const char* inputFile, std::string legendName, int lineColor, int lineStyle) {
   if (openedsummaryfile) {
     std::cout << "Can't load a root file after opening the summary file!" << std::endl;
     assert(0);
   }
   sourceList.push_back(new TkOfflineVariables(inputFile, treeBaseDir, legendName, lineColor, lineStyle));
 }
 
 //------------------------------------------------------------------------------
 /*! \fn useFitForDMRplots
  *  \brief Store the selected boolean in one of the private members of the class
  */
 void PlotAlignmentValidation::useFitForDMRplots(bool usefit) { useFit_ = usefit; }
 
 //------------------------------------------------------------------------------
 /*! \fn numberOfLayers
  *  \brief Select the number of layers associated to a subdetector.
  */
 //TODO Possible improvement: reduce the number of switches in the code by implementing a map
 int PlotAlignmentValidation::numberOfLayers(int phase, int subdetector) {
   switch (phase) {
     case 0:
       switch (subdetector) {
         case 1:
           return 3;
         case 2:
           return 2;
         case 3:
           return 4;
         case 4:
           return 3;
         case 5:
           return 6;
         case 6:
           return 9;
         default:
           assert(false);
       }
     case 1:
       switch (subdetector) {
         case 1:
           return 4;
         case 2:
           return 3;
         case 3:
           return 4;
         case 4:
           return 3;
         case 5:
           return 6;
         case 6:
           return 9;
         default:
           assert(false);
       }
     default:
       assert(false);
   }
   return 0;
 }
 
 //------------------------------------------------------------------------------
 /*! \fn maxNumberOfLayers
  *  \brief Return the number of layers of a subdetector
  */
 int PlotAlignmentValidation::maxNumberOfLayers(int subdetector) {
   int result = 0;
   for (auto it = sourceList.begin(); it != sourceList.end(); ++it) {
     result = std::max(result, numberOfLayers((*it)->getPhase(), subdetector));
   }
   return result;
 }
 
 //------------------------------------------------------------------------------
 /*! \fn legendOptions
  *  \brief Assign legend options to members of the class
  */
 void PlotAlignmentValidation::legendOptions(TString options) {
   showMean_ = false;
   showRMS_ = false;
   showMeanError_ = false;
   showRMSError_ = false;
   showModules_ = false;
   showUnderOverFlow_ = false;
   options.ReplaceAll(" ", "").ToLower();
   if (options.Contains("mean") || options.Contains("all"))
     showMean_ = true;
   if (options.Contains("meanerror") || options.Contains("all"))
     showMeanError_ = true;
   if (options.Contains("rms") || options.Contains("all"))
     showRMS_ = true;
   if (options.Contains("rmserror") || options.Contains("all"))
     showRMSError_ = true;
   if (options.Contains("modules") || options.Contains("all"))
     showModules_ = true;
   if (options.Contains("under") || options.Contains("over") || options.Contains("outside") || options.Contains("all"))
     showUnderOverFlow_ = true;
 
   twolines_ = (showUnderOverFlow_ && (showMean_ + showMeanError_ + showRMS_ + showRMSError_ >= 1) && bigtext_);
 }
 
 //------------------------------------------------------------------------------
 /*! \fn setOutputDir
  *  \brief Set the output direcotry
  */
 void PlotAlignmentValidation::setOutputDir(std::string dir) {
   if (openedsummaryfile) {
     std::cout << "Can't set the output dir after opening the summary file!" << std::endl;
     assert(0);
   }
   outputDir = dir;
   gSystem->mkdir(outputDir.data(), true);
 }
 
 //------------------------------------------------------------------------------
 /*! \fn plotSubDetResiduals
  *  \brief Function used to plot residuals for a subdetector
  */
 void PlotAlignmentValidation::plotSubDetResiduals(bool plotNormHisto, unsigned int subDetId) {
   gStyle->SetOptStat(11111);
   gStyle->SetOptFit(0000);
 
   TCanvas* c = new TCanvas("c", "c");
   c->SetTopMargin(0.15);
   TString histoName = "";
   if (plotNormHisto) {
     histoName = "h_NormXprime";
   } else
     histoName = "h_Xprime_";
   switch (subDetId) {
     case 1:
       histoName += "TPBBarrel_0";
       break;
     case 2:
       histoName += "TPEendcap_1";
       break;
     case 3:
       histoName += "TPEendcap_2";
       break;
     case 4:
       histoName += "TIBBarrel_0";
       break;
     case 5:
       histoName += "TIDEndcap_1";
       break;
     case 6:
       histoName += "TIDEndcap_2";
       break;
     case 7:
       histoName += "TOBBarrel_3";
       break;
     case 8:
       histoName += "TECEndcap_4";
       break;
     case 9:
       histoName += "TECEndcap_5";
       break;
   }
   int tmpcounter = 0;
   TH1* sumHisto = nullptr;
   for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
     if (tmpcounter == 0) {
       TFile* f = (*it)->getFile();
       sumHisto = (TH1*)f->FindKeyAny(histoName)->ReadObj();  //FindObjectAny(histoName.Data());
       sumHisto->SetLineColor(tmpcounter + 1);
       sumHisto->SetLineStyle(tmpcounter + 1);
       sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
       sumHisto->Draw();
 
       //get statistic box coordinate to plot all boxes one below the other
       //gStyle->SetStatY(0.91);
       //gStyle->SetStatW(0.15);
       //gStyle->SetStatBorderSize(1);
       //gStyle->SetStatH(0.10);
 
       tmpcounter++;
     } else {
       sumHisto = (TH1*)(*it)->getFile()->FindObjectAny(histoName);
       sumHisto->SetLineColor(tmpcounter + 1);
       sumHisto->SetLineStyle(tmpcounter + 1);
       sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
       //hstack->Add(sumHisto);
 
       c->Update();
       tmpcounter++;
     }
     TObject* statObj = sumHisto->GetListOfFunctions()->FindObject("stats");
     if (statObj && statObj->InheritsFrom(TPaveStats::Class())) {
       TPaveStats* stats = static_cast<TPaveStats*>(statObj);
       stats->SetLineColor(tmpcounter + 1);
       stats->SetTextColor(tmpcounter + 1);
       stats->SetFillColor(10);
       stats->SetX1NDC(0.91 - tmpcounter * 0.1);
       stats->SetX2NDC(0.15);
       stats->SetY1NDC(1);
       stats->SetY2NDC(0.10);
       sumHisto->Draw("sames");
     }
   }
   //hstack->Draw("nostack");
   char PlotName[1000];
   sprintf(PlotName, "%s/%s.png", outputDir.c_str(), histoName.Data());
   c->Print(PlotName);
   sprintf(PlotName, "%s/%s.eps", outputDir.c_str(), histoName.Data());
   c->Print(PlotName);
   sprintf(PlotName, "%s/%s.pdf", outputDir.c_str(), histoName.Data());
   c->Print(PlotName);
   sprintf(PlotName, "%s/%s.root", outputDir.c_str(), histoName.Data());
   c->Print(PlotName);
   //delete c;
   //c=0;
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::plotHitMaps() {
   //gStyle->SetOptStat(0);
 
   TCanvas* c = new TCanvas("c", "c", 1200, 400);
   c->Divide(3, 1);
   //ps->NewPage();
 
   //-------------------------------------------------
   //plot Hit map
   //-------------------------------------------------
   std::string histName_ = "Entriesprofile";
   c->cd(1);
   TTree* tree = (*sourceList.begin())->getTree();
   tree->Draw("entries:posR:posZ", "", "COLZ2Prof");
   c->cd(2);
   tree->Draw("entries:posY:posX", "", "COLZ2Prof");
   c->cd(3);
   tree->Draw("entries:posR:posPhi", "", "COLZ2Prof");
 
   char PlotName[1000];
   sprintf(PlotName, "%s/%s.png", outputDir.c_str(), histName_.c_str());
   c->Print(PlotName);
   sprintf(PlotName, "%s/%s.eps", outputDir.c_str(), histName_.c_str());
   c->Print(PlotName);
   sprintf(PlotName, "%s/%s.pdf", outputDir.c_str(), histName_.c_str());
   c->Print(PlotName);
   sprintf(PlotName, "%s/%s.root", outputDir.c_str(), histName_.c_str());
   c->Print(PlotName);
   //   //c->Update();
   c->Close();
   //----------------------------------------------------
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::plotOutlierModules(const char* outputFileName,
                                                  std::string plotVariable,
                                                  float plotVariable_cut,
                                                  int unsigned minHits) {
   Int_t counter = 0;
 
   gStyle->SetOptStat(111111);
   gStyle->SetStatY(0.9);
   //TList* treelist=getTreeList();
 
   TCanvas* c1 = new TCanvas("canv", "canv", 800, 500);
   outputFile = outputDir + '/' + outputFileName;
   c1->Print((outputFile + '[').Data());
 
   c1->Divide(2, 1);
 
   TTree* tree = (*sourceList.begin())->getTree();
   TkOffTreeVariables* treeMem = nullptr;  // ROOT will initilise
   tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
 
   Long64_t nentries = tree->GetEntriesFast();
 
   for (Long64_t i = 0; i < nentries; i++) {
     tree->GetEntry(i);
     float var = 0;
     if (plotVariable == "chi2PerDofX")
       var = treeMem->chi2PerDofX;
     else if (plotVariable == "chi2PerDofY")
       var = treeMem->chi2PerDofY;
     else if (plotVariable == "fitMeanX")
       var = treeMem->fitMeanX;
     else if (plotVariable == "fitMeanY")
       var = treeMem->fitMeanY;
     else if (plotVariable == "fitSigmaX")
       var = treeMem->fitSigmaX;
     else if (plotVariable == "fitSigmaY")
       var = treeMem->fitSigmaY;
     else {
       std::cout << "There is no variable " << plotVariable << " included in the tree." << std::endl;
       break;
     }
     //   std::cout<<"treeMem->entries  "<<treeMem->entries<<std::endl;
     //  std::cout<<"var                  "<<var<<std::endl;
     //  std::cout<<"plotVariable_cut     "<<plotVariable_cut<<std::endl;
 
     if (var > plotVariable_cut && treeMem->entries > minHits) {
       TFile* f = (*sourceList.begin())->getFile();  //(TFile*)sourcelist->First();
 
       if (f->FindKeyAny(treeMem->histNameX.c_str()) != nullptr) {
         TH1* h =
             (TH1*)f->FindKeyAny(treeMem->histNameX.c_str())->ReadObj();  //f->FindObjectAny(treeMem->histNameX.c_str());
         gStyle->SetOptFit(0111);
         std::cout << "hist name " << h->GetName() << std::endl;
 
         TString path = (char*)strstr(gDirectory->GetPath(), "TrackerOfflineValidation");
         //std::cout<<"hist path "<<path<<std::endl;
         //std::cout<<"wrote text "<<std::endl;
         if (h)
           std::cout << h->GetEntries() << std::endl;
 
         //modules' location as title
         c1->cd(0);
         TPaveText* text = new TPaveText(0, 0.95, 0.99, 0.99);
         text->AddText(path);
         text->SetFillColor(0);
         text->SetShadowColor(0);
         text->SetBorderSize(0);
         text->Draw();
 
         //residual histogram
         c1->cd(1);
         TPad* subpad = (TPad*)c1->GetPad(1);
         subpad->SetPad(0, 0, 0.5, 0.94);
         h->Draw();
 
         //norm. residual histogram
         h = (TH1*)f->FindObjectAny(treeMem->histNameNormX.c_str());
         if (h)
           std::cout << h->GetEntries() << std::endl;
         c1->cd(2);
         TPad* subpad2 = (TPad*)c1->GetPad(2);
         subpad2->SetPad(0.5, 0, 0.99, 0.94);
         h->Draw();
 
         c1->Print(outputFile);
         counter++;
       } else {
         std::cout << "There are no residual histograms on module level stored!" << std::endl;
         std::cout << "Please make sure that moduleLevelHistsTransient = cms.bool(False) in the validation job!"
                   << std::endl;
         break;
       }
     }
   }
   c1->Print((outputFile + "]").Data());
   if (counter == 0)
     std::cout << "no bad modules found" << std::endl;
 
   //read the number of entries in the t3
   //TTree* tree=0;
   //tree=(TTree*)treeList->At(0);
 
   //c1->Close();
 }
 
 //------------------------------------------------------------------------------
 /*! \fn getTreeList
  *  \brief Extract from the rootfiles stored in the sourcelist the corresponding trees.
  */
 TList* PlotAlignmentValidation::getTreeList() {
   TList* treeList = new TList();
   TFile* first_source = (TFile*)sourcelist->First();
   std::cout << first_source->GetName() << std::endl;
   TDirectoryFile* d = (TDirectoryFile*)first_source->Get(treeBaseDir.c_str());
   treeList->Add((TTree*)(*d).Get("TkOffVal"));
 
   if (moreThanOneSource == true) {
     TFile* nextsource = (TFile*)sourcelist->After(first_source);
     while (nextsource) {
       std::cout << nextsource->GetName() << std::endl;
       d = (TDirectoryFile*)nextsource->Get("TrackerOfflineValidation");
 
       treeList->Add((TTree*)(*d).Get("TkOffVal"));
 
       nextsource = (TFile*)sourcelist->After(nextsource);
     }
   }
   return treeList;
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::setTreeBaseDir(std::string dir) { treeBaseDir = dir; }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::plotSurfaceShapes(const std::string& options, const std::string& residType) {
   std::cout << "-------- plotSurfaceShapes called with " << options << std::endl;
   if (options == "none")
     return;
   else if (options == "coarse") {
     plotSS("subdet=1");
     plotSS("subdet=2");
     plotSS("subdet=3");
     plotSS("subdet=4");
     plotSS("subdet=5");
     plotSS("subdet=6");
   }
   // else if (options == "fine") ...
   else
     plotSS(options, residType);
 
   return;
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::plotSS(const std::string& options, const std::string& residType) {
   if (residType.empty()) {
     plotSS(options, "ResXvsXProfile");
     plotSS(options, "ResXvsYProfile");
     return;
   }
 
   int bkperrorx = gStyle->GetErrorX();
   gStyle->SetErrorX(1);  //regardless of style settings, we want x error bars here
 
   int plotLayerN = 0;
   //  int plotRingN  = 0;
   //  bool plotPlain = false;
   bool plotLayers = false;  // overrides plotLayerN
   //  bool plotRings  = false;  // Todo: implement this?
   int plotSubDetN = 0;  // if zero, plot all
 
   TRegexp layer_re("layer=[0-9]+");
   Ssiz_t index, len;
   if (options.find("layers") != std::string::npos) {
     plotLayers = true;
   }
   if ((index = layer_re.Index(options, &len)) != -1) {
     if (plotLayers) {
       std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
     } else {
       std::string substr = options.substr(index + 6, len - 6);
       plotLayerN = atoi(substr.c_str());
     }
   }
 
   TRegexp subdet_re("subdet=[1-6]+");
   if ((index = subdet_re.Index(options, &len)) != -1) {
     std::string substr = options.substr(index + 7, len - 7);
     plotSubDetN = atoi(substr.c_str());
   }
 
   gStyle->SetOptStat(0);
 
   TCanvas c("canv", "canv");
 
   // todo: title, min/max, nbins?
 
   // Loop over detectors
   for (int iSubDet = 1; iSubDet <= 6; ++iSubDet) {
     // TEC requires special care since rings 1-4 and 5-7 are plotted separately
     bool isTEC = (iSubDet == 6);
 
     // if subdet is specified, skip other subdets
     if (plotSubDetN != 0 && iSubDet != plotSubDetN)
       continue;
 
     // Skips plotting too high layers
     // if it's a mixture of phase 0 and 1, the phase 0 files will be skipped
     //  when plotting the higher layers of BPIX and FPIX
     if (plotLayerN > maxNumberOfLayers(iSubDet)) {
       continue;
     }
 
     int minlayer = plotLayers ? 1 : plotLayerN;
     int maxlayer = plotLayers ? maxNumberOfLayers(iSubDet) : plotLayerN;
     // see later where this is used
     int maxlayerphase0 = plotLayers ? numberOfLayers(0, iSubDet) : plotLayerN;
 
     for (int layer = minlayer; layer <= maxlayer; layer++) {
       // two plots for TEC, skip first
       for (int iTEC = 0; iTEC < 2; iTEC++) {
         if (!isTEC && iTEC == 0)
           continue;
 
         char selection[1000];
         if (!isTEC) {
           if (layer == 0)
             sprintf(selection, "subDetId==%d", iSubDet);
           else
             sprintf(selection, "subDetId==%d && layer == %d", iSubDet, layer);
         } else {          // TEC
           if (iTEC == 0)  // rings
             sprintf(selection, "subDetId==%d && ring <= 4", iSubDet);
           else
             sprintf(selection, "subDetId==%d && ring > 4", iSubDet);
         }
 
         // Title for plot and name for the file
 
         TString subDetName;
         switch (iSubDet) {
           case 1:
             subDetName = "BPIX";
             break;
           case 2:
             subDetName = "FPIX";
             break;
           case 3:
             subDetName = "TIB";
             break;
           case 4:
             subDetName = "TID";
             break;
           case 5:
             subDetName = "TOB";
             break;
           case 6:
             subDetName = "TEC";
             break;
         }
 
         TString secondline = "";
         if (layer != 0) {
           // TEC and TID have discs, the rest have layers
           if (iSubDet == 4 || iSubDet == 6)
             secondline = "disc ";
           else {
             secondline = "layer ";
           }
           secondline += Form("%d", layer);
           secondline += " ";
         }
         if (isTEC && iTEC == 0)
           secondline += TString("R1-4");
         if (isTEC && iTEC > 0)
           secondline += TString("R5-7");
 
         // Generate histograms with selection
         TLegend* legend = nullptr;
         // Any file from phase 0 will be skipped if the last argument is false
         THStack* hs = addHists(selection, residType, &legend, false, /*validforphase0 = */ layer <= maxlayerphase0);
         if (!hs || hs->GetHists() == nullptr || hs->GetHists()->GetSize() == 0) {
           std::cout << "No histogram for " << subDetName << ", perhaps not enough data? Creating default histogram."
                     << std::endl;
           if (hs == nullptr)
             hs = new THStack("hstack", "");
 
           TProfile* defhist = new TProfile("defhist", "Empty default histogram", 100, -1, 1, -1, 1);
           hs->Add(defhist);
           hs->Draw();
         } else {
           hs->Draw("nostack PE");
           modifySSHistAndLegend(hs, legend);
           legend->Draw();
           setTitleStyle(*hs, "", "", iSubDet, true, secondline);
 
           // Adjust Labels
           TH1* firstHisto = (TH1*)hs->GetHists()->First();
           TString xName = firstHisto->GetXaxis()->GetTitle();
           TString yName = firstHisto->GetYaxis()->GetTitle();
           hs->GetHistogram()->GetXaxis()->SetTitleColor(kBlack);
           hs->GetHistogram()->GetXaxis()->SetTitle(xName);
           hs->GetHistogram()->GetYaxis()->SetTitleColor(kBlack);
           // micrometers:
           yName.ReplaceAll("cm", "#mum");
           hs->GetHistogram()->GetYaxis()->SetTitle(yName);
         }
 
         // Save plot to file
         std::ostringstream plotName;
         plotName << outputDir << "/SurfaceShape_" << subDetName << "_";
         plotName << residType;
         if (layer != 0) {
           plotName << "_";
           // TEC and TID have discs, the rest have layers
           if (iSubDet == 4 || iSubDet == 6)
             plotName << "disc";
           else {
             plotName << "layer";
           }
           plotName << layer;
         }
         if (isTEC && iTEC == 0)
           plotName << "_"
                    << "R1-4";
         if (isTEC && iTEC > 0)
           plotName << "_"
                    << "R5-7";
 
         // PNG,EPS,PDF files
         c.Update();
         c.Print((plotName.str() + ".png").c_str());
         c.Print((plotName.str() + ".eps").c_str());
         c.Print((plotName.str() + ".pdf").c_str());
 
         // ROOT file
         TFile f((plotName.str() + ".root").c_str(), "recreate");
         c.Write();
         f.Close();
 
         delete legend;
         delete hs;
       }
     }
   }
   gStyle->SetErrorX(bkperrorx);
 
   return;
 }
 
 //------------------------------------------------------------------------------
 /*! \fn plotDMR
  *  \brief Main function used to plot DMRs for a single IOV printing the canvases in the output directory and saving histograms and fit funtions in a root file.
  */
 void PlotAlignmentValidation::plotDMR(const std::string& variable,
                                       Int_t minHits,
                                       const std::string& options,
                                       const std::string& filterName) {
   // If several, comma-separated values are given in 'variable',
   // call plotDMR with each value separately.
   // If a comma is found, the string is divided to two.
   // (no space allowed)
   std::size_t findres = variable.find(',');
   if (findres != std::string::npos) {
     std::string substring1 = variable.substr(0, findres);
     std::string substring2 = variable.substr(findres + 1, std::string::npos);
     plotDMR(substring1, minHits, options, filterName);
     plotDMR(substring2, minHits, options, filterName);
     return;
   }
 
   // Variable name should end with X or Y. If it doesn't, recursively calls plotDMR twice with
   // X and Y added, respectively
   if (variable == "mean" || variable == "median" || variable == "meanNorm" || variable == "rms" ||
       variable == "rmsNorm") {
     plotDMR(variable + "X", minHits, options, filterName);
     plotDMR(variable + "Y", minHits, options, filterName);
     return;
   }
 
   // options:
   // -plain (default, the whole distribution)
   // -split (distribution splitted to two)
   // -layers (plain db for each layer/disc superimposed in one plot)
   // -layersSeparate (plain db for each layer/disc in separate plots)
   // -layersSplit (splitted db for each layers/disc in one plot)
   // -layersSplitSeparate (splitted db, for each layers/disc in separate plots)
 
   TRegexp layer_re("layer=[0-9]+");
   bool plotPlain = false, plotSplits = false, plotLayers = false;
   int plotLayerN = 0;
   Ssiz_t index, len;
   if (options.find("plain") != std::string::npos) {
     plotPlain = true;
   }
   if (options.find("split") != std::string::npos) {
     plotSplits = true;
   }
   if (options.find("layers") != std::string::npos) {
     plotLayers = true;
   }
   if ((index = layer_re.Index(options, &len)) != -1) {
     if (plotLayers) {
       std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
     } else {
       std::string substr = options.substr(index + 6, len - 6);
       plotLayerN = atoi(substr.c_str());
     }
   }
 
   // Defaults to plotting only plain plot if empty (or invalid)
   // option string is given
   if (!plotPlain && !plotSplits) {
     plotPlain = true;
   }
 
   // This boolean array tells for which detector modules to plot split DMR plots
   // They are plotted for BPIX, FPIX, TIB and TOB
   static bool plotSplitsFor[6] = {true, true, true, false, true, false};
 
   DMRPlotInfo plotinfo;
 
   gStyle->SetOptStat(0);
 
   TCanvas c("canv", "canv");
 
   plotinfo.variable = variable;
   plotinfo.minHits = minHits;
   plotinfo.plotPlain = plotPlain;
   plotinfo.plotLayers = plotLayers;
   plotinfo.filterName = filterName;
 
   // width in cm
   // for DMRS, use 100 bins in range +-10 um, bin width 0.2um
   // if modified, check also TrackerOfflineValidationSummary_cfi.py and TrackerOfflineValidation_Standalone_cff.py
   if (variable == "meanX") {
     plotinfo.nbins = 50;
     plotinfo.min = -0.001;
     plotinfo.max = 0.001;
   } else if (variable == "meanY") {
     plotinfo.nbins = 50;
     plotinfo.min = -0.005;
     plotinfo.max = 0.005;
   } else if (variable == "medianX")
     if (plotSplits) {
       plotinfo.nbins = 50;
       plotinfo.min = -0.0005;
       plotinfo.max = 0.0005;
     } else {
       plotinfo.nbins = 100;
       plotinfo.min = -0.001;
       plotinfo.max = 0.001;
     }
   else if (variable == "medianY")
     if (plotSplits) {
       plotinfo.nbins = 50;
       plotinfo.min = -0.0005;
       plotinfo.max = 0.0005;
     } else {
       plotinfo.nbins = 100;
       plotinfo.min = -0.001;
       plotinfo.max = 0.001;
     }
   else if (variable == "meanNormX") {
     plotinfo.nbins = 100;
     plotinfo.min = -2.0;
     plotinfo.max = 2.0;
   } else if (variable == "meanNormY") {
     plotinfo.nbins = 100;
     plotinfo.min = -2.0;
     plotinfo.max = 2.0;
   } else if (variable == "rmsX") {
     plotinfo.nbins = 100;
     plotinfo.min = 0.0;
     plotinfo.max = 0.1;
   } else if (variable == "rmsY") {
     plotinfo.nbins = 100;
     plotinfo.min = 0.0;
     plotinfo.max = 0.1;
   } else if (variable == "rmsNormX") {
     plotinfo.nbins = 100;
     plotinfo.min = 0.3;
     plotinfo.max = 1.8;
   } else if (variable == "rmsNormY") {
     plotinfo.nbins = 100;
     plotinfo.min = 0.3;
     plotinfo.max = 1.8;
   } else {
     std::cerr << "Unknown variable " << variable << std::endl;
     plotinfo.nbins = 100;
     plotinfo.min = -0.1;
     plotinfo.max = 0.1;
   }
   //Begin loop on structures
   for (int i = 1; i <= 6; ++i) {
     // Skip strip detectors if plotting any "Y" variable
     if (i != 1 && i != 2 && variable.length() > 0 && variable[variable.length() - 1] == 'Y') {
       continue;
     }
 
     // Skips plotting too high layers
     if (plotLayerN > maxNumberOfLayers(i)) {
       continue;
     }
 
     plotinfo.plotSplits = plotSplits && plotSplitsFor[i - 1];
     if (!plotinfo.plotPlain && !plotinfo.plotSplits) {
       continue;
     }
 
     // Sets dimension of legend according to the number of plots
 
     bool hasheader = (TkAlStyle::legendheader != "");
 
     int nPlots = 1;
     if (plotinfo.plotSplits) {
       nPlots = 3;
     }
     // This will make the legend a bit bigger than necessary if there is a mixture of phase 0 and phase 1.
     // Not worth it to implement more complicated logic.
     if (plotinfo.plotLayers) {
       nPlots *= maxNumberOfLayers(i);
     }
     nPlots *= sourceList.size();
     if (twolines_) {
       nPlots *= 2;
     }
     nPlots += hasheader;
 
     double legendY = 0.80;
     if (nPlots > 3) {
       legendY -= 0.01 * (nPlots - 3);
     }
     if (bigtext_) {
       legendY -= 0.05;
     }
     if (legendY < 0.6) {
       std::cerr << "Warning: Huge legend!" << std::endl;
       legendY = 0.6;
     }
 
     THStack hstack("hstack", "hstack");
     plotinfo.maxY = 0;
     plotinfo.subDetId = i;
     plotinfo.legend = new TLegend(0.17, legendY, 0.85, 0.88);
     plotinfo.legend->SetNColumns(2);
     if (hasheader)
       plotinfo.legend->SetHeader(TkAlStyle::legendheader);
     if (bigtext_)
       plotinfo.legend->SetTextSize(TkAlStyle::textSize);
     plotinfo.legend->SetFillStyle(0);
     plotinfo.hstack = &hstack;
     plotinfo.h = plotinfo.h1 = plotinfo.h2 = nullptr;
     plotinfo.firsthisto = true;
 
     openSummaryFile();
     vmean.clear();
     vrms.clear();
     vdeltamean.clear();
     vmeanerror.clear();
     vPValueEqualSplitMeans.clear(), vAlignmentUncertainty.clear();
     vPValueMeanEqualIdeal.clear();
     vPValueRMSEqualIdeal.clear();
 
     std::string stringsubdet;
     switch (i) {
       case 1:
         stringsubdet = "BPIX";
         break;
       case 2:
         stringsubdet = "FPIX";
         break;
       case 3:
         stringsubdet = "TIB";
         break;
       case 4:
         stringsubdet = "TID";
         break;
       case 5:
         stringsubdet = "TOB";
         break;
       case 6:
         stringsubdet = "TEC";
         break;
     }
 
     for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
       plotinfo.vars = *it;
       plotinfo.h1 = plotinfo.h2 = plotinfo.h = nullptr;
 
       int minlayer = plotLayers ? 1 : plotLayerN;
       //Layer 0 is associated to the entire structure, this check ensures that even when both the plotLayers and the plotPlain options are active, also the histogram for the entire structure is made.
       if (plotinfo.plotPlain)
         minlayer = 0;
       int maxlayer = plotLayers ? numberOfLayers(plotinfo.vars->getPhase(), plotinfo.subDetId) : plotLayerN;
 
       for (int layer = minlayer; layer <= maxlayer; layer++) {
         if (plotinfo.plotPlain) {
           plotDMRHistogram(plotinfo, 0, layer, stringsubdet);
         }
 
         if (plotinfo.plotSplits) {
           plotDMRHistogram(plotinfo, -1, layer, stringsubdet);
           plotDMRHistogram(plotinfo, 1, layer, stringsubdet);
         }
 
         if (plotinfo.plotPlain) {
           if (plotinfo.h) {
             setDMRHistStyleAndLegend(plotinfo.h, plotinfo, 0, layer);
           } else {
             if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") &&
                 /*!plotinfo.plotLayers && */ layer == 0) {
               vmean.push_back(nan(""));
               vrms.push_back(nan(""));
               vmeanerror.push_back(nan(""));
               vAlignmentUncertainty.push_back(nan(""));
               vPValueMeanEqualIdeal.push_back(nan(""));
               vPValueRMSEqualIdeal.push_back(nan(""));
               if (plotinfo.plotSplits && plotinfo.plotPlain) {
                 vdeltamean.push_back(nan(""));
                 vPValueEqualSplitMeans.push_back(nan(""));
               }
             }
           }
         }
 
         if (plotinfo.plotSplits) {
           // Add delta mu to the histogram
           if (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && !plotinfo.plotPlain) {
             std::ostringstream legend;
             std::string unit = " #mum";
             legend.precision(3);
             legend << std::fixed;  // to always show 3 decimals
             float factor = 10000.0f;
             if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" ||
                 plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY") {
               factor = 1.0f;
               unit = "";
             }
             float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
             legend << plotinfo.vars->getName();
             if (layer > 0) {
               // TEC and TID have discs, the rest have layers
               if (i == 4 || i == 6)
                 legend << ", disc ";
               else
                 legend << ", layer ";
               legend << layer;
             }
             plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
             legend.str("");
             legend << "#Delta#mu = " << deltamu << unit;
             plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
 
             if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && !plotLayers && layer == 0) {
               vdeltamean.push_back(deltamu);
             }
           }
           if (plotinfo.h1) {
             setDMRHistStyleAndLegend(plotinfo.h1, plotinfo, -1, layer);
           }
           if (plotinfo.h2) {
             setDMRHistStyleAndLegend(plotinfo.h2, plotinfo, 1, layer);
           }
         }
       }
     }
 
     if (hstack.GetHists() != nullptr && hstack.GetHists()->GetSize() != 0) {
       hstack.Draw("nostack");
       hstack.SetMaximum(plotinfo.maxY * 1.3);
       setTitleStyle(hstack, variable.c_str(), "#modules", plotinfo.subDetId);
       setHistStyle(*hstack.GetHistogram(), variable.c_str(), "#modules", 1);
 
       plotinfo.legend->Draw();
     } else {
       // Draw an empty default histogram
       plotinfo.h = new TH1F("defhist", "Empty default histogram", plotinfo.nbins, plotinfo.min, plotinfo.max);
       plotinfo.h->SetMaximum(10);
       if (plotinfo.variable.find("Norm") == std::string::npos)
         scaleXaxis(plotinfo.h, 10000);
       setTitleStyle(*plotinfo.h, variable.c_str(), "#modules", plotinfo.subDetId);
       setHistStyle(*plotinfo.h, variable.c_str(), "#modules", 1);
       plotinfo.h->Draw();
     }
 
     std::ostringstream plotName;
     plotName << outputDir << "/D";
 
     if (variable == "medianX")
       plotName << "medianR_";
     else if (variable == "medianY")
       plotName << "medianYR_";
     else if (variable == "meanX")
       plotName << "meanR_";
     else if (variable == "meanY")
       plotName << "meanYR_";
     else if (variable == "meanNormX")
       plotName << "meanNR_";
     else if (variable == "meanNormY")
       plotName << "meanNYR_";
     else if (variable == "rmsX")
       plotName << "rmsR_";
     else if (variable == "rmsY")
       plotName << "rmsYR_";
     else if (variable == "rmsNormX")
       plotName << "rmsNR_";
     else if (variable == "rmsNormY")
       plotName << "rmsNYR_";
 
     TString subdet;
     switch (i) {
       case 1:
         subdet = "BPIX";
         break;
       case 2:
         subdet = "FPIX";
         break;
       case 3:
         subdet = "TIB";
         break;
       case 4:
         subdet = "TID";
         break;
       case 5:
         subdet = "TOB";
         break;
       case 6:
         subdet = "TEC";
         break;
     }
 
     plotName << subdet;
 
     if (plotPlain && !plotSplits) {
       plotName << "_plain";
     } else if (!plotPlain && plotSplits) {
       plotName << "_split";
     }
     if (plotLayers) {
       // TEC and TID have discs, the rest have layers
       if (i == 4 || i == 6)
         plotName << "_discs";
       else
         plotName << "_layers";
     }
     if (plotLayerN > 0) {
       // TEC and TID have discs, the rest have layers
       if (i == 4 || i == 6)
         plotName << "_disc";
       else
         plotName << "_layer";
       plotName << plotLayerN;
     }
 
     // PNG,EPS,PDF files
     c.Update();
     c.Print((plotName.str() + ".png").c_str());
     c.Print((plotName.str() + ".eps").c_str());
     c.Print((plotName.str() + ".pdf").c_str());
 
     // ROOT file
     TFile f((plotName.str() + ".root").c_str(), "recreate");
     c.Write();
     f.Close();
 
     // Free allocated memory.
     delete plotinfo.h;
     delete plotinfo.h1;
     delete plotinfo.h2;
 
     if (!vmean.empty()) {
       summaryfile << "   mu_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << " (um)\t"
                   << "latexname=$\\mu_\\text{" << subdet << "}";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$ ($\\mu$m)\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto mu : vmean)
         summaryfile << "\t" << mu;
       summaryfile << "\n";
     }
     if (!vrms.empty()) {
       summaryfile << "sigma_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << " (um)\t"
                   << "latexname=$\\sigma_\\text{" << subdet << "}";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$ ($\\mu$m)\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto sigma : vrms)
         summaryfile << "\t" << sigma;
       summaryfile << "\n";
     }
     if (!vdeltamean.empty()) {
       summaryfile << "  dmu_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << " (um)\t"
                   << "latexname=$\\Delta\\mu_\\text{" << subdet << "}";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$ ($\\mu$m)\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto dmu : vdeltamean)
         summaryfile << "\t" << dmu;
       summaryfile << "\n";
     }
     if (!vmeanerror.empty()) {
       summaryfile << "  sigma_mu_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << " (um)\t"
                   << "latexname=$\\sigma\\mu_\\text{" << subdet << "}";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$ ($\\mu$m)\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto dmu : vmeanerror)
         summaryfile << "\t" << dmu;
       summaryfile << "\n";
     }
     if (!vPValueEqualSplitMeans.empty()) {
       summaryfile << "  p_delta_mu_equal_zero_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << "\t"
                   << "latexname=$P(\\delta\\mu_\\text{" << subdet << "}=0)";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto dmu : vPValueEqualSplitMeans)
         summaryfile << "\t" << dmu;
       summaryfile << "\n";
     }
     if (!vAlignmentUncertainty.empty()) {
       summaryfile << "  alignment_uncertainty_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << " (um)\t"
                   << "latexname=$\\sigma_\\text{align}_\\text{" << subdet << "}";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$ ($\\mu$m)\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto dmu : vAlignmentUncertainty)
         summaryfile << "\t" << dmu;
       summaryfile << "\n";
     }
     if (!vPValueMeanEqualIdeal.empty()) {
       summaryfile << "  p_mean_equals_ideal_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << "\t"
                   << "latexname=$P(\\mu_\\text{" << subdet << "}=\\mu_\\text{ideal})";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto dmu : vPValueMeanEqualIdeal)
         summaryfile << "\t" << dmu;
       summaryfile << "\n";
     }
     if (!vPValueRMSEqualIdeal.empty()) {
       summaryfile << "  p_RMS_equals_ideal_" << subdet;
       if (plotinfo.variable == "medianY")
         summaryfile << "_y";
       summaryfile << "\t"
                   << "latexname=$P(\\sigma_\\text{" << subdet << "}=\\sigma_\\text{ideal})";
       if (plotinfo.variable == "medianY")
         summaryfile << "^{y}";
       summaryfile << "$\t"
                   << "format={:.3g}\t"
                   << "latexformat=${:.3g}$";
       for (auto dmu : vPValueRMSEqualIdeal)
         summaryfile << "\t" << dmu;
       summaryfile << "\n";
     }
   }
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::plotChi2(const char* inputFile) {
   // Opens the file (it should be OfflineValidation(Parallel)_result.root)
   // and reads and plots the norm_chi^2 and h_chi2Prob -distributions.
 
   Bool_t errorflag = kFALSE;
   TFile* fi1 = TFile::Open(inputFile, "read");
   if (fi1 != nullptr) {
     if (fi1->GetDirectory("TrackerOfflineValidation/GlobalTrackVariables") == nullptr) {
       errorflag = kTRUE;
     }
   } else {
     errorflag = kTRUE;
   }
   if (errorflag) {
     std::cout << "PlotAlignmentValidation::plotChi2: Can't find GlobalTrackVariables given file,"
               << " no chi^2-plots produced" << std::endl;
     return;
   }
 
   auto normchi = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_normchi2");
   // please remove following line once bug is fixed with ROOT Version: >6.27/01
   normchi->GetFrame()->ResetBit(kCanDelete);
 
   auto chiprob = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_chi2Prob");
   // please remove following line once bug is fixed with ROOT Version: >6.27/01
   chiprob->GetFrame()->ResetBit(kCanDelete);
 
   if (normchi == nullptr || chiprob == nullptr) {
     errorflag = kTRUE;
   }
   if (errorflag) {
     std::cout << "PlotAlignmentValidation::plotChi2: Can't find data from given file,"
               << " no chi^2-plots produced" << std::endl;
     return;
   }
 
   TLegend* legend = nullptr;
   for (auto primitive : *normchi->GetListOfPrimitives()) {
     legend = dynamic_cast<TLegend*>(primitive);
     if (legend)
       break;
   }
   if (legend) {
     openSummaryFile();
     summaryfile << "ntracks";
     for (auto alignment : sourceList) {
       summaryfile << "\t";
       TString title = alignment->getName();
       int color = alignment->getLineColor();
       int style = alignment->getLineStyle();
       bool foundit = false;
       for (auto entry : *legend->GetListOfPrimitives()) {
         TLegendEntry* legendentry = dynamic_cast<TLegendEntry*>(entry);
         assert(legendentry);
         TH1* h = dynamic_cast<TH1*>(legendentry->GetObject());
         if (!h)
           continue;
         if (legendentry->GetLabel() == title && h->GetLineColor() == color && h->GetLineStyle() == style) {
           foundit = true;
           summaryfile << h->GetEntries();
           break;
         }
       }
       if (!foundit) {
         summaryfile << 0;
       }
     }
     summaryfile << "\n";
   }
 
   chiprob->Draw();
   normchi->Draw();
 
   // PNG,EPS,PDF files
   normchi->Print((outputDir + "/h_normchi2.png").c_str());
   chiprob->Print((outputDir + "/h_chi2Prob.png").c_str());
   normchi->Print((outputDir + "/h_normchi2.eps").c_str());
   chiprob->Print((outputDir + "/h_chi2Prob.eps").c_str());
   normchi->Print((outputDir + "/h_normchi2.pdf").c_str());
   chiprob->Print((outputDir + "/h_chi2Prob.pdf").c_str());
 
   // ROOT files
   TFile fi2((outputDir + "/h_normchi2.root").c_str(), "recreate");
   normchi->Write();
   fi2.Close();
 
   TFile fi3((outputDir + "/h_chi2Prob.root").c_str(), "recreate");
   chiprob->Write();
   fi3.Close();
 
   fi1->Close();
   delete fi1;
 }
 
 //------------------------------------------------------------------------------
 THStack* PlotAlignmentValidation::addHists(
     const TString& selection, const TString& residType, TLegend** myLegend, bool printModuleIds, bool validforphase0) {
   enum ResidType {
     xPrimeRes,
     yPrimeRes,
     xPrimeNormRes,
     yPrimeNormRes,
     xRes,
     yRes,
     xNormRes, /*yResNorm*/
     ResXvsXProfile,
     ResXvsYProfile,
     ResYvsXProfile,
     ResYvsYProfile
   };
   ResidType rType = xPrimeRes;
   if (residType == "xPrime")
     rType = xPrimeRes;
   else if (residType == "yPrime")
     rType = yPrimeRes;
   else if (residType == "xPrimeNorm")
     rType = xPrimeNormRes;
   else if (residType == "yPrimeNorm")
     rType = yPrimeNormRes;
   else if (residType == "x")
     rType = xRes;
   else if (residType == "y")
     rType = yRes;
   else if (residType == "xNorm")
     rType = xNormRes;
   // else if (residType == "yNorm") rType = yResNorm;
   else if (residType == "ResXvsXProfile")
     rType = ResXvsXProfile;
   else if (residType == "ResYvsXProfile")
     rType = ResYvsXProfile;
   else if (residType == "ResXvsYProfile")
     rType = ResXvsYProfile;
   else if (residType == "ResYvsYProfile")
     rType = ResYvsYProfile;
   else {
     std::cout << "PlotAlignmentValidation::addHists: Unknown residual type " << residType << std::endl;
     return nullptr;
   }
 
   std::cout << "PlotAlignmentValidation::addHists: using selection " << selection << std::endl;
   THStack* retHistoStack = new THStack("hstack", "");
   if (myLegend != nullptr)
     if (*myLegend == nullptr) {
       *myLegend = new TLegend(0.17, 0.80, 0.85, 0.88);
     }
 
   for (std::vector<TkOfflineVariables*>::iterator itSourceFile = sourceList.begin(); itSourceFile != sourceList.end();
        ++itSourceFile) {
     std::vector<TString> histnames;
 
     TFile* f = (*itSourceFile)->getFile();
     TTree* tree = (*itSourceFile)->getTree();
     int myLineColor = (*itSourceFile)->getLineColor();
     int myLineStyle = (*itSourceFile)->getLineStyle();
     TString myLegendName = (*itSourceFile)->getName();
     TH1* h = nullptr;   // becomes result
     UInt_t nEmpty = 0;  // selected, but empty hists
     Long64_t nentries = tree->GetEntriesFast();
     if (!f || !tree) {
       std::cout << "PlotAlignmentValidation::addHists: no tree or no file" << std::endl;
       return nullptr;
     }
 
     bool histnamesfilled = false;
     int phase = (bool)(f->Get("TrackerOfflineValidation/Pixel/P1PXBBarrel_1"));
     if (residType.Contains("Res") && residType.Contains("Profile")) {
       TString basename = TString(residType)
                              .ReplaceAll("Res", "p_res")
                              .ReplaceAll("vs", "")
                              .ReplaceAll("Profile", "_");  //gives e.g.: p_resXX_
       if (selection == "subDetId==1") {
         if (phase == 1)
           histnames.push_back(TString(basename) += "P1PXBBarrel_1");
         else
           histnames.push_back(TString(basename) += "TPBBarrel_1");
         histnamesfilled = true;
       } else if (selection == "subDetId==2") {
         if (phase == 1) {
           histnames.push_back(TString(basename) += "P1PXECEndcap_2");
           histnames.push_back(TString(basename) += "P1PXECEndcap_3");
         } else {
           histnames.push_back(TString(basename) += "TPEEndcap_2");
           histnames.push_back(TString(basename) += "TPEEndcap_3");
         }
         histnamesfilled = true;
       } else if (selection == "subDetId==3") {
         histnames.push_back(TString(basename) += "TIBBarrel_1");
         histnamesfilled = true;
       } else if (selection == "subDetId==4") {
         histnames.push_back(TString(basename) += "TIDEndcap_2");
         histnames.push_back(TString(basename) += "TIDEndcap_3");
         histnamesfilled = true;
       } else if (selection == "subDetId==5") {
         histnames.push_back(TString(basename) += "TOBBarrel_4");
         histnamesfilled = true;
       } else if (selection == "subDetId==6") {  //whole TEC - doesn't happen by default but easy enough to account for
         histnames.push_back(TString(basename) += "TECEndcap_5");
         histnames.push_back(TString(basename) += "TECEndcap_6");
         histnamesfilled = true;
       } else if (selection == "subDetId==6 && ring <= 4") {
         //There are multiple with the same name and all are needed, so give the full path.  For these TFile::Get is used later instead of FindKeyAny.
         for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
           for (int iDisk = 1; iDisk <= 9; iDisk++)
             for (int iSide = 1; iSide <= 2; iSide++)
               for (int iPetal = 1; iPetal <= 8; iPetal++)
                 for (int iRing = 1; iRing <= 4 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9); iRing++)
                 //in the higher disks, the inner rings go away.  But the numbering in the file structure removes the higher numbers
                 // so the numbers there do not correspond to the actual ring numbers
                 {
                   std::stringstream s;
                   s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
                     << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
                   histnames.push_back(TString(s.str()));
                 }
         histnamesfilled = true;
       } else if (selection == "subDetId==6 && ring > 4") {
         //There are multiple with the same name and all are needed, so give the full path.  For these TFile::Get is used later instead of FindKeyAny.
         for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
           for (int iDisk = 1; iDisk <= 9; iDisk++)
             for (int iSide = 1; iSide <= 2; iSide++)
               for (int iPetal = 1; iPetal <= 8; iPetal++)
                 for (int iRing = 5 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
                      iRing <= 7 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
                      iRing++)
                 //in the higher disks, the inner rings go away.  But the numbering in the file structure removes the higher numbers
                 // so the numbers there do not correspond to the actual ring numbers
                 {
                   std::stringstream s;
                   s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
                     << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
                   histnames.push_back(TString(s.str()));
                 }
         histnamesfilled = true;
       }
     }
 
     Long64_t nSel = 0;
     if (histnamesfilled && !histnames.empty()) {
       nSel = (Long64_t)histnames.size();
     }
     if (!histnamesfilled) {
       // first loop on tree to find out which entries (i.e. modules) fulfill the selection
       // 'Entry$' gives the entry number in the tree
       nSel = tree->Draw("Entry$", selection, "goff");
       if (nSel == -1)
         return nullptr;  // error in selection
       if (nSel == 0) {
         std::cout << "PlotAlignmentValidation::addHists: no selected module." << std::endl;
         return nullptr;
       }
       // copy entry numbers that fulfil the selection
       const std::vector<double> selected(tree->GetV1(), tree->GetV1() + nSel);
 
       std::vector<double>::const_iterator iterEnt = selected.begin();
 
       // second loop on tree:
       // for each selected entry get the hist from the file and merge
       TkOffTreeVariables* treeMem = nullptr;  // ROOT will initialise
       tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
       for (Long64_t i = 0; i < nentries; i++) {
         if (i < *iterEnt - 0.1               // smaller index (with tolerance): skip
             || iterEnt == selected.end()) {  // at the end: skip
           continue;
         } else if (TMath::Abs(i - *iterEnt) < 0.11) {
           ++iterEnt;  // take this entry!
         } else
           std::cout << "Must not happen: " << i << " " << *iterEnt << std::endl;
 
         tree->GetEntry(i);
         if (printModuleIds) {
           std::cout << treeMem->moduleId << ": " << treeMem->entries << " entries" << std::endl;
         }
         if (treeMem->entries <= 0) {  // little speed up: skip empty hists
           ++nEmpty;
           continue;
         }
         TString hName;
         switch (rType) {
           case xPrimeRes:
             hName = treeMem->histNameX.c_str();
             break;
           case yPrimeRes:
             hName = treeMem->histNameY.c_str();
             break;
           case xPrimeNormRes:
             hName = treeMem->histNameNormX.c_str();
             break;
           case yPrimeNormRes:
             hName = treeMem->histNameNormY.c_str();
             break;
           case xRes:
             hName = treeMem->histNameLocalX.c_str();
             break;
           case yRes:
             hName = treeMem->histNameLocalY.c_str();
             break;
           case xNormRes:
             hName = treeMem->histNameNormLocalX.c_str();
             break;
             /*case yResNorm:      hName = treeMem->histNameNormLocalY.c_str(); break;*/
           case ResXvsXProfile:
             hName = treeMem->profileNameResXvsX.c_str();
             break;
           case ResXvsYProfile:
             hName = treeMem->profileNameResXvsY.c_str();
             break;
           case ResYvsXProfile:
             hName = treeMem->profileNameResYvsX.c_str();
             break;
           case ResYvsYProfile:
             hName = treeMem->profileNameResYvsY.c_str();
             break;
         }
         histnames.push_back(hName);
       }
     }
 
     for (std::vector<TString>::iterator ithistname = histnames.begin(); ithistname != histnames.end(); ++ithistname) {
       if (phase == 0 && !validforphase0)
         break;
       TH1* newHist;
       if (ithistname->Contains("/")) {
         newHist = (TH1*)f->Get(*ithistname);
       } else {
         TKey* histKey = f->FindKeyAny(*ithistname);
         newHist = (histKey ? static_cast<TH1*>(histKey->ReadObj()) : nullptr);
       }
       if (!newHist) {
         std::cout << "Hist " << *ithistname << " not found in file, break loop." << std::endl;
         break;
       }
       if (newHist->GetEntries() == 0) {
         nEmpty++;
         continue;
       }
       newHist->SetLineColor(myLineColor);
       newHist->SetLineStyle(myLineStyle);
       if (!h) {  // first hist: clone, but rename keeping only first part of name
         TString name(newHist->GetName());
         Ssiz_t pos_ = 0;
         for (UInt_t i2 = 0; i2 < 3; ++i2)
           pos_ = name.Index("_", pos_ + 1);
         name = name(0, pos_);  // only up to three '_'
         h = static_cast<TH1*>(newHist->Clone("summed_" + name));
         //      TString myTitle = Form("%s: %lld modules", selection, nSel);
         //  h->SetTitle( myTitle );
       } else {  // otherwise just add
         h->Add(newHist);
       }
       delete newHist;
     }
 
     std::cout << "PlotAlignmentValidation::addHists"
               << "Result is merged from " << nSel - nEmpty << " hists, " << nEmpty << " hists were empty." << std::endl;
 
     if (nSel - nEmpty == 0)
       continue;
 
     if (myLegend != nullptr)
       (*myLegend)->AddEntry(h, myLegendName, "L");
 
     retHistoStack->Add(h);
   }
 
   return retHistoStack;
 }
 
 //------------------------------------------------------------------------------
 /*! \fn fitGauss
  *  \brief Operate a Gaussian fit to the given histogram
  */
 TF1* PlotAlignmentValidation::fitGauss(TH1* hist, int color) {
   //1. fits a Gauss function to the inner range of abs(2 rms)
   //2. repeates the Gauss fit in a 2 sigma range around mean of first fit
   //returns mean and sigma from fit in micron
   if (!hist || hist->GetEntries() < 20)
     return nullptr;
 
   float mean = hist->GetMean();
   float sigma = hist->GetRMS();
   std::string functionname = "gaussian_";
   functionname += hist->GetName();
   TF1* func = new TF1(functionname.c_str(), "gaus", mean - 2. * sigma, mean + 2. * sigma);
 
   func->SetLineColor(color);
   func->SetLineStyle(2);
   if (0 == hist->Fit(func, "QNR")) {  // N: do not blow up file by storing fit!
     mean = func->GetParameter(1);
     sigma = func->GetParameter(2);
     // second fit: three sigma of first fit around mean of first fit
     func->SetRange(mean - 3. * sigma, mean + 3. * sigma);
     // I: integral gives more correct results if binning is too wide
     // L: Likelihood can treat empty bins correctly (if hist not weighted...)
     if (0 == hist->Fit(func, "Q0ILR")) {
       if (hist->GetFunction(func->GetName())) {  // Take care that it is later on drawn:
                                                  //hist->GetFunction(func->GetName())->ResetBit(TF1::kNotDraw);
       }
     }
   }
   return func;
 }
 
 //------------------------------------------------------------------------------
 /*! \fn storeHistogramInRootfile
  *  \brief Store the histogram and the gaussian function resulting from the fitGauss function into a root file
  */
 void PlotAlignmentValidation::storeHistogramInRootfile(TH1* hist) {
   //Store histogram and fit function in the root summary file
   rootsummaryfile->cd();
   hist->Write();
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::scaleXaxis(TH1* hist, Int_t scale) {
   Double_t xmin = hist->GetXaxis()->GetXmin();
   Double_t xmax = hist->GetXaxis()->GetXmax();
   hist->GetXaxis()->SetLimits(xmin * scale, xmax * scale);
 }
 
 //------------------------------------------------------------------------------
 TObject* PlotAlignmentValidation::findObjectFromCanvas(TCanvas* canv, const char* className, Int_t n) {
   // Finds the n-th instance of the given class from the canvas
   TIter next(canv->GetListOfPrimitives());
   TObject* obj = nullptr;
   Int_t found = 0;
   while ((obj = next())) {
     if (strncmp(obj->ClassName(), className, 10) == 0) {
       if (++found == n)
         return obj;
     }
   }
 
   return nullptr;
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::setTitleStyle(
     TNamed& hist, const char* titleX, const char* titleY, int subDetId, bool isSurfaceDeformation, TString secondline) {
   std::stringstream title_Xaxis;
   std::stringstream title_Yaxis;
   TString titleXAxis = titleX;
   TString titleYAxis = titleY;
   if (titleXAxis != "" && titleYAxis != "")
     std::cout << "plot " << titleXAxis << " vs " << titleYAxis << std::endl;
 
   hist.SetTitle("");
   TkAlStyle::drawStandardTitle();
 
   //Thanks Candice!
   TString subD;
   switch (subDetId) {
     case 1:
       subD = "BPIX";
       break;
     case 2:
       subD = "FPIX";
       break;
     case 3:
       subD = "TIB";
       break;
     case 4:
       subD = "TID";
       break;
     case 5:
       subD = "TOB";
       break;
     case 6:
       subD = "TEC";
       break;
   }
 
   TPaveText* text2;
   if (!isSurfaceDeformation) {
     text2 = new TPaveText(0.7, 0.3, 0.9, 0.6, "brNDC");
   } else {
     std::cout << "Surface Deformation" << std::endl;
     text2 = new TPaveText(0.8, 0.75, 0.9, 0.9, "brNDC");
   }
   text2->SetTextSize(0.06);
   text2->SetTextFont(42);
   text2->SetFillStyle(0);
   text2->SetBorderSize(0);
   text2->SetMargin(0.01);
   text2->SetTextAlign(12);  // align left
   text2->AddText(0.01, 0.75, subD);
   if (secondline != "") {
     text2->AddText(0.01, 0.25, secondline);
   }
   text2->Draw();
 }
 
 //------------------------------------------------------------------------------
 /*! \fn 
  *  \brief 
  */
 void PlotAlignmentValidation::setHistStyle(TH1& hist, const char* titleX, const char* titleY, int color) {
   std::stringstream title_Xaxis;
   std::stringstream title_Yaxis;
   TString titleXAxis = titleX;
   TString titleYAxis = titleY;
 
   if (titleXAxis.Contains("Phi"))
     title_Xaxis << titleX << "[rad]";
   else if (titleXAxis.Contains("meanX"))
     title_Xaxis << "#LTx'_{pred}-x'_{hit}#GT[#mum]";
   else if (titleXAxis.Contains("meanY"))
     title_Xaxis << "#LTy'_{pred}-y'_{hit}#GT[#mum]";
   else if (titleXAxis.Contains("rmsX"))
     title_Xaxis << "RMS(x'_{pred}-x'_{hit})[#mum]";
   else if (titleXAxis.Contains("rmsY"))
     title_Xaxis << "RMS(y'_{pred}-y'_{hit})[#mum]";
   else if (titleXAxis.Contains("meanNormX"))
     title_Xaxis << "#LTx'_{pred}-x'_{hit}/#sigma#GT";
   else if (titleXAxis.Contains("meanNormY"))
     title_Xaxis << "#LTy'_{pred}-y'_{hit}/#sigma#GT";
   else if (titleXAxis.Contains("rmsNormX"))
     title_Xaxis << "RMS(x'_{pred}-x'_{hit}/#sigma)";
   else if (titleXAxis.Contains("rmsNormY"))
     title_Xaxis << "RMS(y'_{pred}-y'_{hit}/#sigma)";
   else if (titleXAxis.Contains("meanLocalX"))
     title_Xaxis << "#LTx_{pred}-x_{hit}#GT[#mum]";
   else if (titleXAxis.Contains("rmsLocalX"))
     title_Xaxis << "RMS(x_{pred}-x_{hit})[#mum]";
   else if (titleXAxis.Contains("meanNormLocalX"))
     title_Xaxis << "#LTx_{pred}-x_{hit}/#sigma#GT[#mum]";
   else if (titleXAxis.Contains("rmsNormLocalX"))
     title_Xaxis << "RMS(x_{pred}-x_{hit}/#sigma)[#mum]";
   else if (titleXAxis.Contains("medianX"))
     title_Xaxis << "median(x'_{pred}-x'_{hit})[#mum]";
   else if (titleXAxis.Contains("medianY"))
     title_Xaxis << "median(y'_{pred}-y'_{hit})[#mum]";
   else
     title_Xaxis << titleX << "[cm]";
 
   if (hist.IsA()->InheritsFrom(TH1F::Class()))
     hist.SetLineColor(color);
   if (hist.IsA()->InheritsFrom(TProfile::Class())) {
     hist.SetMarkerStyle(20);
     hist.SetMarkerSize(0.8);
     hist.SetMarkerColor(color);
   }
 
   hist.GetXaxis()->SetTitle((title_Xaxis.str()).c_str());
 
   double binning = (hist.GetXaxis()->GetXmax() - hist.GetXaxis()->GetXmin()) / hist.GetNbinsX();
   title_Yaxis.precision(2);
 
   if (((titleYAxis.Contains("layer") || titleYAxis.Contains("ring")) && titleYAxis.Contains("subDetId")) ||
       titleYAxis.Contains("#modules")) {
     title_Yaxis << "number of modules";
     if (TString(title_Xaxis.str()).Contains("[#mum]"))
       title_Yaxis << " / " << binning << " #mum";
     else if (TString(title_Xaxis.str()).Contains("[cm]"))
       title_Yaxis << " / " << binning << " cm";
     else
       title_Yaxis << " / " << binning;
   } else
     title_Yaxis << titleY << "[cm]";
 
   hist.GetYaxis()->SetTitle((title_Yaxis.str()).c_str());
 
   hist.GetXaxis()->SetTitleFont(42);
   hist.GetYaxis()->SetTitleFont(42);
 }
 
 //------------------------------------------------------------------------------
 std::string PlotAlignmentValidation::getSelectionForDMRPlot(int minHits, int subDetId, int direction, int layer) {
   std::ostringstream builder;
   builder << "entries >= " << minHits;
   builder << " && subDetId == " << subDetId;
   if (direction != 0) {
     if (subDetId == 2) {  // FPIX is split by zDirection
       builder << " && zDirection == " << direction;
     } else {
       builder << " && rDirection == " << direction;
     }
   }
   if (layer > 0) {
     builder << " && layer == " << layer;
   }
   return builder.str();
 }
 
 std::string PlotAlignmentValidation::getVariableForDMRPlot(
     const std::string& histoname, const std::string& variable, int nbins, double min, double max) {
   std::ostringstream builder;
   builder << variable << ">>" << histoname << "(" << nbins << "," << min << "," << max << ")";
   return builder.str();
 }
 
 //------------------------------------------------------------------------------
 void PlotAlignmentValidation::setDMRHistStyleAndLegend(TH1F* h,
                                                        PlotAlignmentValidation::DMRPlotInfo& plotinfo,
                                                        int direction,
                                                        int layer) {
   TF1* fitResults = nullptr;
 
   h->SetDirectory(nullptr);
 
   // The whole DMR plot is plotted with wider line than the split plots
   // If only split plots are plotted, they will be stronger too, though
   h->SetLineWidth((direction == 0 || (plotinfo.plotSplits && !plotinfo.plotPlain)) ? 2 : 1);
 
   // These lines determine the style of the plots according to rules:
   // -If the plot is for direction != 0, +1 or +2 is added to the given style for distinction
   // -However if only direction split plots are to be plotted, the additions should be 0 and +1 respectively
   // -Modulo 4 arithmetic, because the styles run from 1..4
   int linestyle = plotinfo.vars->getLineStyle() - 1, linestyleplus = 0;
   if (direction == 1) {
     linestyleplus = 1;
   }
   if (direction == -1) {
     linestyleplus = 2;
   }
   if (direction != 0 && plotinfo.plotSplits && !plotinfo.plotPlain) {
     linestyleplus--;
   }
   linestyle = (linestyle + linestyleplus) % 4 + 1;
 
   int linecolor = plotinfo.vars->getLineColor();
   if (plotinfo.plotLayers && layer > 0) {
     linecolor += layer - 1;
   }
 
   if (plotinfo.firsthisto) {
     setHistStyle(*h, plotinfo.variable.c_str(), "#modules", 1);  //set color later
     plotinfo.firsthisto = false;
   }
 
   h->SetLineColor(linecolor);
   h->SetLineStyle(linestyle);
 
   if (plotinfo.maxY < h->GetMaximum()) {
     plotinfo.maxY = h->GetMaximum();
   }
 
   //fit histogram for median and mean
   if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
       plotinfo.variable == "meanY") {
     fitResults = fitGauss(h, linecolor);
   }
 
   plotinfo.hstack->Add(h);
 
   std::ostringstream legend;
   legend.precision(3);
   legend << std::fixed;  // to always show 3 decimals
 
   // Legend: header part
   if (direction == -1 && plotinfo.subDetId != 2) {
     legend << "rDirection < 0";
   } else if (direction == 1 && plotinfo.subDetId != 2) {
     legend << "rDirection > 0";
   } else if (direction == -1 && plotinfo.subDetId == 2) {
     legend << "zDirection < 0";
   } else if (direction == 1 && plotinfo.subDetId == 2) {
     legend << "zDirection > 0";
   } else {
     legend << plotinfo.vars->getName();
     if (layer > 0) {
       // TEC and TID have discs, the rest have layers
       if (plotinfo.subDetId == 4 || plotinfo.subDetId == 6)
         legend << ", disc ";
       else
         legend << ", layer ";
       legend << layer << "";
     }
   }
 
   plotinfo.legend->AddEntry(h, legend.str().c_str(), "l");
   legend.str("");
 
   // Legend: Statistics
   double mean = 0.0, meanerror = 0.0, rms = 0.0, rmserror = 0.0;
   TString rmsname, units;
   bool showdeltamu =
       (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && plotinfo.plotSplits && plotinfo.plotPlain && direction == 0);
   if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
       plotinfo.variable == "meanY" || plotinfo.variable == "rmsX" || plotinfo.variable == "rmsY") {
     if (useFit_ && fitResults) {
       mean = fitResults->GetParameter(1) * 10000;
       meanerror = fitResults->GetParError(1) * 10000;
       rms = fitResults->GetParameter(2) * 10000;
       rmserror = fitResults->GetParError(2) * 10000;
       rmsname = "#sigma";
       delete fitResults;
     } else {
       mean = h->GetMean(1) * 10000;
       meanerror = h->GetMeanError(1) * 10000;
       rms = h->GetRMS(1) * 10000;
       rmserror = h->GetRMSError(1) * 10000;
       rmsname = "rms";
     }
     units = " #mum";
   } else if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
              plotinfo.variable == "rmsNormY") {
     mean = h->GetMean(1);
     meanerror = h->GetMeanError(1);
     rms = h->GetRMS(1);
     rmserror = h->GetRMSError(1);
     rmsname = "rms";
     units = "";
   }
   if (showMean_) {
     legend << " #mu = " << mean;
     if (showMeanError_)
       legend << " #pm " << meanerror;
     legend << units;
     if (showRMS_ || showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
       legend << ", ";
   }
   if (showRMS_) {
     legend << " " << rmsname << " = " << rms;
     if (showRMSError_)
       legend << " #pm " << rmserror;
     legend << units;
     if (showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
       legend << ", ";
   }
 
   if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && /*!plotinfo.plotLayers && */ layer == 0 &&
       direction == 0) {
     vmean.push_back(mean);
     vrms.push_back(rms);
     vmeanerror.push_back(meanerror);
     TH1F* ideal = (TH1F*)plotinfo.hstack->GetHists()->At(0);
     TH1F* h = plotinfo.h;
     if (h->GetRMS() >= ideal->GetRMS()) {
       vAlignmentUncertainty.push_back(sqrt(pow(h->GetRMS(), 2) - pow(ideal->GetRMS(), 2)));
     } else {
       vAlignmentUncertainty.push_back(nan(""));
     }
     float p = (float)resampleTestOfEqualMeans(ideal, h, 10000);
     vPValueMeanEqualIdeal.push_back(p);
     p = resampleTestOfEqualRMS(ideal, h, 10000);
     vPValueRMSEqualIdeal.push_back(p);
   }
 
   // Legend: Delta mu for split plots
   if (showdeltamu) {
     float factor = 10000.0f;
     if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
         plotinfo.variable == "rmsNormY") {
       factor = 1.0f;
     }
     float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
     legend << "#Delta#mu = " << deltamu << units;
     if ((showModules_ || showUnderOverFlow_) && !twolines_)
       legend << ", ";
 
     if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && /*!plotinfo.plotLayers && */ layer == 0 &&
         direction == 0) {
       vdeltamean.push_back(deltamu);
       if (plotinfo.h1->GetEntries() && plotinfo.h2->GetEntries()) {
         float p = (float)resampleTestOfEqualMeans(plotinfo.h1, plotinfo.h2, 10000);
         vPValueEqualSplitMeans.push_back(p);
       }
     }
   }
 
   if (twolines_) {
     plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
     plotinfo.legend->AddEntry((TObject*)nullptr, "", "");
     legend.str("");
   }
 
   if (!showUnderOverFlow_ && showModules_) {
     legend << (int)h->GetEntries() << " modules";
   }
   if (showUnderOverFlow_) {
     if (showModules_) {
       legend << (int)h->GetEntries() << " modules ("
              << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " outside range)";
     } else {
       legend << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " modules outside range";
     }
   }
   plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
 
   // Scale the x-axis (cm to um), if needed
   if (plotinfo.variable.find("Norm") == std::string::npos)
     scaleXaxis(h, 10000);
 }
 
 /*!
  * \fn plotDMRHistogram 
  * \brief Create the DMR histrogram using data stored in trees and store them in the plotinfo structure.
  */
 
 void PlotAlignmentValidation::plotDMRHistogram(PlotAlignmentValidation::DMRPlotInfo& plotinfo,
                                                int direction,
                                                int layer,
                                                std::string subdet) {
   TH1F* h = nullptr;
   //Create a name for the histogram that summarize all relevant information: name of the geometry, variable plotted, structure, layer, and whether the modules considered point inward or outward.
 
   TString histoname = "";
   if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY")
     histoname = "median";
   else if (plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY")
     histoname = "DrmsNR";
   histoname += "_";
   histoname += plotinfo.vars->getName();
   histoname.ReplaceAll(" ", "_");
   histoname += "_";
   histoname += subdet.c_str();
   if (plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormY")
     histoname += "_y";
   if (layer != 0) {
     if (subdet == "TID" || subdet == "TEC")
       histoname += "_disc";
     else
       histoname += "_layer";
     histoname += std::to_string(layer);
   }
   if (direction == -1) {
     histoname += "_minus";
   } else if (direction == 1) {
     histoname += "_plus";
   } else {
     histoname += "";
   }
 
   //Plotting
   std::string plotVariable =
       getVariableForDMRPlot(histoname.Data(), plotinfo.variable, plotinfo.nbins, plotinfo.min, plotinfo.max);
   std::string selection = "";
   if (plotinfo.filterName.empty()) {
     //Use only default selection and no filter
     selection = getSelectionForDMRPlot(plotinfo.minHits, plotinfo.subDetId, direction, layer);
     plotinfo.vars->getTree()->Draw(plotVariable.c_str(), selection.c_str(), "goff");
     if (gDirectory)
       gDirectory->GetObject(histoname.Data(), h);
     if (h && h->GetEntries() > 0) {
       if (direction == -1) {
         plotinfo.h1 = h;
       } else if (direction == 1) {
         plotinfo.h2 = h;
       } else {
         plotinfo.h = h;
       }
     }
   } else {
     TTreeReader reader(plotinfo.vars->getTree());
     TTreeReaderValue<Float_t> varToPlot(reader, plotinfo.variable.c_str());
     TTreeReaderValue<unsigned int> _entries(reader, "entries");
     TTreeReaderValue<unsigned int> _subDetId(reader, "subDetId");
     TTreeReaderValue<unsigned int> _moduleId(reader, "moduleId");
     TTreeReaderValue<Float_t> _zDirection(reader, "zDirection");
     TTreeReaderValue<Float_t> _rDirection(reader, "rDirection");
     TTreeReaderValue<unsigned int> _layer(reader, "layer");
     std::string badModulesFile_ = plotinfo.filterName;
     TFile* fBadModules = new TFile(badModulesFile_.c_str(), "READ");
     TTree* tBadModules = (TTree*)fBadModules->Get("alignTree");
     TTreeReader readerBad(tBadModules);
     TTreeReaderValue<int> _valid(readerBad, "valid");
     TTreeReaderValue<int> _bad_id(readerBad, "id");
     TTreeReaderValue<double> _bad_lumi(readerBad, "lumi");
 
     //Record which modules were used
     std::ofstream fUsedModules;
     fUsedModules.open("usedModules.txt", std::ios::out | std::ios::app);
 
     //Filter on modules by hand together with base selection
     for (uint i = 0; i < plotinfo.vars->getTree()->GetEntries(); i++) {
       reader.SetEntry(i);
       if (*_entries < uint(plotinfo.minHits))
         continue;
       if (*_subDetId != uint(plotinfo.subDetId))
         continue;
       if (direction != 0) {
         if (plotinfo.subDetId == 2) {  // FPIX is split by zDirection
           if (*_zDirection != direction)
             continue;
         } else {
           if (*_rDirection != direction)
             continue;
         }
       }
       if (layer > 0) {
         if (*_layer != uint(layer))
           continue;
       }
       bool isBadModule = false;
       for (int ibad = 0; ibad < tBadModules->GetEntries(); ibad++) {
         readerBad.SetEntry(ibad);
         //if (*_valid == 0) {continue;} //only modules that failed 0 times are OK = very strict
         if (subdet == "BPIX" || subdet == "FPIX") {
           if (*_bad_lumi <= 2.0)
             continue;
         } else {
           if (*_bad_lumi <= 7.0)
             continue;
         }
         //modules that misbehave for less than 2/fb are OK = mild strict
         if (*_moduleId == uint(*_bad_id))
           isBadModule = true;
       }
 
       if (isBadModule)
         continue;
       fUsedModules << *_moduleId << "\n";
       if (h) {
         h->Fill(*varToPlot);
       }
     }
 
     //Finalize
     fUsedModules.close();
     fBadModules->Close();
     if (h) {
       h->SetName(histoname.Data());
     }
   }
 
   //Store histogram
   if (h && h->GetEntries() > 0) {
     if (direction == -1) {
       plotinfo.h1 = h;
     } else if (direction == 1) {
       plotinfo.h2 = h;
     } else {
       plotinfo.h = h;
     }
   }
   if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormX" ||
       plotinfo.variable == "rmsNormY") {
     storeHistogramInRootfile(h);
   }
 }
 
 void PlotAlignmentValidation::modifySSHistAndLegend(THStack* hs, TLegend* legend) {
   // Add mean-y-values to the legend and scale the histograms.
 
   Double_t legendY = 0.80;
   bool hasheader = (TkAlStyle::legendheader != "");
   if (hasheader)
     legend->SetHeader(TkAlStyle::legendheader);
   legend->SetFillStyle(0);
   int legendsize = hs->GetHists()->GetSize() + hasheader;
 
   if (legendsize > 3)
     legendY -= 0.01 * (legendsize - 3);
   if (bigtext_) {
     legendY -= 0.05;
   }
   if (legendY < 0.6) {
     std::cerr << "Warning: Huge legend!" << std::endl;
     legendY = 0.6;
   }
   legend->SetY1(legendY);
   if (bigtext_)
     legend->SetTextSize(TkAlStyle::textSize);
 
   // Loop over all profiles
   TProfile* prof = nullptr;
   TIter next(hs->GetHists());
   Int_t index = hasheader;  //if hasheader, first entry is the header itself
   while ((prof = (TProfile*)next())) {
     //Scaling: from cm to um
     Double_t scale = 10000;
     prof->Scale(scale);
 
     Double_t stats[6] = {0};
     prof->GetStats(stats);
 
     std::ostringstream legendtext;
     legendtext.precision(3);
     legendtext << std::fixed;  // to always show 3 decimals
     legendtext << ": y mean = " << stats[4] / stats[0] * scale << " #mum";
 
     TLegendEntry* entry = (TLegendEntry*)legend->GetListOfPrimitives()->At(index);
     if (entry == nullptr)
       std::cout << "PlotAlignmentValidation::PlotAlignmentValidation::modifySSLegend: Bad legend!" << std::endl;
     else
       entry->SetLabel((entry->GetLabel() + legendtext.str()).c_str());
     index++;
   }
 
   // Make some room for the legend
   hs->SetMaximum(hs->GetMaximum("nostack PE") * 1.3);
 }
 
 //random variable: \sigma_{X_1}-\sigma_{X_2}-\delta_{RMS}
 //is centered approx around 0
 //null hypothesis: \delta_{RMS}=0
 //so \delta_\sigma is a realization of this random variable
 //how probable is it to get our value of \delta_\sigma?
 //->p-value
 double PlotAlignmentValidation::resampleTestOfEqualRMS(TH1F* h1, TH1F* h2, int numSamples) {
   //vector to store realizations of random variable
   std::vector<double> diff;
   diff.clear();
   //"true" (in bootstrap terms) difference of the samples' RMS
   double rmsdiff = abs(h1->GetRMS() - h2->GetRMS());
   //means of the samples to calculate RMS
   double m1 = h1->GetMean();
   double m2 = h2->GetMean();
   //realization of random variable
   double d1 = 0;
   double d2 = 0;
   //mean of random variable
   double test_mean = 0;
   for (int i = 0; i < numSamples; i++) {
     d1 = 0;
     d2 = 0;
     for (int i = 0; i < h1->GetEntries(); i++) {
       d1 += h1->GetRandom() - m1;
     }
     for (int i = 0; i < h2->GetEntries(); i++) {
       d2 += h2->GetRandom() + m2;
     }
     d1 /= h1->GetEntries();
     d2 /= h2->GetEntries();
     diff.push_back(abs(d1 - d2 - rmsdiff));
     test_mean += abs(d1 - d2 - rmsdiff);
   }
   test_mean /= numSamples;
   edm::LogPrint("") << "test mean:" << test_mean;
   //p value
   double p = 0;
   for (double d : diff) {
     if (d > rmsdiff) {
       p += 1;
     }
   }
 
   p /= numSamples;
   return p;
 }
 
 //random variable: (\overline{X_1}-\mu_1)-(\overline{X_2}-\mu_2)
 //is centered approx around 0
 //null hypothesis: \mu_1-\mu_2=0
 //so \delta_\mu is a realization of this random variable
 //how probable is it to get our value of \delta_\mu?
 //->p-value
 double PlotAlignmentValidation::resampleTestOfEqualMeans(TH1F* h1, TH1F* h2, int numSamples) {
   //vector to store realization of random variable
   std::vector<double> diff;
   diff.clear();
   //"true" (in bootstrap terms) difference of the samples' means
   double meandiff = abs(h1->GetMean() - h2->GetMean());
   //realization of random variable
   double d1 = 0;
   double d2 = 0;
   //mean of random variable
   double test_mean = 0;
   for (int i = 0; i < numSamples; i++) {
     d1 = 0;
     d2 = 0;
     for (int i = 0; i < h1->GetEntries(); i++) {
       d1 += h1->GetRandom();
     }
     for (int i = 0; i < h2->GetEntries(); i++) {
       d2 += h2->GetRandom();
     }
     d1 /= h1->GetEntries();
     d2 /= h2->GetEntries();
     diff.push_back(abs(d1 - d2 - meandiff));
     test_mean += abs(d1 - d2 - meandiff);
   }
   test_mean /= numSamples;
   edm::LogPrint("") << "test mean:" << test_mean;
   //p-value
   double p = 0;
   for (double d : diff) {
     if (d > meandiff) {
       p += 1;
     }
   }
 
   p /= numSamples;
   return p;
 }
 
 float PlotAlignmentValidation::twotailedStudentTTestEqualMean(float t, float v) {
   return 2 * (1 - ROOT::Math::tdistribution_cdf(abs(t), v));
 }
 
 const TString PlotAlignmentValidation::summaryfilename = "OfflineValidationSummary";
 
 std::vector<TH1*> PlotAlignmentValidation::findmodule(TFile* f, unsigned int moduleid) {
   //TFile *f = TFile::Open(filename, "READ");
   TString histnamex;
   TString histnamey;
   //read necessary branch/folder
   auto t = (TTree*)f->Get("TrackerOfflineValidation/TkOffVal");
 
   TkOffTreeVariables* variables = nullptr;
   t->SetBranchAddress("TkOffTreeVariables", &variables);
   unsigned int number_of_entries = t->GetEntries();
   for (unsigned int i = 0; i < number_of_entries; i++) {
     t->GetEntry(i);
     if (variables->moduleId == moduleid) {
       histnamex = variables->histNameX;
       histnamey = variables->histNameY;
       break;
     }
   }
 
   std::vector<TH1*> h;
 
   auto h1 = (TH1*)f->FindObjectAny(histnamex);
   auto h2 = (TH1*)f->FindObjectAny(histnamey);
 
   h1->SetDirectory(nullptr);
   h2->SetDirectory(nullptr);
 
   h.push_back(h1);
   h.push_back(h2);
 
   return h;
 }
 
 void PlotAlignmentValidation::residual_by_moduleID(unsigned int moduleid) {
   TCanvas* cx = new TCanvas("x_residual");
   TCanvas* cy = new TCanvas("y_residual");
   TLegend* legendx = new TLegend(0.55, 0.7, 1, 0.9);
   TLegend* legendy = new TLegend(0.55, 0.7, 1, 0.9);
 
   legendx->SetTextSize(0.016);
   legendx->SetTextAlign(12);
   legendy->SetTextSize(0.016);
   legendy->SetTextAlign(12);
 
   for (auto it : sourceList) {
     TFile* file = it->getFile();
     int color = it->getLineColor();
     int linestyle = it->getLineStyle();  //this you set by doing h->SetLineStyle(linestyle)
     TString legendname = it->getName();  //this goes in the legend
     std::vector<TH1*> hist = findmodule(file, moduleid);
 
     TString histnamex = legendname + " NEntries: " + std::to_string(int(hist[0]->GetEntries()));
     hist[0]->SetTitle(histnamex);
     hist[0]->SetStats(false);
     hist[0]->Rebin(50);
     hist[0]->SetBit(TH1::kNoTitle);
     hist[0]->SetLineColor(color);
     hist[0]->SetLineStyle(linestyle);
     cx->cd();
     hist[0]->Draw("Same");
     legendx->AddEntry(hist[0], histnamex, "l");
 
     TString histnamey = legendname + " NEntries: " + std::to_string(int(hist[1]->GetEntries()));
     hist[1]->SetTitle(histnamey);
     hist[1]->SetStats(false);
     hist[1]->Rebin(50);
     hist[1]->SetBit(TH1::kNoTitle);
     hist[1]->SetLineColor(color);
     hist[1]->SetLineStyle(linestyle);
     cy->cd();
     hist[1]->Draw("Same");
     legendy->AddEntry(hist[1], histnamey, "l");
   }
 
   TString filenamex = "x_residual_" + std::to_string(moduleid);
   TString filenamey = "y_residual_" + std::to_string(moduleid);
   cx->cd();
   legendx->Draw();
   cx->SaveAs(outputDir + "/" + filenamex + ".root");
   cx->SaveAs(outputDir + "/" + filenamex + ".pdf");
   cx->SaveAs(outputDir + "/" + filenamex + ".png");
   cx->SaveAs(outputDir + "/" + filenamex + ".eps");
 
   cy->cd();
   legendy->Draw();
   cy->SaveAs(outputDir + "/" + filenamey + ".root");
   cy->SaveAs(outputDir + "/" + filenamey + ".pdf");
   cy->SaveAs(outputDir + "/" + filenamey + ".png");
   cy->SaveAs(outputDir + "/" + filenamey + ".eps");
 }