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