Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:55:57 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ApeEstimatorSummary
0004 // Class:      ApeEstimatorSummary
0005 //
0006 /**\class ApeEstimatorSummary ApeEstimatorSummary.cc Alignment/APEEstimation/src/ApeEstimatorSummary.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Johannes Hauk,6 2-039,+41227673512,
0015 //         Created:  Mon Oct 11 13:44:03 CEST 2010
0016 //         Modified by: Christian Schomakers (RWTH Aachen)
0017 // $Id: ApeEstimatorSummary.cc,v 1.14 2012/01/26 00:06:33 hauk Exp $
0018 //
0019 //
0020 
0021 // system include files
0022 #include <memory>
0023 #include <fstream>
0024 #include <sstream>
0025 
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 #include "FWCore/Framework/interface/EventSetup.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/Utilities/interface/EDMException.h"
0037 #include "FWCore/Utilities/interface/isFinite.h"
0038 
0039 #include "CLHEP/Matrix/SymMatrix.h"
0040 
0041 //#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0042 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0043 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0044 
0045 //...............
0046 #include "Alignment/APEEstimation/interface/TrackerSectorStruct.h"
0047 
0048 #include "TH1.h"
0049 #include "TString.h"
0050 #include "TFile.h"
0051 #include "TDirectory.h"
0052 #include "TTree.h"
0053 #include "TF1.h"
0054 #include "TMath.h"
0055 //
0056 // class declaration
0057 //
0058 
0059 class ApeEstimatorSummary : public edm::one::EDAnalyzer<> {
0060 public:
0061   explicit ApeEstimatorSummary(const edm::ParameterSet&);
0062   ~ApeEstimatorSummary() override;
0063 
0064 private:
0065   void beginJob() override;
0066   void analyze(const edm::Event&, const edm::EventSetup&) override;
0067   void endJob() override;
0068 
0069   void openInputFile();
0070   void getTrackerSectorStructs();
0071   void bookHists();
0072   std::vector<double> residualErrorBinning();
0073   void writeHists();
0074 
0075   void calculateApe();
0076 
0077   enum ApeWeight { wInvalid, wUnity, wEntries, wEntriesOverSigmaX2 };
0078 
0079   // ----------member data ---------------------------
0080   const edm::ParameterSet parameterSet_;
0081   const edm::ESGetToken<AlignmentErrorsExtended, TrackerAlignmentErrorExtendedRcd> alignmentErrorToken_;
0082 
0083   bool firstEvent = true;
0084 
0085   TFile* inputFile_;
0086 
0087   std::map<unsigned int, TrackerSectorStruct> m_tkSector_;
0088   unsigned int noSectors_;
0089 };
0090 
0091 //
0092 // constants, enums and typedefs
0093 //
0094 
0095 //
0096 // static data member definitions
0097 //
0098 
0099 //
0100 // constructors and destructor
0101 //
0102 ApeEstimatorSummary::ApeEstimatorSummary(const edm::ParameterSet& iConfig)
0103     : parameterSet_(iConfig), alignmentErrorToken_(esConsumes()), inputFile_(nullptr) {}
0104 
0105 ApeEstimatorSummary::~ApeEstimatorSummary() {}
0106 
0107 //
0108 // member functions
0109 //
0110 
0111 void ApeEstimatorSummary::openInputFile() {
0112   const std::string inputFileName(parameterSet_.getParameter<std::string>("InputFile"));
0113   std::ifstream inputFileStream;
0114   // Check if baseline file exists
0115   inputFileStream.open(inputFileName.c_str());
0116   if (inputFileStream.is_open()) {
0117     inputFile_ = new TFile(inputFileName.c_str(), "READ");
0118   }
0119   if (inputFile_) {
0120     edm::LogInfo("CalculateAPE") << "Input file from loop over tracks and corresponding hits sucessfully opened";
0121   } else {
0122     edm::LogError("CalculateAPE") << "There is NO input file\n"
0123                                   << "...APE calculation stopped. Please check path of input file name in config:\n"
0124                                   << "\t" << inputFileName;
0125     throw edm::Exception(edm::errors::Configuration, "Bad input file name");
0126   }
0127 }
0128 
0129 void ApeEstimatorSummary::getTrackerSectorStructs() {
0130   // At present first of the plugins registered in TFileService needs to be the one containing the normalized residual histos per sector per error bin
0131   TString pluginName(inputFile_->GetListOfKeys()->At(0)->GetName());
0132 
0133   pluginName += "/";
0134   TDirectory *sectorDir(nullptr), *intervalDir(nullptr);
0135   bool sectorBool(true);
0136   unsigned int iSector(1);
0137   for (; sectorBool; ++iSector) {
0138     std::stringstream sectorName, fullSectorName;
0139     sectorName << "Sector_" << iSector << "/";
0140     fullSectorName << pluginName << sectorName.str();
0141     TString fullName(fullSectorName.str().c_str());
0142     inputFile_->cd();
0143     sectorDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
0144     if (sectorDir) {
0145       TrackerSectorStruct tkSector;
0146       sectorDir->GetObject("z_name;1", tkSector.Name);
0147       tkSector.name = tkSector.Name->GetTitle();
0148       bool intervalBool(true);
0149       for (unsigned int iInterval(1); intervalBool; ++iInterval) {
0150         std::stringstream intervalName, fullIntervalName;
0151         intervalName << "Interval_" << iInterval << "/";
0152         fullIntervalName << fullSectorName.str() << intervalName.str();
0153         fullName = fullIntervalName.str().c_str();
0154         intervalDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
0155         if (intervalDir) {
0156           intervalDir->GetObject("h_sigmaX;1", tkSector.m_binnedHists[iInterval]["sigmaX"]);
0157           intervalDir->GetObject("h_norResX;1", tkSector.m_binnedHists[iInterval]["norResX"]);
0158           intervalDir->GetObject("h_sigmaY;1", tkSector.m_binnedHists[iInterval]["sigmaY"]);
0159           intervalDir->GetObject("h_norResY;1", tkSector.m_binnedHists[iInterval]["norResY"]);
0160           if (tkSector.m_binnedHists[iInterval]["sigmaY"] && tkSector.m_binnedHists[iInterval]["norResY"]) {
0161             tkSector.isPixel = true;
0162           }
0163         } else {
0164           intervalBool = false;
0165           if (iSector == 1)
0166             edm::LogInfo("CalculateAPE") << "There are " << iInterval - 1
0167                                          << " intervals per sector defined in input file";
0168         }
0169       }
0170       TDirectory* resultsDir(nullptr);
0171       std::stringstream fullResultName;
0172       fullResultName << fullSectorName.str() << "Results/";
0173       fullName = fullResultName.str().c_str();
0174       resultsDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
0175       if (resultsDir) {
0176         resultsDir->GetObject("h_entriesX;1", tkSector.EntriesX);
0177         if (tkSector.isPixel)
0178           resultsDir->GetObject("h_entriesY;1", tkSector.EntriesY);
0179         TTree* rawIdTree(nullptr);
0180         resultsDir->GetObject("rawIdTree", rawIdTree);
0181         unsigned int rawId(0);
0182         rawIdTree->SetBranchAddress("RawId", &rawId);
0183         for (int entry = 0; entry < rawIdTree->GetEntries(); ++entry) {
0184           rawIdTree->GetEntry(entry);
0185           // command "hadd" adds entries in TTree, so rawId are existing as often as number of files are added
0186           bool alreadyAdded(false);
0187           for (auto const& i_rawId : tkSector.v_rawId) {
0188             if (rawId == i_rawId)
0189               alreadyAdded = true;
0190           }
0191           if (alreadyAdded)
0192             break;
0193           tkSector.v_rawId.push_back(rawId);
0194         }
0195       }
0196       m_tkSector_[iSector] = tkSector;
0197     } else {
0198       sectorBool = false;
0199       edm::LogInfo("CalculateAPE") << "There are " << iSector - 1 << " sectors defined in input file";
0200     }
0201   }
0202   noSectors_ = iSector + 1;
0203 }
0204 
0205 void ApeEstimatorSummary::bookHists() {
0206   const std::vector<double> v_binX(this->residualErrorBinning());
0207   for (auto& i_sector : m_tkSector_) {
0208     i_sector.second.WeightX = new TH1F("h_weightX",
0209                                        "relative weight w_{x}/w_{tot,x};#sigma_{x}  [#mum];w_{x}/w_{tot,x}",
0210                                        v_binX.size() - 1,
0211                                        &(v_binX[0]));
0212     i_sector.second.MeanX = new TH1F("h_meanX",
0213                                      "residual mean <r_{x}/#sigma_{r,x}>;#sigma_{x}  [#mum];<r_{x}/#sigma_{r,x}>",
0214                                      v_binX.size() - 1,
0215                                      &(v_binX[0]));
0216     i_sector.second.RmsX = new TH1F("h_rmsX",
0217                                     "residual rms RMS(r_{x}/#sigma_{r,x});#sigma_{x}  [#mum];RMS(r_{x}/#sigma_{r,x})",
0218                                     v_binX.size() - 1,
0219                                     &(v_binX[0]));
0220     i_sector.second.FitMeanX1 = new TH1F(
0221         "h_fitMeanX1", "fitted residual mean #mu_{x};#sigma_{x}  [#mum];#mu_{x}", v_binX.size() - 1, &(v_binX[0]));
0222     i_sector.second.ResidualWidthX1 = new TH1F("h_residualWidthX1",
0223                                                "residual width #Delta_{x};#sigma_{x}  [#mum];#Delta_{x}",
0224                                                v_binX.size() - 1,
0225                                                &(v_binX[0]));
0226     i_sector.second.CorrectionX1 = new TH1F("h_correctionX1",
0227                                             "correction to APE_{x};#sigma_{x}  [#mum];#Delta#sigma_{align,x}  [#mum]",
0228                                             v_binX.size() - 1,
0229                                             &(v_binX[0]));
0230     i_sector.second.FitMeanX2 = new TH1F(
0231         "h_fitMeanX2", "fitted residual mean #mu_{x};#sigma_{x}  [#mum];#mu_{x}", v_binX.size() - 1, &(v_binX[0]));
0232     i_sector.second.ResidualWidthX2 = new TH1F("h_residualWidthX2",
0233                                                "residual width #Delta_{x};#sigma_{x}  [#mum];#Delta_{x}",
0234                                                v_binX.size() - 1,
0235                                                &(v_binX[0]));
0236     i_sector.second.CorrectionX2 = new TH1F("h_correctionX2",
0237                                             "correction to APE_{x};#sigma_{x}  [#mum];#Delta#sigma_{align,x}  [#mum]",
0238                                             v_binX.size() - 1,
0239                                             &(v_binX[0]));
0240 
0241     if (i_sector.second.isPixel) {
0242       i_sector.second.WeightY = new TH1F("h_weightY",
0243                                          "relative weight w_{y}/w_{tot,y};#sigma_{y}  [#mum];w_{y}/w_{tot,y}",
0244                                          v_binX.size() - 1,
0245                                          &(v_binX[0]));
0246       i_sector.second.MeanY = new TH1F("h_meanY",
0247                                        "residual mean <r_{y}/#sigma_{r,y}>;#sigma_{y}  [#mum];<r_{y}/#sigma_{r,y}>",
0248                                        v_binX.size() - 1,
0249                                        &(v_binX[0]));
0250       i_sector.second.RmsY = new TH1F("h_rmsY",
0251                                       "residual rms RMS(r_{y}/#sigma_{r,y});#sigma_{y}  [#mum];RMS(r_{y}/#sigma_{r,y})",
0252                                       v_binX.size() - 1,
0253                                       &(v_binX[0]));
0254       i_sector.second.FitMeanY1 = new TH1F(
0255           "h_fitMeanY1", "fitted residual mean #mu_{y};#sigma_{y}  [#mum];#mu_{y}", v_binX.size() - 1, &(v_binX[0]));
0256       i_sector.second.ResidualWidthY1 = new TH1F("h_residualWidthY1",
0257                                                  "residual width #Delta_{y};#sigma_{y}  [#mum];#Delta_{y}",
0258                                                  v_binX.size() - 1,
0259                                                  &(v_binX[0]));
0260       i_sector.second.CorrectionY1 = new TH1F("h_correctionY1",
0261                                               "correction to APE_{y};#sigma_{y}  [#mum];#Delta#sigma_{align,y}  [#mum]",
0262                                               v_binX.size() - 1,
0263                                               &(v_binX[0]));
0264       i_sector.second.FitMeanY2 = new TH1F(
0265           "h_fitMeanY2", "fitted residual mean #mu_{y};#sigma_{y}  [#mum];#mu_{y}", v_binX.size() - 1, &(v_binX[0]));
0266       i_sector.second.ResidualWidthY2 = new TH1F("h_residualWidthY2",
0267                                                  "residual width #Delta_{y};#sigma_{y}  [#mum];#Delta_{y}",
0268                                                  v_binX.size() - 1,
0269                                                  &(v_binX[0]));
0270       i_sector.second.CorrectionY2 = new TH1F("h_correctionY2",
0271                                               "correction to APE_{y};#sigma_{y}  [#mum];#Delta#sigma_{align,y}  [#mum]",
0272                                               v_binX.size() - 1,
0273                                               &(v_binX[0]));
0274     }
0275   }
0276 }
0277 
0278 std::vector<double> ApeEstimatorSummary::residualErrorBinning() {
0279   std::vector<double> v_binX;
0280   TH1* EntriesX(m_tkSector_[1].EntriesX);
0281   for (int iBin = 1; iBin <= EntriesX->GetNbinsX() + 1; ++iBin) {
0282     v_binX.push_back(EntriesX->GetBinLowEdge(iBin));
0283   }
0284   return v_binX;
0285 }
0286 
0287 void ApeEstimatorSummary::writeHists() {
0288   TFile* resultsFile = new TFile(parameterSet_.getParameter<std::string>("ResultsFile").c_str(), "RECREATE");
0289   TDirectory* baseDir = resultsFile->mkdir("ApeEstimatorSummary");
0290   for (auto const& i_sector : m_tkSector_) {
0291     std::stringstream dirName;
0292     dirName << "Sector_" << i_sector.first;
0293     TDirectory* dir = baseDir->mkdir(dirName.str().c_str());
0294     dir->cd();
0295 
0296     i_sector.second.Name->Write();
0297 
0298     i_sector.second.WeightX->Write();
0299     i_sector.second.MeanX->Write();
0300     i_sector.second.RmsX->Write();
0301     i_sector.second.FitMeanX1->Write();
0302     i_sector.second.ResidualWidthX1->Write();
0303     i_sector.second.CorrectionX1->Write();
0304     i_sector.second.FitMeanX2->Write();
0305     i_sector.second.ResidualWidthX2->Write();
0306     i_sector.second.CorrectionX2->Write();
0307 
0308     if (i_sector.second.isPixel) {
0309       i_sector.second.WeightY->Write();
0310       i_sector.second.MeanY->Write();
0311       i_sector.second.RmsY->Write();
0312       i_sector.second.FitMeanY1->Write();
0313       i_sector.second.ResidualWidthY1->Write();
0314       i_sector.second.CorrectionY1->Write();
0315       i_sector.second.FitMeanY2->Write();
0316       i_sector.second.ResidualWidthY2->Write();
0317       i_sector.second.CorrectionY2->Write();
0318     }
0319   }
0320   resultsFile->Close();
0321 }
0322 
0323 void ApeEstimatorSummary::calculateApe() {
0324   // Set baseline or calculate APE value?
0325   const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
0326 
0327   // Read in baseline file for calculation of APE value (if not setting baseline)
0328   // Has same format as iterationFile
0329   const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
0330   TFile* baselineFile(nullptr);
0331   TTree* baselineTreeX(nullptr);
0332   TTree* baselineTreeY(nullptr);
0333   TTree* sectorNameBaselineTree(nullptr);
0334   if (!setBaseline) {
0335     std::ifstream baselineFileStream;
0336     // Check if baseline file exists
0337     baselineFileStream.open(baselineFileName.c_str());
0338     if (baselineFileStream.is_open()) {
0339       baselineFileStream.close();
0340       baselineFile = new TFile(baselineFileName.c_str(), "READ");
0341     }
0342     if (baselineFile) {
0343       edm::LogInfo("CalculateAPE") << "Baseline file for APE values sucessfully opened";
0344       baselineFile->GetObject("iterTreeX;1", baselineTreeX);
0345       baselineFile->GetObject("iterTreeY;1", baselineTreeY);
0346       baselineFile->GetObject("nameTree;1", sectorNameBaselineTree);
0347     } else {
0348       edm::LogWarning("CalculateAPE") << "There is NO baseline file for APE values, so normalized residual width =1 "
0349                                          "for ideal conditions is assumed";
0350     }
0351   }
0352 
0353   // Set up root file for iterations on APE value (or for setting baseline in setBaseline mode)
0354   const std::string iterationFileName(setBaseline ? baselineFileName
0355                                                   : parameterSet_.getParameter<std::string>("IterationFile"));
0356   // For iterations, updates are needed to not overwrite the iterations before
0357   TFile* iterationFile = new TFile(iterationFileName.c_str(), setBaseline ? "RECREATE" : "UPDATE");
0358 
0359   // Set up TTree for iterative APE values on first pass (first iteration) or read from file (further iterations)
0360   TTree* iterationTreeX(nullptr);
0361   TTree* iterationTreeY(nullptr);
0362   iterationFile->GetObject("iterTreeX;1", iterationTreeX);
0363   iterationFile->GetObject("iterTreeY;1", iterationTreeY);
0364   // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as iterationTree)
0365   TTree* sectorNameTree(nullptr);
0366   iterationFile->GetObject("nameTree;1", sectorNameTree);
0367 
0368   const bool firstIter(!iterationTreeX);
0369   if (firstIter) {  // should be always true in setBaseline mode, since file is recreated
0370     if (!setBaseline) {
0371       iterationTreeX = new TTree("iterTreeX", "Tree for APE x values of all iterations");
0372       iterationTreeY = new TTree("iterTreeY", "Tree for APE y values of all iterations");
0373       edm::LogInfo("CalculateAPE") << "First APE iteration (number 0.), create iteration file with TTree";
0374       sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
0375     } else {
0376       iterationTreeX = new TTree("iterTreeX", "Tree for baseline x values of normalized residual width");
0377       iterationTreeY = new TTree("iterTreeY", "Tree for baseline y values of normalized residual width");
0378       edm::LogInfo("CalculateAPE") << "Set baseline, create baseline file with TTree";
0379       sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
0380     }
0381   } else {
0382     const unsigned int iteration(iterationTreeX->GetEntries());
0383     edm::LogWarning("CalculateAPE") << "NOT the first APE iteration (number 0.) but the " << iteration
0384                                     << ". one, is this wanted or forgot to delete old iteration file with TTree?";
0385   }
0386 
0387   // Assign the information stored in the trees to arrays
0388   double a_apeSectorX[noSectors_];
0389   double a_apeSectorY[noSectors_];
0390   double a_baselineSectorX[noSectors_];
0391   double a_baselineSectorY[noSectors_];
0392 
0393   std::string* a_sectorName[noSectors_];
0394   std::string* a_sectorBaselineName[noSectors_];
0395   for (auto const& i_sector : m_tkSector_) {
0396     const unsigned int iSector(i_sector.first);
0397     const bool pixelSector(i_sector.second.isPixel);
0398     a_apeSectorX[iSector] = 99.;
0399     a_apeSectorY[iSector] = 99.;
0400     a_baselineSectorX[iSector] = -99.;
0401     a_baselineSectorY[iSector] = -99.;
0402 
0403     a_sectorName[iSector] = nullptr;
0404     a_sectorBaselineName[iSector] = nullptr;
0405     std::stringstream ss_sector, ss_sectorSuffixed;
0406     ss_sector << "Ape_Sector_" << iSector;
0407     if (!setBaseline && baselineTreeX) {
0408       baselineTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorX[iSector]);
0409       baselineTreeX->GetEntry(0);
0410       if (pixelSector) {
0411         baselineTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorY[iSector]);
0412         baselineTreeY->GetEntry(0);
0413       }
0414       sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
0415       sectorNameBaselineTree->GetEntry(0);
0416     } else {
0417       // Set default ideal normalized residual width to 1
0418       a_baselineSectorX[iSector] = 1.;
0419       a_baselineSectorY[iSector] = 1.;
0420     }
0421     if (firstIter) {  // should be always true in setBaseline mode, since file is recreated
0422       ss_sectorSuffixed << ss_sector.str() << "/D";
0423       iterationTreeX->Branch(ss_sector.str().c_str(), &a_apeSectorX[iSector], ss_sectorSuffixed.str().c_str());
0424       if (pixelSector) {
0425         iterationTreeY->Branch(ss_sector.str().c_str(), &a_apeSectorY[iSector], ss_sectorSuffixed.str().c_str());
0426       }
0427       sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
0428     } else {
0429       iterationTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorX[iSector]);
0430       iterationTreeX->GetEntry(iterationTreeX->GetEntries() - 1);
0431       if (pixelSector) {
0432         iterationTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorY[iSector]);
0433         iterationTreeY->GetEntry(iterationTreeY->GetEntries() - 1);
0434       }
0435       sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
0436       sectorNameTree->GetEntry(0);
0437     }
0438   }
0439 
0440   // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
0441   bool failed(false);
0442   for (auto& i_sector : m_tkSector_) {
0443     const std::string& name(i_sector.second.name);
0444     if (firstIter) {
0445       a_sectorName[i_sector.first] = new std::string(name);
0446     } else {
0447       const std::string& nameLastIter(*a_sectorName[i_sector.first]);
0448       if (name != nameLastIter) {
0449         edm::LogError("CalculateAPE") << "Inconsistent sector definition in iterationFile for sector " << i_sector.first
0450                                       << ",\n"
0451                                       << "Recent iteration has name \"" << name << "\", while previous had \""
0452                                       << nameLastIter << "\"\n"
0453                                       << "...APE calculation stopped. Please check sector definitions in config!\n";
0454         failed = true;
0455       }
0456     }
0457     if (!setBaseline && baselineFile) {
0458       const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
0459       if (name != nameBaseline) {
0460         failed = true;
0461         edm::LogError("CalculateAPE") << "Inconsistent sector definition in baselineFile for sector " << i_sector.first
0462                                       << ",\n"
0463                                       << "Recent iteration has name \"" << name << "\", while baseline had \""
0464                                       << nameBaseline << "\"\n"
0465                                       << "...APE calculation stopped. Please check sector definitions in config!\n";
0466       }
0467     }
0468   }
0469 
0470   if (failed) {
0471     if (firstIter) {
0472       for (auto& i_sector : m_tkSector_) {
0473         delete a_sectorName[i_sector.first];
0474       }
0475     }
0476     return;
0477   }
0478 
0479   // Set up text file for writing out APE values per module
0480   std::ofstream apeOutputFile;
0481   if (!setBaseline) {
0482     const std::string apeOutputFileName(parameterSet_.getParameter<std::string>("ApeOutputFile"));
0483     apeOutputFile.open(apeOutputFileName.c_str());
0484     if (apeOutputFile.is_open()) {
0485       edm::LogInfo("CalculateAPE") << "Text file for writing APE values successfully opened";
0486     } else {
0487       edm::LogError("CalculateAPE") << "Text file for writing APE values NOT opened,\n"
0488                                     << "...APE calculation stopped. Please check path of text file name in config:\n"
0489                                     << "\t" << apeOutputFileName;
0490       return;
0491     }
0492   }
0493 
0494   // Loop over sectors for calculating APE
0495   const std::string apeWeightName(parameterSet_.getParameter<std::string>("apeWeight"));
0496   ApeWeight apeWeight(wInvalid);
0497   if (apeWeightName == "unity")
0498     apeWeight = wUnity;
0499   else if (apeWeightName == "entries")
0500     apeWeight = wEntries;
0501   else if (apeWeightName == "entriesOverSigmaX2")
0502     apeWeight = wEntriesOverSigmaX2;
0503   if (apeWeight == wInvalid) {
0504     edm::LogError("CalculateAPE") << "Invalid parameter 'apeWeight' in cfg file: \"" << apeWeightName
0505                                   << "\"\nimplemented apeWeights are \"unity\", \"entries\", \"entriesOverSigmaX2\""
0506                                   << "\n...APE calculation stopped.";
0507     return;
0508   }
0509   const double minHitsPerInterval(parameterSet_.getParameter<double>("minHitsPerInterval"));
0510   const double sigmaFactorFit(parameterSet_.getParameter<double>("sigmaFactorFit"));
0511   const double correctionScaling(parameterSet_.getParameter<double>("correctionScaling"));
0512   const bool smoothIteration(parameterSet_.getParameter<bool>("smoothIteration"));
0513   const double smoothFraction(parameterSet_.getParameter<double>("smoothFraction"));
0514   if (smoothFraction <= 0. || smoothFraction > 1.) {
0515     edm::LogError("CalculateAPE") << "Incorrect parameter in cfg file,"
0516                                   << "\nsmoothFraction has to be in [0,1], but is " << smoothFraction
0517                                   << "\n...APE calculation stopped.";
0518     return;
0519   }
0520   for (auto& i_sector : m_tkSector_) {
0521     typedef std::pair<double, double> Error2AndResidualWidth2PerBin;
0522     typedef std::pair<double, Error2AndResidualWidth2PerBin> WeightAndResultsPerBin;
0523     std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinX;
0524     std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinY;
0525 
0526     double baselineWidthX2(a_baselineSectorX[i_sector.first]);
0527     double baselineWidthY2(a_baselineSectorY[i_sector.first]);
0528 
0529     // Loop over residual error bins to calculate APE for every bin
0530     for (auto const& i_errBins : i_sector.second.m_binnedHists) {
0531       std::map<std::string, TH1*> mHists = i_errBins.second;
0532 
0533       double entriesX = mHists["sigmaX"]->GetEntries();
0534       double meanSigmaX = mHists["sigmaX"]->GetMean();
0535 
0536       // Fitting Parameters
0537       double xMin = mHists["norResX"]->GetXaxis()->GetXmin();
0538       double xMax = mHists["norResX"]->GetXaxis()->GetXmax();
0539       double integralX = mHists["norResX"]->Integral();
0540       double meanX = mHists["norResX"]->GetMean();
0541       double rmsX = mHists["norResX"]->GetRMS();
0542       double maximumX = mHists["norResX"]->GetBinContent(mHists["norResX"]->GetMaximumBin());
0543 
0544       // First Gaus Fit
0545       TF1 funcX_1("mygausX", "gaus", xMin, xMax);
0546       funcX_1.SetParameters(maximumX, meanX, rmsX);
0547       TString fitOpt("ILERQ");  //TString fitOpt("IMR"); ("IRLL"); ("IRQ");
0548       if (integralX > minHitsPerInterval) {
0549         if (mHists["norResX"]->Fit(&funcX_1, fitOpt)) {
0550           edm::LogWarning("CalculateAPE") << "Fit1 did not work: " << mHists["norResX"]->Fit(&funcX_1, fitOpt);
0551           continue;
0552         }
0553         LogDebug("CalculateAPE") << "FitResultX1\t" << mHists["norResX"]->Fit(&funcX_1, fitOpt) << "\n";
0554       }
0555       double meanX_1 = funcX_1.GetParameter(1);
0556       double sigmaX_1 = funcX_1.GetParameter(2);
0557 
0558       // Second gaus fit
0559       TF1 funcX_2("mygausX2",
0560                   "gaus",
0561                   meanX_1 - sigmaFactorFit * TMath::Abs(sigmaX_1),
0562                   meanX_1 + sigmaFactorFit * TMath::Abs(sigmaX_1));
0563       funcX_2.SetParameters(funcX_1.GetParameter(0), meanX_1, sigmaX_1);
0564       if (integralX > minHitsPerInterval) {
0565         if (mHists["norResX"]->Fit(&funcX_2, fitOpt)) {
0566           edm::LogWarning("CalculateAPE") << "Fit2 did not work for x : " << mHists["norResX"]->Fit(&funcX_2, fitOpt);
0567           continue;
0568         }
0569         LogDebug("CalculateAPE") << "FitResultX2\t" << mHists["norResX"]->Fit(&funcX_2, fitOpt) << "\n";
0570       }
0571       double meanX_2 = funcX_2.GetParameter(1);
0572       double sigmaX_2 = funcX_2.GetParameter(2);
0573 
0574       // Now the same for y coordinate
0575       double entriesY(0.);
0576       double meanSigmaY(0.);
0577       if (i_sector.second.isPixel) {
0578         entriesY = mHists["sigmaY"]->GetEntries();
0579         meanSigmaY = mHists["sigmaY"]->GetMean();
0580       }
0581 
0582       double meanY_1(0.);
0583       double sigmaY_1(0.);
0584       double meanY_2(0.);
0585       double sigmaY_2(0.);
0586       double meanY(0.);
0587       double rmsY(0.);
0588       if (i_sector.second.isPixel) {
0589         // Fitting Parameters
0590         double yMin = mHists["norResY"]->GetXaxis()->GetXmin();
0591         double yMax = mHists["norResY"]->GetXaxis()->GetXmax();
0592         double integralY = mHists["norResY"]->Integral();
0593         meanY = mHists["norResY"]->GetMean();
0594         rmsY = mHists["norResY"]->GetRMS();
0595         double maximumY = mHists["norResY"]->GetBinContent(mHists["norResY"]->GetMaximumBin());
0596 
0597         // First Gaus Fit
0598         TF1 funcY_1("mygausY", "gaus", yMin, yMax);
0599         funcY_1.SetParameters(maximumY, meanY, rmsY);
0600         if (integralY > minHitsPerInterval) {
0601           if (mHists["norResY"]->Fit(&funcY_1, fitOpt)) {
0602             edm::LogWarning("CalculateAPE") << "Fit1 did not work: " << mHists["norResY"]->Fit(&funcY_1, fitOpt);
0603             continue;
0604           }
0605           LogDebug("CalculateAPE") << "FitResultY1\t" << mHists["norResY"]->Fit(&funcY_1, fitOpt) << "\n";
0606         }
0607         meanY_1 = funcY_1.GetParameter(1);
0608         sigmaY_1 = funcY_1.GetParameter(2);
0609 
0610         // Second gaus fit
0611         TF1 funcY_2("mygausY2",
0612                     "gaus",
0613                     meanY_1 - sigmaFactorFit * TMath::Abs(sigmaY_1),
0614                     meanY_1 + sigmaFactorFit * TMath::Abs(sigmaY_1));
0615         funcY_2.SetParameters(funcY_1.GetParameter(0), meanY_1, sigmaY_1);
0616         if (integralY > minHitsPerInterval) {
0617           if (mHists["norResY"]->Fit(&funcY_2, fitOpt)) {
0618             edm::LogWarning("CalculateAPE") << "Fit2 did not work for y : " << mHists["norResY"]->Fit(&funcY_2, fitOpt);
0619             continue;
0620           }
0621           LogDebug("CalculateAPE") << "FitResultY2\t" << mHists["norResY"]->Fit(&funcY_2, fitOpt) << "\n";
0622         }
0623         meanY_2 = funcY_2.GetParameter(1);
0624         sigmaY_2 = funcY_2.GetParameter(2);
0625       }
0626 
0627       // Fill histograms
0628       double fitMeanX_1(meanX_1), fitMeanX_2(meanX_2);
0629       double residualWidthX_1(sigmaX_1), residualWidthX_2(sigmaX_2);
0630       double fitMeanY_1(meanY_1), fitMeanY_2(meanY_2);
0631       double residualWidthY_1(sigmaY_1), residualWidthY_2(sigmaY_2);
0632 
0633       double correctionX2_1(-0.0010), correctionX2_2(-0.0010);
0634       double correctionY2_1(-0.0010), correctionY2_2(-0.0010);
0635       correctionX2_1 = meanSigmaX * meanSigmaX * (residualWidthX_1 * residualWidthX_1 - baselineWidthX2);
0636       correctionX2_2 = meanSigmaX * meanSigmaX * (residualWidthX_2 * residualWidthX_2 - baselineWidthX2);
0637       if (i_sector.second.isPixel) {
0638         correctionY2_1 = meanSigmaY * meanSigmaY * (residualWidthY_1 * residualWidthY_1 - baselineWidthY2);
0639         correctionY2_2 = meanSigmaY * meanSigmaY * (residualWidthY_2 * residualWidthY_2 - baselineWidthY2);
0640       }
0641 
0642       double correctionX_1 = correctionX2_1 >= 0. ? std::sqrt(correctionX2_1) : -std::sqrt(-correctionX2_1);
0643       double correctionX_2 = correctionX2_2 >= 0. ? std::sqrt(correctionX2_2) : -std::sqrt(-correctionX2_2);
0644       double correctionY_1 = correctionY2_1 >= 0. ? std::sqrt(correctionY2_1) : -std::sqrt(-correctionY2_1);
0645       double correctionY_2 = correctionY2_2 >= 0. ? std::sqrt(correctionY2_2) : -std::sqrt(-correctionY2_2);
0646       // Meanwhile, this got very bad default values, or? (negative corrections allowed)
0647       if (edm::isNotFinite(correctionX_1))
0648         correctionX_1 = -0.0010;
0649       if (edm::isNotFinite(correctionX_2))
0650         correctionX_2 = -0.0010;
0651       if (edm::isNotFinite(correctionY_1))
0652         correctionY_1 = -0.0010;
0653       if (edm::isNotFinite(correctionY_2))
0654         correctionY_2 = -0.0010;
0655 
0656       if (entriesX < minHitsPerInterval) {
0657         meanX = 0.;
0658         rmsX = -0.0010;
0659         fitMeanX_1 = 0.;
0660         correctionX_1 = residualWidthX_1 = -0.0010;
0661         fitMeanX_2 = 0.;
0662         correctionX_2 = residualWidthX_2 = -0.0010;
0663       }
0664 
0665       if (entriesY < minHitsPerInterval) {
0666         meanY = 0.;
0667         rmsY = -0.0010;
0668         fitMeanY_1 = 0.;
0669         correctionY_1 = residualWidthY_1 = -0.0010;
0670         fitMeanY_2 = 0.;
0671         correctionY_2 = residualWidthY_2 = -0.0010;
0672       }
0673 
0674       i_sector.second.MeanX->SetBinContent(i_errBins.first, meanX);
0675       i_sector.second.RmsX->SetBinContent(i_errBins.first, rmsX);
0676 
0677       i_sector.second.FitMeanX1->SetBinContent(i_errBins.first, fitMeanX_1);
0678       i_sector.second.ResidualWidthX1->SetBinContent(i_errBins.first, residualWidthX_1);
0679       i_sector.second.CorrectionX1->SetBinContent(i_errBins.first, correctionX_1 * 10000.);
0680 
0681       i_sector.second.FitMeanX2->SetBinContent(i_errBins.first, fitMeanX_2);
0682       i_sector.second.ResidualWidthX2->SetBinContent(i_errBins.first, residualWidthX_2);
0683       i_sector.second.CorrectionX2->SetBinContent(i_errBins.first, correctionX_2 * 10000.);
0684 
0685       if (i_sector.second.isPixel) {
0686         i_sector.second.MeanY->SetBinContent(i_errBins.first, meanY);
0687         i_sector.second.RmsY->SetBinContent(i_errBins.first, rmsY);
0688 
0689         i_sector.second.FitMeanY1->SetBinContent(i_errBins.first, fitMeanY_1);
0690         i_sector.second.ResidualWidthY1->SetBinContent(i_errBins.first, residualWidthY_1);
0691         i_sector.second.CorrectionY1->SetBinContent(i_errBins.first, correctionY_1 * 10000.);
0692 
0693         i_sector.second.FitMeanY2->SetBinContent(i_errBins.first, fitMeanY_2);
0694         i_sector.second.ResidualWidthY2->SetBinContent(i_errBins.first, residualWidthY_2);
0695         i_sector.second.CorrectionY2->SetBinContent(i_errBins.first, correctionY_2 * 10000.);
0696       }
0697 
0698       // Use result for bin only when entries>=minHitsPerInterval
0699       if (entriesX < minHitsPerInterval && entriesY < minHitsPerInterval)
0700         continue;
0701 
0702       double weightX(0.);
0703       double weightY(0.);
0704       if (apeWeight == wUnity) {
0705         weightX = 1.;
0706         weightY = 1.;
0707       } else if (apeWeight == wEntries) {
0708         weightX = entriesX;
0709         weightY = entriesY;
0710       } else if (apeWeight == wEntriesOverSigmaX2) {
0711         weightX = entriesX / (meanSigmaX * meanSigmaX);
0712         weightY = entriesY / (meanSigmaY * meanSigmaY);
0713       }
0714 
0715       const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinX(meanSigmaX * meanSigmaX,
0716                                                                          residualWidthX_1 * residualWidthX_1);
0717       const WeightAndResultsPerBin weightAndResultsPerBinX(weightX, error2AndResidualWidth2PerBinX);
0718       if (!(entriesX < minHitsPerInterval)) {
0719         //Fill absolute weights
0720         i_sector.second.WeightX->SetBinContent(i_errBins.first, weightX);
0721         v_weightAndResultsPerBinX.push_back(weightAndResultsPerBinX);
0722       }
0723 
0724       const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinY(meanSigmaY * meanSigmaY,
0725                                                                          residualWidthY_1 * residualWidthY_1);
0726       const WeightAndResultsPerBin weightAndResultsPerBinY(weightY, error2AndResidualWidth2PerBinY);
0727       if (!(entriesY < minHitsPerInterval)) {
0728         //Fill absolute weights
0729         i_sector.second.WeightY->SetBinContent(i_errBins.first, weightY);
0730         v_weightAndResultsPerBinY.push_back(weightAndResultsPerBinY);
0731       }
0732     }
0733 
0734     // Do the final calculations
0735 
0736     if (v_weightAndResultsPerBinX.empty()) {
0737       edm::LogError("CalculateAPE") << "NO error interval of sector " << i_sector.first
0738                                     << " has a valid x APE calculated,\n...so cannot set APE";
0739       continue;
0740     }
0741 
0742     if (i_sector.second.isPixel && v_weightAndResultsPerBinY.empty()) {
0743       edm::LogError("CalculateAPE") << "NO error interval of sector " << i_sector.first
0744                                     << " has a valid y APE calculated,\n...so cannot set APE";
0745       continue;
0746     }
0747 
0748     double correctionX2(999.);
0749     double correctionY2(999.);
0750 
0751     // Get sum of all weights
0752     double weightSumX(0.);
0753     for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
0754       weightSumX += i_apeBin.first;
0755     }
0756     i_sector.second.WeightX->Scale(1 / weightSumX);
0757     double weightSumY(0.);
0758     if (i_sector.second.isPixel) {
0759       for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
0760         weightSumY += i_apeBin.first;
0761       }
0762       i_sector.second.WeightY->Scale(1 / weightSumY);
0763     }
0764 
0765     if (!setBaseline) {
0766       // Calculate weighted mean
0767       bool firstIntervalX(true);
0768       for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
0769         if (firstIntervalX) {
0770           correctionX2 = i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthX2);
0771           firstIntervalX = false;
0772         } else {
0773           correctionX2 += i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthX2);
0774         }
0775       }
0776       correctionX2 = correctionX2 / weightSumX;
0777     } else {
0778       double numeratorX(0.), denominatorX(0.);
0779       for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
0780         numeratorX += i_apeBin.first * i_apeBin.second.first * i_apeBin.second.second;
0781         denominatorX += i_apeBin.first * i_apeBin.second.first;
0782       }
0783       correctionX2 = numeratorX / denominatorX;
0784     }
0785 
0786     if (i_sector.second.isPixel) {
0787       if (!setBaseline) {
0788         // Calculate weighted mean
0789         bool firstIntervalY(true);
0790         for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
0791           if (firstIntervalY) {
0792             correctionY2 = i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthY2);
0793             firstIntervalY = false;
0794           } else {
0795             correctionY2 += i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthY2);
0796           }
0797         }
0798         correctionY2 = correctionY2 / weightSumY;
0799       } else {
0800         double numeratorY(0.), denominatorY(0.);
0801         for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
0802           numeratorY += i_apeBin.first * i_apeBin.second.first * i_apeBin.second.second;
0803           denominatorY += i_apeBin.first * i_apeBin.second.first;
0804         }
0805         correctionY2 = numeratorY / denominatorY;
0806       }
0807     }
0808 
0809     if (!setBaseline) {
0810       // Calculate updated squared APE of current iteration
0811       double apeX2(999.);
0812       double apeY2(999.);
0813 
0814       // old APE value from last iteration
0815       if (firstIter) {
0816         apeX2 = 0.;
0817         apeY2 = 0.;
0818       } else {
0819         apeX2 = a_apeSectorX[i_sector.first];
0820         apeY2 = a_apeSectorY[i_sector.first];
0821       }
0822       const double apeX2old = apeX2;
0823       const double apeY2old = apeY2;
0824 
0825       // scale APE Correction with value given in cfg (not if smoothed iteration is used)
0826       edm::LogInfo("CalculateAPE") << "Unscaled correction x for sector " << i_sector.first << " is "
0827                                    << (correctionX2 > 0. ? +1. : -1.) * std::sqrt(std::fabs(correctionX2));
0828       if (!smoothIteration || firstIter)
0829         correctionX2 = correctionX2 * correctionScaling * correctionScaling;
0830       if (i_sector.second.isPixel) {
0831         edm::LogInfo("CalculateAPE") << "Unscaled correction y for sector " << i_sector.first << " is "
0832                                      << (correctionY2 > 0. ? +1. : -1.) * std::sqrt(std::fabs(correctionY2));
0833         if (!smoothIteration || firstIter)
0834           correctionY2 = correctionY2 * correctionScaling * correctionScaling;
0835       }
0836 
0837       // new APE value
0838       // smooth iteration or not?
0839       if (apeX2 + correctionX2 < 0.)
0840         correctionX2 = -apeX2;
0841       if (apeY2 + correctionY2 < 0.)
0842         correctionY2 = -apeY2;
0843       const double apeX2new(apeX2old + correctionX2);
0844       const double apeY2new(apeY2old + correctionY2);
0845       if (!smoothIteration || firstIter) {
0846         apeX2 = apeX2new;
0847         apeY2 = apeY2new;
0848       } else {
0849         const double apeXtmp = smoothFraction * std::sqrt(apeX2new) + (1 - smoothFraction) * std::sqrt(apeX2old);
0850         const double apeYtmp = smoothFraction * std::sqrt(apeY2new) + (1 - smoothFraction) * std::sqrt(apeY2old);
0851         apeX2 = apeXtmp * apeXtmp;
0852         apeY2 = apeYtmp * apeYtmp;
0853       }
0854       if (apeX2 < 0. || apeY2 < 0.)
0855         edm::LogError("CalculateAPE") << "\n\n\tBad APE, but why???\n\n\n";
0856       a_apeSectorX[i_sector.first] = apeX2;
0857       a_apeSectorY[i_sector.first] = apeY2;
0858 
0859       // Set the calculated APE spherical for all modules of strip sectors
0860       const double apeX(std::sqrt(apeX2));
0861       const double apeY(std::sqrt(apeY2));
0862       const double apeZ(std::sqrt(0.5 * (apeX2 + apeY2)));
0863       for (auto const& i_rawId : i_sector.second.v_rawId) {
0864         if (i_sector.second.isPixel) {
0865           apeOutputFile << i_rawId << " " << std::fixed << std::setprecision(5) << apeX << " " << apeY << " " << apeZ
0866                         << "\n";
0867         } else {
0868           apeOutputFile << i_rawId << " " << std::fixed << std::setprecision(5) << apeX << " " << apeX << " " << apeX
0869                         << "\n";
0870         }
0871       }
0872     } else {  // In setBaseline mode, just fill estimated mean value of residual width
0873       a_apeSectorX[i_sector.first] = correctionX2;
0874       a_apeSectorY[i_sector.first] = correctionY2;
0875     }
0876   }
0877   if (!setBaseline)
0878     apeOutputFile.close();
0879 
0880   iterationTreeX->Fill();
0881   iterationTreeX->Write("iterTreeX",
0882                         TObject::kOverwrite);  // TObject::kOverwrite needed to not produce another iterTreeX;2
0883   iterationTreeY->Fill();
0884   iterationTreeY->Write("iterTreeY",
0885                         TObject::kOverwrite);  // TObject::kOverwrite needed to not produce another iterTreeY;2
0886   if (firstIter) {
0887     sectorNameTree->Fill();
0888     sectorNameTree->Write("nameTree");
0889     for (auto& i_sector : m_tkSector_) {
0890       delete a_sectorName[i_sector.first];
0891     }
0892   }
0893   iterationFile->Close();
0894   delete iterationFile;
0895 
0896   if (baselineFile) {
0897     baselineFile->Close();
0898     delete baselineFile;
0899   }
0900 }
0901 
0902 // ------------ method called to for each event  ------------
0903 void ApeEstimatorSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0904   // Load APEs from the GT and write them to root files similar to the ones from calculateAPE()
0905 
0906   if (firstEvent) {
0907     // Set baseline or calculate APE value?
0908     const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
0909 
0910     const AlignmentErrorsExtended* alignmentErrors = &iSetup.getData(alignmentErrorToken_);
0911 
0912     // Read in baseline file for calculation of APE value (if not setting baseline)
0913     // Has same format as iterationFile
0914     const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
0915     TFile* baselineFile(nullptr);
0916     TTree* sectorNameBaselineTree(nullptr);
0917     if (!setBaseline) {
0918       std::ifstream baselineFileStream;
0919       // Check if baseline file exists
0920       baselineFileStream.open(baselineFileName.c_str());
0921       if (baselineFileStream.is_open()) {
0922         baselineFileStream.close();
0923         baselineFile = new TFile(baselineFileName.c_str(), "READ");
0924       }
0925       if (baselineFile) {
0926         edm::LogInfo("CalculateAPE") << "Baseline file for APE values sucessfully opened";
0927         baselineFile->GetObject("nameTree;1", sectorNameBaselineTree);
0928       } else {
0929         edm::LogWarning("CalculateAPE") << "There is NO baseline file for APE values, so normalized residual width =1 "
0930                                            "for ideal conditions is assumed";
0931       }
0932     }
0933 
0934     // Set up root file for default APE values
0935     const std::string defaultFileName(parameterSet_.getParameter<std::string>("DefaultFile"));
0936     TFile* defaultFile = new TFile(defaultFileName.c_str(), "RECREATE");
0937 
0938     // Naming in the root files has to be iterTreeX to be consistent for the plotting tool
0939     TTree* defaultTreeX(nullptr);
0940     TTree* defaultTreeY(nullptr);
0941     defaultFile->GetObject("iterTreeX;1", defaultTreeX);
0942     defaultFile->GetObject("iterTreeY;1", defaultTreeY);
0943     // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as defaultTree)
0944     TTree* sectorNameTree(nullptr);
0945     defaultFile->GetObject("nameTree;1", sectorNameTree);
0946 
0947     const bool firstIter(!defaultTreeX);
0948     if (firstIter) {  // should be always true in setBaseline mode, since file is recreated
0949       defaultTreeX = new TTree("iterTreeX", "Tree for default APE x values from GT");
0950       defaultTreeY = new TTree("iterTreeY", "Tree for default APE y values from GT");
0951       edm::LogInfo("CalculateAPE") << "First APE iteration (number 0.), create default file with TTree";
0952       sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
0953     } else {
0954       edm::LogWarning("CalculateAPE") << "NOT the first APE iteration (number 0.), is this wanted or forgot to delete "
0955                                          "old iteration file with TTree?";
0956     }
0957 
0958     // Assign the information stored in the trees to arrays
0959     double a_defaultSectorX[noSectors_];
0960     double a_defaultSectorY[noSectors_];
0961 
0962     std::string* a_sectorName[noSectors_];
0963     std::string* a_sectorBaselineName[noSectors_];
0964     for (auto const& i_sector : m_tkSector_) {
0965       const unsigned int iSector(i_sector.first);
0966       const bool pixelSector(i_sector.second.isPixel);
0967 
0968       a_defaultSectorX[iSector] = -99.;
0969       a_defaultSectorY[iSector] = -99.;
0970 
0971       a_sectorName[iSector] = nullptr;
0972       a_sectorBaselineName[iSector] = nullptr;
0973       std::stringstream ss_sector, ss_sectorSuffixed;
0974       ss_sector << "Ape_Sector_" << iSector;
0975       if (!setBaseline && sectorNameBaselineTree) {
0976         sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
0977         sectorNameBaselineTree->GetEntry(0);
0978       }
0979 
0980       if (firstIter) {  // should be always true in setBaseline mode, since file is recreated
0981         ss_sectorSuffixed << ss_sector.str() << "/D";
0982         defaultTreeX->Branch(ss_sector.str().c_str(), &a_defaultSectorX[iSector], ss_sectorSuffixed.str().c_str());
0983         if (pixelSector) {
0984           defaultTreeY->Branch(ss_sector.str().c_str(), &a_defaultSectorY[iSector], ss_sectorSuffixed.str().c_str());
0985         }
0986         sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
0987       } else {
0988         defaultTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorX[iSector]);
0989         defaultTreeX->GetEntry(0);
0990         if (pixelSector) {
0991           defaultTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorY[iSector]);
0992           defaultTreeY->GetEntry(0);
0993         }
0994         sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
0995         sectorNameTree->GetEntry(0);
0996       }
0997     }
0998 
0999     // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
1000     bool failed(false);
1001     for (auto& i_sector : m_tkSector_) {
1002       const std::string& name(i_sector.second.name);
1003       if (firstIter) {
1004         a_sectorName[i_sector.first] = new std::string(name);
1005       } else {
1006         const std::string& nameLastIter(*a_sectorName[i_sector.first]);
1007         if (name != nameLastIter) {
1008           edm::LogError("CalculateAPE") << "Inconsistent sector definition in iterationFile for sector "
1009                                         << i_sector.first << ",\n"
1010                                         << "Recent iteration has name \"" << name << "\", while previous had \""
1011                                         << nameLastIter << "\"\n"
1012                                         << "...APE calculation stopped. Please check sector definitions in config!\n";
1013           failed = true;
1014         }
1015       }
1016       if (!setBaseline && baselineFile) {
1017         const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
1018         if (name != nameBaseline) {
1019           edm::LogError("CalculateAPE") << "Inconsistent sector definition in baselineFile for sector "
1020                                         << i_sector.first << ",\n"
1021                                         << "Recent iteration has name \"" << name << "\", while baseline had \""
1022                                         << nameBaseline << "\"\n"
1023                                         << "...APE calculation stopped. Please check sector definitions in config!\n";
1024           failed = true;
1025         }
1026       }
1027     }
1028 
1029     if (failed) {
1030       if (firstIter) {
1031         delete defaultTreeX;
1032         delete defaultTreeY;
1033         delete sectorNameTree;
1034         for (auto& i_sector : m_tkSector_) {
1035           delete a_sectorName[i_sector.first];
1036         }
1037       }
1038       return;
1039     }
1040 
1041     // Loop over sectors for calculating getting default APE
1042     for (auto& i_sector : m_tkSector_) {
1043       double defaultApeX(0.);
1044       double defaultApeY(0.);
1045       unsigned int nModules(0);
1046       for (auto const& i_rawId : i_sector.second.v_rawId) {
1047         std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
1048         for (std::vector<AlignTransformErrorExtended>::const_iterator i_alignError = alignErrors.begin();
1049              i_alignError != alignErrors.end();
1050              ++i_alignError) {
1051           if (i_rawId == i_alignError->rawId()) {
1052             CLHEP::HepSymMatrix errMatrix = i_alignError->matrix();
1053             defaultApeX += errMatrix[0][0];
1054             defaultApeY += errMatrix[1][1];
1055             nModules++;
1056           }
1057         }
1058       }
1059       a_defaultSectorX[i_sector.first] = defaultApeX / nModules;
1060       a_defaultSectorY[i_sector.first] = defaultApeY / nModules;
1061     }
1062     sectorNameTree->Fill();
1063     sectorNameTree->Write("nameTree");
1064     defaultTreeX->Fill();
1065     defaultTreeX->Write("iterTreeX");
1066     defaultTreeY->Fill();
1067     defaultTreeY->Write("iterTreeY");
1068 
1069     defaultFile->Close();
1070     if (baselineFile) {
1071       baselineFile->Close();
1072       delete baselineFile;
1073     }
1074 
1075     delete defaultFile;
1076     if (firstIter) {
1077       for (auto& i_sector : m_tkSector_) {
1078         delete a_sectorName[i_sector.first];
1079       }
1080     }
1081 
1082     firstEvent = false;
1083   }
1084 }
1085 
1086 // ------------ method called once each job just before starting event loop  ------------
1087 void ApeEstimatorSummary::beginJob() {
1088   this->openInputFile();
1089 
1090   this->getTrackerSectorStructs();
1091 
1092   this->bookHists();
1093 }
1094 
1095 // ------------ method called once each job just after ending the event loop  ------------
1096 void ApeEstimatorSummary::endJob() {
1097   this->calculateApe();
1098 
1099   this->writeHists();
1100 }
1101 
1102 //define this as a plug-in
1103 DEFINE_FWK_MODULE(ApeEstimatorSummary);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!