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::one::EDAnalyzer<> edm::one::EDAnalyzerBase edm::EDConsumerBase

Public Member Functions

 ApeEstimatorSummary (const edm::ParameterSet &)
 
 ~ApeEstimatorSummary () override
 
- Public Member Functions inherited from edm::one::EDAnalyzer<>
 EDAnalyzer ()=default
 
SerialTaskQueue * globalLuminosityBlocksQueue () final
 
SerialTaskQueue * globalRunsQueue () final
 
bool wantsGlobalLuminosityBlocks () const final
 
bool wantsGlobalRuns () const final
 
- Public Member Functions inherited from edm::one::EDAnalyzerBase
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzerBase ()
 
ModuleDescription const & moduleDescription () const
 
bool wantsStreamLuminosityBlocks () const
 
bool wantsStreamRuns () const
 
 ~EDAnalyzerBase () 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
 
ProductResolverIndexAndSkipBit uncheckedIndexFrom (EDGetToken) 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_
 
unsigned int noSectors_
 
const edm::ParameterSet parameterSet_
 

Additional Inherited Members

- Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase ModuleType
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::one::EDAnalyzerBase
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &descriptions)
 
- 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 106 of file ApeEstimatorSummary.cc.

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

Definition at line 113 of file ApeEstimatorSummary.cc.

114 {
115 }

Member Function Documentation

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

Definition at line 837 of file ApeEstimatorSummary.cc.

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

838 {
839  // Load APEs from the GT and write them to root files similar to the ones from calculateAPE()
840 
841  if(firstEvent){
842  // Set baseline or calculate APE value?
843  const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
844 
845  edm::ESHandle<AlignmentErrorsExtended> alignmentErrors;
846  iSetup.get<TrackerAlignmentErrorExtendedRcd>().get(alignmentErrors);
847 
848  // Read in baseline file for calculation of APE value (if not setting baseline)
849  // Has same format as iterationFile
850  const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
851  TFile* baselineFile(nullptr);
852  TTree* sectorNameBaselineTree(nullptr);
853  if(!setBaseline){
854  std::ifstream baselineFileStream;
855  // Check if baseline file exists
856  baselineFileStream.open(baselineFileName.c_str());
857  if(baselineFileStream.is_open()){
858  baselineFileStream.close();
859  baselineFile = new TFile(baselineFileName.c_str(),"READ");
860  }
861  if(baselineFile){
862  edm::LogInfo("CalculateAPE")<<"Baseline file for APE values sucessfully opened";
863  baselineFile->GetObject("nameTree;1",sectorNameBaselineTree);
864  }else{
865  edm::LogWarning("CalculateAPE")<<"There is NO baseline file for APE values, so normalized residual width =1 for ideal conditions is assumed";
866  }
867  }
868 
869  // Set up root file for default APE values
870  const std::string defaultFileName(parameterSet_.getParameter<std::string>("DefaultFile"));
871  TFile* defaultFile = new TFile(defaultFileName.c_str(),"RECREATE");
872 
873  // Naming in the root files has to be iterTreeX to be consistent for the plotting tool
874  TTree* defaultTreeX(nullptr);
875  TTree* defaultTreeY(nullptr);
876  defaultFile->GetObject("iterTreeX;1",defaultTreeX);
877  defaultFile->GetObject("iterTreeY;1",defaultTreeY);
878  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as defaultTree)
879  TTree* sectorNameTree(nullptr);
880  defaultFile->GetObject("nameTree;1",sectorNameTree);
881 
882  bool firstIter(false);
883  if(!defaultTreeX){ // should be always true in setBaseline mode, since file is recreated
884  firstIter = true;
885  defaultTreeX = new TTree("iterTreeX","Tree for default APE x values from GT");
886  defaultTreeY = new TTree("iterTreeY","Tree for default APE y values from GT");
887  edm::LogInfo("CalculateAPE")<<"First APE iteration (number 0.), create default file with TTree";
888  sectorNameTree = new TTree("nameTree","Tree with names of sectors");
889  }else{
890  edm::LogWarning("CalculateAPE")<<"NOT the first APE iteration (number 0.), is this wanted or forgot to delete old iteration file with TTree?";
891  }
892 
893  // Assign the information stored in the trees to arrays
894  double a_defaultSectorX[noSectors_];
895  double a_defaultSectorY[noSectors_];
896 
897  std::string* a_sectorName[noSectors_];
898  std::string* a_sectorBaselineName[noSectors_];
899  for(auto const & i_sector : m_tkSector_){
900  const unsigned int iSector(i_sector.first);
901  const bool pixelSector(i_sector.second.isPixel);
902 
903  a_defaultSectorX[iSector] = -99.;
904  a_defaultSectorY[iSector] = -99.;
905 
906  a_sectorName[iSector] = nullptr;
907  a_sectorBaselineName[iSector] = nullptr;
908  std::stringstream ss_sector, ss_sectorSuffixed;
909  ss_sector << "Ape_Sector_" << iSector;
910  if(!setBaseline && sectorNameBaselineTree){
911  sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
912  sectorNameBaselineTree->GetEntry(0);
913  }
914 
915  if(firstIter){ // should be always true in setBaseline mode, since file is recreated
916  ss_sectorSuffixed << ss_sector.str() << "/D";
917  defaultTreeX->Branch(ss_sector.str().c_str(), &a_defaultSectorX[iSector], ss_sectorSuffixed.str().c_str());
918  if(pixelSector){
919  defaultTreeY->Branch(ss_sector.str().c_str(), &a_defaultSectorY[iSector], ss_sectorSuffixed.str().c_str());
920  }
921  sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
922  }else{
923  defaultTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorX[iSector]);
924  defaultTreeX->GetEntry(0);
925  if(pixelSector){
926  defaultTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorY[iSector]);
927  defaultTreeY->GetEntry(0);
928  }
929  sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
930  sectorNameTree->GetEntry(0);
931  }
932  }
933 
934 
935  // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
936  for(auto & i_sector : m_tkSector_){
937  const std::string& name(i_sector.second.name);
938  if(!firstIter){
939  const std::string& nameLastIter(*a_sectorName[i_sector.first]);
940  if(name!=nameLastIter){
941  edm::LogError("CalculateAPE")<<"Inconsistent sector definition in iterationFile for sector "<<i_sector.first<<",\n"
942  <<"Recent iteration has name \""<<name<<"\", while previous had \""<<nameLastIter<<"\"\n"
943  <<"...APE calculation stopped. Please check sector definitions in config!\n";
944  return;
945  }
946  }
947  else a_sectorName[i_sector.first] = new std::string(name);
948  if(!setBaseline && baselineFile){
949  const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
950  if(name!=nameBaseline){
951  edm::LogError("CalculateAPE")<<"Inconsistent sector definition in baselineFile for sector "<<i_sector.first<<",\n"
952  <<"Recent iteration has name \""<<name<<"\", while baseline had \""<<nameBaseline<<"\"\n"
953  <<"...APE calculation stopped. Please check sector definitions in config!\n";
954  return;
955  }
956  }
957  }
958 
959  // Loop over sectors for calculating getting default APE
960  for(auto & i_sector : m_tkSector_){
961  double defaultApeX(0.);
962  double defaultApeY(0.);
963  unsigned int nModules(0);
964  for(auto const & i_rawId : i_sector.second.v_rawId){
965  std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
966  for(std::vector<AlignTransformErrorExtended>::const_iterator i_alignError = alignErrors.begin(); i_alignError != alignErrors.end(); ++i_alignError){
967  if(i_rawId == i_alignError->rawId()){
968  CLHEP::HepSymMatrix errMatrix = i_alignError->matrix();
969  defaultApeX += errMatrix[0][0];
970  defaultApeY += errMatrix[1][1];
971  nModules++;
972  }
973  }
974  }
975  a_defaultSectorX[i_sector.first] = defaultApeX/nModules;
976  a_defaultSectorY[i_sector.first] = defaultApeY/nModules;
977  }
978  sectorNameTree->Fill();
979  sectorNameTree->Write("nameTree");
980  defaultTreeX->Fill();
981  defaultTreeX->Write("iterTreeX");
982  defaultTreeY->Fill();
983  defaultTreeY->Write("iterTreeY");
984 
985  defaultFile->Close();
986  delete defaultFile;
987 
988  if(baselineFile){
989  baselineFile->Close();
990  delete baselineFile;
991  }
992 
993  firstEvent = false;
994  }
995 }
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
ApeEstimatorSummary::parameterSet_
const edm::ParameterSet parameterSet_
Definition: ApeEstimatorSummary.cc:84
edm::EventSetup::get
T get() const
Definition: EventSetup.h:63
void ApeEstimatorSummary::beginJob ( )
overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 1000 of file ApeEstimatorSummary.cc.

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

1001 {
1002  this->openInputFile();
1003 
1004  this->getTrackerSectorStructs();
1005 
1006  this->bookHists();
1007 }
void ApeEstimatorSummary::bookHists ( )
private

Definition at line 222 of file ApeEstimatorSummary.cc.

References m_tkSector_, and residualErrorBinning().

Referenced by beginJob().

222  {
223  const std::vector<double> v_binX(this->residualErrorBinning());
224  for(auto & i_sector : m_tkSector_){
225  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]));
226  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]));
227  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]));
228  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]));
229  i_sector.second.ResidualWidthX1 = new TH1F("h_residualWidthX1","residual width #Delta_{x};#sigma_{x} [#mum];#Delta_{x}",v_binX.size()-1,&(v_binX[0]));
230  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]));
231  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]));
232  i_sector.second.ResidualWidthX2 = new TH1F("h_residualWidthX2","residual width #Delta_{x};#sigma_{x} [#mum];#Delta_{x}",v_binX.size()-1,&(v_binX[0]));
233  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]));
234 
235  if(i_sector.second.isPixel){
236  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]));
237  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]));
238  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]));
239  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]));
240  i_sector.second.ResidualWidthY1 = new TH1F("h_residualWidthY1","residual width #Delta_{y};#sigma_{y} [#mum];#Delta_{y}",v_binX.size()-1,&(v_binX[0]));
241  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]));
242  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]));
243  i_sector.second.ResidualWidthY2 = new TH1F("h_residualWidthY2","residual width #Delta_{y};#sigma_{y} [#mum];#Delta_{y}",v_binX.size()-1,&(v_binX[0]));
244  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]));
245  }
246  }
247 }
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
ApeEstimatorSummary::residualErrorBinning
std::vector< double > residualErrorBinning()
Definition: ApeEstimatorSummary.cc:252
void ApeEstimatorSummary::calculateApe ( )
private

Definition at line 302 of file ApeEstimatorSummary.cc.

References Abs(), alignBH_cfg::fixed, edm::ParameterSet::getParameter(), edm::detail::isnan(), align_cfg::iteration, HIPAlignmentAlgorithm_cfi::iterationFile, LogDebug, m_tkSector_, dataset::name, noSectors_, parameterSet_, 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().

302  {
303  // Set baseline or calculate APE value?
304  const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));
305 
306  // Read in baseline file for calculation of APE value (if not setting baseline)
307  // Has same format as iterationFile
308  const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
309  TFile* baselineFile(nullptr);
310  TTree* baselineTreeX(nullptr);
311  TTree* baselineTreeY(nullptr);
312  TTree* sectorNameBaselineTree(nullptr);
313  if(!setBaseline){
314  std::ifstream baselineFileStream;
315  // Check if baseline file exists
316  baselineFileStream.open(baselineFileName.c_str());
317  if(baselineFileStream.is_open()){
318  baselineFileStream.close();
319  baselineFile = new TFile(baselineFileName.c_str(),"READ");
320  }
321  if(baselineFile){
322  edm::LogInfo("CalculateAPE")<<"Baseline file for APE values sucessfully opened";
323  baselineFile->GetObject("iterTreeX;1",baselineTreeX);
324  baselineFile->GetObject("iterTreeY;1",baselineTreeY);
325  baselineFile->GetObject("nameTree;1",sectorNameBaselineTree);
326  }else{
327  edm::LogWarning("CalculateAPE")<<"There is NO baseline file for APE values, so normalized residual width =1 for ideal conditions is assumed";
328  }
329  }
330 
331  // Set up root file for iterations on APE value (or for setting baseline in setBaseline mode)
332  const std::string iterationFileName(setBaseline ? baselineFileName : parameterSet_.getParameter<std::string>("IterationFile"));
333  // For iterations, updates are needed to not overwrite the iterations before
334  TFile* iterationFile = new TFile(iterationFileName.c_str(),setBaseline ? "RECREATE" : "UPDATE");
335 
336 
337  // Set up TTree for iterative APE values on first pass (first iteration) or read from file (further iterations)
338  TTree* iterationTreeX(nullptr);
339  TTree* iterationTreeY(nullptr);
340  iterationFile->GetObject("iterTreeX;1",iterationTreeX);
341  iterationFile->GetObject("iterTreeY;1",iterationTreeY);
342  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as iterationTree)
343  TTree* sectorNameTree(nullptr);
344  iterationFile->GetObject("nameTree;1",sectorNameTree);
345 
346  bool firstIter(false);
347  if(!iterationTreeX){ // should be always true in setBaseline mode, since file is recreated
348  firstIter = true;
349  if(!setBaseline){
350  iterationTreeX = new TTree("iterTreeX","Tree for APE x values of all iterations");
351  iterationTreeY = new TTree("iterTreeY","Tree for APE y values of all iterations");
352  edm::LogInfo("CalculateAPE")<<"First APE iteration (number 0.), create iteration file with TTree";
353  sectorNameTree = new TTree("nameTree","Tree with names of sectors");
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  }else{
361  const unsigned int iteration(iterationTreeX->GetEntries());
362  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?";
363  }
364 
365 
366  // Assign the information stored in the trees to arrays
367  double a_apeSectorX[noSectors_];
368  double a_apeSectorY[noSectors_];
369  double a_baselineSectorX[noSectors_];
370  double a_baselineSectorY[noSectors_];
371 
372  std::string* a_sectorName[noSectors_];
373  std::string* a_sectorBaselineName[noSectors_];
374  for(auto const & i_sector : m_tkSector_){
375  const unsigned int iSector(i_sector.first);
376  const bool pixelSector(i_sector.second.isPixel);
377  a_apeSectorX[iSector] = 99.;
378  a_apeSectorY[iSector] = 99.;
379  a_baselineSectorX[iSector] = -99.;
380  a_baselineSectorY[iSector] = -99.;
381 
382  a_sectorName[iSector] = nullptr;
383  a_sectorBaselineName[iSector] = nullptr;
384  std::stringstream ss_sector, ss_sectorSuffixed;
385  ss_sector << "Ape_Sector_" << iSector;
386  if(!setBaseline && baselineTreeX){
387  baselineTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorX[iSector]);
388  baselineTreeX->GetEntry(0);
389  if(pixelSector){
390  baselineTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorY[iSector]);
391  baselineTreeY->GetEntry(0);
392  }
393  sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
394  sectorNameBaselineTree->GetEntry(0);
395  }else{
396  // Set default ideal normalized residual width to 1
397  a_baselineSectorX[iSector] = 1.;
398  a_baselineSectorY[iSector] = 1.;
399  }
400  if(firstIter){ // should be always true in setBaseline mode, since file is recreated
401  ss_sectorSuffixed << ss_sector.str() << "/D";
402  iterationTreeX->Branch(ss_sector.str().c_str(), &a_apeSectorX[iSector], ss_sectorSuffixed.str().c_str());
403  if(pixelSector){
404  iterationTreeY->Branch(ss_sector.str().c_str(), &a_apeSectorY[iSector], ss_sectorSuffixed.str().c_str());
405  }
406  sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
407  }else{
408  iterationTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorX[iSector]);
409  iterationTreeX->GetEntry(iterationTreeX->GetEntries()-1);
410  if(pixelSector){
411  iterationTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorY[iSector]);
412  iterationTreeY->GetEntry(iterationTreeY->GetEntries()-1);
413  }
414  sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
415  sectorNameTree->GetEntry(0);
416  }
417  }
418 
419 
420  // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
421  for(auto & i_sector : m_tkSector_){
422  const std::string& name(i_sector.second.name);
423  if(!firstIter){
424  const std::string& nameLastIter(*a_sectorName[i_sector.first]);
425  if(name!=nameLastIter){
426  edm::LogError("CalculateAPE")<<"Inconsistent sector definition in iterationFile for sector "<<i_sector.first<<",\n"
427  <<"Recent iteration has name \""<<name<<"\", while previous had \""<<nameLastIter<<"\"\n"
428  <<"...APE calculation stopped. Please check sector definitions in config!\n";
429  return;
430  }
431  }else{
432  a_sectorName[i_sector.first] = new std::string(name);
433  }
434  if(!setBaseline && baselineFile){
435  const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
436  if(name!=nameBaseline){
437  edm::LogError("CalculateAPE")<<"Inconsistent sector definition in baselineFile for sector "<<i_sector.first<<",\n"
438  <<"Recent iteration has name \""<<name<<"\", while baseline had \""<<nameBaseline<<"\"\n"
439  <<"...APE calculation stopped. Please check sector definitions in config!\n";
440  return;
441  }
442  }
443  }
444 
445 
446  // Set up text file for writing out APE values per module
447  std::ofstream apeOutputFile;
448  if(!setBaseline){
449  const std::string apeOutputFileName(parameterSet_.getParameter<std::string>("ApeOutputFile"));
450  apeOutputFile.open(apeOutputFileName.c_str());
451  if(apeOutputFile.is_open()){
452  edm::LogInfo("CalculateAPE")<<"Text file for writing APE values successfully opened";
453  }else{
454  edm::LogError("CalculateAPE")<<"Text file for writing APE values NOT opened,\n"
455  <<"...APE calculation stopped. Please check path of text file name in config:\n"
456  <<"\t"<<apeOutputFileName;
457  return;
458  }
459  }
460 
461  // Loop over sectors for calculating APE
462  const std::string apeWeightName(parameterSet_.getParameter<std::string>("apeWeight"));
463  ApeWeight apeWeight(wInvalid);
464  if(apeWeightName=="unity") apeWeight = wUnity;
465  else if(apeWeightName=="entries")apeWeight = wEntries;
466  else if(apeWeightName=="entriesOverSigmaX2")apeWeight = wEntriesOverSigmaX2;
467  if(apeWeight==wInvalid){
468  edm::LogError("CalculateAPE")<<"Invalid parameter 'apeWeight' in cfg file: \""<<apeWeightName
469  <<"\"\nimplemented apeWeights are \"unity\", \"entries\", \"entriesOverSigmaX2\""
470  <<"\n...APE calculation stopped.";
471  return;
472  }
473  const double minHitsPerInterval(parameterSet_.getParameter<double>("minHitsPerInterval"));
474  const double sigmaFactorFit(parameterSet_.getParameter<double>("sigmaFactorFit"));
475  const double correctionScaling(parameterSet_.getParameter<double>("correctionScaling"));
476  const bool smoothIteration(parameterSet_.getParameter<bool>("smoothIteration"));
477  const double smoothFraction(parameterSet_.getParameter<double>("smoothFraction"));
478  if(smoothFraction<=0. || smoothFraction>1.){
479  edm::LogError("CalculateAPE")<<"Incorrect parameter in cfg file,"
480  <<"\nsmoothFraction has to be in [0,1], but is "<<smoothFraction
481  <<"\n...APE calculation stopped.";
482  return;
483  }
484  for(auto & i_sector : m_tkSector_){
485  typedef std::pair<double,double> Error2AndResidualWidth2PerBin;
486  typedef std::pair<double,Error2AndResidualWidth2PerBin> WeightAndResultsPerBin;
487  std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinX;
488  std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinY;
489 
490  double baselineWidthX2(a_baselineSectorX[i_sector.first]);
491  double baselineWidthY2(a_baselineSectorY[i_sector.first]);
492 
493  // Loop over residual error bins to calculate APE for every bin
494  for(auto const & i_errBins : i_sector.second.m_binnedHists){
495  std::map<std::string,TH1*> mHists = i_errBins.second;
496 
497  double entriesX = mHists["sigmaX"]->GetEntries();
498  double meanSigmaX = mHists["sigmaX"]->GetMean();
499 
500  // Fitting Parameters
501  double xMin = mHists["norResX"]->GetXaxis()->GetXmin();
502  double xMax = mHists["norResX"]->GetXaxis()->GetXmax();
503  double integralX = mHists["norResX"]->Integral();
504  double meanX = mHists["norResX"]->GetMean();
505  double rmsX = mHists["norResX"]->GetRMS();
506  double maximumX = mHists["norResX"]->GetBinContent(mHists["norResX"]->GetMaximumBin());
507 
508  // First Gaus Fit
509  TF1 funcX_1("mygausX", "gaus", xMin, xMax);
510  funcX_1.SetParameters(maximumX, meanX, rmsX);
511  TString fitOpt("ILERQ"); //TString fitOpt("IMR"); ("IRLL"); ("IRQ");
512  if(integralX>minHitsPerInterval){
513  if(mHists["norResX"]->Fit(&funcX_1, fitOpt)){
514  edm::LogWarning("CalculateAPE")<<"Fit1 did not work: "<<mHists["norResX"]->Fit(&funcX_1, fitOpt);
515  continue;
516  }
517  LogDebug("CalculateAPE")<<"FitResultX1\t"<<mHists["norResX"]->Fit(&funcX_1, fitOpt)<<"\n";
518  }
519  double meanX_1 = funcX_1.GetParameter(1);
520  double sigmaX_1 = funcX_1.GetParameter(2);
521 
522  // Second gaus fit
523  TF1 funcX_2("mygausX2","gaus",meanX_1 - sigmaFactorFit*TMath::Abs(sigmaX_1), meanX_1 + sigmaFactorFit*TMath::Abs(sigmaX_1));
524  funcX_2.SetParameters(funcX_1.GetParameter(0),meanX_1,sigmaX_1);
525  if(integralX>minHitsPerInterval){
526  if(mHists["norResX"]->Fit(&funcX_2, fitOpt)){
527  edm::LogWarning("CalculateAPE")<<"Fit2 did not work for x : "<<mHists["norResX"]->Fit(&funcX_2, fitOpt);
528  continue;
529  }
530  LogDebug("CalculateAPE")<<"FitResultX2\t"<<mHists["norResX"]->Fit(&funcX_2, fitOpt)<<"\n";
531  }
532  double meanX_2 = funcX_2.GetParameter(1);
533  double sigmaX_2 = funcX_2.GetParameter(2);
534 
535  // Now the same for y coordinate
536  double entriesY(0.);
537  double meanSigmaY(0.);
538  if(i_sector.second.isPixel){
539  entriesY = mHists["sigmaY"]->GetEntries();
540  meanSigmaY = mHists["sigmaY"]->GetMean();
541  }
542 
543  double meanY_1(0.);
544  double sigmaY_1(0.);
545  double meanY_2(0.);
546  double sigmaY_2(0.);
547  double meanY(0.);
548  double rmsY(0.);
549  if(i_sector.second.isPixel){
550  // Fitting Parameters
551  double yMin = mHists["norResY"]->GetXaxis()->GetXmin();
552  double yMax = mHists["norResY"]->GetXaxis()->GetXmax();
553  double integralY = mHists["norResY"]->Integral();
554  meanY = mHists["norResY"]->GetMean();
555  rmsY = mHists["norResY"]->GetRMS();
556  double maximumY = mHists["norResY"]->GetBinContent(mHists["norResY"]->GetMaximumBin());
557 
558  // First Gaus Fit
559  TF1 funcY_1("mygausY", "gaus", yMin, yMax);
560  funcY_1.SetParameters(maximumY, meanY, rmsY);
561  if(integralY>minHitsPerInterval){
562  if(mHists["norResY"]->Fit(&funcY_1, fitOpt)){
563  edm::LogWarning("CalculateAPE")<<"Fit1 did not work: "<<mHists["norResY"]->Fit(&funcY_1, fitOpt);
564  continue;
565  }
566  LogDebug("CalculateAPE")<<"FitResultY1\t"<<mHists["norResY"]->Fit(&funcY_1, fitOpt)<<"\n";
567  }
568  meanY_1 = funcY_1.GetParameter(1);
569  sigmaY_1 = funcY_1.GetParameter(2);
570 
571  // Second gaus fit
572  TF1 funcY_2("mygausY2","gaus",meanY_1 - sigmaFactorFit*TMath::Abs(sigmaY_1), meanY_1 + sigmaFactorFit*TMath::Abs(sigmaY_1));
573  funcY_2.SetParameters(funcY_1.GetParameter(0),meanY_1,sigmaY_1);
574  if(integralY>minHitsPerInterval){
575  if(mHists["norResY"]->Fit(&funcY_2, fitOpt)){
576  edm::LogWarning("CalculateAPE")<<"Fit2 did not work for y : "<<mHists["norResY"]->Fit(&funcY_2, fitOpt);
577  continue;
578  }
579  LogDebug("CalculateAPE")<<"FitResultY2\t"<<mHists["norResY"]->Fit(&funcY_2, fitOpt)<<"\n";
580  }
581  meanY_2 = funcY_2.GetParameter(1);
582  sigmaY_2 = funcY_2.GetParameter(2);
583  }
584 
585  // Fill histograms
586  double fitMeanX_1(meanX_1), fitMeanX_2(meanX_2);
587  double residualWidthX_1(sigmaX_1), residualWidthX_2(sigmaX_2);
588  double fitMeanY_1(meanY_1), fitMeanY_2(meanY_2);
589  double residualWidthY_1(sigmaY_1), residualWidthY_2(sigmaY_2);
590 
591  double correctionX2_1(-0.0010), correctionX2_2(-0.0010);
592  double correctionY2_1(-0.0010), correctionY2_2(-0.0010);
593  correctionX2_1 = meanSigmaX*meanSigmaX*(residualWidthX_1*residualWidthX_1 -baselineWidthX2);
594  correctionX2_2 = meanSigmaX*meanSigmaX*(residualWidthX_2*residualWidthX_2 -baselineWidthX2);
595  if(i_sector.second.isPixel){
596  correctionY2_1 = meanSigmaY*meanSigmaY*(residualWidthY_1*residualWidthY_1 -baselineWidthY2);
597  correctionY2_2 = meanSigmaY*meanSigmaY*(residualWidthY_2*residualWidthY_2 -baselineWidthY2);
598  }
599 
600  double correctionX_1 = correctionX2_1>=0. ? std::sqrt(correctionX2_1) : -std::sqrt(-correctionX2_1);
601  double correctionX_2 = correctionX2_2>=0. ? std::sqrt(correctionX2_2) : -std::sqrt(-correctionX2_2);
602  double correctionY_1 = correctionY2_1>=0. ? std::sqrt(correctionY2_1) : -std::sqrt(-correctionY2_1);
603  double correctionY_2 = correctionY2_2>=0. ? std::sqrt(correctionY2_2) : -std::sqrt(-correctionY2_2);
604  // Meanwhile, this got very bad default values, or? (negative corrections allowed)
605  if(isnan(correctionX_1))correctionX_1 = -0.0010;
606  if(isnan(correctionX_2))correctionX_2 = -0.0010;
607  if(isnan(correctionY_1))correctionY_1 = -0.0010;
608  if(isnan(correctionY_2))correctionY_2 = -0.0010;
609 
610  if(entriesX<minHitsPerInterval){
611  meanX = 0.; rmsX = -0.0010;
612  fitMeanX_1 = 0.; correctionX_1 = residualWidthX_1 = -0.0010;
613  fitMeanX_2 = 0.; correctionX_2 = residualWidthX_2 = -0.0010;
614  }
615 
616  if(entriesY<minHitsPerInterval){
617  meanY = 0.; rmsY = -0.0010;
618  fitMeanY_1 = 0.; correctionY_1 = residualWidthY_1 = -0.0010;
619  fitMeanY_2 = 0.; correctionY_2 = residualWidthY_2 = -0.0010;
620  }
621 
622  i_sector.second.MeanX ->SetBinContent(i_errBins.first,meanX);
623  i_sector.second.RmsX ->SetBinContent(i_errBins.first,rmsX);
624 
625  i_sector.second.FitMeanX1 ->SetBinContent(i_errBins.first,fitMeanX_1);
626  i_sector.second.ResidualWidthX1->SetBinContent(i_errBins.first,residualWidthX_1);
627  i_sector.second.CorrectionX1 ->SetBinContent(i_errBins.first,correctionX_1*10000.);
628 
629  i_sector.second.FitMeanX2 ->SetBinContent(i_errBins.first,fitMeanX_2);
630  i_sector.second.ResidualWidthX2->SetBinContent(i_errBins.first,residualWidthX_2);
631  i_sector.second.CorrectionX2 ->SetBinContent(i_errBins.first,correctionX_2*10000.);
632 
633  if(i_sector.second.isPixel){
634  i_sector.second.MeanY ->SetBinContent(i_errBins.first,meanY);
635  i_sector.second.RmsY ->SetBinContent(i_errBins.first,rmsY);
636 
637  i_sector.second.FitMeanY1 ->SetBinContent(i_errBins.first,fitMeanY_1);
638  i_sector.second.ResidualWidthY1->SetBinContent(i_errBins.first,residualWidthY_1);
639  i_sector.second.CorrectionY1 ->SetBinContent(i_errBins.first,correctionY_1*10000.);
640 
641  i_sector.second.FitMeanY2 ->SetBinContent(i_errBins.first,fitMeanY_2);
642  i_sector.second.ResidualWidthY2->SetBinContent(i_errBins.first,residualWidthY_2);
643  i_sector.second.CorrectionY2 ->SetBinContent(i_errBins.first,correctionY_2*10000.);
644  }
645 
646 
647  // Use result for bin only when entries>=minHitsPerInterval
648  if(entriesX<minHitsPerInterval && entriesY<minHitsPerInterval)continue;
649 
650  double weightX(0.);
651  double weightY(0.);
652  if(apeWeight==wUnity){
653  weightX = 1.;
654  weightY = 1.;
655  }else if(apeWeight==wEntries){
656  weightX = entriesX;
657  weightY = entriesY;
658  }else if(apeWeight==wEntriesOverSigmaX2){
659  weightX = entriesX/(meanSigmaX*meanSigmaX);
660  weightY = entriesY/(meanSigmaY*meanSigmaY);
661  }
662 
663  const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinX(meanSigmaX*meanSigmaX, residualWidthX_1*residualWidthX_1);
664  const WeightAndResultsPerBin weightAndResultsPerBinX(weightX, error2AndResidualWidth2PerBinX);
665  if(!(entriesX<minHitsPerInterval)){
666  //Fill absolute weights
667  i_sector.second.WeightX->SetBinContent(i_errBins.first,weightX);
668  v_weightAndResultsPerBinX.push_back(weightAndResultsPerBinX);
669  }
670 
671  const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinY(meanSigmaY*meanSigmaY, residualWidthY_1*residualWidthY_1);
672  const WeightAndResultsPerBin weightAndResultsPerBinY(weightY, error2AndResidualWidth2PerBinY);
673  if(!(entriesY<minHitsPerInterval)){
674  //Fill absolute weights
675  i_sector.second.WeightY->SetBinContent(i_errBins.first,weightY);
676  v_weightAndResultsPerBinY.push_back(weightAndResultsPerBinY);
677  }
678  }
679 
680 
681  // Do the final calculations
682 
683  if(v_weightAndResultsPerBinX.empty()){
684  edm::LogError("CalculateAPE")<<"NO error interval of sector "<<i_sector.first<<" has a valid x APE calculated,\n...so cannot set APE";
685  continue;
686  }
687 
688  if(i_sector.second.isPixel && v_weightAndResultsPerBinY.empty()){
689  edm::LogError("CalculateAPE")<<"NO error interval of sector "<<i_sector.first<<" has a valid y APE calculated,\n...so cannot set APE";
690  continue;
691  }
692 
693  double correctionX2(999.);
694  double correctionY2(999.);
695 
696  // Get sum of all weights
697  double weightSumX(0.);
698  for(auto const & i_apeBin : v_weightAndResultsPerBinX){
699  weightSumX += i_apeBin.first;
700  }
701  i_sector.second.WeightX->Scale(1/weightSumX);
702  double weightSumY(0.);
703  if(i_sector.second.isPixel){
704  for(auto const & i_apeBin : v_weightAndResultsPerBinY){
705  weightSumY += i_apeBin.first;
706  }
707  i_sector.second.WeightY->Scale(1/weightSumY);
708  }
709 
710  if(!setBaseline){
711  // Calculate weighted mean
712  bool firstIntervalX(true);
713  for(auto const & i_apeBin : v_weightAndResultsPerBinX){
714  if(firstIntervalX){
715  correctionX2 = i_apeBin.first*i_apeBin.second.first*(i_apeBin.second.second - baselineWidthX2);
716  firstIntervalX = false;
717  }else{
718  correctionX2 += i_apeBin.first*i_apeBin.second.first*(i_apeBin.second.second - baselineWidthX2);
719  }
720  }
721  correctionX2 = correctionX2/weightSumX;
722  }else{
723  double numeratorX(0.), denominatorX(0.);
724  for(auto const & i_apeBin : v_weightAndResultsPerBinX){
725  numeratorX += i_apeBin.first*i_apeBin.second.first*i_apeBin.second.second;
726  denominatorX += i_apeBin.first*i_apeBin.second.first;
727  }
728  correctionX2 = numeratorX/denominatorX;
729  }
730 
731  if(i_sector.second.isPixel){
732  if(!setBaseline){
733  // Calculate weighted mean
734  bool firstIntervalY(true);
735  for(auto const & i_apeBin : v_weightAndResultsPerBinY){
736  if(firstIntervalY){
737  correctionY2 = i_apeBin.first*i_apeBin.second.first*(i_apeBin.second.second - baselineWidthY2);
738  firstIntervalY = false;
739  }else{
740  correctionY2 += i_apeBin.first*i_apeBin.second.first*(i_apeBin.second.second - baselineWidthY2);
741  }
742  }
743  correctionY2 = correctionY2/weightSumY;
744  }else{
745  double numeratorY(0.), denominatorY(0.);
746  for(auto const & i_apeBin : v_weightAndResultsPerBinY){
747  numeratorY += i_apeBin.first*i_apeBin.second.first*i_apeBin.second.second;
748  denominatorY += i_apeBin.first*i_apeBin.second.first;
749  }
750  correctionY2 = numeratorY/denominatorY;
751  }
752  }
753 
754  if(!setBaseline){
755  // Calculate updated squared APE of current iteration
756  double apeX2(999.);
757  double apeY2(999.);
758 
759  // old APE value from last iteration
760  if(firstIter){
761  apeX2 = 0.;
762  apeY2 = 0.;
763  }
764  else{
765  apeX2 = a_apeSectorX[i_sector.first];
766  apeY2 = a_apeSectorY[i_sector.first];
767  }
768  const double apeX2old = apeX2;
769  const double apeY2old = apeY2;
770 
771  // scale APE Correction with value given in cfg (not if smoothed iteration is used)
772  edm::LogInfo("CalculateAPE")<<"Unscaled correction x for sector "<<i_sector.first<<" is "<<(correctionX2>0. ? +1. : -1.)*std::sqrt(std::fabs(correctionX2));
773  if(!smoothIteration || firstIter)correctionX2 = correctionX2*correctionScaling*correctionScaling;
774  if(i_sector.second.isPixel){
775  edm::LogInfo("CalculateAPE")<<"Unscaled correction y for sector "<<i_sector.first<<" is "<<(correctionY2>0. ? +1. : -1.)*std::sqrt(std::fabs(correctionY2));
776  if(!smoothIteration || firstIter)correctionY2 = correctionY2*correctionScaling*correctionScaling;
777  }
778 
779  // new APE value
780  // smooth iteration or not?
781  if(apeX2 + correctionX2 < 0.) correctionX2 = -apeX2;
782  if(apeY2 + correctionY2 < 0.) correctionY2 = -apeY2;
783  const double apeX2new(apeX2old + correctionX2);
784  const double apeY2new(apeY2old + correctionY2);
785  if(!smoothIteration || firstIter){
786  apeX2 = apeX2new;
787  apeY2 = apeY2new;
788  }else{
789  const double apeXtmp = smoothFraction*std::sqrt(apeX2new) + (1-smoothFraction)*std::sqrt(apeX2old);
790  const double apeYtmp = smoothFraction*std::sqrt(apeY2new) + (1-smoothFraction)*std::sqrt(apeY2old);
791  apeX2 = apeXtmp*apeXtmp;
792  apeY2 = apeYtmp*apeYtmp;
793  }
794  if(apeX2<0. || apeY2<0.)edm::LogError("CalculateAPE")<<"\n\n\tBad APE, but why???\n\n\n";
795  a_apeSectorX[i_sector.first] = apeX2;
796  a_apeSectorY[i_sector.first] = apeY2;
797 
798  // Set the calculated APE spherical for all modules of strip sectors
799  const double apeX(std::sqrt(apeX2));
800  const double apeY(std::sqrt(apeY2));
801  const double apeZ(std::sqrt(0.5*(apeX2+apeY2)));
802  for(auto const & i_rawId : i_sector.second.v_rawId){
803  if(i_sector.second.isPixel){
804  apeOutputFile<<i_rawId<<" "<<std::fixed<<std::setprecision(5)<<apeX<<" "<<apeY<<" "<<apeZ<<"\n";
805  }else{
806  apeOutputFile<<i_rawId<<" "<<std::fixed<<std::setprecision(5)<<apeX<<" "<<apeX<<" "<<apeX<<"\n";
807  }
808  }
809  }else{ // In setBaseline mode, just fill estimated mean value of residual width
810  a_apeSectorX[i_sector.first] = correctionX2;
811  a_apeSectorY[i_sector.first] = correctionY2;
812  }
813  }
814  if(!setBaseline)apeOutputFile.close();
815 
816  iterationTreeX->Fill();
817  iterationTreeX->Write("iterTreeX", TObject::kOverwrite); // TObject::kOverwrite needed to not produce another iterTreeX;2
818  iterationTreeY->Fill();
819  iterationTreeY->Write("iterTreeY", TObject::kOverwrite); // TObject::kOverwrite needed to not produce another iterTreeY;2
820  if(firstIter){
821  sectorNameTree->Fill();
822  sectorNameTree->Write("nameTree");
823  }
824  iterationFile->Close();
825  delete iterationFile;
826 
827  if(baselineFile){
828  baselineFile->Close();
829  delete baselineFile;
830  }
831 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:601
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
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
void ApeEstimatorSummary::endJob ( )
overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 1011 of file ApeEstimatorSummary.cc.

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

Referenced by o2olib.O2ORunMgr::executeJob().

1011  {
1012 
1013  this->calculateApe();
1014 
1015  this->writeHists();
1016 
1017 
1018 }
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, noSectors_, 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  unsigned int iSector(1);
154  for(;sectorBool;++iSector){
155  std::stringstream sectorName, fullSectorName;
156  sectorName << "Sector_" << iSector << "/";
157  fullSectorName << pluginName << sectorName.str();
158  TString fullName(fullSectorName.str().c_str());
159  inputFile_->cd();
160  sectorDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
161  if(sectorDir){
162  TrackerSectorStruct tkSector;
163  sectorDir->GetObject("z_name;1", tkSector.Name);
164  tkSector.name = tkSector.Name->GetTitle();
165  bool intervalBool(true);
166  for(unsigned int iInterval(1);intervalBool;++iInterval){
167  std::stringstream intervalName, fullIntervalName;
168  intervalName << "Interval_" << iInterval <<"/";
169  fullIntervalName << fullSectorName.str() << intervalName.str();
170  fullName = fullIntervalName.str().c_str();
171  intervalDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
172  if(intervalDir){
173  intervalDir->GetObject("h_sigmaX;1", tkSector.m_binnedHists[iInterval]["sigmaX"]);
174  intervalDir->GetObject("h_norResX;1", tkSector.m_binnedHists[iInterval]["norResX"]);
175  intervalDir->GetObject("h_sigmaY;1", tkSector.m_binnedHists[iInterval]["sigmaY"]);
176  intervalDir->GetObject("h_norResY;1", tkSector.m_binnedHists[iInterval]["norResY"]);
177  if(tkSector.m_binnedHists[iInterval]["sigmaY"] && tkSector.m_binnedHists[iInterval]["norResY"]){
178  tkSector.isPixel = true;
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 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(auto const & i_rawId : tkSector.v_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  }else{
210  sectorBool = false;
211  edm::LogInfo("CalculateAPE")<<"There are "<<iSector-1<<" sectors defined in input file";
212  }
213  }
214  noSectors_ = iSector+1;
215 
216 
217 }
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 124 of file ApeEstimatorSummary.cc.

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

Referenced by beginJob().

124  {
125  const std::string inputFileName(parameterSet_.getParameter<std::string>("InputFile"));
126  std::ifstream inputFileStream;
127  // Check if baseline file exists
128  inputFileStream.open(inputFileName.c_str());
129  if(inputFileStream.is_open()){
130  inputFile_ = new TFile(inputFileName.c_str(),"READ");
131  }
132  if(inputFile_){
133  edm::LogInfo("CalculateAPE")<<"Input file from loop over tracks and corresponding hits sucessfully opened";
134  }
135  else{
136  edm::LogError("CalculateAPE")<<"There is NO input file\n"
137  <<"...APE calculation stopped. Please check path of input file name in config:\n"
138  <<"\t"<<inputFileName;
139  throw edm::Exception( edm::errors::Configuration, "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 252 of file ApeEstimatorSummary.cc.

References m_tkSector_.

Referenced by bookHists().

252  {
253  std::vector<double> v_binX;
254  TH1* EntriesX(m_tkSector_[1].EntriesX);
255  for(int iBin=1; iBin<=EntriesX->GetNbinsX()+1; ++iBin){
256  v_binX.push_back(EntriesX->GetBinLowEdge(iBin));
257  }
258  return v_binX;
259 }
ApeEstimatorSummary::m_tkSector_
std::map< unsigned int, TrackerSectorStruct > m_tkSector_
Definition: ApeEstimatorSummary.cc:91
void ApeEstimatorSummary::writeHists ( )
private

Definition at line 263 of file ApeEstimatorSummary.cc.

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

Referenced by endJob().

263  {
264  TFile* resultsFile = new TFile(parameterSet_.getParameter<std::string>("ResultsFile").c_str(), "RECREATE");
265  TDirectory* baseDir = resultsFile->mkdir("ApeEstimatorSummary");
266  for(auto const & i_sector : m_tkSector_){
267  std::stringstream dirName;
268  dirName<<"Sector_" << i_sector.first;
269  TDirectory* dir = baseDir->mkdir(dirName.str().c_str());
270  dir->cd();
271 
272  i_sector.second.Name->Write();
273 
274  i_sector.second.WeightX->Write();
275  i_sector.second.MeanX->Write();
276  i_sector.second.RmsX->Write();
277  i_sector.second.FitMeanX1->Write();
278  i_sector.second.ResidualWidthX1 ->Write();
279  i_sector.second.CorrectionX1->Write();
280  i_sector.second.FitMeanX2->Write();
281  i_sector.second.ResidualWidthX2->Write();
282  i_sector.second.CorrectionX2->Write();
283 
284  if(i_sector.second.isPixel){
285  i_sector.second.WeightY->Write();
286  i_sector.second.MeanY->Write();
287  i_sector.second.RmsY->Write();
288  i_sector.second.FitMeanY1->Write();
289  i_sector.second.ResidualWidthY1 ->Write();
290  i_sector.second.CorrectionY1->Write();
291  i_sector.second.FitMeanY2->Write();
292  i_sector.second.ResidualWidthY2->Write();
293  i_sector.second.CorrectionY2->Write();
294  }
295  }
296  resultsFile->Close();
297 }
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().

unsigned int ApeEstimatorSummary::noSectors_
private

Definition at line 92 of file ApeEstimatorSummary.cc.

Referenced by analyze(), calculateApe(), and getTrackerSectorStructs().

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