test
CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros 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 
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  Int_t fitResultX(0);
524  if(integralX>minHitsPerInterval){
525  if(mHists["norResX"]->Fit(&funcX_1, fitOpt)){
526  edm::LogWarning("CalculateAPE")<<"Fit1 did not work: "<<mHists["norResX"]->Fit(&funcX_1, fitOpt);
527  continue;
528  }
529  fitResultX = mHists["norResX"]->Fit(&funcX_1, fitOpt);
530  std::cout<<"FitResult1\t"<<fitResultX<<"\n";
531  }
532  Double_t meanX_1 = funcX_1.GetParameter(1);
533  Double_t sigmaX_1 = funcX_1.GetParameter(2);
534 
535 
536  // Second gaus fit
537  TF1 funcX_2("mygausX2","gaus",meanX_1 - sigmaFactorFit*TMath::Abs(sigmaX_1), meanX_1 + sigmaFactorFit*TMath::Abs(sigmaX_1));
538  funcX_2.SetParameters(funcX_1.GetParameter(0),meanX_1,sigmaX_1);
539  if(integralX>minHitsPerInterval){
540  if(mHists["norResX"]->Fit(&funcX_2, fitOpt)){
541  edm::LogWarning("CalculateAPE")<<"Fit2 did not work for x : "<<mHists["norResX"]->Fit(&funcX_2, fitOpt);
542  continue;
543  }
544  fitResultX = mHists["norResX"]->Fit(&funcX_2, fitOpt);
545  }
546  Double_t meanX_2 = funcX_2.GetParameter(1);
547  Double_t sigmaX_2 = funcX_2.GetParameter(2);
548 
549 
550  // Now the same for y coordinate
551  double entriesY(0.);
552  double meanSigmaY(0.);
553  if((*i_sector).second.isPixel){
554  entriesY = mHists["sigmaY"]->GetEntries();
555  meanSigmaY = mHists["sigmaY"]->GetMean();
556  }
557 
558  Double_t meanY_1(0.);
559  Double_t sigmaY_1(0.);
560  Double_t meanY_2(0.);
561  Double_t sigmaY_2(0.);
562  double meanY(0.);
563  double rmsY(0.);
564  if((*i_sector).second.isPixel){
565  // Fitting Parameters
566  double yMin = mHists["norResY"]->GetXaxis()->GetXmin();
567  double yMax = mHists["norResY"]->GetXaxis()->GetXmax();
568  double integralY = mHists["norResY"]->Integral();
569  meanY = mHists["norResY"]->GetMean();
570  rmsY = mHists["norResY"]->GetRMS();
571  double maximumY = mHists["norResY"]->GetBinContent(mHists["norResY"]->GetMaximumBin());
572 
573  // First Gaus Fit
574  TF1 funcY_1("mygausY", "gaus", yMin, yMax);
575  funcY_1.SetParameters(maximumY, meanY, rmsY);
576  Int_t fitResultY(0);
577  if(integralY>minHitsPerInterval){
578  if(mHists["norResY"]->Fit(&funcY_1, fitOpt)){
579  edm::LogWarning("CalculateAPE")<<"Fit1 did not work: "<<mHists["norResY"]->Fit(&funcY_1, fitOpt);
580  continue;
581  }
582  fitResultY = mHists["norResY"]->Fit(&funcY_1, fitOpt);
583  std::cout<<"FitResult2\t"<<fitResultY<<"\n";
584  }
585  meanY_1 = funcY_1.GetParameter(1);
586  sigmaY_1 = funcY_1.GetParameter(2);
587 
588  // Second gaus fit
589  TF1 funcY_2("mygausY2","gaus",meanY_1 - sigmaFactorFit*TMath::Abs(sigmaY_1), meanY_1 + sigmaFactorFit*TMath::Abs(sigmaY_1));
590  funcY_2.SetParameters(funcY_1.GetParameter(0),meanY_1,sigmaY_1);
591  if(integralY>minHitsPerInterval){
592  if(mHists["norResY"]->Fit(&funcY_2, fitOpt)){
593  edm::LogWarning("CalculateAPE")<<"Fit2 did not work for y : "<<mHists["norResY"]->Fit(&funcY_2, fitOpt);
594  continue;
595  }
596  fitResultY = mHists["norResY"]->Fit(&funcY_2, fitOpt);
597  }
598  meanY_2 = funcY_2.GetParameter(1);
599  sigmaY_2 = funcY_2.GetParameter(2);
600  }
601 
602 
603 
604  // Fill histograms
605  double fitMeanX_1(meanX_1), fitMeanX_2(meanX_2);
606  double residualWidthX_1(sigmaX_1), residualWidthX_2(sigmaX_2);
607  double fitMeanY_1(meanY_1), fitMeanY_2(meanY_2);
608  double residualWidthY_1(sigmaY_1), residualWidthY_2(sigmaY_2);
609 
610  double correctionX2_1(-0.0010), correctionX2_2(-0.0010);
611  double correctionY2_1(-0.0010), correctionY2_2(-0.0010);
612  correctionX2_1 = meanSigmaX*meanSigmaX*(residualWidthX_1*residualWidthX_1 -baselineWidthX2);
613  correctionX2_2 = meanSigmaX*meanSigmaX*(residualWidthX_2*residualWidthX_2 -baselineWidthX2);
614  if((*i_sector).second.isPixel){
615  correctionY2_1 = meanSigmaY*meanSigmaY*(residualWidthY_1*residualWidthY_1 -baselineWidthY2);
616  correctionY2_2 = meanSigmaY*meanSigmaY*(residualWidthY_2*residualWidthY_2 -baselineWidthY2);
617  }
618 
619  double correctionX_1 = correctionX2_1>=0. ? std::sqrt(correctionX2_1) : -std::sqrt(-correctionX2_1);
620  double correctionX_2 = correctionX2_2>=0. ? std::sqrt(correctionX2_2) : -std::sqrt(-correctionX2_2);
621  double correctionY_1 = correctionY2_1>=0. ? std::sqrt(correctionY2_1) : -std::sqrt(-correctionY2_1);
622  double correctionY_2 = correctionY2_2>=0. ? std::sqrt(correctionY2_2) : -std::sqrt(-correctionY2_2);
623  // Meanwhile, this got very bad default values, or? (negative corrections allowed)
624  if(isnan(correctionX_1))correctionX_1 = -0.0010;
625  if(isnan(correctionX_2))correctionX_2 = -0.0010;
626  if(isnan(correctionY_1))correctionY_1 = -0.0010;
627  if(isnan(correctionY_2))correctionY_2 = -0.0010;
628 
629  if(entriesX<minHitsPerInterval){
630  meanX = 0.; rmsX = -0.0010;
631  fitMeanX_1 = 0.; correctionX_1 = residualWidthX_1 = -0.0010;
632  fitMeanX_2 = 0.; correctionX_2 = residualWidthX_2 = -0.0010;
633  }
634 
635  if(entriesY<minHitsPerInterval){
636  meanY = 0.; rmsY = -0.0010;
637  fitMeanY_1 = 0.; correctionY_1 = residualWidthY_1 = -0.0010;
638  fitMeanY_2 = 0.; correctionY_2 = residualWidthY_2 = -0.0010;
639  }
640 
641  (*i_sector).second.MeanX ->SetBinContent((*i_errBins).first,meanX);
642  (*i_sector).second.RmsX ->SetBinContent((*i_errBins).first,rmsX);
643 
644  (*i_sector).second.FitMeanX1 ->SetBinContent((*i_errBins).first,fitMeanX_1);
645  (*i_sector).second.ResidualWidthX1->SetBinContent((*i_errBins).first,residualWidthX_1);
646  (*i_sector).second.CorrectionX1 ->SetBinContent((*i_errBins).first,correctionX_1*10000.);
647 
648  (*i_sector).second.FitMeanX2 ->SetBinContent((*i_errBins).first,fitMeanX_2);
649  (*i_sector).second.ResidualWidthX2->SetBinContent((*i_errBins).first,residualWidthX_2);
650  (*i_sector).second.CorrectionX2 ->SetBinContent((*i_errBins).first,correctionX_2*10000.);
651 
652  if((*i_sector).second.isPixel){
653  (*i_sector).second.MeanY ->SetBinContent((*i_errBins).first,meanY);
654  (*i_sector).second.RmsY ->SetBinContent((*i_errBins).first,rmsY);
655 
656  (*i_sector).second.FitMeanY1 ->SetBinContent((*i_errBins).first,fitMeanY_1);
657  (*i_sector).second.ResidualWidthY1->SetBinContent((*i_errBins).first,residualWidthY_1);
658  (*i_sector).second.CorrectionY1 ->SetBinContent((*i_errBins).first,correctionY_1*10000.);
659 
660  (*i_sector).second.FitMeanY2 ->SetBinContent((*i_errBins).first,fitMeanY_2);
661  (*i_sector).second.ResidualWidthY2->SetBinContent((*i_errBins).first,residualWidthY_2);
662  (*i_sector).second.CorrectionY2 ->SetBinContent((*i_errBins).first,correctionY_2*10000.);
663  }
664 
665 
666  // Use result for bin only when entries>=minHitsPerInterval
667  if(entriesX<minHitsPerInterval && entriesY<minHitsPerInterval)continue;
668 
669  double weightX(0.);
670  double weightY(0.);
671  if(apeWeight==wUnity){
672  weightX = 1.;
673  weightY = 1.;
674  }
675  else if(apeWeight==wEntries){
676  weightX = entriesX;
677  weightY = entriesY;
678  }
679  else if(apeWeight==wEntriesOverSigmaX2){
680  weightX = entriesX/(meanSigmaX*meanSigmaX);
681  weightY = entriesY/(meanSigmaY*meanSigmaY);
682  }
683 
684  const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinX(meanSigmaX*meanSigmaX, residualWidthX_1*residualWidthX_1);
685  const WeightAndResultsPerBin weightAndResultsPerBinX(weightX, error2AndResidualWidth2PerBinX);
686  if(!(entriesX<minHitsPerInterval)){
687  //Fill absolute weights
688  (*i_sector).second.WeightX->SetBinContent((*i_errBins).first,weightX);
689  v_weightAndResultsPerBinX.push_back(weightAndResultsPerBinX);
690  }
691 
692  const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinY(meanSigmaY*meanSigmaY, residualWidthY_1*residualWidthY_1);
693  const WeightAndResultsPerBin weightAndResultsPerBinY(weightY, error2AndResidualWidth2PerBinY);
694  if(!(entriesY<minHitsPerInterval)){
695  //Fill absolute weights
696  (*i_sector).second.WeightY->SetBinContent((*i_errBins).first,weightY);
697  v_weightAndResultsPerBinY.push_back(weightAndResultsPerBinY);
698  }
699  }
700 
701 
702  // Do the final calculations
703 
704  if(v_weightAndResultsPerBinX.size()==0){
705  edm::LogError("CalculateAPE")<<"NO error interval of sector "<<(*i_sector).first<<" has a valid x APE calculated,\n...so cannot set APE";
706  continue;
707  }
708 
709  if((*i_sector).second.isPixel && v_weightAndResultsPerBinY.size()==0){
710  edm::LogError("CalculateAPE")<<"NO error interval of sector "<<(*i_sector).first<<" has a valid y APE calculated,\n...so cannot set APE";
711  continue;
712  }
713 
714  double correctionX2(999.);
715  double correctionY2(999.);
716 
717  // Get sum of all weights
718  double weightSumX(0.);
719  std::vector<WeightAndResultsPerBin>::const_iterator i_apeBin;
720  for(i_apeBin=v_weightAndResultsPerBinX.begin(); i_apeBin!=v_weightAndResultsPerBinX.end(); ++i_apeBin){
721  weightSumX += i_apeBin->first;
722  }
723  (*i_sector).second.WeightX->Scale(1/weightSumX);
724  double weightSumY(0.);
725  if((*i_sector).second.isPixel){
726  std::vector<WeightAndResultsPerBin>::const_iterator i_apeBin;
727  for(i_apeBin=v_weightAndResultsPerBinY.begin(); i_apeBin!=v_weightAndResultsPerBinY.end(); ++i_apeBin){
728  weightSumY += i_apeBin->first;
729  }
730  (*i_sector).second.WeightY->Scale(1/weightSumY);
731  }
732 
733 
734  if(!setBaseline){
735  // Calculate weighted mean
736  bool firstIntervalX(true);
737  for(i_apeBin=v_weightAndResultsPerBinX.begin(); i_apeBin!=v_weightAndResultsPerBinX.end(); ++i_apeBin){
738  if(firstIntervalX){
739  correctionX2 = i_apeBin->first*i_apeBin->second.first*(i_apeBin->second.second - baselineWidthX2);
740  firstIntervalX = false;
741  }
742  else{
743  correctionX2 += i_apeBin->first*i_apeBin->second.first*(i_apeBin->second.second - baselineWidthX2);
744  }
745  }
746  correctionX2 = correctionX2/weightSumX;
747  }
748  else{
749  double numeratorX(0.), denominatorX(0.);
750  std::vector<WeightAndResultsPerBin>::const_iterator i_apeBin;
751  for(i_apeBin=v_weightAndResultsPerBinX.begin(); i_apeBin!=v_weightAndResultsPerBinX.end(); ++i_apeBin){
752  numeratorX += i_apeBin->first*i_apeBin->second.first*i_apeBin->second.second;
753  denominatorX += i_apeBin->first*i_apeBin->second.first;
754  }
755  correctionX2 = numeratorX/denominatorX;
756  }
757 
758  if((*i_sector).second.isPixel){
759  if(!setBaseline){
760  // Calculate weighted mean
761  bool firstIntervalY(true);
762  for(i_apeBin=v_weightAndResultsPerBinY.begin(); i_apeBin!=v_weightAndResultsPerBinY.end(); ++i_apeBin){
763  if(firstIntervalY){
764  correctionY2 = i_apeBin->first*i_apeBin->second.first*(i_apeBin->second.second - baselineWidthY2);
765  firstIntervalY = false;
766  }
767  else{
768  correctionY2 += i_apeBin->first*i_apeBin->second.first*(i_apeBin->second.second - baselineWidthY2);
769  }
770  }
771  correctionY2 = correctionY2/weightSumY;
772  }
773  else{
774  double numeratorY(0.), denominatorY(0.);
775  std::vector<WeightAndResultsPerBin>::const_iterator i_apeBin;
776  for(i_apeBin=v_weightAndResultsPerBinY.begin(); i_apeBin!=v_weightAndResultsPerBinY.end(); ++i_apeBin){
777  numeratorY += i_apeBin->first*i_apeBin->second.first*i_apeBin->second.second;
778  denominatorY += i_apeBin->first*i_apeBin->second.first;
779  }
780  correctionY2 = numeratorY/denominatorY;
781  }
782  }
783 
784 
785 
786  if(!setBaseline){
787  // Calculate updated squared APE of current iteration
788  double apeX2(999.);
789  double apeY2(999.);
790 
791  // old APE value from last iteration
792  if(firstIter){
793  apeX2 = 0.;
794  apeY2 = 0.;
795  }
796  else{
797  apeX2 = a_apeSectorX[(*i_sector).first];
798  apeY2 = a_apeSectorY[(*i_sector).first];
799  }
800  const double apeX2old = apeX2;
801  const double apeY2old = apeY2;
802 
803  // scale APE Correction with value given in cfg (not if smoothed iteration is used)
804  edm::LogInfo("CalculateAPE")<<"Unscaled correction x for sector "<<(*i_sector).first<<" is "<<(correctionX2>0. ? +1. : -1.)*std::sqrt(std::fabs(correctionX2));
805  if(!smoothIteration || firstIter)correctionX2 = correctionX2*correctionScaling*correctionScaling;
806  if((*i_sector).second.isPixel){
807  edm::LogInfo("CalculateAPE")<<"Unscaled correction y for sector "<<(*i_sector).first<<" is "<<(correctionY2>0. ? +1. : -1.)*std::sqrt(std::fabs(correctionY2));
808  if(!smoothIteration || firstIter)correctionY2 = correctionY2*correctionScaling*correctionScaling;
809  }
810 
811  // new APE value
812  // smooth iteration or not?
813  if(apeX2 + correctionX2 < 0.) correctionX2 = -apeX2;
814  if(apeY2 + correctionY2 < 0.) correctionY2 = -apeY2;
815  const double apeX2new(apeX2old + correctionX2);
816  const double apeY2new(apeY2old + correctionY2);
817  if(!smoothIteration || firstIter){
818  apeX2 = apeX2new;
819  apeY2 = apeY2new;
820  }
821  else{
822  apeX2 = std::pow(smoothFraction*std::sqrt(apeX2new) + (1-smoothFraction)*std::sqrt(apeX2old), 2);
823  apeY2 = std::pow(smoothFraction*std::sqrt(apeY2new) + (1-smoothFraction)*std::sqrt(apeY2old), 2);
824  }
825  if(apeX2<0. || apeY2<0.)edm::LogError("CalculateAPE")<<"\n\n\tBad APE, but why???\n\n\n";
826  a_apeSectorX[(*i_sector).first] = apeX2;
827  a_apeSectorY[(*i_sector).first] = apeY2;
828 
829  // Set the calculated APE spherical for all modules of strip sectors
830  const double apeX(std::sqrt(apeX2));
831  const double apeY(std::sqrt(apeY2));
832  const double apeZ(std::sqrt(0.5*(apeX2+apeY2)));
833  std::vector<unsigned int>::const_iterator i_rawId;
834  for(i_rawId = (*i_sector).second.v_rawId.begin(); i_rawId != (*i_sector).second.v_rawId.end(); ++i_rawId){
835  if((*i_sector).second.isPixel){
836  apeOutputFile<<*i_rawId<<" "<<std::fixed<<std::setprecision(5)<<apeX<<" "<<apeY<<" "<<apeZ<<"\n";
837  }
838  else{
839  apeOutputFile<<*i_rawId<<" "<<std::fixed<<std::setprecision(5)<<apeX<<" "<<apeX<<" "<<apeX<<"\n";
840  }
841  }
842  }
843  // In setBaseline mode, just fill estimated mean value of residual width
844  else{
845  a_apeSectorX[(*i_sector).first] = correctionX2;
846  a_apeSectorY[(*i_sector).first] = correctionY2;
847  }
848  }
849  if(!setBaseline)apeOutputFile.close();
850 
851  iterationTreeX->Fill();
852  iterationTreeX->Write("iterTreeX", TObject::kOverwrite); // TObject::kOverwrite needed to not produce another iterTreeX;2
853  iterationTreeY->Fill();
854  iterationTreeY->Write("iterTreeY", TObject::kOverwrite); // TObject::kOverwrite needed to not produce another iterTreeY;2
855  if(firstIter){
856  sectorNameTree->Fill();
857  sectorNameTree->Write("nameTree");
858  }
859  iterationFile->Close();
860 
861  if(baselineFile)baselineFile->Close();
862 }
863 
864 
865 
866 // ------------ method called to for each event ------------
867 void
868 ApeEstimatorSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
869 {
870  // Load APEs from the GT and write them to root files similar to the ones from calculateAPE()
871 
872  if(firstEvent){
873  // Set baseline or calculate APE value?
874  const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
875 
876  edm::ESHandle<AlignmentErrorsExtended> alignmentErrors;
877  iSetup.get<TrackerAlignmentErrorExtendedRcd>().get(alignmentErrors);
878 
879  // Read in baseline file for calculation of APE value (if not setting baseline)
880  // Has same format as iterationFile
881  const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
882  TFile* baselineFile(0);
883  TTree* sectorNameBaselineTree(0);
884  if(!setBaseline){
885  std::ifstream baselineFileStream;
886  // Check if baseline file exists
887  baselineFileStream.open(baselineFileName.c_str());
888  if(baselineFileStream.is_open()){
889  baselineFileStream.close();
890  baselineFile = new TFile(baselineFileName.c_str(),"READ");
891  }
892  if(baselineFile){
893  edm::LogInfo("CalculateAPE")<<"Baseline file for APE values sucessfully opened";
894  baselineFile->GetObject("nameTree;1",sectorNameBaselineTree);
895  }
896  else{
897  edm::LogWarning("CalculateAPE")<<"There is NO baseline file for APE values, so normalized residual width =1 for ideal conditions is assumed";
898  }
899  }
900 
901  // Set up root file for default APE values
902  const std::string defaultFileName(parameterSet_.getParameter<std::string>("DefaultFile"));
903  TFile* defaultFile = new TFile(defaultFileName.c_str(),"RECREATE");
904 
905  // Naming in the root files has to be iterTreeX to be consistent for the plotting tool
906  TTree* defaultTreeX(0);
907  TTree* defaultTreeY(0);
908  defaultFile->GetObject("iterTreeX;1",defaultTreeX);
909  defaultFile->GetObject("iterTreeY;1",defaultTreeY);
910  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as defaultTree)
911  TTree* sectorNameTree(0);
912  defaultFile->GetObject("nameTree;1",sectorNameTree);
913 
914  bool firstIter(false);
915  if(!defaultTreeX){ // should be always true in setBaseline mode, since file is recreated
916  firstIter = true;
917  defaultTreeX = new TTree("iterTreeX","Tree for default APE x values from GT");
918  defaultTreeY = new TTree("iterTreeY","Tree for default APE y values from GT");
919  edm::LogInfo("CalculateAPE")<<"First APE iteration (number 0.), create default file with TTree";
920  sectorNameTree = new TTree("nameTree","Tree with names of sectors");
921  }
922  else{
923  edm::LogWarning("CalculateAPE")<<"NOT the first APE iteration (number 0.), is this wanted or forgot to delete old iteration file with TTree?";
924  }
925 
926  // Assign the information stored in the trees to arrays
927  double a_defaultSectorX[16589];
928  double a_defaultSectorY[16589];
929 
930  std::string* a_sectorName[16589];
931  std::string* a_sectorBaselineName[16589];
932  std::map<unsigned int, TrackerSectorStruct>::const_iterator i_sector;
933  for(i_sector = m_tkSector_.begin(); i_sector != m_tkSector_.end(); ++i_sector){
934  const unsigned int iSector(i_sector->first);
935  const bool pixelSector(i_sector->second.isPixel);
936 
937  a_defaultSectorX[iSector] = -99.;
938  a_defaultSectorY[iSector] = -99.;
939 
940  a_sectorName[iSector] = 0;
941  a_sectorBaselineName[iSector] = 0;
942  std::stringstream ss_sector, ss_sectorSuffixed;
943  ss_sector << "Ape_Sector_" << iSector;
944  if(!setBaseline && sectorNameBaselineTree){
945  sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
946  sectorNameBaselineTree->GetEntry(0);
947  }
948 
949  if(firstIter){ // should be always true in setBaseline mode, since file is recreated
950  ss_sectorSuffixed << ss_sector.str() << "/D";
951  defaultTreeX->Branch(ss_sector.str().c_str(), &a_defaultSectorX[iSector], ss_sectorSuffixed.str().c_str());
952  if(pixelSector){
953  defaultTreeY->Branch(ss_sector.str().c_str(), &a_defaultSectorY[iSector], ss_sectorSuffixed.str().c_str());
954  }
955  sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
956  }
957  else{
958  defaultTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorX[iSector]);
959  defaultTreeX->GetEntry(0);
960  if(pixelSector){
961  defaultTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorY[iSector]);
962  defaultTreeY->GetEntry(0);
963  }
964  sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
965  sectorNameTree->GetEntry(0);
966  }
967  }
968 
969 
970  // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
971  for(std::map<unsigned int,TrackerSectorStruct>::iterator i_sector = m_tkSector_.begin(); i_sector != m_tkSector_.end(); ++i_sector){
972  const std::string& name(i_sector->second.name);
973  if(!firstIter){
974  const std::string& nameLastIter(*a_sectorName[(*i_sector).first]);
975  if(name!=nameLastIter){
976  edm::LogError("CalculateAPE")<<"Inconsistent sector definition in iterationFile for sector "<<i_sector->first<<",\n"
977  <<"Recent iteration has name \""<<name<<"\", while previous had \""<<nameLastIter<<"\"\n"
978  <<"...APE calculation stopped. Please check sector definitions in config!\n";
979  return;
980  }
981  }
982  else a_sectorName[(*i_sector).first] = new std::string(name);
983  if(!setBaseline && baselineFile){
984  const std::string& nameBaseline(*a_sectorBaselineName[(*i_sector).first]);
985  if(name!=nameBaseline){
986  edm::LogError("CalculateAPE")<<"Inconsistent sector definition in baselineFile for sector "<<i_sector->first<<",\n"
987  <<"Recent iteration has name \""<<name<<"\", while baseline had \""<<nameBaseline<<"\"\n"
988  <<"...APE calculation stopped. Please check sector definitions in config!\n";
989  return;
990  }
991  }
992  }
993 
994 
995  // Loop over sectors for calculating getting default APE
996 
997  for(std::map<unsigned int,TrackerSectorStruct>::iterator i_sector = m_tkSector_.begin(); i_sector != m_tkSector_.end(); ++i_sector){
998 
999  double defaultApeX(0.);
1000  double defaultApeY(0.);
1001  unsigned int nModules(0);
1002  std::vector<unsigned int>::const_iterator i_rawId;
1003  for(i_rawId = (*i_sector).second.v_rawId.begin(); i_rawId != (*i_sector).second.v_rawId.end(); ++i_rawId){
1004  std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
1005  for(std::vector<AlignTransformErrorExtended>::const_iterator i_alignError = alignErrors.begin(); i_alignError != alignErrors.end(); ++i_alignError){
1006  if(*i_rawId == i_alignError->rawId()){
1007  CLHEP::HepSymMatrix errMatrix = i_alignError->matrix();
1008  defaultApeX += errMatrix[0][0];
1009  defaultApeY += errMatrix[1][1];
1010  nModules++;
1011  }
1012  }
1013  }
1014  a_defaultSectorX[(*i_sector).first] = defaultApeX/nModules;
1015  a_defaultSectorY[(*i_sector).first] = defaultApeY/nModules;
1016 
1017  }
1018  sectorNameTree->Fill();
1019  sectorNameTree->Write("nameTree");
1020  defaultTreeX->Fill();
1021  defaultTreeX->Write("iterTreeX");
1022  defaultTreeY->Fill();
1023  defaultTreeY->Write("iterTreeY");
1024 
1025  defaultFile->Close();
1026  if(baselineFile)baselineFile->Close();
1027 
1028 
1029  firstEvent = false;
1030  }
1031 
1032 }
1033 
1034 
1035 // ------------ method called once each job just before starting event loop ------------
1036 void
1037 ApeEstimatorSummary::beginJob()
1038 {
1039  this->openInputFile();
1040 
1041  this->getTrackerSectorStructs();
1042 
1043  this->bookHists();
1044 }
1045 
1046 // ------------ method called once each job just after ending the event loop ------------
1047 void
1048 ApeEstimatorSummary::endJob() {
1049 
1050  this->calculateApe();
1051 
1052  this->writeHists();
1053 
1054 
1055 }
1056 
1057 //define this as a plug-in
1058 DEFINE_FWK_MODULE(ApeEstimatorSummary);
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
writeAntiElectronDiscrMVA_cfi.inputFileName
tuple inputFileName
Definition: writeAntiElectronDiscrMVA_cfi.py:43
ApeEstimatorSummary::ApeEstimatorSummary
ApeEstimatorSummary(const edm::ParameterSet &)
Definition: ApeEstimatorSummary.cc:105
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
ApeEstimatorSummary::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: ApeEstimatorSummary.cc:868
iEvent
int iEvent
Definition: GenABIO.cc:230
align_cfg.iteration
tuple iteration
Definition: align_cfg.py:5
Fit
Definition: Fit.h:34
edm::detail::isnan
bool isnan(float x)
Definition: math.h:13
TrackerOfflineValidation_Dqm_cff.dirName
tuple dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:56
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
HLT_25ns14e33_v1_cff.yMin
tuple yMin
Definition: HLT_25ns14e33_v1_cff.py:2488
Abs
T Abs(T a)
Definition: MathUtil.h:49
newFWLiteAna.fullName
string fullName
Definition: newFWLiteAna.py:122
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
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
HLT_25ns14e33_v1_cff.yMax
tuple yMax
Definition: HLT_25ns14e33_v1_cff.py:2491
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
dir
dbl *** dir
Definition: mlp_gen.cc:35
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.5