CMS 3D CMS Logo
List of all members | Public Member Functions | Private Types | Private Member Functions | Private Attributes
ApeEstimatorSummary Class Reference

#include <Alignment/APEEstimation/src/ApeEstimatorSummary.cc>

Inheritance diagram for ApeEstimatorSummary:
edm::EDAnalyzer edm::EDConsumerBase

Public Member Functions

 ApeEstimatorSummary (const edm::ParameterSet &)
 
 ~ApeEstimatorSummary () override
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
 ~EDAnalyzer () override
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo > consumesInfo () const
 
void convertCurrentProcessAlias (std::string const &processName)
 Convert "@currentProcess" in InputTag process names to the actual current process name. More...
 
 EDConsumerBase ()
 
 EDConsumerBase (EDConsumerBase const &)=delete
 
 EDConsumerBase (EDConsumerBase &&)=default
 
ProductResolverIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
std::vector< ProductResolverIndexAndSkipBit > const & itemsToGetFrom (BranchType iType) const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
EDConsumerBase const & operator= (EDConsumerBase const &)=delete
 
EDConsumerBase & operator= (EDConsumerBase &&)=default
 
bool registeredToConsume (ProductResolverIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
 
virtual ~EDConsumerBase () noexcept(false)
 

Private Types

enum  ApeWeight { wInvalid, wUnity, wEntries, wEntriesOverSigmaX2 }
 

Private Member Functions

void analyze (const edm::Event &, const edm::EventSetup &) override
 
void beginJob () override
 
void bookHists ()
 
void calculateApe ()
 
void endJob () override
 
void getTrackerSectorStructs ()
 
void openInputFile ()
 
std::vector< double > residualErrorBinning ()
 
void writeHists ()
 

Private Attributes

bool firstEvent =true
 
TFile * inputFile_
 
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
 
const edm::ParameterSet parameterSet_
 

Additional Inherited Members

- Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer ModuleType
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Description: [one line class summary]

Implementation: [Notes on implementation]

Definition at line 62 of file ApeEstimatorSummary.cc.

Member Enumeration Documentation

enum ApeEstimatorSummary::ApeWeight
private
Enumerator
wInvalid 
wUnity 
wEntries 
wEntriesOverSigmaX2 

Definition at line 81 of file ApeEstimatorSummary.cc.

81 {wInvalid, wUnity, wEntries, wEntriesOverSigmaX2};

Constructor & Destructor Documentation

ApeEstimatorSummary::ApeEstimatorSummary ( const edm::ParameterSet &  iConfig)
explicit

Definition at line 105 of file ApeEstimatorSummary.cc.

105  :
106 parameterSet_(iConfig),
107 inputFile_(nullptr)
108 {
109 }
ApeEstimatorSummary::parameterSet_
const edm::ParameterSet parameterSet_
Definition: ApeEstimatorSummary.cc:84
ApeEstimatorSummary::~ApeEstimatorSummary ( )
override

Definition at line 112 of file ApeEstimatorSummary.cc.

113 {
114 }

Member Function Documentation

void ApeEstimatorSummary::analyze ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)
overrideprivate

Definition at line 864 of file ApeEstimatorSummary.cc.

References firstEvent, edm::EventSetup::get(), edm::ParameterSet::getParameter(), AlignmentErrorsExtended::m_alignError, m_tkSector_, dataset::name, parameterSet_, and AlCaHLTBitMon_QueryRunRegistry::string.

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(nullptr);
879  TTree* sectorNameBaselineTree(nullptr);
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(nullptr);
903  TTree* defaultTreeY(nullptr);
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(nullptr);
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] = nullptr;
937  a_sectorBaselineName[iSector] = nullptr;
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 }
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
AlignmentErrorsExtended::m_alignError
std::vector< AlignTransformErrorExtended > m_alignError
Definition: AlignmentErrorsExtended.h:19
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
ApeEstimatorSummary::parameterSet_
const edm::ParameterSet parameterSet_
Definition: ApeEstimatorSummary.cc:84
void ApeEstimatorSummary::beginJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 1033 of file ApeEstimatorSummary.cc.

References bookHists(), getTrackerSectorStructs(), and openInputFile().

1034 {
1035  this->openInputFile();
1036 
1037  this->getTrackerSectorStructs();
1038 
1039  this->bookHists();
1040 }
void ApeEstimatorSummary::bookHists ( )
private

Definition at line 220 of file ApeEstimatorSummary.cc.

References m_tkSector_, and residualErrorBinning().

Referenced by beginJob().

220  {
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 }
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
ApeEstimatorSummary::residualErrorBinning
std::vector< double > residualErrorBinning()
Definition: ApeEstimatorSummary.cc:250
void ApeEstimatorSummary::calculateApe ( )
private

Definition at line 300 of file ApeEstimatorSummary.cc.

References Abs(), alignBH_cfg::fixed, edm::ParameterSet::getParameter(), edm::detail::isnan(), align_cfg::iteration, HIPAlignmentAlgorithm_cfi::iterationFile, LogDebug, m_tkSector_, genParticles_cff::map, dataset::name, parameterSet_, funct::pow(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, wEntries, wEntriesOverSigmaX2, wInvalid, wUnity, anotherprimaryvertexanalyzer_cfi::xMax, anotherprimaryvertexanalyzer_cfi::xMin, CMSBoostedTauSeedingParameters_cfi::yMax, and CMSBoostedTauSeedingParameters_cfi::yMin.

Referenced by endJob().

300  {
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(nullptr);
308  TTree* baselineTreeX(nullptr);
309  TTree* baselineTreeY(nullptr);
310  TTree* sectorNameBaselineTree(nullptr);
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(nullptr);
338  TTree* iterationTreeY(nullptr);
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(nullptr);
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] = nullptr;
385  a_sectorBaselineName[iSector] = nullptr;
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.empty()){
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.empty()){
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 }
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
CMSBoostedTauSeedingParameters_cfi.yMin
yMin
Definition: CMSBoostedTauSeedingParameters_cfi.py:10
anotherprimaryvertexanalyzer_cfi.xMin
xMin
Definition: anotherprimaryvertexanalyzer_cfi.py:12
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
genParticles_cff.map
map
Definition: genParticles_cff.py:11
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
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
void ApeEstimatorSummary::endJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 1044 of file ApeEstimatorSummary.cc.

References calculateApe(), DEFINE_FWK_MODULE, and writeHists().

Referenced by o2olib.O2ORunMgr::executeJob().

1044  {
1045 
1046  this->calculateApe();
1047 
1048  this->writeHists();
1049 
1050 
1051 }
void ApeEstimatorSummary::getTrackerSectorStructs ( )
private

Definition at line 146 of file ApeEstimatorSummary.cc.

References TrackerSectorStruct::EntriesX, TrackerSectorStruct::EntriesY, mps_splice::entry, inputFile_, TrackerSectorStruct::isPixel, TrackerSectorStruct::m_binnedHists, m_tkSector_, TrackerSectorStruct::name, TrackerSectorStruct::Name, and TrackerSectorStruct::v_rawId.

Referenced by beginJob().

146  {
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(nullptr), *intervalDir(nullptr);
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(nullptr);
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(nullptr);
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 }
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
mps_splice.entry
entry
Definition: mps_splice.py:66
TrackerSectorStruct::m_binnedHists
std::map< unsigned int, std::map< std::string, TH1 * > > m_binnedHists
Definition: TrackerSectorStruct.h:105
TrackerSectorStruct::v_rawId
std::vector< unsigned int > v_rawId
Definition: TrackerSectorStruct.h:86
void ApeEstimatorSummary::openInputFile ( )
private

Definition at line 123 of file ApeEstimatorSummary.cc.

References edm::errors::Configuration, Exception, edm::ParameterSet::getParameter(), inputFile_, writeAntiElectronDiscrMVA_cfi::inputFileName, parameterSet_, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by beginJob().

123  {
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 }
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
writeAntiElectronDiscrMVA_cfi.inputFileName
inputFileName
Definition: writeAntiElectronDiscrMVA_cfi.py:43
ApeEstimatorSummary::parameterSet_
const edm::ParameterSet parameterSet_
Definition: ApeEstimatorSummary.cc:84
std::vector< double > ApeEstimatorSummary::residualErrorBinning ( )
private

Definition at line 250 of file ApeEstimatorSummary.cc.

References m_tkSector_.

Referenced by bookHists().

250  {
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 }
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
void ApeEstimatorSummary::writeHists ( )
private

Definition at line 261 of file ApeEstimatorSummary.cc.

References dir, TrackerOfflineValidation_Dqm_cff::dirName, edm::ParameterSet::getParameter(), m_tkSector_, parameterSet_, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by endJob().

261  {
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 }
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
TrackerOfflineValidation_Dqm_cff.dirName
dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:56
ApeEstimatorSummary::parameterSet_
const edm::ParameterSet parameterSet_
Definition: ApeEstimatorSummary.cc:84
dir
dbl *** dir
Definition: mlp_gen.cc:35

Member Data Documentation

bool ApeEstimatorSummary::firstEvent =true
private

Definition at line 86 of file ApeEstimatorSummary.cc.

Referenced by analyze(), looper.Looper::loop(), and Printer.Printer::process().

TFile* ApeEstimatorSummary::inputFile_
private

Definition at line 89 of file ApeEstimatorSummary.cc.

Referenced by getTrackerSectorStructs(), and openInputFile().

std::map<unsigned int, TrackerSectorStruct> ApeEstimatorSummary::m_tkSector_
private

Definition at line 91 of file ApeEstimatorSummary.cc.

Referenced by analyze(), bookHists(), calculateApe(), getTrackerSectorStructs(), residualErrorBinning(), and writeHists().

const edm::ParameterSet ApeEstimatorSummary::parameterSet_
private

Definition at line 84 of file ApeEstimatorSummary.cc.

Referenced by analyze(), calculateApe(), openInputFile(), and writeHists().

Generated for CMSSW Reference Manual by doxygen 1.8.11