CMS 3D CMS Logo
ApeEstimatorSummary.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: ApeEstimatorSummary
4 // Class: ApeEstimatorSummary
5 //
13 //
14 // Original Author: Johannes Hauk,6 2-039,+41227673512,
15 // Created: Mon Oct 11 13:44:03 CEST 2010
16 // Modified by: Christian Schomakers (RWTH Aachen)
17 // $Id: ApeEstimatorSummary.cc,v 1.14 2012/01/26 00:06:33 hauk Exp $
18 //
19 //
20 
21 // system include files
22 #include <memory>
23 #include <fstream>
24 #include <sstream>
25 
26 // user include files
27 #include "FWCore/Framework/interface/Frameworkfwd.h"
28 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
29 #include "FWCore/Framework/interface/ESHandle.h"
30 #include "FWCore/Framework/interface/EventSetup.h"
31 #include "FWCore/Framework/interface/Event.h"
32 #include "FWCore/Framework/interface/MakerMacros.h"
33 
34 #include "FWCore/ParameterSet/interface/ParameterSet.h"
35 
36 #include "FWCore/MessageLogger/interface/MessageLogger.h"
37 #include "FWCore/Utilities/interface/EDMException.h"
38 
39 #include "CLHEP/Matrix/SymMatrix.h"
40 
41 //#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
42 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
43 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
44 
45 //...............
46 #include "Alignment/APEEstimation/interface/TrackerSectorStruct.h"
47 
48 #include "TH1.h"
49 #include "TString.h"
50 #include "TFile.h"
51 #include "TDirectory.h"
52 #include "TTree.h"
53 #include "TF1.h"
54 #include "TMath.h"
55 //
56 // class declaration
57 //
58 
59 class ApeEstimatorSummary : public edm::one::EDAnalyzer<> {
60 public:
61  explicit ApeEstimatorSummary(const edm::ParameterSet&);
62  ~ApeEstimatorSummary() override;
63 
64 private:
65  void beginJob() override;
66  void analyze(const edm::Event&, const edm::EventSetup&) override;
67  void endJob() override;
68 
69  void openInputFile();
70  void getTrackerSectorStructs();
71  void bookHists();
72  std::vector<double> residualErrorBinning();
73  void writeHists();
74 
75  void calculateApe();
76 
77  enum ApeWeight { wInvalid, wUnity, wEntries, wEntriesOverSigmaX2 };
78 
79  // ----------member data ---------------------------
80  const edm::ParameterSet parameterSet_;
81 
82  bool firstEvent = true;
83 
84  TFile* inputFile_;
85 
86  std::map<unsigned int, TrackerSectorStruct> m_tkSector_;
87  unsigned int noSectors_;
88 };
89 
90 //
91 // constants, enums and typedefs
92 //
93 
94 //
95 // static data member definitions
96 //
97 
98 //
99 // constructors and destructor
100 //
101 ApeEstimatorSummary::ApeEstimatorSummary(const edm::ParameterSet& iConfig)
102  : parameterSet_(iConfig), inputFile_(nullptr) {}
103 
104 ApeEstimatorSummary::~ApeEstimatorSummary() {}
105 
106 //
107 // member functions
108 //
109 
110 void ApeEstimatorSummary::openInputFile() {
111  const std::string inputFileName(parameterSet_.getParameter<std::string>("InputFile"));
112  std::ifstream inputFileStream;
113  // Check if baseline file exists
114  inputFileStream.open(inputFileName.c_str());
115  if (inputFileStream.is_open()) {
116  inputFile_ = new TFile(inputFileName.c_str(), "READ");
117  }
118  if (inputFile_) {
119  edm::LogInfo("CalculateAPE") << "Input file from loop over tracks and corresponding hits sucessfully opened";
120  } else {
121  edm::LogError("CalculateAPE") << "There is NO input file\n"
122  << "...APE calculation stopped. Please check path of input file name in config:\n"
123  << "\t" << inputFileName;
124  throw edm::Exception(edm::errors::Configuration, "Bad input file name");
125  }
126 }
127 
128 void ApeEstimatorSummary::getTrackerSectorStructs() {
129  // At present first of the plugins registered in TFileService needs to be the one containing the normalized residual histos per sector per error bin
130  TString pluginName(inputFile_->GetListOfKeys()->At(0)->GetName());
131 
132  pluginName += "/";
133  TDirectory *sectorDir(nullptr), *intervalDir(nullptr);
134  bool sectorBool(true);
135  unsigned int iSector(1);
136  for (; sectorBool; ++iSector) {
137  std::stringstream sectorName, fullSectorName;
138  sectorName << "Sector_" << iSector << "/";
139  fullSectorName << pluginName << sectorName.str();
140  TString fullName(fullSectorName.str().c_str());
141  inputFile_->cd();
142  sectorDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
143  if (sectorDir) {
144  TrackerSectorStruct tkSector;
145  sectorDir->GetObject("z_name;1", tkSector.Name);
146  tkSector.name = tkSector.Name->GetTitle();
147  bool intervalBool(true);
148  for (unsigned int iInterval(1); intervalBool; ++iInterval) {
149  std::stringstream intervalName, fullIntervalName;
150  intervalName << "Interval_" << iInterval << "/";
151  fullIntervalName << fullSectorName.str() << intervalName.str();
152  fullName = fullIntervalName.str().c_str();
153  intervalDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
154  if (intervalDir) {
155  intervalDir->GetObject("h_sigmaX;1", tkSector.m_binnedHists[iInterval]["sigmaX"]);
156  intervalDir->GetObject("h_norResX;1", tkSector.m_binnedHists[iInterval]["norResX"]);
157  intervalDir->GetObject("h_sigmaY;1", tkSector.m_binnedHists[iInterval]["sigmaY"]);
158  intervalDir->GetObject("h_norResY;1", tkSector.m_binnedHists[iInterval]["norResY"]);
159  if (tkSector.m_binnedHists[iInterval]["sigmaY"] && tkSector.m_binnedHists[iInterval]["norResY"]) {
160  tkSector.isPixel = true;
161  }
162  } else {
163  intervalBool = false;
164  if (iSector == 1)
165  edm::LogInfo("CalculateAPE") << "There are " << iInterval - 1
166  << " intervals per sector defined in input file";
167  }
168  }
169  TDirectory* resultsDir(nullptr);
170  std::stringstream fullResultName;
171  fullResultName << fullSectorName.str() << "Results/";
172  fullName = fullResultName.str().c_str();
173  resultsDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
174  if (resultsDir) {
175  resultsDir->GetObject("h_entriesX;1", tkSector.EntriesX);
176  if (tkSector.isPixel)
177  resultsDir->GetObject("h_entriesY;1", tkSector.EntriesY);
178  TTree* rawIdTree(nullptr);
179  resultsDir->GetObject("rawIdTree", rawIdTree);
180  unsigned int rawId(0);
181  rawIdTree->SetBranchAddress("RawId", &rawId);
182  for (int entry = 0; entry < rawIdTree->GetEntries(); ++entry) {
183  rawIdTree->GetEntry(entry);
184  // command "hadd" adds entries in TTree, so rawId are existing as often as number of files are added
185  bool alreadyAdded(false);
186  for (auto const& i_rawId : tkSector.v_rawId) {
187  if (rawId == i_rawId)
188  alreadyAdded = true;
189  }
190  if (alreadyAdded)
191  break;
192  tkSector.v_rawId.push_back(rawId);
193  }
194  }
195  m_tkSector_[iSector] = tkSector;
196  } else {
197  sectorBool = false;
198  edm::LogInfo("CalculateAPE") << "There are " << iSector - 1 << " sectors defined in input file";
199  }
200  }
201  noSectors_ = iSector + 1;
202 }
203 
204 void ApeEstimatorSummary::bookHists() {
205  const std::vector<double> v_binX(this->residualErrorBinning());
206  for (auto& i_sector : m_tkSector_) {
207  i_sector.second.WeightX = new TH1F("h_weightX",
208  "relative weight w_{x}/w_{tot,x};#sigma_{x} [#mum];w_{x}/w_{tot,x}",
209  v_binX.size() - 1,
210  &(v_binX[0]));
211  i_sector.second.MeanX = new TH1F("h_meanX",
212  "residual mean <r_{x}/#sigma_{r,x}>;#sigma_{x} [#mum];<r_{x}/#sigma_{r,x}>",
213  v_binX.size() - 1,
214  &(v_binX[0]));
215  i_sector.second.RmsX = new TH1F("h_rmsX",
216  "residual rms RMS(r_{x}/#sigma_{r,x});#sigma_{x} [#mum];RMS(r_{x}/#sigma_{r,x})",
217  v_binX.size() - 1,
218  &(v_binX[0]));
219  i_sector.second.FitMeanX1 = new TH1F(
220  "h_fitMeanX1", "fitted residual mean #mu_{x};#sigma_{x} [#mum];#mu_{x}", v_binX.size() - 1, &(v_binX[0]));
221  i_sector.second.ResidualWidthX1 = new TH1F("h_residualWidthX1",
222  "residual width #Delta_{x};#sigma_{x} [#mum];#Delta_{x}",
223  v_binX.size() - 1,
224  &(v_binX[0]));
225  i_sector.second.CorrectionX1 = new TH1F("h_correctionX1",
226  "correction to APE_{x};#sigma_{x} [#mum];#Delta#sigma_{align,x} [#mum]",
227  v_binX.size() - 1,
228  &(v_binX[0]));
229  i_sector.second.FitMeanX2 = new TH1F(
230  "h_fitMeanX2", "fitted residual mean #mu_{x};#sigma_{x} [#mum];#mu_{x}", v_binX.size() - 1, &(v_binX[0]));
231  i_sector.second.ResidualWidthX2 = new TH1F("h_residualWidthX2",
232  "residual width #Delta_{x};#sigma_{x} [#mum];#Delta_{x}",
233  v_binX.size() - 1,
234  &(v_binX[0]));
235  i_sector.second.CorrectionX2 = new TH1F("h_correctionX2",
236  "correction to APE_{x};#sigma_{x} [#mum];#Delta#sigma_{align,x} [#mum]",
237  v_binX.size() - 1,
238  &(v_binX[0]));
239 
240  if (i_sector.second.isPixel) {
241  i_sector.second.WeightY = new TH1F("h_weightY",
242  "relative weight w_{y}/w_{tot,y};#sigma_{y} [#mum];w_{y}/w_{tot,y}",
243  v_binX.size() - 1,
244  &(v_binX[0]));
245  i_sector.second.MeanY = new TH1F("h_meanY",
246  "residual mean <r_{y}/#sigma_{r,y}>;#sigma_{y} [#mum];<r_{y}/#sigma_{r,y}>",
247  v_binX.size() - 1,
248  &(v_binX[0]));
249  i_sector.second.RmsY = new TH1F("h_rmsY",
250  "residual rms RMS(r_{y}/#sigma_{r,y});#sigma_{y} [#mum];RMS(r_{y}/#sigma_{r,y})",
251  v_binX.size() - 1,
252  &(v_binX[0]));
253  i_sector.second.FitMeanY1 = new TH1F(
254  "h_fitMeanY1", "fitted residual mean #mu_{y};#sigma_{y} [#mum];#mu_{y}", v_binX.size() - 1, &(v_binX[0]));
255  i_sector.second.ResidualWidthY1 = new TH1F("h_residualWidthY1",
256  "residual width #Delta_{y};#sigma_{y} [#mum];#Delta_{y}",
257  v_binX.size() - 1,
258  &(v_binX[0]));
259  i_sector.second.CorrectionY1 = new TH1F("h_correctionY1",
260  "correction to APE_{y};#sigma_{y} [#mum];#Delta#sigma_{align,y} [#mum]",
261  v_binX.size() - 1,
262  &(v_binX[0]));
263  i_sector.second.FitMeanY2 = new TH1F(
264  "h_fitMeanY2", "fitted residual mean #mu_{y};#sigma_{y} [#mum];#mu_{y}", v_binX.size() - 1, &(v_binX[0]));
265  i_sector.second.ResidualWidthY2 = new TH1F("h_residualWidthY2",
266  "residual width #Delta_{y};#sigma_{y} [#mum];#Delta_{y}",
267  v_binX.size() - 1,
268  &(v_binX[0]));
269  i_sector.second.CorrectionY2 = new TH1F("h_correctionY2",
270  "correction to APE_{y};#sigma_{y} [#mum];#Delta#sigma_{align,y} [#mum]",
271  v_binX.size() - 1,
272  &(v_binX[0]));
273  }
274  }
275 }
276 
277 std::vector<double> ApeEstimatorSummary::residualErrorBinning() {
278  std::vector<double> v_binX;
279  TH1* EntriesX(m_tkSector_[1].EntriesX);
280  for (int iBin = 1; iBin <= EntriesX->GetNbinsX() + 1; ++iBin) {
281  v_binX.push_back(EntriesX->GetBinLowEdge(iBin));
282  }
283  return v_binX;
284 }
285 
286 void ApeEstimatorSummary::writeHists() {
287  TFile* resultsFile = new TFile(parameterSet_.getParameter<std::string>("ResultsFile").c_str(), "RECREATE");
288  TDirectory* baseDir = resultsFile->mkdir("ApeEstimatorSummary");
289  for (auto const& i_sector : m_tkSector_) {
290  std::stringstream dirName;
291  dirName << "Sector_" << i_sector.first;
292  TDirectory* dir = baseDir->mkdir(dirName.str().c_str());
293  dir->cd();
294 
295  i_sector.second.Name->Write();
296 
297  i_sector.second.WeightX->Write();
298  i_sector.second.MeanX->Write();
299  i_sector.second.RmsX->Write();
300  i_sector.second.FitMeanX1->Write();
301  i_sector.second.ResidualWidthX1->Write();
302  i_sector.second.CorrectionX1->Write();
303  i_sector.second.FitMeanX2->Write();
304  i_sector.second.ResidualWidthX2->Write();
305  i_sector.second.CorrectionX2->Write();
306 
307  if (i_sector.second.isPixel) {
308  i_sector.second.WeightY->Write();
309  i_sector.second.MeanY->Write();
310  i_sector.second.RmsY->Write();
311  i_sector.second.FitMeanY1->Write();
312  i_sector.second.ResidualWidthY1->Write();
313  i_sector.second.CorrectionY1->Write();
314  i_sector.second.FitMeanY2->Write();
315  i_sector.second.ResidualWidthY2->Write();
316  i_sector.second.CorrectionY2->Write();
317  }
318  }
319  resultsFile->Close();
320 }
321 
322 void ApeEstimatorSummary::calculateApe() {
323  // Set baseline or calculate APE value?
324  const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
325 
326  // Read in baseline file for calculation of APE value (if not setting baseline)
327  // Has same format as iterationFile
328  const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
329  TFile* baselineFile(nullptr);
330  TTree* baselineTreeX(nullptr);
331  TTree* baselineTreeY(nullptr);
332  TTree* sectorNameBaselineTree(nullptr);
333  if (!setBaseline) {
334  std::ifstream baselineFileStream;
335  // Check if baseline file exists
336  baselineFileStream.open(baselineFileName.c_str());
337  if (baselineFileStream.is_open()) {
338  baselineFileStream.close();
339  baselineFile = new TFile(baselineFileName.c_str(), "READ");
340  }
341  if (baselineFile) {
342  edm::LogInfo("CalculateAPE") << "Baseline file for APE values sucessfully opened";
343  baselineFile->GetObject("iterTreeX;1", baselineTreeX);
344  baselineFile->GetObject("iterTreeY;1", baselineTreeY);
345  baselineFile->GetObject("nameTree;1", sectorNameBaselineTree);
346  } else {
347  edm::LogWarning("CalculateAPE") << "There is NO baseline file for APE values, so normalized residual width =1 "
348  "for ideal conditions is assumed";
349  }
350  }
351 
352  // Set up root file for iterations on APE value (or for setting baseline in setBaseline mode)
353  const std::string iterationFileName(setBaseline ? baselineFileName
354  : parameterSet_.getParameter<std::string>("IterationFile"));
355  // For iterations, updates are needed to not overwrite the iterations before
356  TFile* iterationFile = new TFile(iterationFileName.c_str(), setBaseline ? "RECREATE" : "UPDATE");
357 
358  // Set up TTree for iterative APE values on first pass (first iteration) or read from file (further iterations)
359  TTree* iterationTreeX(nullptr);
360  TTree* iterationTreeY(nullptr);
361  iterationFile->GetObject("iterTreeX;1", iterationTreeX);
362  iterationFile->GetObject("iterTreeY;1", iterationTreeY);
363  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as iterationTree)
364  TTree* sectorNameTree(nullptr);
365  iterationFile->GetObject("nameTree;1", sectorNameTree);
366 
367  bool firstIter(false);
368  if (!iterationTreeX) { // should be always true in setBaseline mode, since file is recreated
369  firstIter = true;
370  if (!setBaseline) {
371  iterationTreeX = new TTree("iterTreeX", "Tree for APE x values of all iterations");
372  iterationTreeY = new TTree("iterTreeY", "Tree for APE y values of all iterations");
373  edm::LogInfo("CalculateAPE") << "First APE iteration (number 0.), create iteration file with TTree";
374  sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
375  } else {
376  iterationTreeX = new TTree("iterTreeX", "Tree for baseline x values of normalized residual width");
377  iterationTreeY = new TTree("iterTreeY", "Tree for baseline y values of normalized residual width");
378  edm::LogInfo("CalculateAPE") << "Set baseline, create baseline file with TTree";
379  sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
380  }
381  } else {
382  const unsigned int iteration(iterationTreeX->GetEntries());
383  edm::LogWarning("CalculateAPE") << "NOT the first APE iteration (number 0.) but the " << iteration
384  << ". one, is this wanted or forgot to delete old iteration file with TTree?";
385  }
386 
387  // Assign the information stored in the trees to arrays
388  double a_apeSectorX[noSectors_];
389  double a_apeSectorY[noSectors_];
390  double a_baselineSectorX[noSectors_];
391  double a_baselineSectorY[noSectors_];
392 
393  std::string* a_sectorName[noSectors_];
394  std::string* a_sectorBaselineName[noSectors_];
395  for (auto const& i_sector : m_tkSector_) {
396  const unsigned int iSector(i_sector.first);
397  const bool pixelSector(i_sector.second.isPixel);
398  a_apeSectorX[iSector] = 99.;
399  a_apeSectorY[iSector] = 99.;
400  a_baselineSectorX[iSector] = -99.;
401  a_baselineSectorY[iSector] = -99.;
402 
403  a_sectorName[iSector] = nullptr;
404  a_sectorBaselineName[iSector] = nullptr;
405  std::stringstream ss_sector, ss_sectorSuffixed;
406  ss_sector << "Ape_Sector_" << iSector;
407  if (!setBaseline && baselineTreeX) {
408  baselineTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorX[iSector]);
409  baselineTreeX->GetEntry(0);
410  if (pixelSector) {
411  baselineTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorY[iSector]);
412  baselineTreeY->GetEntry(0);
413  }
414  sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
415  sectorNameBaselineTree->GetEntry(0);
416  } else {
417  // Set default ideal normalized residual width to 1
418  a_baselineSectorX[iSector] = 1.;
419  a_baselineSectorY[iSector] = 1.;
420  }
421  if (firstIter) { // should be always true in setBaseline mode, since file is recreated
422  ss_sectorSuffixed << ss_sector.str() << "/D";
423  iterationTreeX->Branch(ss_sector.str().c_str(), &a_apeSectorX[iSector], ss_sectorSuffixed.str().c_str());
424  if (pixelSector) {
425  iterationTreeY->Branch(ss_sector.str().c_str(), &a_apeSectorY[iSector], ss_sectorSuffixed.str().c_str());
426  }
427  sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
428  } else {
429  iterationTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorX[iSector]);
430  iterationTreeX->GetEntry(iterationTreeX->GetEntries() - 1);
431  if (pixelSector) {
432  iterationTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorY[iSector]);
433  iterationTreeY->GetEntry(iterationTreeY->GetEntries() - 1);
434  }
435  sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
436  sectorNameTree->GetEntry(0);
437  }
438  }
439 
440  // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
441  for (auto& i_sector : m_tkSector_) {
442  const std::string& name(i_sector.second.name);
443  if (!firstIter) {
444  const std::string& nameLastIter(*a_sectorName[i_sector.first]);
445  if (name != nameLastIter) {
446  edm::LogError("CalculateAPE") << "Inconsistent sector definition in iterationFile for sector " << i_sector.first
447  << ",\n"
448  << "Recent iteration has name \"" << name << "\", while previous had \""
449  << nameLastIter << "\"\n"
450  << "...APE calculation stopped. Please check sector definitions in config!\n";
451  return;
452  }
453  } else {
454  a_sectorName[i_sector.first] = new std::string(name);
455  }
456  if (!setBaseline && baselineFile) {
457  const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
458  if (name != nameBaseline) {
459  edm::LogError("CalculateAPE") << "Inconsistent sector definition in baselineFile for sector " << i_sector.first
460  << ",\n"
461  << "Recent iteration has name \"" << name << "\", while baseline had \""
462  << nameBaseline << "\"\n"
463  << "...APE calculation stopped. Please check sector definitions in config!\n";
464  return;
465  }
466  }
467  }
468 
469  // Set up text file for writing out APE values per module
470  std::ofstream apeOutputFile;
471  if (!setBaseline) {
472  const std::string apeOutputFileName(parameterSet_.getParameter<std::string>("ApeOutputFile"));
473  apeOutputFile.open(apeOutputFileName.c_str());
474  if (apeOutputFile.is_open()) {
475  edm::LogInfo("CalculateAPE") << "Text file for writing APE values successfully opened";
476  } else {
477  edm::LogError("CalculateAPE") << "Text file for writing APE values NOT opened,\n"
478  << "...APE calculation stopped. Please check path of text file name in config:\n"
479  << "\t" << apeOutputFileName;
480  return;
481  }
482  }
483 
484  // Loop over sectors for calculating APE
485  const std::string apeWeightName(parameterSet_.getParameter<std::string>("apeWeight"));
486  ApeWeight apeWeight(wInvalid);
487  if (apeWeightName == "unity")
488  apeWeight = wUnity;
489  else if (apeWeightName == "entries")
490  apeWeight = wEntries;
491  else if (apeWeightName == "entriesOverSigmaX2")
492  apeWeight = wEntriesOverSigmaX2;
493  if (apeWeight == wInvalid) {
494  edm::LogError("CalculateAPE") << "Invalid parameter 'apeWeight' in cfg file: \"" << apeWeightName
495  << "\"\nimplemented apeWeights are \"unity\", \"entries\", \"entriesOverSigmaX2\""
496  << "\n...APE calculation stopped.";
497  return;
498  }
499  const double minHitsPerInterval(parameterSet_.getParameter<double>("minHitsPerInterval"));
500  const double sigmaFactorFit(parameterSet_.getParameter<double>("sigmaFactorFit"));
501  const double correctionScaling(parameterSet_.getParameter<double>("correctionScaling"));
502  const bool smoothIteration(parameterSet_.getParameter<bool>("smoothIteration"));
503  const double smoothFraction(parameterSet_.getParameter<double>("smoothFraction"));
504  if (smoothFraction <= 0. || smoothFraction > 1.) {
505  edm::LogError("CalculateAPE") << "Incorrect parameter in cfg file,"
506  << "\nsmoothFraction has to be in [0,1], but is " << smoothFraction
507  << "\n...APE calculation stopped.";
508  return;
509  }
510  for (auto& i_sector : m_tkSector_) {
511  typedef std::pair<double, double> Error2AndResidualWidth2PerBin;
512  typedef std::pair<double, Error2AndResidualWidth2PerBin> WeightAndResultsPerBin;
513  std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinX;
514  std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinY;
515 
516  double baselineWidthX2(a_baselineSectorX[i_sector.first]);
517  double baselineWidthY2(a_baselineSectorY[i_sector.first]);
518 
519  // Loop over residual error bins to calculate APE for every bin
520  for (auto const& i_errBins : i_sector.second.m_binnedHists) {
521  std::map<std::string, TH1*> mHists = i_errBins.second;
522 
523  double entriesX = mHists["sigmaX"]->GetEntries();
524  double meanSigmaX = mHists["sigmaX"]->GetMean();
525 
526  // Fitting Parameters
527  double xMin = mHists["norResX"]->GetXaxis()->GetXmin();
528  double xMax = mHists["norResX"]->GetXaxis()->GetXmax();
529  double integralX = mHists["norResX"]->Integral();
530  double meanX = mHists["norResX"]->GetMean();
531  double rmsX = mHists["norResX"]->GetRMS();
532  double maximumX = mHists["norResX"]->GetBinContent(mHists["norResX"]->GetMaximumBin());
533 
534  // First Gaus Fit
535  TF1 funcX_1("mygausX", "gaus", xMin, xMax);
536  funcX_1.SetParameters(maximumX, meanX, rmsX);
537  TString fitOpt("ILERQ"); //TString fitOpt("IMR"); ("IRLL"); ("IRQ");
538  if (integralX > minHitsPerInterval) {
539  if (mHists["norResX"]->Fit(&funcX_1, fitOpt)) {
540  edm::LogWarning("CalculateAPE") << "Fit1 did not work: " << mHists["norResX"]->Fit(&funcX_1, fitOpt);
541  continue;
542  }
543  LogDebug("CalculateAPE") << "FitResultX1\t" << mHists["norResX"]->Fit(&funcX_1, fitOpt) << "\n";
544  }
545  double meanX_1 = funcX_1.GetParameter(1);
546  double sigmaX_1 = funcX_1.GetParameter(2);
547 
548  // Second gaus fit
549  TF1 funcX_2("mygausX2",
550  "gaus",
551  meanX_1 - sigmaFactorFit * TMath::Abs(sigmaX_1),
552  meanX_1 + sigmaFactorFit * TMath::Abs(sigmaX_1));
553  funcX_2.SetParameters(funcX_1.GetParameter(0), meanX_1, sigmaX_1);
554  if (integralX > minHitsPerInterval) {
555  if (mHists["norResX"]->Fit(&funcX_2, fitOpt)) {
556  edm::LogWarning("CalculateAPE") << "Fit2 did not work for x : " << mHists["norResX"]->Fit(&funcX_2, fitOpt);
557  continue;
558  }
559  LogDebug("CalculateAPE") << "FitResultX2\t" << mHists["norResX"]->Fit(&funcX_2, fitOpt) << "\n";
560  }
561  double meanX_2 = funcX_2.GetParameter(1);
562  double sigmaX_2 = funcX_2.GetParameter(2);
563 
564  // Now the same for y coordinate
565  double entriesY(0.);
566  double meanSigmaY(0.);
567  if (i_sector.second.isPixel) {
568  entriesY = mHists["sigmaY"]->GetEntries();
569  meanSigmaY = mHists["sigmaY"]->GetMean();
570  }
571 
572  double meanY_1(0.);
573  double sigmaY_1(0.);
574  double meanY_2(0.);
575  double sigmaY_2(0.);
576  double meanY(0.);
577  double rmsY(0.);
578  if (i_sector.second.isPixel) {
579  // Fitting Parameters
580  double yMin = mHists["norResY"]->GetXaxis()->GetXmin();
581  double yMax = mHists["norResY"]->GetXaxis()->GetXmax();
582  double integralY = mHists["norResY"]->Integral();
583  meanY = mHists["norResY"]->GetMean();
584  rmsY = mHists["norResY"]->GetRMS();
585  double maximumY = mHists["norResY"]->GetBinContent(mHists["norResY"]->GetMaximumBin());
586 
587  // First Gaus Fit
588  TF1 funcY_1("mygausY", "gaus", yMin, yMax);
589  funcY_1.SetParameters(maximumY, meanY, rmsY);
590  if (integralY > minHitsPerInterval) {
591  if (mHists["norResY"]->Fit(&funcY_1, fitOpt)) {
592  edm::LogWarning("CalculateAPE") << "Fit1 did not work: " << mHists["norResY"]->Fit(&funcY_1, fitOpt);
593  continue;
594  }
595  LogDebug("CalculateAPE") << "FitResultY1\t" << mHists["norResY"]->Fit(&funcY_1, fitOpt) << "\n";
596  }
597  meanY_1 = funcY_1.GetParameter(1);
598  sigmaY_1 = funcY_1.GetParameter(2);
599 
600  // Second gaus fit
601  TF1 funcY_2("mygausY2",
602  "gaus",
603  meanY_1 - sigmaFactorFit * TMath::Abs(sigmaY_1),
604  meanY_1 + sigmaFactorFit * TMath::Abs(sigmaY_1));
605  funcY_2.SetParameters(funcY_1.GetParameter(0), meanY_1, sigmaY_1);
606  if (integralY > minHitsPerInterval) {
607  if (mHists["norResY"]->Fit(&funcY_2, fitOpt)) {
608  edm::LogWarning("CalculateAPE") << "Fit2 did not work for y : " << mHists["norResY"]->Fit(&funcY_2, fitOpt);
609  continue;
610  }
611  LogDebug("CalculateAPE") << "FitResultY2\t" << mHists["norResY"]->Fit(&funcY_2, fitOpt) << "\n";
612  }
613  meanY_2 = funcY_2.GetParameter(1);
614  sigmaY_2 = funcY_2.GetParameter(2);
615  }
616 
617  // Fill histograms
618  double fitMeanX_1(meanX_1), fitMeanX_2(meanX_2);
619  double residualWidthX_1(sigmaX_1), residualWidthX_2(sigmaX_2);
620  double fitMeanY_1(meanY_1), fitMeanY_2(meanY_2);
621  double residualWidthY_1(sigmaY_1), residualWidthY_2(sigmaY_2);
622 
623  double correctionX2_1(-0.0010), correctionX2_2(-0.0010);
624  double correctionY2_1(-0.0010), correctionY2_2(-0.0010);
625  correctionX2_1 = meanSigmaX * meanSigmaX * (residualWidthX_1 * residualWidthX_1 - baselineWidthX2);
626  correctionX2_2 = meanSigmaX * meanSigmaX * (residualWidthX_2 * residualWidthX_2 - baselineWidthX2);
627  if (i_sector.second.isPixel) {
628  correctionY2_1 = meanSigmaY * meanSigmaY * (residualWidthY_1 * residualWidthY_1 - baselineWidthY2);
629  correctionY2_2 = meanSigmaY * meanSigmaY * (residualWidthY_2 * residualWidthY_2 - baselineWidthY2);
630  }
631 
632  double correctionX_1 = correctionX2_1 >= 0. ? std::sqrt(correctionX2_1) : -std::sqrt(-correctionX2_1);
633  double correctionX_2 = correctionX2_2 >= 0. ? std::sqrt(correctionX2_2) : -std::sqrt(-correctionX2_2);
634  double correctionY_1 = correctionY2_1 >= 0. ? std::sqrt(correctionY2_1) : -std::sqrt(-correctionY2_1);
635  double correctionY_2 = correctionY2_2 >= 0. ? std::sqrt(correctionY2_2) : -std::sqrt(-correctionY2_2);
636  // Meanwhile, this got very bad default values, or? (negative corrections allowed)
637  if (isnan(correctionX_1))
638  correctionX_1 = -0.0010;
639  if (isnan(correctionX_2))
640  correctionX_2 = -0.0010;
641  if (isnan(correctionY_1))
642  correctionY_1 = -0.0010;
643  if (isnan(correctionY_2))
644  correctionY_2 = -0.0010;
645 
646  if (entriesX < minHitsPerInterval) {
647  meanX = 0.;
648  rmsX = -0.0010;
649  fitMeanX_1 = 0.;
650  correctionX_1 = residualWidthX_1 = -0.0010;
651  fitMeanX_2 = 0.;
652  correctionX_2 = residualWidthX_2 = -0.0010;
653  }
654 
655  if (entriesY < minHitsPerInterval) {
656  meanY = 0.;
657  rmsY = -0.0010;
658  fitMeanY_1 = 0.;
659  correctionY_1 = residualWidthY_1 = -0.0010;
660  fitMeanY_2 = 0.;
661  correctionY_2 = residualWidthY_2 = -0.0010;
662  }
663 
664  i_sector.second.MeanX->SetBinContent(i_errBins.first, meanX);
665  i_sector.second.RmsX->SetBinContent(i_errBins.first, rmsX);
666 
667  i_sector.second.FitMeanX1->SetBinContent(i_errBins.first, fitMeanX_1);
668  i_sector.second.ResidualWidthX1->SetBinContent(i_errBins.first, residualWidthX_1);
669  i_sector.second.CorrectionX1->SetBinContent(i_errBins.first, correctionX_1 * 10000.);
670 
671  i_sector.second.FitMeanX2->SetBinContent(i_errBins.first, fitMeanX_2);
672  i_sector.second.ResidualWidthX2->SetBinContent(i_errBins.first, residualWidthX_2);
673  i_sector.second.CorrectionX2->SetBinContent(i_errBins.first, correctionX_2 * 10000.);
674 
675  if (i_sector.second.isPixel) {
676  i_sector.second.MeanY->SetBinContent(i_errBins.first, meanY);
677  i_sector.second.RmsY->SetBinContent(i_errBins.first, rmsY);
678 
679  i_sector.second.FitMeanY1->SetBinContent(i_errBins.first, fitMeanY_1);
680  i_sector.second.ResidualWidthY1->SetBinContent(i_errBins.first, residualWidthY_1);
681  i_sector.second.CorrectionY1->SetBinContent(i_errBins.first, correctionY_1 * 10000.);
682 
683  i_sector.second.FitMeanY2->SetBinContent(i_errBins.first, fitMeanY_2);
684  i_sector.second.ResidualWidthY2->SetBinContent(i_errBins.first, residualWidthY_2);
685  i_sector.second.CorrectionY2->SetBinContent(i_errBins.first, correctionY_2 * 10000.);
686  }
687 
688  // Use result for bin only when entries>=minHitsPerInterval
689  if (entriesX < minHitsPerInterval && entriesY < minHitsPerInterval)
690  continue;
691 
692  double weightX(0.);
693  double weightY(0.);
694  if (apeWeight == wUnity) {
695  weightX = 1.;
696  weightY = 1.;
697  } else if (apeWeight == wEntries) {
698  weightX = entriesX;
699  weightY = entriesY;
700  } else if (apeWeight == wEntriesOverSigmaX2) {
701  weightX = entriesX / (meanSigmaX * meanSigmaX);
702  weightY = entriesY / (meanSigmaY * meanSigmaY);
703  }
704 
705  const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinX(meanSigmaX * meanSigmaX,
706  residualWidthX_1 * residualWidthX_1);
707  const WeightAndResultsPerBin weightAndResultsPerBinX(weightX, error2AndResidualWidth2PerBinX);
708  if (!(entriesX < minHitsPerInterval)) {
709  //Fill absolute weights
710  i_sector.second.WeightX->SetBinContent(i_errBins.first, weightX);
711  v_weightAndResultsPerBinX.push_back(weightAndResultsPerBinX);
712  }
713 
714  const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinY(meanSigmaY * meanSigmaY,
715  residualWidthY_1 * residualWidthY_1);
716  const WeightAndResultsPerBin weightAndResultsPerBinY(weightY, error2AndResidualWidth2PerBinY);
717  if (!(entriesY < minHitsPerInterval)) {
718  //Fill absolute weights
719  i_sector.second.WeightY->SetBinContent(i_errBins.first, weightY);
720  v_weightAndResultsPerBinY.push_back(weightAndResultsPerBinY);
721  }
722  }
723 
724  // Do the final calculations
725 
726  if (v_weightAndResultsPerBinX.empty()) {
727  edm::LogError("CalculateAPE") << "NO error interval of sector " << i_sector.first
728  << " has a valid x APE calculated,\n...so cannot set APE";
729  continue;
730  }
731 
732  if (i_sector.second.isPixel && v_weightAndResultsPerBinY.empty()) {
733  edm::LogError("CalculateAPE") << "NO error interval of sector " << i_sector.first
734  << " has a valid y APE calculated,\n...so cannot set APE";
735  continue;
736  }
737 
738  double correctionX2(999.);
739  double correctionY2(999.);
740 
741  // Get sum of all weights
742  double weightSumX(0.);
743  for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
744  weightSumX += i_apeBin.first;
745  }
746  i_sector.second.WeightX->Scale(1 / weightSumX);
747  double weightSumY(0.);
748  if (i_sector.second.isPixel) {
749  for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
750  weightSumY += i_apeBin.first;
751  }
752  i_sector.second.WeightY->Scale(1 / weightSumY);
753  }
754 
755  if (!setBaseline) {
756  // Calculate weighted mean
757  bool firstIntervalX(true);
758  for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
759  if (firstIntervalX) {
760  correctionX2 = i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthX2);
761  firstIntervalX = false;
762  } else {
763  correctionX2 += i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthX2);
764  }
765  }
766  correctionX2 = correctionX2 / weightSumX;
767  } else {
768  double numeratorX(0.), denominatorX(0.);
769  for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
770  numeratorX += i_apeBin.first * i_apeBin.second.first * i_apeBin.second.second;
771  denominatorX += i_apeBin.first * i_apeBin.second.first;
772  }
773  correctionX2 = numeratorX / denominatorX;
774  }
775 
776  if (i_sector.second.isPixel) {
777  if (!setBaseline) {
778  // Calculate weighted mean
779  bool firstIntervalY(true);
780  for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
781  if (firstIntervalY) {
782  correctionY2 = i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthY2);
783  firstIntervalY = false;
784  } else {
785  correctionY2 += i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthY2);
786  }
787  }
788  correctionY2 = correctionY2 / weightSumY;
789  } else {
790  double numeratorY(0.), denominatorY(0.);
791  for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
792  numeratorY += i_apeBin.first * i_apeBin.second.first * i_apeBin.second.second;
793  denominatorY += i_apeBin.first * i_apeBin.second.first;
794  }
795  correctionY2 = numeratorY / denominatorY;
796  }
797  }
798 
799  if (!setBaseline) {
800  // Calculate updated squared APE of current iteration
801  double apeX2(999.);
802  double apeY2(999.);
803 
804  // old APE value from last iteration
805  if (firstIter) {
806  apeX2 = 0.;
807  apeY2 = 0.;
808  } else {
809  apeX2 = a_apeSectorX[i_sector.first];
810  apeY2 = a_apeSectorY[i_sector.first];
811  }
812  const double apeX2old = apeX2;
813  const double apeY2old = apeY2;
814 
815  // scale APE Correction with value given in cfg (not if smoothed iteration is used)
816  edm::LogInfo("CalculateAPE") << "Unscaled correction x for sector " << i_sector.first << " is "
817  << (correctionX2 > 0. ? +1. : -1.) * std::sqrt(std::fabs(correctionX2));
818  if (!smoothIteration || firstIter)
819  correctionX2 = correctionX2 * correctionScaling * correctionScaling;
820  if (i_sector.second.isPixel) {
821  edm::LogInfo("CalculateAPE") << "Unscaled correction y for sector " << i_sector.first << " is "
822  << (correctionY2 > 0. ? +1. : -1.) * std::sqrt(std::fabs(correctionY2));
823  if (!smoothIteration || firstIter)
824  correctionY2 = correctionY2 * correctionScaling * correctionScaling;
825  }
826 
827  // new APE value
828  // smooth iteration or not?
829  if (apeX2 + correctionX2 < 0.)
830  correctionX2 = -apeX2;
831  if (apeY2 + correctionY2 < 0.)
832  correctionY2 = -apeY2;
833  const double apeX2new(apeX2old + correctionX2);
834  const double apeY2new(apeY2old + correctionY2);
835  if (!smoothIteration || firstIter) {
836  apeX2 = apeX2new;
837  apeY2 = apeY2new;
838  } else {
839  const double apeXtmp = smoothFraction * std::sqrt(apeX2new) + (1 - smoothFraction) * std::sqrt(apeX2old);
840  const double apeYtmp = smoothFraction * std::sqrt(apeY2new) + (1 - smoothFraction) * std::sqrt(apeY2old);
841  apeX2 = apeXtmp * apeXtmp;
842  apeY2 = apeYtmp * apeYtmp;
843  }
844  if (apeX2 < 0. || apeY2 < 0.)
845  edm::LogError("CalculateAPE") << "\n\n\tBad APE, but why???\n\n\n";
846  a_apeSectorX[i_sector.first] = apeX2;
847  a_apeSectorY[i_sector.first] = apeY2;
848 
849  // Set the calculated APE spherical for all modules of strip sectors
850  const double apeX(std::sqrt(apeX2));
851  const double apeY(std::sqrt(apeY2));
852  const double apeZ(std::sqrt(0.5 * (apeX2 + apeY2)));
853  for (auto const& i_rawId : i_sector.second.v_rawId) {
854  if (i_sector.second.isPixel) {
855  apeOutputFile << i_rawId << " " << std::fixed << std::setprecision(5) << apeX << " " << apeY << " " << apeZ
856  << "\n";
857  } else {
858  apeOutputFile << i_rawId << " " << std::fixed << std::setprecision(5) << apeX << " " << apeX << " " << apeX
859  << "\n";
860  }
861  }
862  } else { // In setBaseline mode, just fill estimated mean value of residual width
863  a_apeSectorX[i_sector.first] = correctionX2;
864  a_apeSectorY[i_sector.first] = correctionY2;
865  }
866  }
867  if (!setBaseline)
868  apeOutputFile.close();
869 
870  iterationTreeX->Fill();
871  iterationTreeX->Write("iterTreeX",
872  TObject::kOverwrite); // TObject::kOverwrite needed to not produce another iterTreeX;2
873  iterationTreeY->Fill();
874  iterationTreeY->Write("iterTreeY",
875  TObject::kOverwrite); // TObject::kOverwrite needed to not produce another iterTreeY;2
876  if (firstIter) {
877  sectorNameTree->Fill();
878  sectorNameTree->Write("nameTree");
879  }
880  iterationFile->Close();
881  delete iterationFile;
882 
883  if (baselineFile) {
884  baselineFile->Close();
885  delete baselineFile;
886  }
887 }
888 
889 // ------------ method called to for each event ------------
890 void ApeEstimatorSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
891  // Load APEs from the GT and write them to root files similar to the ones from calculateAPE()
892 
893  if (firstEvent) {
894  // Set baseline or calculate APE value?
895  const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
896 
897  edm::ESHandle<AlignmentErrorsExtended> alignmentErrors;
898  iSetup.get<TrackerAlignmentErrorExtendedRcd>().get(alignmentErrors);
899 
900  // Read in baseline file for calculation of APE value (if not setting baseline)
901  // Has same format as iterationFile
902  const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
903  TFile* baselineFile(nullptr);
904  TTree* sectorNameBaselineTree(nullptr);
905  if (!setBaseline) {
906  std::ifstream baselineFileStream;
907  // Check if baseline file exists
908  baselineFileStream.open(baselineFileName.c_str());
909  if (baselineFileStream.is_open()) {
910  baselineFileStream.close();
911  baselineFile = new TFile(baselineFileName.c_str(), "READ");
912  }
913  if (baselineFile) {
914  edm::LogInfo("CalculateAPE") << "Baseline file for APE values sucessfully opened";
915  baselineFile->GetObject("nameTree;1", sectorNameBaselineTree);
916  } else {
917  edm::LogWarning("CalculateAPE") << "There is NO baseline file for APE values, so normalized residual width =1 "
918  "for ideal conditions is assumed";
919  }
920  }
921 
922  // Set up root file for default APE values
923  const std::string defaultFileName(parameterSet_.getParameter<std::string>("DefaultFile"));
924  TFile* defaultFile = new TFile(defaultFileName.c_str(), "RECREATE");
925 
926  // Naming in the root files has to be iterTreeX to be consistent for the plotting tool
927  TTree* defaultTreeX(nullptr);
928  TTree* defaultTreeY(nullptr);
929  defaultFile->GetObject("iterTreeX;1", defaultTreeX);
930  defaultFile->GetObject("iterTreeY;1", defaultTreeY);
931  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as defaultTree)
932  TTree* sectorNameTree(nullptr);
933  defaultFile->GetObject("nameTree;1", sectorNameTree);
934 
935  bool firstIter(false);
936  if (!defaultTreeX) { // should be always true in setBaseline mode, since file is recreated
937  firstIter = true;
938  defaultTreeX = new TTree("iterTreeX", "Tree for default APE x values from GT");
939  defaultTreeY = new TTree("iterTreeY", "Tree for default APE y values from GT");
940  edm::LogInfo("CalculateAPE") << "First APE iteration (number 0.), create default file with TTree";
941  sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
942  } else {
943  edm::LogWarning("CalculateAPE") << "NOT the first APE iteration (number 0.), is this wanted or forgot to delete "
944  "old iteration file with TTree?";
945  }
946 
947  // Assign the information stored in the trees to arrays
948  double a_defaultSectorX[noSectors_];
949  double a_defaultSectorY[noSectors_];
950 
951  std::string* a_sectorName[noSectors_];
952  std::string* a_sectorBaselineName[noSectors_];
953  for (auto const& i_sector : m_tkSector_) {
954  const unsigned int iSector(i_sector.first);
955  const bool pixelSector(i_sector.second.isPixel);
956 
957  a_defaultSectorX[iSector] = -99.;
958  a_defaultSectorY[iSector] = -99.;
959 
960  a_sectorName[iSector] = nullptr;
961  a_sectorBaselineName[iSector] = nullptr;
962  std::stringstream ss_sector, ss_sectorSuffixed;
963  ss_sector << "Ape_Sector_" << iSector;
964  if (!setBaseline && sectorNameBaselineTree) {
965  sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
966  sectorNameBaselineTree->GetEntry(0);
967  }
968 
969  if (firstIter) { // should be always true in setBaseline mode, since file is recreated
970  ss_sectorSuffixed << ss_sector.str() << "/D";
971  defaultTreeX->Branch(ss_sector.str().c_str(), &a_defaultSectorX[iSector], ss_sectorSuffixed.str().c_str());
972  if (pixelSector) {
973  defaultTreeY->Branch(ss_sector.str().c_str(), &a_defaultSectorY[iSector], ss_sectorSuffixed.str().c_str());
974  }
975  sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
976  } else {
977  defaultTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorX[iSector]);
978  defaultTreeX->GetEntry(0);
979  if (pixelSector) {
980  defaultTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorY[iSector]);
981  defaultTreeY->GetEntry(0);
982  }
983  sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
984  sectorNameTree->GetEntry(0);
985  }
986  }
987 
988  // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
989  for (auto& i_sector : m_tkSector_) {
990  const std::string& name(i_sector.second.name);
991  if (!firstIter) {
992  const std::string& nameLastIter(*a_sectorName[i_sector.first]);
993  if (name != nameLastIter) {
994  edm::LogError("CalculateAPE") << "Inconsistent sector definition in iterationFile for sector "
995  << i_sector.first << ",\n"
996  << "Recent iteration has name \"" << name << "\", while previous had \""
997  << nameLastIter << "\"\n"
998  << "...APE calculation stopped. Please check sector definitions in config!\n";
999  return;
1000  }
1001  } else
1002  a_sectorName[i_sector.first] = new std::string(name);
1003  if (!setBaseline && baselineFile) {
1004  const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
1005  if (name != nameBaseline) {
1006  edm::LogError("CalculateAPE") << "Inconsistent sector definition in baselineFile for sector "
1007  << i_sector.first << ",\n"
1008  << "Recent iteration has name \"" << name << "\", while baseline had \""
1009  << nameBaseline << "\"\n"
1010  << "...APE calculation stopped. Please check sector definitions in config!\n";
1011  return;
1012  }
1013  }
1014  }
1015 
1016  // Loop over sectors for calculating getting default APE
1017  for (auto& i_sector : m_tkSector_) {
1018  double defaultApeX(0.);
1019  double defaultApeY(0.);
1020  unsigned int nModules(0);
1021  for (auto const& i_rawId : i_sector.second.v_rawId) {
1022  std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
1023  for (std::vector<AlignTransformErrorExtended>::const_iterator i_alignError = alignErrors.begin();
1024  i_alignError != alignErrors.end();
1025  ++i_alignError) {
1026  if (i_rawId == i_alignError->rawId()) {
1027  CLHEP::HepSymMatrix errMatrix = i_alignError->matrix();
1028  defaultApeX += errMatrix[0][0];
1029  defaultApeY += errMatrix[1][1];
1030  nModules++;
1031  }
1032  }
1033  }
1034  a_defaultSectorX[i_sector.first] = defaultApeX / nModules;
1035  a_defaultSectorY[i_sector.first] = defaultApeY / nModules;
1036  }
1037  sectorNameTree->Fill();
1038  sectorNameTree->Write("nameTree");
1039  defaultTreeX->Fill();
1040  defaultTreeX->Write("iterTreeX");
1041  defaultTreeY->Fill();
1042  defaultTreeY->Write("iterTreeY");
1043 
1044  defaultFile->Close();
1045  delete defaultFile;
1046 
1047  if (baselineFile) {
1048  baselineFile->Close();
1049  delete baselineFile;
1050  }
1051 
1052  firstEvent = false;
1053  }
1054 }
1055 
1056 // ------------ method called once each job just before starting event loop ------------
1057 void ApeEstimatorSummary::beginJob() {
1058  this->openInputFile();
1059 
1060  this->getTrackerSectorStructs();
1061 
1062  this->bookHists();
1063 }
1064 
1065 // ------------ method called once each job just after ending the event loop ------------
1066 void ApeEstimatorSummary::endJob() {
1067  this->calculateApe();
1068 
1069  this->writeHists();
1070 }
1071 
1072 //define this as a plug-in
1073 DEFINE_FWK_MODULE(ApeEstimatorSummary);
AlignmentErrorsExtended.h
ApeEstimatorSummary::beginJob
void beginJob() override
Definition: ApeEstimatorSummary.cc:1057
ApeEstimatorSummary::residualErrorBinning
std::vector< double > residualErrorBinning()
Definition: ApeEstimatorSummary.cc:277
alignBH_cfg.fixed
fixed
Definition: alignBH_cfg.py:54
EDAnalyzer.h
ApeEstimatorSummary::inputFile_
TFile * inputFile_
Definition: ApeEstimatorSummary.cc:84
ApeEstimatorSummary::parameterSet_
const edm::ParameterSet parameterSet_
Definition: ApeEstimatorSummary.cc:80
MessageLogger.h
ApeEstimatorSummary::firstEvent
bool firstEvent
Definition: ApeEstimatorSummary.cc:82
ESHandle.h
ApeEstimatorSummary::wEntriesOverSigmaX2
Definition: ApeEstimatorSummary.cc:77
mps_splice.entry
entry
Definition: mps_splice.py:68
ApeEstimatorSummary::writeHists
void writeHists()
Definition: ApeEstimatorSummary.cc:286
ApeEstimatorSummary::calculateApe
void calculateApe()
Definition: ApeEstimatorSummary.cc:322
ALCARECODTCalibSynchDQM_cff.baseDir
baseDir
Definition: ALCARECODTCalibSynchDQM_cff.py:14
photonAnalyzer_cfi.yMin
yMin
Definition: photonAnalyzer_cfi.py:86
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
edm::one::EDAnalyzer
Definition: EDAnalyzer.h:30
tools.TF1
TF1
Definition: tools.py:23
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
TrackerSectorStruct::EntriesY
TH1 * EntriesY
Definition: TrackerSectorStruct.h:98
TrackerSectorStruct::EntriesX
TH1 * EntriesX
Definition: TrackerSectorStruct.h:96
TrackerAlignmentErrorExtendedRcd.h
ApeEstimatorSummary::wInvalid
Definition: ApeEstimatorSummary.cc:77
photonAnalyzer_cfi.xMin
xMin
Definition: photonAnalyzer_cfi.py:82
EDMException.h
MakerMacros.h
edm::EventSetup::get
T get() const
Definition: EventSetup.h:80
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
ApeEstimatorSummary
Definition: ApeEstimatorSummary.cc:59
Abs
T Abs(T a)
Definition: MathUtil.h:49
InefficientDoubleROC.inputFileName
inputFileName
Definition: InefficientDoubleROC.py:437
Fit
Definition: Fit.h:32
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HIPAlignmentAlgorithm_cfi.iterationFile
iterationFile
Definition: HIPAlignmentAlgorithm_cfi.py:84
edm::ESHandle
Definition: DTSurvey.h:22
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:86
CommonMethods.isnan
def isnan(num)
Definition: CommonMethods.py:98
TrackerSectorStruct::m_binnedHists
std::map< unsigned int, std::map< std::string, TH1 * > > m_binnedHists
Definition: TrackerSectorStruct.h:92
multiplicitycorr_cfi.yMax
yMax
Definition: multiplicitycorr_cfi.py:6
ApeEstimatorSummary::endJob
void endJob() override
Definition: ApeEstimatorSummary.cc:1066
ApeEstimatorSummary::wEntries
Definition: ApeEstimatorSummary.cc:77
ApeEstimatorSummary::~ApeEstimatorSummary
~ApeEstimatorSummary() override
Definition: ApeEstimatorSummary.cc:104
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
newFWLiteAna.fullName
fullName
Definition: newFWLiteAna.py:122
ApeEstimatorSummary::openInputFile
void openInputFile()
Definition: ApeEstimatorSummary.cc:110
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:223
edm::ParameterSet
Definition: ParameterSet.h:47
Event.h
ApeEstimatorSummary::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: ApeEstimatorSummary.cc:890
iEvent
int iEvent
Definition: GenABIO.cc:224
ApeEstimatorSummary::getTrackerSectorStructs
void getTrackerSectorStructs()
Definition: ApeEstimatorSummary.cc:128
ApeEstimatorSummary::noSectors_
unsigned int noSectors_
Definition: ApeEstimatorSummary.cc:87
edm::EventSetup
Definition: EventSetup.h:57
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
get
#define get
Frameworkfwd.h
TrackerSectorStruct::v_rawId
std::vector< unsigned int > v_rawId
Definition: TrackerSectorStruct.h:78
Exception
Definition: hltDiff.cc:246
ApeEstimatorSummary::wUnity
Definition: ApeEstimatorSummary.cc:77
TrackerSectorStruct.h
TrackerOfflineValidation_Dqm_cff.dirName
dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:55
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
multiplicitycorr_cfi.xMax
xMax
Definition: multiplicitycorr_cfi.py:5
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
AlignmentErrorsExtended::m_alignError
std::vector< AlignTransformErrorExtended > m_alignError
Definition: AlignmentErrorsExtended.h:19
ParameterSet.h
ApeEstimatorSummary::bookHists
void bookHists()
Definition: ApeEstimatorSummary.cc:204
TrackerSectorStruct::Name
TH1 * Name
Definition: TrackerSectorStruct.h:75
edm::Event
Definition: Event.h:73
TrackerSectorStruct::name
std::string name
Definition: TrackerSectorStruct.h:74
align_cfg.iteration
iteration
Definition: align_cfg.py:5
ApeEstimatorSummary::ApeEstimatorSummary
ApeEstimatorSummary(const edm::ParameterSet &)
Definition: ApeEstimatorSummary.cc:101
TrackerAlignmentErrorExtendedRcd
Definition: TrackerAlignmentErrorExtendedRcd.h:6
edm::errors::Configuration
Definition: EDMException.h:36
ApeEstimatorSummary::ApeWeight
ApeWeight
Definition: ApeEstimatorSummary.cc:77
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
TrackerSectorStruct::isPixel
bool isPixel
Definition: TrackerSectorStruct.h:102
TrackerSectorStruct
Definition: TrackerSectorStruct.h:25
Generated for CMSSW Reference Manual by doxygen 1.8.16