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