Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "Alignment/OfflineValidation/interface/PlotAlignmentValidation.h"
0002 
0003 /*! \class PlotAlignmentValidation
0004  *  \brief Class PlotAlignmentValidation
0005  *         Class used as the last step for Offline Track Validation tool.
0006  *         The main goal of this class is creating the plots regarding DMRs and Surface Deformations for modules and substructures.
0007  */
0008 
0009 //------------------------------------------------------------------------------
0010 /*! \fn PlotAlignmentValidation
0011  *  \brief Constructor for the class
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   // Force ROOT to use scientific notation even with smaller datasets
0023   TGaxis::SetMaxDigits(4);
0024   // (This sets a static variable: correct in .eps images but must be set
0025   // again manually when viewing the .root files)
0026 
0027   // Make ROOT calculate histogram statistics using all data,
0028   // regardless of displayed range
0029   TH1::StatOverflows(kTRUE);
0030 
0031   //show all information in the legend by default
0032   legendOptions(TkAlStyle::legendoptions);
0033 }
0034 
0035 //------------------------------------------------------------------------------
0036 /*! \fn PlotAlignmentValidation
0037  *  \brief Constructor for the class. This function also retrieves the list of root files used to produce DMRs and Surface Deformations
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 /*! \fn ~PlotAlignmentValidation
0047  *  \brief Default destructor
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  * \fn openSummaryFile
0060  * \brief Create/open the root and txt summary files, where the DMR histograms and the associtated mean/sigma are stored respectively
0061  */
0062 
0063 void PlotAlignmentValidation::openSummaryFile() {
0064   if (!openedsummaryfile) {
0065     openedsummaryfile = true;
0066     summaryfile.open(outputDir + "/" + summaryfilename + ".txt");
0067     //Rootfile introduced to store the DMR histograms
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     //Check for the rootfile to be open, and open it in case it is not already.
0076     if (!rootsummaryfile->IsOpen())
0077       rootsummaryfile->Open(outputDir + "/" + summaryfilename + ".root", "UPDATE");
0078   }
0079 }
0080 
0081 //------------------------------------------------------------------------------
0082 /*! \fn loadFileList
0083  *  \brief Add to the list of sources the rootfile associated to a particular geometry
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 /*! \fn useFitForDMRplots
0095  *  \brief Store the selected boolean in one of the private members of the class
0096  */
0097 void PlotAlignmentValidation::useFitForDMRplots(bool usefit) { useFit_ = usefit; }
0098 
0099 //------------------------------------------------------------------------------
0100 /*! \fn numberOfLayers
0101  *  \brief Select the number of layers associated to a subdetector.
0102  */
0103 //TODO Possible improvement: reduce the number of switches in the code by implementing a map
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 /*! \fn maxNumberOfLayers
0148  *  \brief Return the number of layers of a subdetector
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 /*! \fn legendOptions
0160  *  \brief Assign legend options to members of the class
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 /*! \fn setOutputDir
0188  *  \brief Set the output direcotry
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 /*! \fn plotSubDetResiduals
0201  *  \brief Function used to plot residuals for a subdetector
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();  //FindObjectAny(histoName.Data());
0249       sumHisto->SetLineColor(tmpcounter + 1);
0250       sumHisto->SetLineStyle(tmpcounter + 1);
0251       sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
0252       sumHisto->Draw();
0253 
0254       //get statistic box coordinate to plot all boxes one below the other
0255       //gStyle->SetStatY(0.91);
0256       //gStyle->SetStatW(0.15);
0257       //gStyle->SetStatBorderSize(1);
0258       //gStyle->SetStatH(0.10);
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       //hstack->Add(sumHisto);
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   //hstack->Draw("nostack");
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   //delete c;
0295   //c=0;
0296 }
0297 
0298 //------------------------------------------------------------------------------
0299 void PlotAlignmentValidation::plotHitMaps() {
0300   //gStyle->SetOptStat(0);
0301 
0302   TCanvas* c = new TCanvas("c", "c", 1200, 400);
0303   c->Divide(3, 1);
0304   //ps->NewPage();
0305 
0306   //-------------------------------------------------
0307   //plot Hit map
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   //   //c->Update();
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   //TList* treelist=getTreeList();
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;  // ROOT will initilise
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     //   std::cout<<"treeMem->entries  "<<treeMem->entries<<std::endl;
0375     //  std::cout<<"var                  "<<var<<std::endl;
0376     //  std::cout<<"plotVariable_cut     "<<plotVariable_cut<<std::endl;
0377 
0378     if (var > plotVariable_cut && treeMem->entries > minHits) {
0379       TFile* f = (*sourceList.begin())->getFile();  //(TFile*)sourcelist->First();
0380 
0381       if (f->FindKeyAny(treeMem->histNameX.c_str()) != nullptr) {
0382         TH1* h =
0383             (TH1*)f->FindKeyAny(treeMem->histNameX.c_str())->ReadObj();  //f->FindObjectAny(treeMem->histNameX.c_str());
0384         gStyle->SetOptFit(0111);
0385         std::cout << "hist name " << h->GetName() << std::endl;
0386 
0387         TString path = (char*)strstr(gDirectory->GetPath(), "TrackerOfflineValidation");
0388         //std::cout<<"hist path "<<path<<std::endl;
0389         //std::cout<<"wrote text "<<std::endl;
0390         if (h)
0391           std::cout << h->GetEntries() << std::endl;
0392 
0393         //modules' location as title
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         //residual histogram
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         //norm. residual histogram
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   //read the number of entries in the t3
0432   //TTree* tree=0;
0433   //tree=(TTree*)treeList->At(0);
0434 
0435   //c1->Close();
0436 }
0437 
0438 //------------------------------------------------------------------------------
0439 /*! \fn getTreeList
0440  *  \brief Extract from the rootfiles stored in the sourcelist the corresponding trees.
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   // else if (options == "fine") ...
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);  //regardless of style settings, we want x error bars here
0496 
0497   int plotLayerN = 0;
0498   //  int plotRingN  = 0;
0499   //  bool plotPlain = false;
0500   bool plotLayers = false;  // overrides plotLayerN
0501   //  bool plotRings  = false;  // Todo: implement this?
0502   int plotSubDetN = 0;  // if zero, plot all
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   // todo: title, min/max, nbins?
0529 
0530   // Loop over detectors
0531   for (int iSubDet = 1; iSubDet <= 6; ++iSubDet) {
0532     // TEC requires special care since rings 1-4 and 5-7 are plotted separately
0533     bool isTEC = (iSubDet == 6);
0534 
0535     // if subdet is specified, skip other subdets
0536     if (plotSubDetN != 0 && iSubDet != plotSubDetN)
0537       continue;
0538 
0539     // Skips plotting too high layers
0540     // if it's a mixture of phase 0 and 1, the phase 0 files will be skipped
0541     //  when plotting the higher layers of BPIX and FPIX
0542     if (plotLayerN > maxNumberOfLayers(iSubDet)) {
0543       continue;
0544     }
0545 
0546     int minlayer = plotLayers ? 1 : plotLayerN;
0547     int maxlayer = plotLayers ? maxNumberOfLayers(iSubDet) : plotLayerN;
0548     // see later where this is used
0549     int maxlayerphase0 = plotLayers ? numberOfLayers(0, iSubDet) : plotLayerN;
0550 
0551     for (int layer = minlayer; layer <= maxlayer; layer++) {
0552       // two plots for TEC, skip first
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 {          // TEC
0564           if (iTEC == 0)  // rings
0565             sprintf(selection, "subDetId==%d && ring <= 4", iSubDet);
0566           else
0567             sprintf(selection, "subDetId==%d && ring > 4", iSubDet);
0568         }
0569 
0570         // Title for plot and name for the file
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           // TEC and TID have discs, the rest have layers
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         // Generate histograms with selection
0611         TLegend* legend = nullptr;
0612         // Any file from phase 0 will be skipped if the last argument is false
0613         THStack* hs = addHists(selection, residType, &legend, false, /*validforphase0 = */ 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           // Adjust Labels
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           // micrometers:
0637           yName.ReplaceAll("cm", "#mum");
0638           hs->GetHistogram()->GetYaxis()->SetTitle(yName);
0639         }
0640 
0641         // Save plot to file
0642         std::ostringstream plotName;
0643         plotName << outputDir << "/SurfaceShape_" << subDetName << "_";
0644         plotName << residType;
0645         if (layer != 0) {
0646           plotName << "_";
0647           // TEC and TID have discs, the rest have layers
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         // PNG,EPS,PDF files
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         // ROOT file
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 /*! \fn plotDMR
0685  *  \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.
0686  */
0687 void PlotAlignmentValidation::plotDMR(const std::string& variable,
0688                                       Int_t minHits,
0689                                       const std::string& options,
0690                                       const std::string& filterName) {
0691   // If several, comma-separated values are given in 'variable',
0692   // call plotDMR with each value separately.
0693   // If a comma is found, the string is divided to two.
0694   // (no space allowed)
0695   std::size_t findres = variable.find(',');
0696   if (findres != std::string::npos) {
0697     std::string substring1 = variable.substr(0, findres);
0698     std::string substring2 = variable.substr(findres + 1, std::string::npos);
0699     plotDMR(substring1, minHits, options, filterName);
0700     plotDMR(substring2, minHits, options, filterName);
0701     return;
0702   }
0703 
0704   // Variable name should end with X or Y. If it doesn't, recursively calls plotDMR twice with
0705   // X and Y added, respectively
0706   if (variable == "mean" || variable == "median" || variable == "meanNorm" || variable == "rms" ||
0707       variable == "rmsNorm") {
0708     plotDMR(variable + "X", minHits, options, filterName);
0709     plotDMR(variable + "Y", minHits, options, filterName);
0710     return;
0711   }
0712 
0713   // options:
0714   // -plain (default, the whole distribution)
0715   // -split (distribution splitted to two)
0716   // -layers (plain db for each layer/disc superimposed in one plot)
0717   // -layersSeparate (plain db for each layer/disc in separate plots)
0718   // -layersSplit (splitted db for each layers/disc in one plot)
0719   // -layersSplitSeparate (splitted db, for each layers/disc in separate plots)
0720 
0721   TRegexp layer_re("layer=[0-9]+");
0722   bool plotPlain = false, plotSplits = false, plotLayers = false;
0723   int plotLayerN = 0;
0724   Ssiz_t index, len;
0725   if (options.find("plain") != std::string::npos) {
0726     plotPlain = true;
0727   }
0728   if (options.find("split") != std::string::npos) {
0729     plotSplits = true;
0730   }
0731   if (options.find("layers") != std::string::npos) {
0732     plotLayers = true;
0733   }
0734   if ((index = layer_re.Index(options, &len)) != -1) {
0735     if (plotLayers) {
0736       std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
0737     } else {
0738       std::string substr = options.substr(index + 6, len - 6);
0739       plotLayerN = atoi(substr.c_str());
0740     }
0741   }
0742 
0743   // Defaults to plotting only plain plot if empty (or invalid)
0744   // option string is given
0745   if (!plotPlain && !plotSplits) {
0746     plotPlain = true;
0747   }
0748 
0749   // This boolean array tells for which detector modules to plot split DMR plots
0750   // They are plotted for BPIX, FPIX, TIB and TOB
0751   static bool plotSplitsFor[6] = {true, true, true, false, true, false};
0752 
0753   DMRPlotInfo plotinfo;
0754 
0755   gStyle->SetOptStat(0);
0756 
0757   TCanvas c("canv", "canv");
0758 
0759   plotinfo.variable = variable;
0760   plotinfo.minHits = minHits;
0761   plotinfo.plotPlain = plotPlain;
0762   plotinfo.plotLayers = plotLayers;
0763   plotinfo.filterName = filterName;
0764 
0765   // width in cm
0766   // for DMRS, use 100 bins in range +-10 um, bin width 0.2um
0767   // if modified, check also TrackerOfflineValidationSummary_cfi.py and TrackerOfflineValidation_Standalone_cff.py
0768   if (variable == "meanX") {
0769     plotinfo.nbins = 50;
0770     plotinfo.min = -0.001;
0771     plotinfo.max = 0.001;
0772   } else if (variable == "meanY") {
0773     plotinfo.nbins = 50;
0774     plotinfo.min = -0.005;
0775     plotinfo.max = 0.005;
0776   } else if (variable == "medianX")
0777     if (plotSplits) {
0778       plotinfo.nbins = 50;
0779       plotinfo.min = -0.0005;
0780       plotinfo.max = 0.0005;
0781     } else {
0782       plotinfo.nbins = 100;
0783       plotinfo.min = -0.001;
0784       plotinfo.max = 0.001;
0785     }
0786   else if (variable == "medianY")
0787     if (plotSplits) {
0788       plotinfo.nbins = 50;
0789       plotinfo.min = -0.0005;
0790       plotinfo.max = 0.0005;
0791     } else {
0792       plotinfo.nbins = 100;
0793       plotinfo.min = -0.001;
0794       plotinfo.max = 0.001;
0795     }
0796   else if (variable == "meanNormX") {
0797     plotinfo.nbins = 100;
0798     plotinfo.min = -2.0;
0799     plotinfo.max = 2.0;
0800   } else if (variable == "meanNormY") {
0801     plotinfo.nbins = 100;
0802     plotinfo.min = -2.0;
0803     plotinfo.max = 2.0;
0804   } else if (variable == "rmsX") {
0805     plotinfo.nbins = 100;
0806     plotinfo.min = 0.0;
0807     plotinfo.max = 0.1;
0808   } else if (variable == "rmsY") {
0809     plotinfo.nbins = 100;
0810     plotinfo.min = 0.0;
0811     plotinfo.max = 0.1;
0812   } else if (variable == "rmsNormX") {
0813     plotinfo.nbins = 100;
0814     plotinfo.min = 0.3;
0815     plotinfo.max = 1.8;
0816   } else if (variable == "rmsNormY") {
0817     plotinfo.nbins = 100;
0818     plotinfo.min = 0.3;
0819     plotinfo.max = 1.8;
0820   } else {
0821     std::cerr << "Unknown variable " << variable << std::endl;
0822     plotinfo.nbins = 100;
0823     plotinfo.min = -0.1;
0824     plotinfo.max = 0.1;
0825   }
0826   //Begin loop on structures
0827   for (int i = 1; i <= 6; ++i) {
0828     // Skip strip detectors if plotting any "Y" variable
0829     if (i != 1 && i != 2 && variable.length() > 0 && variable[variable.length() - 1] == 'Y') {
0830       continue;
0831     }
0832 
0833     // Skips plotting too high layers
0834     if (plotLayerN > maxNumberOfLayers(i)) {
0835       continue;
0836     }
0837 
0838     plotinfo.plotSplits = plotSplits && plotSplitsFor[i - 1];
0839     if (!plotinfo.plotPlain && !plotinfo.plotSplits) {
0840       continue;
0841     }
0842 
0843     // Sets dimension of legend according to the number of plots
0844 
0845     bool hasheader = (TkAlStyle::legendheader != "");
0846 
0847     int nPlots = 1;
0848     if (plotinfo.plotSplits) {
0849       nPlots = 3;
0850     }
0851     // This will make the legend a bit bigger than necessary if there is a mixture of phase 0 and phase 1.
0852     // Not worth it to implement more complicated logic.
0853     if (plotinfo.plotLayers) {
0854       nPlots *= maxNumberOfLayers(i);
0855     }
0856     nPlots *= sourceList.size();
0857     if (twolines_) {
0858       nPlots *= 2;
0859     }
0860     nPlots += hasheader;
0861 
0862     double legendY = 0.80;
0863     if (nPlots > 3) {
0864       legendY -= 0.01 * (nPlots - 3);
0865     }
0866     if (bigtext_) {
0867       legendY -= 0.05;
0868     }
0869     if (legendY < 0.6) {
0870       std::cerr << "Warning: Huge legend!" << std::endl;
0871       legendY = 0.6;
0872     }
0873 
0874     THStack hstack("hstack", "hstack");
0875     plotinfo.maxY = 0;
0876     plotinfo.subDetId = i;
0877     plotinfo.legend = new TLegend(0.17, legendY, 0.85, 0.88);
0878     plotinfo.legend->SetNColumns(2);
0879     if (hasheader)
0880       plotinfo.legend->SetHeader(TkAlStyle::legendheader);
0881     if (bigtext_)
0882       plotinfo.legend->SetTextSize(TkAlStyle::textSize);
0883     plotinfo.legend->SetFillStyle(0);
0884     plotinfo.hstack = &hstack;
0885     plotinfo.h = plotinfo.h1 = plotinfo.h2 = nullptr;
0886     plotinfo.firsthisto = true;
0887 
0888     openSummaryFile();
0889     vmean.clear();
0890     vrms.clear();
0891     vdeltamean.clear();
0892     vmeanerror.clear();
0893     vPValueEqualSplitMeans.clear(), vAlignmentUncertainty.clear();
0894     vPValueMeanEqualIdeal.clear();
0895     vPValueRMSEqualIdeal.clear();
0896 
0897     std::string stringsubdet;
0898     switch (i) {
0899       case 1:
0900         stringsubdet = "BPIX";
0901         break;
0902       case 2:
0903         stringsubdet = "FPIX";
0904         break;
0905       case 3:
0906         stringsubdet = "TIB";
0907         break;
0908       case 4:
0909         stringsubdet = "TID";
0910         break;
0911       case 5:
0912         stringsubdet = "TOB";
0913         break;
0914       case 6:
0915         stringsubdet = "TEC";
0916         break;
0917     }
0918 
0919     for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0920       plotinfo.vars = *it;
0921       plotinfo.h1 = plotinfo.h2 = plotinfo.h = nullptr;
0922 
0923       int minlayer = plotLayers ? 1 : plotLayerN;
0924       //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.
0925       if (plotinfo.plotPlain)
0926         minlayer = 0;
0927       int maxlayer = plotLayers ? numberOfLayers(plotinfo.vars->getPhase(), plotinfo.subDetId) : plotLayerN;
0928 
0929       for (int layer = minlayer; layer <= maxlayer; layer++) {
0930         if (plotinfo.plotPlain) {
0931           plotDMRHistogram(plotinfo, 0, layer, stringsubdet);
0932         }
0933 
0934         if (plotinfo.plotSplits) {
0935           plotDMRHistogram(plotinfo, -1, layer, stringsubdet);
0936           plotDMRHistogram(plotinfo, 1, layer, stringsubdet);
0937         }
0938 
0939         if (plotinfo.plotPlain) {
0940           if (plotinfo.h) {
0941             setDMRHistStyleAndLegend(plotinfo.h, plotinfo, 0, layer);
0942           } else {
0943             if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") &&
0944                 /*!plotinfo.plotLayers && */ layer == 0) {
0945               vmean.push_back(nan(""));
0946               vrms.push_back(nan(""));
0947               vmeanerror.push_back(nan(""));
0948               vAlignmentUncertainty.push_back(nan(""));
0949               vPValueMeanEqualIdeal.push_back(nan(""));
0950               vPValueRMSEqualIdeal.push_back(nan(""));
0951               if (plotinfo.plotSplits && plotinfo.plotPlain) {
0952                 vdeltamean.push_back(nan(""));
0953                 vPValueEqualSplitMeans.push_back(nan(""));
0954               }
0955             }
0956           }
0957         }
0958 
0959         if (plotinfo.plotSplits) {
0960           // Add delta mu to the histogram
0961           if (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && !plotinfo.plotPlain) {
0962             std::ostringstream legend;
0963             std::string unit = " #mum";
0964             legend.precision(3);
0965             legend << std::fixed;  // to always show 3 decimals
0966             float factor = 10000.0f;
0967             if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" ||
0968                 plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY") {
0969               factor = 1.0f;
0970               unit = "";
0971             }
0972             float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
0973             legend << plotinfo.vars->getName();
0974             if (layer > 0) {
0975               // TEC and TID have discs, the rest have layers
0976               if (i == 4 || i == 6)
0977                 legend << ", disc ";
0978               else
0979                 legend << ", layer ";
0980               legend << layer;
0981             }
0982             plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
0983             legend.str("");
0984             legend << "#Delta#mu = " << deltamu << unit;
0985             plotinfo.legend->AddEntry(static_cast<TObject*>(nullptr), legend.str().c_str(), "");
0986 
0987             if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && !plotLayers && layer == 0) {
0988               vdeltamean.push_back(deltamu);
0989             }
0990           }
0991           if (plotinfo.h1) {
0992             setDMRHistStyleAndLegend(plotinfo.h1, plotinfo, -1, layer);
0993           }
0994           if (plotinfo.h2) {
0995             setDMRHistStyleAndLegend(plotinfo.h2, plotinfo, 1, layer);
0996           }
0997         }
0998       }
0999     }
1000 
1001     if (hstack.GetHists() != nullptr && hstack.GetHists()->GetSize() != 0) {
1002       hstack.Draw("nostack");
1003       hstack.SetMaximum(plotinfo.maxY * 1.3);
1004       setTitleStyle(hstack, variable.c_str(), "#modules", plotinfo.subDetId);
1005       setHistStyle(*hstack.GetHistogram(), variable.c_str(), "#modules", 1);
1006 
1007       plotinfo.legend->Draw();
1008     } else {
1009       // Draw an empty default histogram
1010       plotinfo.h = new TH1F("defhist", "Empty default histogram", plotinfo.nbins, plotinfo.min, plotinfo.max);
1011       plotinfo.h->SetMaximum(10);
1012       if (plotinfo.variable.find("Norm") == std::string::npos)
1013         scaleXaxis(plotinfo.h, 10000);
1014       setTitleStyle(*plotinfo.h, variable.c_str(), "#modules", plotinfo.subDetId);
1015       setHistStyle(*plotinfo.h, variable.c_str(), "#modules", 1);
1016       plotinfo.h->Draw();
1017     }
1018 
1019     std::ostringstream plotName;
1020     plotName << outputDir << "/D";
1021 
1022     if (variable == "medianX")
1023       plotName << "medianR_";
1024     else if (variable == "medianY")
1025       plotName << "medianYR_";
1026     else if (variable == "meanX")
1027       plotName << "meanR_";
1028     else if (variable == "meanY")
1029       plotName << "meanYR_";
1030     else if (variable == "meanNormX")
1031       plotName << "meanNR_";
1032     else if (variable == "meanNormY")
1033       plotName << "meanNYR_";
1034     else if (variable == "rmsX")
1035       plotName << "rmsR_";
1036     else if (variable == "rmsY")
1037       plotName << "rmsYR_";
1038     else if (variable == "rmsNormX")
1039       plotName << "rmsNR_";
1040     else if (variable == "rmsNormY")
1041       plotName << "rmsNYR_";
1042 
1043     TString subdet;
1044     switch (i) {
1045       case 1:
1046         subdet = "BPIX";
1047         break;
1048       case 2:
1049         subdet = "FPIX";
1050         break;
1051       case 3:
1052         subdet = "TIB";
1053         break;
1054       case 4:
1055         subdet = "TID";
1056         break;
1057       case 5:
1058         subdet = "TOB";
1059         break;
1060       case 6:
1061         subdet = "TEC";
1062         break;
1063     }
1064 
1065     plotName << subdet;
1066 
1067     if (plotPlain && !plotSplits) {
1068       plotName << "_plain";
1069     } else if (!plotPlain && plotSplits) {
1070       plotName << "_split";
1071     }
1072     if (plotLayers) {
1073       // TEC and TID have discs, the rest have layers
1074       if (i == 4 || i == 6)
1075         plotName << "_discs";
1076       else
1077         plotName << "_layers";
1078     }
1079     if (plotLayerN > 0) {
1080       // TEC and TID have discs, the rest have layers
1081       if (i == 4 || i == 6)
1082         plotName << "_disc";
1083       else
1084         plotName << "_layer";
1085       plotName << plotLayerN;
1086     }
1087 
1088     // PNG,EPS,PDF files
1089     c.Update();
1090     c.Print((plotName.str() + ".png").c_str());
1091     c.Print((plotName.str() + ".eps").c_str());
1092     c.Print((plotName.str() + ".pdf").c_str());
1093 
1094     // ROOT file
1095     TFile f((plotName.str() + ".root").c_str(), "recreate");
1096     c.Write();
1097     f.Close();
1098 
1099     // Free allocated memory.
1100     delete plotinfo.h;
1101     delete plotinfo.h1;
1102     delete plotinfo.h2;
1103 
1104     if (!vmean.empty()) {
1105       summaryfile << "   mu_" << subdet;
1106       if (plotinfo.variable == "medianY")
1107         summaryfile << "_y";
1108       summaryfile << " (um)\t"
1109                   << "latexname=$\\mu_\\text{" << subdet << "}";
1110       if (plotinfo.variable == "medianY")
1111         summaryfile << "^{y}";
1112       summaryfile << "$ ($\\mu$m)\t"
1113                   << "format={:.3g}\t"
1114                   << "latexformat=${:.3g}$";
1115       for (auto mu : vmean)
1116         summaryfile << "\t" << mu;
1117       summaryfile << "\n";
1118     }
1119     if (!vrms.empty()) {
1120       summaryfile << "sigma_" << subdet;
1121       if (plotinfo.variable == "medianY")
1122         summaryfile << "_y";
1123       summaryfile << " (um)\t"
1124                   << "latexname=$\\sigma_\\text{" << subdet << "}";
1125       if (plotinfo.variable == "medianY")
1126         summaryfile << "^{y}";
1127       summaryfile << "$ ($\\mu$m)\t"
1128                   << "format={:.3g}\t"
1129                   << "latexformat=${:.3g}$";
1130       for (auto sigma : vrms)
1131         summaryfile << "\t" << sigma;
1132       summaryfile << "\n";
1133     }
1134     if (!vdeltamean.empty()) {
1135       summaryfile << "  dmu_" << subdet;
1136       if (plotinfo.variable == "medianY")
1137         summaryfile << "_y";
1138       summaryfile << " (um)\t"
1139                   << "latexname=$\\Delta\\mu_\\text{" << subdet << "}";
1140       if (plotinfo.variable == "medianY")
1141         summaryfile << "^{y}";
1142       summaryfile << "$ ($\\mu$m)\t"
1143                   << "format={:.3g}\t"
1144                   << "latexformat=${:.3g}$";
1145       for (auto dmu : vdeltamean)
1146         summaryfile << "\t" << dmu;
1147       summaryfile << "\n";
1148     }
1149     if (!vmeanerror.empty()) {
1150       summaryfile << "  sigma_mu_" << subdet;
1151       if (plotinfo.variable == "medianY")
1152         summaryfile << "_y";
1153       summaryfile << " (um)\t"
1154                   << "latexname=$\\sigma\\mu_\\text{" << subdet << "}";
1155       if (plotinfo.variable == "medianY")
1156         summaryfile << "^{y}";
1157       summaryfile << "$ ($\\mu$m)\t"
1158                   << "format={:.3g}\t"
1159                   << "latexformat=${:.3g}$";
1160       for (auto dmu : vmeanerror)
1161         summaryfile << "\t" << dmu;
1162       summaryfile << "\n";
1163     }
1164     if (!vPValueEqualSplitMeans.empty()) {
1165       summaryfile << "  p_delta_mu_equal_zero_" << subdet;
1166       if (plotinfo.variable == "medianY")
1167         summaryfile << "_y";
1168       summaryfile << "\t"
1169                   << "latexname=$P(\\delta\\mu_\\text{" << subdet << "}=0)";
1170       if (plotinfo.variable == "medianY")
1171         summaryfile << "^{y}";
1172       summaryfile << "$\t"
1173                   << "format={:.3g}\t"
1174                   << "latexformat=${:.3g}$";
1175       for (auto dmu : vPValueEqualSplitMeans)
1176         summaryfile << "\t" << dmu;
1177       summaryfile << "\n";
1178     }
1179     if (!vAlignmentUncertainty.empty()) {
1180       summaryfile << "  alignment_uncertainty_" << subdet;
1181       if (plotinfo.variable == "medianY")
1182         summaryfile << "_y";
1183       summaryfile << " (um)\t"
1184                   << "latexname=$\\sigma_\\text{align}_\\text{" << subdet << "}";
1185       if (plotinfo.variable == "medianY")
1186         summaryfile << "^{y}";
1187       summaryfile << "$ ($\\mu$m)\t"
1188                   << "format={:.3g}\t"
1189                   << "latexformat=${:.3g}$";
1190       for (auto dmu : vAlignmentUncertainty)
1191         summaryfile << "\t" << dmu;
1192       summaryfile << "\n";
1193     }
1194     if (!vPValueMeanEqualIdeal.empty()) {
1195       summaryfile << "  p_mean_equals_ideal_" << subdet;
1196       if (plotinfo.variable == "medianY")
1197         summaryfile << "_y";
1198       summaryfile << "\t"
1199                   << "latexname=$P(\\mu_\\text{" << subdet << "}=\\mu_\\text{ideal})";
1200       if (plotinfo.variable == "medianY")
1201         summaryfile << "^{y}";
1202       summaryfile << "$\t"
1203                   << "format={:.3g}\t"
1204                   << "latexformat=${:.3g}$";
1205       for (auto dmu : vPValueMeanEqualIdeal)
1206         summaryfile << "\t" << dmu;
1207       summaryfile << "\n";
1208     }
1209     if (!vPValueRMSEqualIdeal.empty()) {
1210       summaryfile << "  p_RMS_equals_ideal_" << subdet;
1211       if (plotinfo.variable == "medianY")
1212         summaryfile << "_y";
1213       summaryfile << "\t"
1214                   << "latexname=$P(\\sigma_\\text{" << subdet << "}=\\sigma_\\text{ideal})";
1215       if (plotinfo.variable == "medianY")
1216         summaryfile << "^{y}";
1217       summaryfile << "$\t"
1218                   << "format={:.3g}\t"
1219                   << "latexformat=${:.3g}$";
1220       for (auto dmu : vPValueRMSEqualIdeal)
1221         summaryfile << "\t" << dmu;
1222       summaryfile << "\n";
1223     }
1224   }
1225 }
1226 
1227 //------------------------------------------------------------------------------
1228 void PlotAlignmentValidation::plotChi2(const char* inputFile) {
1229   // Opens the file (it should be OfflineValidation(Parallel)_result.root)
1230   // and reads and plots the norm_chi^2 and h_chi2Prob -distributions.
1231 
1232   Bool_t errorflag = kFALSE;
1233   TFile* fi1 = TFile::Open(inputFile, "read");
1234   if (fi1 != nullptr) {
1235     if (fi1->GetDirectory("TrackerOfflineValidation/GlobalTrackVariables") == nullptr) {
1236       errorflag = kTRUE;
1237     }
1238   } else {
1239     errorflag = kTRUE;
1240   }
1241   if (errorflag) {
1242     std::cout << "PlotAlignmentValidation::plotChi2: Can't find GlobalTrackVariables given file,"
1243               << " no chi^2-plots produced" << std::endl;
1244     return;
1245   }
1246 
1247   auto normchi = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_normchi2");
1248   // please remove following line once bug is fixed with ROOT Version: >6.27/01
1249   normchi->GetFrame()->ResetBit(kCanDelete);
1250 
1251   auto chiprob = fi1->Get<TCanvas>("TrackerOfflineValidation/GlobalTrackVariables/h_chi2Prob");
1252   // please remove following line once bug is fixed with ROOT Version: >6.27/01
1253   chiprob->GetFrame()->ResetBit(kCanDelete);
1254 
1255   if (normchi == nullptr || chiprob == nullptr) {
1256     errorflag = kTRUE;
1257   }
1258   if (errorflag) {
1259     std::cout << "PlotAlignmentValidation::plotChi2: Can't find data from given file,"
1260               << " no chi^2-plots produced" << std::endl;
1261     return;
1262   }
1263 
1264   TLegend* legend = nullptr;
1265   for (auto primitive : *normchi->GetListOfPrimitives()) {
1266     legend = dynamic_cast<TLegend*>(primitive);
1267     if (legend)
1268       break;
1269   }
1270   if (legend) {
1271     openSummaryFile();
1272     summaryfile << "ntracks";
1273     for (auto alignment : sourceList) {
1274       summaryfile << "\t";
1275       TString title = alignment->getName();
1276       int color = alignment->getLineColor();
1277       int style = alignment->getLineStyle();
1278       bool foundit = false;
1279       for (auto entry : *legend->GetListOfPrimitives()) {
1280         TLegendEntry* legendentry = dynamic_cast<TLegendEntry*>(entry);
1281         assert(legendentry);
1282         TH1* h = dynamic_cast<TH1*>(legendentry->GetObject());
1283         if (!h)
1284           continue;
1285         if (legendentry->GetLabel() == title && h->GetLineColor() == color && h->GetLineStyle() == style) {
1286           foundit = true;
1287           summaryfile << h->GetEntries();
1288           break;
1289         }
1290       }
1291       if (!foundit) {
1292         summaryfile << 0;
1293       }
1294     }
1295     summaryfile << "\n";
1296   }
1297 
1298   chiprob->Draw();
1299   normchi->Draw();
1300 
1301   // PNG,EPS,PDF files
1302   normchi->Print((outputDir + "/h_normchi2.png").c_str());
1303   chiprob->Print((outputDir + "/h_chi2Prob.png").c_str());
1304   normchi->Print((outputDir + "/h_normchi2.eps").c_str());
1305   chiprob->Print((outputDir + "/h_chi2Prob.eps").c_str());
1306   normchi->Print((outputDir + "/h_normchi2.pdf").c_str());
1307   chiprob->Print((outputDir + "/h_chi2Prob.pdf").c_str());
1308 
1309   // ROOT files
1310   TFile fi2((outputDir + "/h_normchi2.root").c_str(), "recreate");
1311   normchi->Write();
1312   fi2.Close();
1313 
1314   TFile fi3((outputDir + "/h_chi2Prob.root").c_str(), "recreate");
1315   chiprob->Write();
1316   fi3.Close();
1317 
1318   fi1->Close();
1319   delete fi1;
1320 }
1321 
1322 //------------------------------------------------------------------------------
1323 THStack* PlotAlignmentValidation::addHists(
1324     const TString& selection, const TString& residType, TLegend** myLegend, bool printModuleIds, bool validforphase0) {
1325   enum ResidType {
1326     xPrimeRes,
1327     yPrimeRes,
1328     xPrimeNormRes,
1329     yPrimeNormRes,
1330     xRes,
1331     yRes,
1332     xNormRes, /*yResNorm*/
1333     ResXvsXProfile,
1334     ResXvsYProfile,
1335     ResYvsXProfile,
1336     ResYvsYProfile
1337   };
1338   ResidType rType = xPrimeRes;
1339   if (residType == "xPrime")
1340     rType = xPrimeRes;
1341   else if (residType == "yPrime")
1342     rType = yPrimeRes;
1343   else if (residType == "xPrimeNorm")
1344     rType = xPrimeNormRes;
1345   else if (residType == "yPrimeNorm")
1346     rType = yPrimeNormRes;
1347   else if (residType == "x")
1348     rType = xRes;
1349   else if (residType == "y")
1350     rType = yRes;
1351   else if (residType == "xNorm")
1352     rType = xNormRes;
1353   // else if (residType == "yNorm") rType = yResNorm;
1354   else if (residType == "ResXvsXProfile")
1355     rType = ResXvsXProfile;
1356   else if (residType == "ResYvsXProfile")
1357     rType = ResYvsXProfile;
1358   else if (residType == "ResXvsYProfile")
1359     rType = ResXvsYProfile;
1360   else if (residType == "ResYvsYProfile")
1361     rType = ResYvsYProfile;
1362   else {
1363     std::cout << "PlotAlignmentValidation::addHists: Unknown residual type " << residType << std::endl;
1364     return nullptr;
1365   }
1366 
1367   std::cout << "PlotAlignmentValidation::addHists: using selection " << selection << std::endl;
1368   THStack* retHistoStack = new THStack("hstack", "");
1369   if (myLegend != nullptr)
1370     if (*myLegend == nullptr) {
1371       *myLegend = new TLegend(0.17, 0.80, 0.85, 0.88);
1372     }
1373 
1374   for (std::vector<TkOfflineVariables*>::iterator itSourceFile = sourceList.begin(); itSourceFile != sourceList.end();
1375        ++itSourceFile) {
1376     std::vector<TString> histnames;
1377 
1378     TFile* f = (*itSourceFile)->getFile();
1379     TTree* tree = (*itSourceFile)->getTree();
1380     int myLineColor = (*itSourceFile)->getLineColor();
1381     int myLineStyle = (*itSourceFile)->getLineStyle();
1382     TString myLegendName = (*itSourceFile)->getName();
1383     TH1* h = nullptr;   // becomes result
1384     UInt_t nEmpty = 0;  // selected, but empty hists
1385     Long64_t nentries = tree->GetEntriesFast();
1386     if (!f || !tree) {
1387       std::cout << "PlotAlignmentValidation::addHists: no tree or no file" << std::endl;
1388       return nullptr;
1389     }
1390 
1391     bool histnamesfilled = false;
1392     int phase = (bool)(f->Get("TrackerOfflineValidation/Pixel/P1PXBBarrel_1"));
1393     if (residType.Contains("Res") && residType.Contains("Profile")) {
1394       TString basename = TString(residType)
1395                              .ReplaceAll("Res", "p_res")
1396                              .ReplaceAll("vs", "")
1397                              .ReplaceAll("Profile", "_");  //gives e.g.: p_resXX_
1398       if (selection == "subDetId==1") {
1399         if (phase == 1)
1400           histnames.push_back(TString(basename) += "P1PXBBarrel_1");
1401         else
1402           histnames.push_back(TString(basename) += "TPBBarrel_1");
1403         histnamesfilled = true;
1404       } else if (selection == "subDetId==2") {
1405         if (phase == 1) {
1406           histnames.push_back(TString(basename) += "P1PXECEndcap_2");
1407           histnames.push_back(TString(basename) += "P1PXECEndcap_3");
1408         } else {
1409           histnames.push_back(TString(basename) += "TPEEndcap_2");
1410           histnames.push_back(TString(basename) += "TPEEndcap_3");
1411         }
1412         histnamesfilled = true;
1413       } else if (selection == "subDetId==3") {
1414         histnames.push_back(TString(basename) += "TIBBarrel_1");
1415         histnamesfilled = true;
1416       } else if (selection == "subDetId==4") {
1417         histnames.push_back(TString(basename) += "TIDEndcap_2");
1418         histnames.push_back(TString(basename) += "TIDEndcap_3");
1419         histnamesfilled = true;
1420       } else if (selection == "subDetId==5") {
1421         histnames.push_back(TString(basename) += "TOBBarrel_4");
1422         histnamesfilled = true;
1423       } else if (selection == "subDetId==6") {  //whole TEC - doesn't happen by default but easy enough to account for
1424         histnames.push_back(TString(basename) += "TECEndcap_5");
1425         histnames.push_back(TString(basename) += "TECEndcap_6");
1426         histnamesfilled = true;
1427       } else if (selection == "subDetId==6 && ring <= 4") {
1428         //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.
1429         for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1430           for (int iDisk = 1; iDisk <= 9; iDisk++)
1431             for (int iSide = 1; iSide <= 2; iSide++)
1432               for (int iPetal = 1; iPetal <= 8; iPetal++)
1433                 for (int iRing = 1; iRing <= 4 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9); iRing++)
1434                 //in the higher disks, the inner rings go away.  But the numbering in the file structure removes the higher numbers
1435                 // so the numbers there do not correspond to the actual ring numbers
1436                 {
1437                   std::stringstream s;
1438                   s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
1439                     << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1440                   histnames.push_back(TString(s.str()));
1441                 }
1442         histnamesfilled = true;
1443       } else if (selection == "subDetId==6 && ring > 4") {
1444         //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.
1445         for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1446           for (int iDisk = 1; iDisk <= 9; iDisk++)
1447             for (int iSide = 1; iSide <= 2; iSide++)
1448               for (int iPetal = 1; iPetal <= 8; iPetal++)
1449                 for (int iRing = 5 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1450                      iRing <= 7 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1451                      iRing++)
1452                 //in the higher disks, the inner rings go away.  But the numbering in the file structure removes the higher numbers
1453                 // so the numbers there do not correspond to the actual ring numbers
1454                 {
1455                   std::stringstream s;
1456                   s << "TrackerOfflineValidation/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk << "/TECSide_"
1457                     << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1458                   histnames.push_back(TString(s.str()));
1459                 }
1460         histnamesfilled = true;
1461       }
1462     }
1463 
1464     Long64_t nSel = 0;
1465     if (histnamesfilled && !histnames.empty()) {
1466       nSel = (Long64_t)histnames.size();
1467     }
1468     if (!histnamesfilled) {
1469       // first loop on tree to find out which entries (i.e. modules) fulfill the selection
1470       // 'Entry$' gives the entry number in the tree
1471       nSel = tree->Draw("Entry$", selection, "goff");
1472       if (nSel == -1)
1473         return nullptr;  // error in selection
1474       if (nSel == 0) {
1475         std::cout << "PlotAlignmentValidation::addHists: no selected module." << std::endl;
1476         return nullptr;
1477       }
1478       // copy entry numbers that fulfil the selection
1479       const std::vector<double> selected(tree->GetV1(), tree->GetV1() + nSel);
1480 
1481       std::vector<double>::const_iterator iterEnt = selected.begin();
1482 
1483       // second loop on tree:
1484       // for each selected entry get the hist from the file and merge
1485       TkOffTreeVariables* treeMem = nullptr;  // ROOT will initialise
1486       tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
1487       for (Long64_t i = 0; i < nentries; i++) {
1488         if (i < *iterEnt - 0.1               // smaller index (with tolerance): skip
1489             || iterEnt == selected.end()) {  // at the end: skip
1490           continue;
1491         } else if (TMath::Abs(i - *iterEnt) < 0.11) {
1492           ++iterEnt;  // take this entry!
1493         } else
1494           std::cout << "Must not happen: " << i << " " << *iterEnt << std::endl;
1495 
1496         tree->GetEntry(i);
1497         if (printModuleIds) {
1498           std::cout << treeMem->moduleId << ": " << treeMem->entries << " entries" << std::endl;
1499         }
1500         if (treeMem->entries <= 0) {  // little speed up: skip empty hists
1501           ++nEmpty;
1502           continue;
1503         }
1504         TString hName;
1505         switch (rType) {
1506           case xPrimeRes:
1507             hName = treeMem->histNameX.c_str();
1508             break;
1509           case yPrimeRes:
1510             hName = treeMem->histNameY.c_str();
1511             break;
1512           case xPrimeNormRes:
1513             hName = treeMem->histNameNormX.c_str();
1514             break;
1515           case yPrimeNormRes:
1516             hName = treeMem->histNameNormY.c_str();
1517             break;
1518           case xRes:
1519             hName = treeMem->histNameLocalX.c_str();
1520             break;
1521           case yRes:
1522             hName = treeMem->histNameLocalY.c_str();
1523             break;
1524           case xNormRes:
1525             hName = treeMem->histNameNormLocalX.c_str();
1526             break;
1527             /*case yResNorm:      hName = treeMem->histNameNormLocalY.c_str(); break;*/
1528           case ResXvsXProfile:
1529             hName = treeMem->profileNameResXvsX.c_str();
1530             break;
1531           case ResXvsYProfile:
1532             hName = treeMem->profileNameResXvsY.c_str();
1533             break;
1534           case ResYvsXProfile:
1535             hName = treeMem->profileNameResYvsX.c_str();
1536             break;
1537           case ResYvsYProfile:
1538             hName = treeMem->profileNameResYvsY.c_str();
1539             break;
1540         }
1541         histnames.push_back(hName);
1542       }
1543     }
1544 
1545     for (std::vector<TString>::iterator ithistname = histnames.begin(); ithistname != histnames.end(); ++ithistname) {
1546       if (phase == 0 && !validforphase0)
1547         break;
1548       TH1* newHist;
1549       if (ithistname->Contains("/")) {
1550         newHist = (TH1*)f->Get(*ithistname);
1551       } else {
1552         TKey* histKey = f->FindKeyAny(*ithistname);
1553         newHist = (histKey ? static_cast<TH1*>(histKey->ReadObj()) : nullptr);
1554       }
1555       if (!newHist) {
1556         std::cout << "Hist " << *ithistname << " not found in file, break loop." << std::endl;
1557         break;
1558       }
1559       if (newHist->GetEntries() == 0) {
1560         nEmpty++;
1561         continue;
1562       }
1563       newHist->SetLineColor(myLineColor);
1564       newHist->SetLineStyle(myLineStyle);
1565       if (!h) {  // first hist: clone, but rename keeping only first part of name
1566         TString name(newHist->GetName());
1567         Ssiz_t pos_ = 0;
1568         for (UInt_t i2 = 0; i2 < 3; ++i2)
1569           pos_ = name.Index("_", pos_ + 1);
1570         name = name(0, pos_);  // only up to three '_'
1571         h = static_cast<TH1*>(newHist->Clone("summed_" + name));
1572         //      TString myTitle = Form("%s: %lld modules", selection, nSel);
1573         //  h->SetTitle( myTitle );
1574       } else {  // otherwise just add
1575         h->Add(newHist);
1576       }
1577       delete newHist;
1578     }
1579 
1580     std::cout << "PlotAlignmentValidation::addHists"
1581               << "Result is merged from " << nSel - nEmpty << " hists, " << nEmpty << " hists were empty." << std::endl;
1582 
1583     if (nSel - nEmpty == 0)
1584       continue;
1585 
1586     if (myLegend != nullptr)
1587       (*myLegend)->AddEntry(h, myLegendName, "L");
1588 
1589     retHistoStack->Add(h);
1590   }
1591 
1592   return retHistoStack;
1593 }
1594 
1595 //------------------------------------------------------------------------------
1596 /*! \fn fitGauss
1597  *  \brief Operate a Gaussian fit to the given histogram
1598  */
1599 TF1* PlotAlignmentValidation::fitGauss(TH1* hist, int color) {
1600   //1. fits a Gauss function to the inner range of abs(2 rms)
1601   //2. repeates the Gauss fit in a 2 sigma range around mean of first fit
1602   //returns mean and sigma from fit in micron
1603   if (!hist || hist->GetEntries() < 20)
1604     return nullptr;
1605 
1606   float mean = hist->GetMean();
1607   float sigma = hist->GetRMS();
1608   std::string functionname = "gaussian_";
1609   functionname += hist->GetName();
1610   TF1* func = new TF1(functionname.c_str(), "gaus", mean - 2. * sigma, mean + 2. * sigma);
1611 
1612   func->SetLineColor(color);
1613   func->SetLineStyle(2);
1614   if (0 == hist->Fit(func, "QNR")) {  // N: do not blow up file by storing fit!
1615     mean = func->GetParameter(1);
1616     sigma = func->GetParameter(2);
1617     // second fit: three sigma of first fit around mean of first fit
1618     func->SetRange(mean - 3. * sigma, mean + 3. * sigma);
1619     // I: integral gives more correct results if binning is too wide
1620     // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1621     if (0 == hist->Fit(func, "Q0ILR")) {
1622       if (hist->GetFunction(func->GetName())) {  // Take care that it is later on drawn:
1623                                                  //hist->GetFunction(func->GetName())->ResetBit(TF1::kNotDraw);
1624       }
1625     }
1626   }
1627   return func;
1628 }
1629 
1630 //------------------------------------------------------------------------------
1631 /*! \fn storeHistogramInRootfile
1632  *  \brief Store the histogram and the gaussian function resulting from the fitGauss function into a root file
1633  */
1634 void PlotAlignmentValidation::storeHistogramInRootfile(TH1* hist) {
1635   //Store histogram and fit function in the root summary file
1636   rootsummaryfile->cd();
1637   hist->Write();
1638 }
1639 
1640 //------------------------------------------------------------------------------
1641 void PlotAlignmentValidation::scaleXaxis(TH1* hist, Int_t scale) {
1642   Double_t xmin = hist->GetXaxis()->GetXmin();
1643   Double_t xmax = hist->GetXaxis()->GetXmax();
1644   hist->GetXaxis()->SetLimits(xmin * scale, xmax * scale);
1645 }
1646 
1647 //------------------------------------------------------------------------------
1648 TObject* PlotAlignmentValidation::findObjectFromCanvas(TCanvas* canv, const char* className, Int_t n) {
1649   // Finds the n-th instance of the given class from the canvas
1650   TIter next(canv->GetListOfPrimitives());
1651   TObject* obj = nullptr;
1652   Int_t found = 0;
1653   while ((obj = next())) {
1654     if (strncmp(obj->ClassName(), className, 10) == 0) {
1655       if (++found == n)
1656         return obj;
1657     }
1658   }
1659 
1660   return nullptr;
1661 }
1662 
1663 //------------------------------------------------------------------------------
1664 void PlotAlignmentValidation::setTitleStyle(
1665     TNamed& hist, const char* titleX, const char* titleY, int subDetId, bool isSurfaceDeformation, TString secondline) {
1666   std::stringstream title_Xaxis;
1667   std::stringstream title_Yaxis;
1668   TString titleXAxis = titleX;
1669   TString titleYAxis = titleY;
1670   if (titleXAxis != "" && titleYAxis != "")
1671     std::cout << "plot " << titleXAxis << " vs " << titleYAxis << std::endl;
1672 
1673   hist.SetTitle("");
1674   TkAlStyle::drawStandardTitle();
1675 
1676   //Thanks Candice!
1677   TString subD;
1678   switch (subDetId) {
1679     case 1:
1680       subD = "BPIX";
1681       break;
1682     case 2:
1683       subD = "FPIX";
1684       break;
1685     case 3:
1686       subD = "TIB";
1687       break;
1688     case 4:
1689       subD = "TID";
1690       break;
1691     case 5:
1692       subD = "TOB";
1693       break;
1694     case 6:
1695       subD = "TEC";
1696       break;
1697   }
1698 
1699   TPaveText* text2;
1700   if (!isSurfaceDeformation) {
1701     text2 = new TPaveText(0.7, 0.3, 0.9, 0.6, "brNDC");
1702   } else {
1703     std::cout << "Surface Deformation" << std::endl;
1704     text2 = new TPaveText(0.8, 0.75, 0.9, 0.9, "brNDC");
1705   }
1706   text2->SetTextSize(0.06);
1707   text2->SetTextFont(42);
1708   text2->SetFillStyle(0);
1709   text2->SetBorderSize(0);
1710   text2->SetMargin(0.01);
1711   text2->SetTextAlign(12);  // align left
1712   text2->AddText(0.01, 0.75, subD);
1713   if (secondline != "") {
1714     text2->AddText(0.01, 0.25, secondline);
1715   }
1716   text2->Draw();
1717 }
1718 
1719 //------------------------------------------------------------------------------
1720 /*! \fn 
1721  *  \brief 
1722  */
1723 void PlotAlignmentValidation::setHistStyle(TH1& hist, const char* titleX, const char* titleY, int color) {
1724   std::stringstream title_Xaxis;
1725   std::stringstream title_Yaxis;
1726   TString titleXAxis = titleX;
1727   TString titleYAxis = titleY;
1728 
1729   if (titleXAxis.Contains("Phi"))
1730     title_Xaxis << titleX << "[rad]";
1731   else if (titleXAxis.Contains("meanX"))
1732     title_Xaxis << "#LTx'_{pred}-x'_{hit}#GT[#mum]";
1733   else if (titleXAxis.Contains("meanY"))
1734     title_Xaxis << "#LTy'_{pred}-y'_{hit}#GT[#mum]";
1735   else if (titleXAxis.Contains("rmsX"))
1736     title_Xaxis << "RMS(x'_{pred}-x'_{hit})[#mum]";
1737   else if (titleXAxis.Contains("rmsY"))
1738     title_Xaxis << "RMS(y'_{pred}-y'_{hit})[#mum]";
1739   else if (titleXAxis.Contains("meanNormX"))
1740     title_Xaxis << "#LTx'_{pred}-x'_{hit}/#sigma#GT";
1741   else if (titleXAxis.Contains("meanNormY"))
1742     title_Xaxis << "#LTy'_{pred}-y'_{hit}/#sigma#GT";
1743   else if (titleXAxis.Contains("rmsNormX"))
1744     title_Xaxis << "RMS(x'_{pred}-x'_{hit}/#sigma)";
1745   else if (titleXAxis.Contains("rmsNormY"))
1746     title_Xaxis << "RMS(y'_{pred}-y'_{hit}/#sigma)";
1747   else if (titleXAxis.Contains("meanLocalX"))
1748     title_Xaxis << "#LTx_{pred}-x_{hit}#GT[#mum]";
1749   else if (titleXAxis.Contains("rmsLocalX"))
1750     title_Xaxis << "RMS(x_{pred}-x_{hit})[#mum]";
1751   else if (titleXAxis.Contains("meanNormLocalX"))
1752     title_Xaxis << "#LTx_{pred}-x_{hit}/#sigma#GT[#mum]";
1753   else if (titleXAxis.Contains("rmsNormLocalX"))
1754     title_Xaxis << "RMS(x_{pred}-x_{hit}/#sigma)[#mum]";
1755   else if (titleXAxis.Contains("medianX"))
1756     title_Xaxis << "median(x'_{pred}-x'_{hit})[#mum]";
1757   else if (titleXAxis.Contains("medianY"))
1758     title_Xaxis << "median(y'_{pred}-y'_{hit})[#mum]";
1759   else
1760     title_Xaxis << titleX << "[cm]";
1761 
1762   if (hist.IsA()->InheritsFrom(TH1F::Class()))
1763     hist.SetLineColor(color);
1764   if (hist.IsA()->InheritsFrom(TProfile::Class())) {
1765     hist.SetMarkerStyle(20);
1766     hist.SetMarkerSize(0.8);
1767     hist.SetMarkerColor(color);
1768   }
1769 
1770   hist.GetXaxis()->SetTitle((title_Xaxis.str()).c_str());
1771 
1772   double binning = (hist.GetXaxis()->GetXmax() - hist.GetXaxis()->GetXmin()) / hist.GetNbinsX();
1773   title_Yaxis.precision(2);
1774 
1775   if (((titleYAxis.Contains("layer") || titleYAxis.Contains("ring")) && titleYAxis.Contains("subDetId")) ||
1776       titleYAxis.Contains("#modules")) {
1777     title_Yaxis << "number of modules";
1778     if (TString(title_Xaxis.str()).Contains("[#mum]"))
1779       title_Yaxis << " / " << binning << " #mum";
1780     else if (TString(title_Xaxis.str()).Contains("[cm]"))
1781       title_Yaxis << " / " << binning << " cm";
1782     else
1783       title_Yaxis << " / " << binning;
1784   } else
1785     title_Yaxis << titleY << "[cm]";
1786 
1787   hist.GetYaxis()->SetTitle((title_Yaxis.str()).c_str());
1788 
1789   hist.GetXaxis()->SetTitleFont(42);
1790   hist.GetYaxis()->SetTitleFont(42);
1791 }
1792 
1793 //------------------------------------------------------------------------------
1794 std::string PlotAlignmentValidation::getSelectionForDMRPlot(int minHits, int subDetId, int direction, int layer) {
1795   std::ostringstream builder;
1796   builder << "entries >= " << minHits;
1797   builder << " && subDetId == " << subDetId;
1798   if (direction != 0) {
1799     if (subDetId == 2) {  // FPIX is split by zDirection
1800       builder << " && zDirection == " << direction;
1801     } else {
1802       builder << " && rDirection == " << direction;
1803     }
1804   }
1805   if (layer > 0) {
1806     builder << " && layer == " << layer;
1807   }
1808   return builder.str();
1809 }
1810 
1811 std::string PlotAlignmentValidation::getVariableForDMRPlot(
1812     const std::string& histoname, const std::string& variable, int nbins, double min, double max) {
1813   std::ostringstream builder;
1814   builder << variable << ">>" << histoname << "(" << nbins << "," << min << "," << max << ")";
1815   return builder.str();
1816 }
1817 
1818 //------------------------------------------------------------------------------
1819 void PlotAlignmentValidation::setDMRHistStyleAndLegend(TH1F* h,
1820                                                        PlotAlignmentValidation::DMRPlotInfo& plotinfo,
1821                                                        int direction,
1822                                                        int layer) {
1823   TF1* fitResults = nullptr;
1824 
1825   h->SetDirectory(nullptr);
1826 
1827   // The whole DMR plot is plotted with wider line than the split plots
1828   // If only split plots are plotted, they will be stronger too, though
1829   h->SetLineWidth((direction == 0 || (plotinfo.plotSplits && !plotinfo.plotPlain)) ? 2 : 1);
1830 
1831   // These lines determine the style of the plots according to rules:
1832   // -If the plot is for direction != 0, +1 or +2 is added to the given style for distinction
1833   // -However if only direction split plots are to be plotted, the additions should be 0 and +1 respectively
1834   // -Modulo 4 arithmetic, because the styles run from 1..4
1835   int linestyle = plotinfo.vars->getLineStyle() - 1, linestyleplus = 0;
1836   if (direction == 1) {
1837     linestyleplus = 1;
1838   }
1839   if (direction == -1) {
1840     linestyleplus = 2;
1841   }
1842   if (direction != 0 && plotinfo.plotSplits && !plotinfo.plotPlain) {
1843     linestyleplus--;
1844   }
1845   linestyle = (linestyle + linestyleplus) % 4 + 1;
1846 
1847   int linecolor = plotinfo.vars->getLineColor();
1848   if (plotinfo.plotLayers && layer > 0) {
1849     linecolor += layer - 1;
1850   }
1851 
1852   if (plotinfo.firsthisto) {
1853     setHistStyle(*h, plotinfo.variable.c_str(), "#modules", 1);  //set color later
1854     plotinfo.firsthisto = false;
1855   }
1856 
1857   h->SetLineColor(linecolor);
1858   h->SetLineStyle(linestyle);
1859 
1860   if (plotinfo.maxY < h->GetMaximum()) {
1861     plotinfo.maxY = h->GetMaximum();
1862   }
1863 
1864   //fit histogram for median and mean
1865   if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1866       plotinfo.variable == "meanY") {
1867     fitResults = fitGauss(h, linecolor);
1868   }
1869 
1870   plotinfo.hstack->Add(h);
1871 
1872   std::ostringstream legend;
1873   legend.precision(3);
1874   legend << std::fixed;  // to always show 3 decimals
1875 
1876   // Legend: header part
1877   if (direction == -1 && plotinfo.subDetId != 2) {
1878     legend << "rDirection < 0";
1879   } else if (direction == 1 && plotinfo.subDetId != 2) {
1880     legend << "rDirection > 0";
1881   } else if (direction == -1 && plotinfo.subDetId == 2) {
1882     legend << "zDirection < 0";
1883   } else if (direction == 1 && plotinfo.subDetId == 2) {
1884     legend << "zDirection > 0";
1885   } else {
1886     legend << plotinfo.vars->getName();
1887     if (layer > 0) {
1888       // TEC and TID have discs, the rest have layers
1889       if (plotinfo.subDetId == 4 || plotinfo.subDetId == 6)
1890         legend << ", disc ";
1891       else
1892         legend << ", layer ";
1893       legend << layer << "";
1894     }
1895   }
1896 
1897   plotinfo.legend->AddEntry(h, legend.str().c_str(), "l");
1898   legend.str("");
1899 
1900   // Legend: Statistics
1901   double mean = 0.0, meanerror = 0.0, rms = 0.0, rmserror = 0.0;
1902   TString rmsname, units;
1903   bool showdeltamu =
1904       (plotinfo.h1 != nullptr && plotinfo.h2 != nullptr && plotinfo.plotSplits && plotinfo.plotPlain && direction == 0);
1905   if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1906       plotinfo.variable == "meanY" || plotinfo.variable == "rmsX" || plotinfo.variable == "rmsY") {
1907     if (useFit_ && fitResults) {
1908       mean = fitResults->GetParameter(1) * 10000;
1909       meanerror = fitResults->GetParError(1) * 10000;
1910       rms = fitResults->GetParameter(2) * 10000;
1911       rmserror = fitResults->GetParError(2) * 10000;
1912       rmsname = "#sigma";
1913       delete fitResults;
1914     } else {
1915       mean = h->GetMean(1) * 10000;
1916       meanerror = h->GetMeanError(1) * 10000;
1917       rms = h->GetRMS(1) * 10000;
1918       rmserror = h->GetRMSError(1) * 10000;
1919       rmsname = "rms";
1920     }
1921     units = " #mum";
1922   } else if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
1923              plotinfo.variable == "rmsNormY") {
1924     mean = h->GetMean(1);
1925     meanerror = h->GetMeanError(1);
1926     rms = h->GetRMS(1);
1927     rmserror = h->GetRMSError(1);
1928     rmsname = "rms";
1929     units = "";
1930   }
1931   if (showMean_) {
1932     legend << " #mu = " << mean;
1933     if (showMeanError_)
1934       legend << " #pm " << meanerror;
1935     legend << units;
1936     if (showRMS_ || showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1937       legend << ", ";
1938   }
1939   if (showRMS_) {
1940     legend << " " << rmsname << " = " << rms;
1941     if (showRMSError_)
1942       legend << " #pm " << rmserror;
1943     legend << units;
1944     if (showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1945       legend << ", ";
1946   }
1947 
1948   if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && /*!plotinfo.plotLayers && */ layer == 0 &&
1949       direction == 0) {
1950     vmean.push_back(mean);
1951     vrms.push_back(rms);
1952     vmeanerror.push_back(meanerror);
1953     TH1F* ideal = (TH1F*)plotinfo.hstack->GetHists()->At(0);
1954     TH1F* h = plotinfo.h;
1955     if (h->GetRMS() >= ideal->GetRMS()) {
1956       vAlignmentUncertainty.push_back(sqrt(pow(h->GetRMS(), 2) - pow(ideal->GetRMS(), 2)));
1957     } else {
1958       vAlignmentUncertainty.push_back(nan(""));
1959     }
1960     float p = (float)resampleTestOfEqualMeans(ideal, h, 10000);
1961     vPValueMeanEqualIdeal.push_back(p);
1962     p = resampleTestOfEqualRMS(ideal, h, 10000);
1963     vPValueRMSEqualIdeal.push_back(p);
1964   }
1965 
1966   // Legend: Delta mu for split plots
1967   if (showdeltamu) {
1968     float factor = 10000.0f;
1969     if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
1970         plotinfo.variable == "rmsNormY") {
1971       factor = 1.0f;
1972     }
1973     float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
1974     legend << "#Delta#mu = " << deltamu << units;
1975     if ((showModules_ || showUnderOverFlow_) && !twolines_)
1976       legend << ", ";
1977 
1978     if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && /*!plotinfo.plotLayers && */ layer == 0 &&
1979         direction == 0) {
1980       vdeltamean.push_back(deltamu);
1981       if (plotinfo.h1->GetEntries() && plotinfo.h2->GetEntries()) {
1982         float p = (float)resampleTestOfEqualMeans(plotinfo.h1, plotinfo.h2, 10000);
1983         vPValueEqualSplitMeans.push_back(p);
1984       }
1985     }
1986   }
1987 
1988   if (twolines_) {
1989     plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
1990     plotinfo.legend->AddEntry((TObject*)nullptr, "", "");
1991     legend.str("");
1992   }
1993 
1994   if (!showUnderOverFlow_ && showModules_) {
1995     legend << (int)h->GetEntries() << " modules";
1996   }
1997   if (showUnderOverFlow_) {
1998     if (showModules_) {
1999       legend << (int)h->GetEntries() << " modules ("
2000              << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " outside range)";
2001     } else {
2002       legend << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " modules outside range";
2003     }
2004   }
2005   plotinfo.legend->AddEntry((TObject*)nullptr, legend.str().c_str(), "");
2006 
2007   // Scale the x-axis (cm to um), if needed
2008   if (plotinfo.variable.find("Norm") == std::string::npos)
2009     scaleXaxis(h, 10000);
2010 }
2011 
2012 /*!
2013  * \fn plotDMRHistogram 
2014  * \brief Create the DMR histrogram using data stored in trees and store them in the plotinfo structure.
2015  */
2016 
2017 void PlotAlignmentValidation::plotDMRHistogram(PlotAlignmentValidation::DMRPlotInfo& plotinfo,
2018                                                int direction,
2019                                                int layer,
2020                                                std::string subdet) {
2021   TH1F* h = nullptr;
2022   //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.
2023 
2024   TString histoname = "";
2025   if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY")
2026     histoname = "median";
2027   else if (plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY")
2028     histoname = "DrmsNR";
2029   histoname += "_";
2030   histoname += plotinfo.vars->getName();
2031   histoname.ReplaceAll(" ", "_");
2032   histoname += "_";
2033   histoname += subdet.c_str();
2034   if (plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormY")
2035     histoname += "_y";
2036   if (layer != 0) {
2037     if (subdet == "TID" || subdet == "TEC")
2038       histoname += "_disc";
2039     else
2040       histoname += "_layer";
2041     histoname += std::to_string(layer);
2042   }
2043   if (direction == -1) {
2044     histoname += "_minus";
2045   } else if (direction == 1) {
2046     histoname += "_plus";
2047   } else {
2048     histoname += "";
2049   }
2050 
2051   //Plotting
2052   std::string plotVariable =
2053       getVariableForDMRPlot(histoname.Data(), plotinfo.variable, plotinfo.nbins, plotinfo.min, plotinfo.max);
2054   std::string selection = "";
2055   if (plotinfo.filterName.empty()) {
2056     //Use only default selection and no filter
2057     selection = getSelectionForDMRPlot(plotinfo.minHits, plotinfo.subDetId, direction, layer);
2058     plotinfo.vars->getTree()->Draw(plotVariable.c_str(), selection.c_str(), "goff");
2059     if (gDirectory)
2060       gDirectory->GetObject(histoname.Data(), h);
2061     if (h && h->GetEntries() > 0) {
2062       if (direction == -1) {
2063         plotinfo.h1 = h;
2064       } else if (direction == 1) {
2065         plotinfo.h2 = h;
2066       } else {
2067         plotinfo.h = h;
2068       }
2069     }
2070   } else {
2071     TTreeReader reader(plotinfo.vars->getTree());
2072     TTreeReaderValue<Float_t> varToPlot(reader, plotinfo.variable.c_str());
2073     TTreeReaderValue<unsigned int> _entries(reader, "entries");
2074     TTreeReaderValue<unsigned int> _subDetId(reader, "subDetId");
2075     TTreeReaderValue<unsigned int> _moduleId(reader, "moduleId");
2076     TTreeReaderValue<Float_t> _zDirection(reader, "zDirection");
2077     TTreeReaderValue<Float_t> _rDirection(reader, "rDirection");
2078     TTreeReaderValue<unsigned int> _layer(reader, "layer");
2079     std::string badModulesFile_ = plotinfo.filterName;
2080     TFile* fBadModules = new TFile(badModulesFile_.c_str(), "READ");
2081     TTree* tBadModules = (TTree*)fBadModules->Get("alignTree");
2082     TTreeReader readerBad(tBadModules);
2083     TTreeReaderValue<int> _valid(readerBad, "valid");
2084     TTreeReaderValue<int> _bad_id(readerBad, "id");
2085     TTreeReaderValue<double> _bad_lumi(readerBad, "lumi");
2086 
2087     //Record which modules were used
2088     std::ofstream fUsedModules;
2089     fUsedModules.open("usedModules.txt", std::ios::out | std::ios::app);
2090 
2091     //Filter on modules by hand together with base selection
2092     for (uint i = 0; i < plotinfo.vars->getTree()->GetEntries(); i++) {
2093       reader.SetEntry(i);
2094       if (*_entries < uint(plotinfo.minHits))
2095         continue;
2096       if (*_subDetId != uint(plotinfo.subDetId))
2097         continue;
2098       if (direction != 0) {
2099         if (plotinfo.subDetId == 2) {  // FPIX is split by zDirection
2100           if (*_zDirection != direction)
2101             continue;
2102         } else {
2103           if (*_rDirection != direction)
2104             continue;
2105         }
2106       }
2107       if (layer > 0) {
2108         if (*_layer != uint(layer))
2109           continue;
2110       }
2111       bool isBadModule = false;
2112       for (int ibad = 0; ibad < tBadModules->GetEntries(); ibad++) {
2113         readerBad.SetEntry(ibad);
2114         //if (*_valid == 0) {continue;} //only modules that failed 0 times are OK = very strict
2115         if (subdet == "BPIX" || subdet == "FPIX") {
2116           if (*_bad_lumi <= 2.0)
2117             continue;
2118         } else {
2119           if (*_bad_lumi <= 7.0)
2120             continue;
2121         }
2122         //modules that misbehave for less than 2/fb are OK = mild strict
2123         if (*_moduleId == uint(*_bad_id))
2124           isBadModule = true;
2125       }
2126 
2127       if (isBadModule)
2128         continue;
2129       fUsedModules << *_moduleId << "\n";
2130       if (h) {
2131         h->Fill(*varToPlot);
2132       }
2133     }
2134 
2135     //Finalize
2136     fUsedModules.close();
2137     fBadModules->Close();
2138     if (h) {
2139       h->SetName(histoname.Data());
2140     }
2141   }
2142 
2143   //Store histogram
2144   if (h && h->GetEntries() > 0) {
2145     if (direction == -1) {
2146       plotinfo.h1 = h;
2147     } else if (direction == 1) {
2148       plotinfo.h2 = h;
2149     } else {
2150       plotinfo.h = h;
2151     }
2152   }
2153   if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormX" ||
2154       plotinfo.variable == "rmsNormY") {
2155     storeHistogramInRootfile(h);
2156   }
2157 }
2158 
2159 void PlotAlignmentValidation::modifySSHistAndLegend(THStack* hs, TLegend* legend) {
2160   // Add mean-y-values to the legend and scale the histograms.
2161 
2162   Double_t legendY = 0.80;
2163   bool hasheader = (TkAlStyle::legendheader != "");
2164   if (hasheader)
2165     legend->SetHeader(TkAlStyle::legendheader);
2166   legend->SetFillStyle(0);
2167   int legendsize = hs->GetHists()->GetSize() + hasheader;
2168 
2169   if (legendsize > 3)
2170     legendY -= 0.01 * (legendsize - 3);
2171   if (bigtext_) {
2172     legendY -= 0.05;
2173   }
2174   if (legendY < 0.6) {
2175     std::cerr << "Warning: Huge legend!" << std::endl;
2176     legendY = 0.6;
2177   }
2178   legend->SetY1(legendY);
2179   if (bigtext_)
2180     legend->SetTextSize(TkAlStyle::textSize);
2181 
2182   // Loop over all profiles
2183   TProfile* prof = nullptr;
2184   TIter next(hs->GetHists());
2185   Int_t index = hasheader;  //if hasheader, first entry is the header itself
2186   while ((prof = (TProfile*)next())) {
2187     //Scaling: from cm to um
2188     Double_t scale = 10000;
2189     prof->Scale(scale);
2190 
2191     Double_t stats[6] = {0};
2192     prof->GetStats(stats);
2193 
2194     std::ostringstream legendtext;
2195     legendtext.precision(3);
2196     legendtext << std::fixed;  // to always show 3 decimals
2197     legendtext << ": y mean = " << stats[4] / stats[0] * scale << " #mum";
2198 
2199     TLegendEntry* entry = (TLegendEntry*)legend->GetListOfPrimitives()->At(index);
2200     if (entry == nullptr)
2201       std::cout << "PlotAlignmentValidation::PlotAlignmentValidation::modifySSLegend: Bad legend!" << std::endl;
2202     else
2203       entry->SetLabel((entry->GetLabel() + legendtext.str()).c_str());
2204     index++;
2205   }
2206 
2207   // Make some room for the legend
2208   hs->SetMaximum(hs->GetMaximum("nostack PE") * 1.3);
2209 }
2210 
2211 //random variable: \sigma_{X_1}-\sigma_{X_2}-\delta_{RMS}
2212 //is centered approx around 0
2213 //null hypothesis: \delta_{RMS}=0
2214 //so \delta_\sigma is a realization of this random variable
2215 //how probable is it to get our value of \delta_\sigma?
2216 //->p-value
2217 double PlotAlignmentValidation::resampleTestOfEqualRMS(TH1F* h1, TH1F* h2, int numSamples) {
2218   //vector to store realizations of random variable
2219   std::vector<double> diff;
2220   diff.clear();
2221   //"true" (in bootstrap terms) difference of the samples' RMS
2222   double rmsdiff = abs(h1->GetRMS() - h2->GetRMS());
2223   //means of the samples to calculate RMS
2224   double m1 = h1->GetMean();
2225   double m2 = h2->GetMean();
2226   //realization of random variable
2227   double d1 = 0;
2228   double d2 = 0;
2229   //mean of random variable
2230   double test_mean = 0;
2231   for (int i = 0; i < numSamples; i++) {
2232     d1 = 0;
2233     d2 = 0;
2234     for (int i = 0; i < h1->GetEntries(); i++) {
2235       d1 += h1->GetRandom() - m1;
2236     }
2237     for (int i = 0; i < h2->GetEntries(); i++) {
2238       d2 += h2->GetRandom() + m2;
2239     }
2240     d1 /= h1->GetEntries();
2241     d2 /= h2->GetEntries();
2242     diff.push_back(abs(d1 - d2 - rmsdiff));
2243     test_mean += abs(d1 - d2 - rmsdiff);
2244   }
2245   test_mean /= numSamples;
2246   edm::LogPrint("") << "test mean:" << test_mean;
2247   //p value
2248   double p = 0;
2249   for (double d : diff) {
2250     if (d > rmsdiff) {
2251       p += 1;
2252     }
2253   }
2254 
2255   p /= numSamples;
2256   return p;
2257 }
2258 
2259 //random variable: (\overline{X_1}-\mu_1)-(\overline{X_2}-\mu_2)
2260 //is centered approx around 0
2261 //null hypothesis: \mu_1-\mu_2=0
2262 //so \delta_\mu is a realization of this random variable
2263 //how probable is it to get our value of \delta_\mu?
2264 //->p-value
2265 double PlotAlignmentValidation::resampleTestOfEqualMeans(TH1F* h1, TH1F* h2, int numSamples) {
2266   //vector to store realization of random variable
2267   std::vector<double> diff;
2268   diff.clear();
2269   //"true" (in bootstrap terms) difference of the samples' means
2270   double meandiff = abs(h1->GetMean() - h2->GetMean());
2271   //realization of random variable
2272   double d1 = 0;
2273   double d2 = 0;
2274   //mean of random variable
2275   double test_mean = 0;
2276   for (int i = 0; i < numSamples; i++) {
2277     d1 = 0;
2278     d2 = 0;
2279     for (int i = 0; i < h1->GetEntries(); i++) {
2280       d1 += h1->GetRandom();
2281     }
2282     for (int i = 0; i < h2->GetEntries(); i++) {
2283       d2 += h2->GetRandom();
2284     }
2285     d1 /= h1->GetEntries();
2286     d2 /= h2->GetEntries();
2287     diff.push_back(abs(d1 - d2 - meandiff));
2288     test_mean += abs(d1 - d2 - meandiff);
2289   }
2290   test_mean /= numSamples;
2291   edm::LogPrint("") << "test mean:" << test_mean;
2292   //p-value
2293   double p = 0;
2294   for (double d : diff) {
2295     if (d > meandiff) {
2296       p += 1;
2297     }
2298   }
2299 
2300   p /= numSamples;
2301   return p;
2302 }
2303 
2304 float PlotAlignmentValidation::twotailedStudentTTestEqualMean(float t, float v) {
2305   return 2 * (1 - ROOT::Math::tdistribution_cdf(abs(t), v));
2306 }
2307 
2308 const TString PlotAlignmentValidation::summaryfilename = "OfflineValidationSummary";
2309 
2310 std::vector<TH1*> PlotAlignmentValidation::findmodule(TFile* f, unsigned int moduleid) {
2311   //TFile *f = TFile::Open(filename, "READ");
2312   TString histnamex;
2313   TString histnamey;
2314   //read necessary branch/folder
2315   auto t = (TTree*)f->Get("TrackerOfflineValidation/TkOffVal");
2316 
2317   TkOffTreeVariables* variables = nullptr;
2318   t->SetBranchAddress("TkOffTreeVariables", &variables);
2319   unsigned int number_of_entries = t->GetEntries();
2320   for (unsigned int i = 0; i < number_of_entries; i++) {
2321     t->GetEntry(i);
2322     if (variables->moduleId == moduleid) {
2323       histnamex = variables->histNameX;
2324       histnamey = variables->histNameY;
2325       break;
2326     }
2327   }
2328 
2329   std::vector<TH1*> h;
2330 
2331   auto h1 = (TH1*)f->FindObjectAny(histnamex);
2332   auto h2 = (TH1*)f->FindObjectAny(histnamey);
2333 
2334   h1->SetDirectory(nullptr);
2335   h2->SetDirectory(nullptr);
2336 
2337   h.push_back(h1);
2338   h.push_back(h2);
2339 
2340   return h;
2341 }
2342 
2343 void PlotAlignmentValidation::residual_by_moduleID(unsigned int moduleid) {
2344   TCanvas* cx = new TCanvas("x_residual");
2345   TCanvas* cy = new TCanvas("y_residual");
2346   TLegend* legendx = new TLegend(0.55, 0.7, 1, 0.9);
2347   TLegend* legendy = new TLegend(0.55, 0.7, 1, 0.9);
2348 
2349   legendx->SetTextSize(0.016);
2350   legendx->SetTextAlign(12);
2351   legendy->SetTextSize(0.016);
2352   legendy->SetTextAlign(12);
2353 
2354   for (auto it : sourceList) {
2355     TFile* file = it->getFile();
2356     int color = it->getLineColor();
2357     int linestyle = it->getLineStyle();  //this you set by doing h->SetLineStyle(linestyle)
2358     TString legendname = it->getName();  //this goes in the legend
2359     std::vector<TH1*> hist = findmodule(file, moduleid);
2360 
2361     TString histnamex = legendname + " NEntries: " + std::to_string(int(hist[0]->GetEntries()));
2362     hist[0]->SetTitle(histnamex);
2363     hist[0]->SetStats(false);
2364     hist[0]->Rebin(50);
2365     hist[0]->SetBit(TH1::kNoTitle);
2366     hist[0]->SetLineColor(color);
2367     hist[0]->SetLineStyle(linestyle);
2368     cx->cd();
2369     hist[0]->Draw("Same");
2370     legendx->AddEntry(hist[0], histnamex, "l");
2371 
2372     TString histnamey = legendname + " NEntries: " + std::to_string(int(hist[1]->GetEntries()));
2373     hist[1]->SetTitle(histnamey);
2374     hist[1]->SetStats(false);
2375     hist[1]->Rebin(50);
2376     hist[1]->SetBit(TH1::kNoTitle);
2377     hist[1]->SetLineColor(color);
2378     hist[1]->SetLineStyle(linestyle);
2379     cy->cd();
2380     hist[1]->Draw("Same");
2381     legendy->AddEntry(hist[1], histnamey, "l");
2382   }
2383 
2384   TString filenamex = "x_residual_" + std::to_string(moduleid);
2385   TString filenamey = "y_residual_" + std::to_string(moduleid);
2386   cx->cd();
2387   legendx->Draw();
2388   cx->SaveAs(outputDir + "/" + filenamex + ".root");
2389   cx->SaveAs(outputDir + "/" + filenamex + ".pdf");
2390   cx->SaveAs(outputDir + "/" + filenamex + ".png");
2391   cx->SaveAs(outputDir + "/" + filenamex + ".eps");
2392 
2393   cy->cd();
2394   legendy->Draw();
2395   cy->SaveAs(outputDir + "/" + filenamey + ".root");
2396   cy->SaveAs(outputDir + "/" + filenamey + ".pdf");
2397   cy->SaveAs(outputDir + "/" + filenamey + ".png");
2398   cy->SaveAs(outputDir + "/" + filenamey + ".eps");
2399 }