CommonAlignmentAlgorithm/plugins/SiStripBackplaneCalibration.cc

0001 /// \class SiStripBackplaneCalibration
0002 ///
0003 /// Calibration of back plane corrections for the strip tracker, integrated in the
0004 /// alignment algorithms. Note that not all algorithms support this...
0005 ///
0006 /// Usally use one instance for deco mode data since peak mode should give the
0007 //  reference - but the strip local reco foresees also calibration parameters
0008 //  for peak data like for the Lorentz angle.
0009 ///
0010 ///  \author    : Gero Flucke
0011 ///  date       : November 2012
0012 ///  $Revision: 1.1.2.13 $
0013 ///  $Date: 2013/05/31 08:37:12 $
0014 ///  (last update by $Author: flucke $)
0015
0016 #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
0017 #include "Alignment/CommonAlignmentAlgorithm/interface/TkModuleGroupSelector.h"
0018 // include 'locally':
0019 #include "SiStripReadoutModeEnums.h"
0020 #include "TreeStruct.h"
0021
0022 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0023 #include "CondFormats/DataRecord/interface/SiStripCondDataRecords.h"
0024 #include "CondFormats/SiStripObjects/interface/SiStripLatency.h"
0025 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
0026 // #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
0027 #include "CondFormats/SiStripObjects/interface/SiStripBackPlaneCorrection.h"
0028
0029 #include "DataFormats/DetId/interface/DetId.h"
0030 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0031
0032 #include "FWCore/Framework/interface/ESWatcher.h"
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "FWCore/ServiceRegistry/interface/Service.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036
0037 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0038 #include "MagneticField/Engine/interface/MagneticField.h"
0039 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0040 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0041 #include "FWCore/Framework/interface/ConsumesCollector.h"
0042
0043 #include "TTree.h"
0044 #include "TFile.h"
0045 #include "TString.h"
0046
0047 // #include <iostream>
0048 #include <vector>
0049 #include <map>
0050 #include <sstream>
0051 #include <cstdio>
0052 #include <functional>
0053
0054 class SiStripBackplaneCalibration : public IntegratedCalibrationBase {
0055 public:
0056   /// Constructor
0057   explicit SiStripBackplaneCalibration(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC);
0058
0059   /// Destructor
0060   ~SiStripBackplaneCalibration() override;
0061
0062   /// How many parameters does this calibration define?
0063   unsigned int numParameters() const override;
0064
0065   // /// Return all derivatives,
0066   // /// default implementation uses other derivatives(..) method,
0067   // /// but can be overwritten in derived class for efficiency.
0068   // virtual std::vector<double> derivatives(const TransientTrackingRecHit &hit,
0069   //                      const TrajectoryStateOnSurface &tsos,
0070   //                      const edm::EventSetup &setup,
0071   //                      const EventInfo &eventInfo) const;
0072
0073   /// Return non-zero derivatives for x- and y-measurements with their indices by reference.
0074   /// Return value is their number.
0075   unsigned int derivatives(std::vector<ValuesIndexPair> &outDerivInds,
0076                            const TransientTrackingRecHit &hit,
0077                            const TrajectoryStateOnSurface &tsos,
0078                            const edm::EventSetup &setup,
0079                            const EventInfo &eventInfo) const override;
0080
0081   /// Setting the determined parameter identified by index,
0082   /// returns false if out-of-bounds, true otherwise.
0083   bool setParameter(unsigned int index, double value) override;
0084
0085   /// Setting the determined parameter uncertainty identified by index,
0086   /// returns false if out-of-bounds, true otherwise.
0087   bool setParameterError(unsigned int index, double error) override;
0088
0089   /// Return current value of parameter identified by index.
0090   /// Returns 0. if index out-of-bounds.
0091   double getParameter(unsigned int index) const override;
0092
0093   /// Return current value of parameter identified by index.
0094   /// Returns 0. if index out-of-bounds or if errors undetermined.
0095   double getParameterError(unsigned int index) const override;
0096
0097   // /// Call at beginning of job:
0098   void beginOfJob(AlignableTracker *tracker, AlignableMuon *muon, AlignableExtras *extras) override;
0099
0100   /// Called at end of a the job of the AlignmentProducer.
0101   /// Write out determined parameters.
0102   void endOfJob() override;
0103
0104 private:
0105   /// If called the first time, fill 'siStripBackPlaneCorrInput_',
0106   /// later check that Backplane has not changed.
0107   bool checkBackPlaneCorrectionInput(const edm::EventSetup &setup, const EventInfo &eventInfo);
0108   /// Input BackPlaneCorrection values:
0109   /// - either from EventSetup of first call to derivatives(..)
0110   /// - or created from files of passed by configuration (i.e. from parallel processing)
0111   const SiStripBackPlaneCorrection *getBackPlaneCorrectionInput();
0112
0113   /// Determined parameter value for this detId (detId not treated => 0.)
0114   /// and the given run.
0115   double getParameterForDetId(unsigned int detId, edm::RunNumber_t run) const;
0116
0117   void writeTree(const SiStripBackPlaneCorrection &backPlaneCorr,
0118                  const std::map<unsigned int, TreeStruct> &treeInfo,
0119                  const char *treeName) const;
0120   SiStripBackPlaneCorrection *createFromTree(const char *fileName, const char *treeName) const;
0121
0122   const std::string readoutModeName_;
0123   int16_t readoutMode_;
0124   const bool saveToDB_;
0125   const std::string recordNameDBwrite_;
0126   const std::string outFileName_;
0127   const std::vector<std::string> mergeFileNames_;
0128
0129   edm::ESWatcher<SiStripBackPlaneCorrectionRcd> watchBackPlaneCorrRcd_;
0130
0131   SiStripBackPlaneCorrection *siStripBackPlaneCorrInput_;
0132
0133   std::vector<double> parameters_;
0134   std::vector<double> paramUncertainties_;
0135
0136   TkModuleGroupSelector *moduleGroupSelector_;
0137   const edm::ParameterSet moduleGroupSelCfg_;
0138   const edm::ESGetToken<SiStripLatency, SiStripLatencyRcd> latencyToken_;
0139   const edm::ESGetToken<SiStripLorentzAngle, SiStripLorentzAngleRcd> lorentzAngleToken_;
0140   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0141   const edm::ESGetToken<SiStripBackPlaneCorrection, SiStripBackPlaneCorrectionRcd> backPlaneCorrToken_;
0142 };
0143
0144 //======================================================================
0145 //======================================================================
0146 //======================================================================
0147
0148 SiStripBackplaneCalibration::SiStripBackplaneCalibration(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC)
0149     : IntegratedCalibrationBase(cfg),
0150       readoutModeName_(cfg.getParameter<std::string>("readoutMode")),
0151       saveToDB_(cfg.getParameter<bool>("saveToDB")),
0152       recordNameDBwrite_(cfg.getParameter<std::string>("recordNameDBwrite")),
0153       outFileName_(cfg.getParameter<std::string>("treeFile")),
0154       mergeFileNames_(cfg.getParameter<std::vector<std::string> >("mergeTreeFiles")),
0155       siStripBackPlaneCorrInput_(nullptr),
0156       moduleGroupSelector_(nullptr),
0157       moduleGroupSelCfg_(cfg.getParameter<edm::ParameterSet>("BackplaneModuleGroups")),
0158       latencyToken_(iC.esConsumes()),
0159       lorentzAngleToken_(iC.esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", readoutModeName_))),
0160       magFieldToken_(iC.esConsumes()),
0161       backPlaneCorrToken_(iC.esConsumes(edm::ESInputTag("", readoutModeName_))) {
0162   // SiStripLatency::singleReadOutMode() returns
0163   // 1: all in peak, 0: all in deco, -1: mixed state
0164   // (in principle one could treat even mixed state APV by APV...)
0165   if (readoutModeName_ == "peak") {
0166     readoutMode_ = kPeakMode;
0167   } else if (readoutModeName_ == "deconvolution") {
0168     readoutMode_ = kDeconvolutionMode;
0169   } else {
0170     throw cms::Exception("BadConfig") << "SiStripBackplaneCalibration:\n"
0171                                       << "Unknown mode '" << readoutModeName_
0172                                       << "', should be 'peak' or 'deconvolution' .\n";
0173   }
0174 }
0175
0176 //======================================================================
0177 SiStripBackplaneCalibration::~SiStripBackplaneCalibration() {
0178   delete moduleGroupSelector_;
0179   //  std::cout << "Destroy SiStripBackplaneCalibration named " << this->name() << std::endl;
0180   delete siStripBackPlaneCorrInput_;
0181 }
0182
0183 //======================================================================
0184 unsigned int SiStripBackplaneCalibration::numParameters() const { return parameters_.size(); }
0185
0186 //======================================================================
0187 unsigned int SiStripBackplaneCalibration::derivatives(std::vector<ValuesIndexPair> &outDerivInds,
0188                                                       const TransientTrackingRecHit &hit,
0189                                                       const TrajectoryStateOnSurface &tsos,
0190                                                       const edm::EventSetup &setup,
0191                                                       const EventInfo &eventInfo) const {
0192   // ugly const-cast:
0193   // But it is either only first initialisation or throwing an exception...
0194   const_cast<SiStripBackplaneCalibration *>(this)->checkBackPlaneCorrectionInput(setup, eventInfo);
0195
0196   outDerivInds.clear();
0197
0198   const SiStripLatency *latency = &setup.getData(latencyToken_);
0199   const int16_t mode = latency->singleReadOutMode();
0200
0201   if (mode == readoutMode_) {
0202     if (hit.det()) {  // otherwise 'constraint hit' or whatever
0203
0204       const int index =
0205           moduleGroupSelector_->getParameterIndexFromDetId(hit.det()->geographicalId(), eventInfo.eventId().run());
0206       if (index >= 0) {  // otherwise not treated
0207         const MagneticField *magneticField = &setup.getData(magFieldToken_);
0208         const GlobalVector bField(magneticField->inTesla(hit.det()->surface().position()));
0209         const LocalVector bFieldLocal(hit.det()->surface().toLocal(bField));
0210         //std::cout << "SiStripBackplaneCalibration derivatives " << readoutModeName_ << std::endl;
0211         const double dZ = hit.det()->surface().bounds().thickness();          // it's a float only...
0212         const double tanPsi = tsos.localParameters().mixedFormatVector()[1];  //float...
0213
0214         const SiStripLorentzAngle *lorentzAngleHandle = &setup.getData(lorentzAngleToken_);
0215         // Yes, mobility (= LA/By) stored in object called LA...
0216         const double mobility = lorentzAngleHandle->getLorentzAngle(hit.det()->geographicalId());
0217         // shift due to dead back plane has two parts:
0218         // 1) Lorentz Angle correction formula gets reduced thickness dz*(1-bp)
0219         // 2) 'Direct' effect is shift of effective module position in local z by bp*dz/2
0220         //   (see GF's presentation in alignment meeting 25.10.2012,
0221         //    https://indico.cern.ch/conferenceDisplay.py?confId=174266#2012-10-25)
0222         //        const double xDerivative = 0.5 * dZ * (mobility * bFieldLocal.y() - tanPsi);
0223         //FIXME: +tanPsi? At least that fits the sign of the dr/dw residual
0224         //       in KarimakiDerivatives...
0225         const double xDerivative = 0.5 * dZ * (mobility * bFieldLocal.y() + tanPsi);
0226         // std::cout << "derivative is " << xDerivative << " for index " << index
0227         //    << std::endl;
0228         //         if (index >= 10) {
0229         //           std::cout << "mobility * y-field " << mobility * bFieldLocal.y()
0230         //                     << ", \n tanPsi " << tanPsi /*<< ", dZ " << dZ */ << std::endl;
0231         //           std::cout << "|tanPsi| - |mobility * y-field| "
0232         //                     << fabs(tanPsi) - fabs(mobility * bFieldLocal.y())
0233         //                     << std::endl;
0234         //         }
0235         const Values derivs(xDerivative, 0.);  // yDerivative = 0.
0236         outDerivInds.push_back(ValuesIndexPair(derivs, index));
0237       }
0238     } else {
0239       edm::LogWarning("Alignment") << "@SUB=SiStripBackplaneCalibration::derivatives1"
0240                                    << "Hit without GeomDet, skip!";
0241     }
0242   } else if (mode != kDeconvolutionMode && mode != kPeakMode) {
0243     // warn only if unknown/mixed mode
0244     edm::LogWarning("Alignment") << "@SUB=SiStripBackplaneCalibration::derivatives2"
0245                                  << "Readout mode is " << mode << ", but looking for " << readoutMode_ << " ("
0246                                  << readoutModeName_ << ").";
0247   }
0248
0249   return outDerivInds.size();
0250 }
0251
0252 //======================================================================
0253 bool SiStripBackplaneCalibration::setParameter(unsigned int index, double value) {
0254   if (index >= parameters_.size()) {
0255     return false;
0256   } else {
0257     parameters_[index] = value;
0258     return true;
0259   }
0260 }
0261
0262 //======================================================================
0263 bool SiStripBackplaneCalibration::setParameterError(unsigned int index, double error) {
0264   if (index >= paramUncertainties_.size()) {
0265     return false;
0266   } else {
0267     paramUncertainties_[index] = error;
0268     return true;
0269   }
0270 }
0271
0272 //======================================================================
0273 double SiStripBackplaneCalibration::getParameter(unsigned int index) const {
0274   return (index >= parameters_.size() ? 0. : parameters_[index]);
0275 }
0276
0277 //======================================================================
0278 double SiStripBackplaneCalibration::getParameterError(unsigned int index) const {
0279   return (index >= paramUncertainties_.size() ? 0. : paramUncertainties_[index]);
0280 }
0281
0282 //======================================================================
0283 void SiStripBackplaneCalibration::beginOfJob(AlignableTracker *aliTracker,
0284                                              AlignableMuon * /*aliMuon*/,
0285                                              AlignableExtras * /*aliExtras*/) {
0286   //specify the sub-detectors for which the back plane correction is determined: all strips
0287   const std::vector<int> sdets = {SiStripDetId::TIB, SiStripDetId::TID, SiStripDetId::TOB, SiStripDetId::TEC};
0288
0289   moduleGroupSelector_ = new TkModuleGroupSelector(aliTracker, moduleGroupSelCfg_, sdets);
0290
0291   parameters_.resize(moduleGroupSelector_->getNumberOfParameters(), 0.);
0292   paramUncertainties_.resize(moduleGroupSelector_->getNumberOfParameters(), 0.);
0293
0294   edm::LogInfo("Alignment") << "@SUB=SiStripBackplaneCalibration"
0295                             << "Created with name " << this->name() << " for readout mode '" << readoutModeName_
0296                             << "',\n"
0297                             << this->numParameters() << " parameters to be determined."
0298                             << "\nsaveToDB = " << saveToDB_ << "\n outFileName = " << outFileName_
0299                             << "\n N(merge files) = " << mergeFileNames_.size()
0300                             << "\n number of IOVs = " << moduleGroupSelector_->numIovs();
0301
0302   if (!mergeFileNames_.empty()) {
0303     edm::LogInfo("Alignment") << "@SUB=SiStripBackplaneCalibration"
0304                               << "First file to merge: " << mergeFileNames_[0];
0305   }
0306 }
0307
0308 //======================================================================
0309 void SiStripBackplaneCalibration::endOfJob() {
0310   // loginfo output
0311   std::ostringstream out;
0312   out << "Parameter results for readout mode '" << readoutModeName_ << "'\n";
0313   for (unsigned int iPar = 0; iPar < parameters_.size(); ++iPar) {
0314     out << iPar << ": " << parameters_[iPar] << " +- " << paramUncertainties_[iPar] << "\n";
0315   }
0316   edm::LogInfo("Alignment") << "@SUB=SiStripBackplaneCalibration::endOfJob" << out.str();
0317
0318   std::map<unsigned int, TreeStruct> treeInfo;  // map of TreeStruct for each detId
0319
0320   // now write 'input' tree
0321   const SiStripBackPlaneCorrection *input = this->getBackPlaneCorrectionInput();  // never NULL
0322   const std::string treeName(this->name() + '_' + readoutModeName_ + '_');
0323   this->writeTree(*input, treeInfo, (treeName + "input").c_str());  // empty treeInfo for input...
0324
0325   if (input->getBackPlaneCorrections().empty()) {
0326     edm::LogError("Alignment") << "@SUB=SiStripBackplaneCalibration::endOfJob"
0327                                << "Input back plane correction map is empty ('" << readoutModeName_
0328                                << "' mode), skip writing output!";
0329     return;
0330   }
0331
0332   const unsigned int nonZeroParamsOrErrors =  // Any determined value?
0333       count_if(parameters_.begin(), parameters_.end(), [](auto c) { return c != 0.; }) +
0334       count_if(paramUncertainties_.begin(), paramUncertainties_.end(), [](auto c) { return c != 0.; });
0335
0336   for (unsigned int iIOV = 0; iIOV < moduleGroupSelector_->numIovs(); ++iIOV) {
0337     cond::Time_t firstRunOfIOV = moduleGroupSelector_->firstRunOfIOV(iIOV);
0338     SiStripBackPlaneCorrection output{};
0339     // Loop on map of values from input and add (possible) parameter results
0340     for (auto iterIdValue = input->getBackPlaneCorrections().begin();
0341          iterIdValue != input->getBackPlaneCorrections().end();
0342          ++iterIdValue) {
0343       // type of (*iterIdValue) is pair<unsigned int, float>
0344       const unsigned int detId = iterIdValue->first;  // key of map is DetId
0345       // Some code one could use to miscalibrate wrt input:
0346       // double param = 0.;
0347       // const DetId id(detId);
0348       // switch (id.subdetId()) {
0349       // case 3:
0350       //   // delta = 0.025;
0351       //   {
0352       //     TIBDetId tibId(detId);
0353       //     param = tibId.layer() * 0.01;
0354       //   }
0355       //   break;
0356       // case 5:
0357       //   //delta = 0.050;
0358       //   {
0359       //     TOBDetId tobId(detId);
0360       //     param = tobId.layer() * 0.015;
0361       //   }
0362       //   break;
0363       // case 4: // TID
0364       //   param = 0.03; break;
0365       // case 6: // TEC
0366       //   param = 0.05; break;
0367       //   break;
0368       // default:
0369       //   std::cout << "unknown subdet " << id.subdetId() << std::endl;
0370       // }
0371       const double param = this->getParameterForDetId(detId, firstRunOfIOV);
0372       // put result in output, i.e. sum of input and determined parameter:
0373       output.putBackPlaneCorrection(detId, iterIdValue->second + param);
0374       const int paramIndex = moduleGroupSelector_->getParameterIndexFromDetId(detId, firstRunOfIOV);
0375       treeInfo[detId] = TreeStruct(param, this->getParameterError(paramIndex), paramIndex);
0376     }
0377
0378     if (saveToDB_ || nonZeroParamsOrErrors != 0) {  // Skip writing mille jobs...
0379       this->writeTree(output, treeInfo, (treeName + Form("result_%lld", firstRunOfIOV)).c_str());
0380     }
0381
0382     if (saveToDB_) {  // If requested, write out to DB
0383       edm::Service<cond::service::PoolDBOutputService> dbService;
0384       if (dbService.isAvailable()) {
0385         dbService->writeOneIOV(output, firstRunOfIOV, recordNameDBwrite_);
0386         // no 'delete output;': writeOne(..) took over ownership
0387       } else {
0388         edm::LogError("BadConfig") << "@SUB=SiStripBackplaneCalibration::endOfJob"
0389                                    << "No PoolDBOutputService available, but saveToDB true!";
0390       }
0391     }
0392   }  // end loop on IOVs
0393 }
0394
0395 //======================================================================
0396 bool SiStripBackplaneCalibration::checkBackPlaneCorrectionInput(const edm::EventSetup &setup,
0397                                                                 const EventInfo &eventInfo) {
0398   const SiStripBackPlaneCorrection *backPlaneCorrHandle = nullptr;
0399   if (!siStripBackPlaneCorrInput_) {
0400     backPlaneCorrHandle = &setup.getData(backPlaneCorrToken_);
0401     siStripBackPlaneCorrInput_ = new SiStripBackPlaneCorrection(*backPlaneCorrHandle);
0402     // FIXME: Should we call 'watchBackPlaneCorrRcd_.check(setup)' as well?
0403     //        Otherwise could be that next check has to check via following 'else', though
0404     //        no new IOV has started... (to be checked)
0405   } else {
0406     if (watchBackPlaneCorrRcd_.check(setup)) {  // new IOV of input - but how to check peak vs deco?
0407       backPlaneCorrHandle = &setup.getData(backPlaneCorrToken_);
0408
0409       if (backPlaneCorrHandle->getBackPlaneCorrections()               // but only bad if non-identical values
0410           != siStripBackPlaneCorrInput_->getBackPlaneCorrections()) {  // (comparing maps)
0411         // Maps are containers sorted by key, but comparison problems may arise from
0412         // 'floating point comparison' problems (FIXME?)
0413         throw cms::Exception("BadInput") << "SiStripBackplaneCalibration::checkBackPlaneCorrectionInput:\n"
0414                                          << "Content of SiStripBackPlaneCorrection changed at run "
0415                                          << eventInfo.eventId().run() << ", but algorithm expects constant input!\n";
0416         return false;  // not reached...
0417       }
0418     }
0419   }
0420
0421   return true;
0422 }
0423
0424 //======================================================================
0425 const SiStripBackPlaneCorrection *SiStripBackplaneCalibration::getBackPlaneCorrectionInput() {
0426   // For parallel processing in Millepede II, create SiStripBackPlaneCorrection
0427   // from info stored in files of parallel jobs and check that they are identical.
0428   // If this job has run on events, still check that back plane corrections are
0429   // identical to the ones from mergeFileNames_.
0430   const std::string treeName(((this->name() + '_') += readoutModeName_) += "_input");
0431   for (auto iFile = mergeFileNames_.begin(); iFile != mergeFileNames_.end(); ++iFile) {
0432     SiStripBackPlaneCorrection *bpCorr = this->createFromTree(iFile->c_str(), treeName.c_str());
0433     // siStripBackPlaneCorrInput_ could be non-null from previous file of this loop
0434     // or from checkBackPlaneCorrectionInput(..) when running on data in this job as well
0435     if (!siStripBackPlaneCorrInput_ || siStripBackPlaneCorrInput_->getBackPlaneCorrections().empty()) {
0436       delete siStripBackPlaneCorrInput_;  // NULL or empty
0437       siStripBackPlaneCorrInput_ = bpCorr;
0438     } else {
0439       // about comparison of maps see comments in checkBackPlaneCorrectionInput
0440       if (bpCorr && !bpCorr->getBackPlaneCorrections().empty() &&  // single job might not have got events
0441           bpCorr->getBackPlaneCorrections() != siStripBackPlaneCorrInput_->getBackPlaneCorrections()) {
0442         // Throw exception instead of error?
0443         edm::LogError("NoInput") << "@SUB=SiStripBackplaneCalibration::getBackPlaneCorrectionInput"
0444                                  << "Different input values from tree " << treeName << " in file " << *iFile << ".";
0445       }
0446       delete bpCorr;
0447     }
0448   }
0449
0450   if (!siStripBackPlaneCorrInput_) {  // no files nor ran on events
0451     siStripBackPlaneCorrInput_ = new SiStripBackPlaneCorrection;
0452     edm::LogError("NoInput") << "@SUB=SiStripBackplaneCalibration::getBackPlaneCorrectionInput"
0453                              << "No input, create an empty one ('" << readoutModeName_ << "' mode)!";
0454   } else if (siStripBackPlaneCorrInput_->getBackPlaneCorrections().empty()) {
0455     edm::LogError("NoInput") << "@SUB=SiStripBackplaneCalibration::getBackPlaneCorrectionInput"
0456                              << "Empty result ('" << readoutModeName_ << "' mode)!";
0457   }
0458
0459   return siStripBackPlaneCorrInput_;
0460 }
0461
0462 //======================================================================
0463 double SiStripBackplaneCalibration::getParameterForDetId(unsigned int detId, edm::RunNumber_t run) const {
0464   const int index = moduleGroupSelector_->getParameterIndexFromDetId(detId, run);
0465
0466   return (index < 0 ? 0. : parameters_[index]);
0467 }
0468
0469 //======================================================================
0470 void SiStripBackplaneCalibration::writeTree(const SiStripBackPlaneCorrection &backPlaneCorrection,
0471                                             const std::map<unsigned int, TreeStruct> &treeInfo,
0472                                             const char *treeName) const {
0473   TFile *file = TFile::Open(outFileName_.c_str(), "UPDATE");
0474   if (!file) {
0475     edm::LogError("BadConfig") << "@SUB=SiStripBackplaneCalibration::writeTree"
0476                                << "Could not open file '" << outFileName_ << "'.";
0477     return;
0478   }
0479
0480   TTree *tree = new TTree(treeName, treeName);
0481   unsigned int id = 0;
0482   float value = 0.;
0483   TreeStruct treeStruct;
0484   tree->Branch("detId", &id, "detId/i");
0485   tree->Branch("value", &value, "value/F");
0486   tree->Branch("treeStruct", &treeStruct, TreeStruct::LeafList());
0487
0488   for (auto iterIdValue = backPlaneCorrection.getBackPlaneCorrections().begin();
0489        iterIdValue != backPlaneCorrection.getBackPlaneCorrections().end();
0490        ++iterIdValue) {
0491     // type of (*iterIdValue) is pair<unsigned int, float>
0492     id = iterIdValue->first;  // key of map is DetId
0493     value = iterIdValue->second;
0494     // type of (*treeStructIter) is pair<unsigned int, TreeStruct>
0495     auto treeStructIter = treeInfo.find(id);  // find info for this id
0496     if (treeStructIter != treeInfo.end()) {
0497       treeStruct = treeStructIter->second;  // info from input map
0498     } else {                                // if none found, fill at least parameter index (using 1st IOV...)
0499       const cond::Time_t run1of1stIov = moduleGroupSelector_->firstRunOfIOV(0);
0500       const int ind = moduleGroupSelector_->getParameterIndexFromDetId(id, run1of1stIov);
0501       treeStruct = TreeStruct(ind);
0502     }
0503     tree->Fill();
0504   }
0505   tree->Write();
0506   delete file;  // tree vanishes with the file...
0507 }
0508
0509 //======================================================================
0510 SiStripBackPlaneCorrection *SiStripBackplaneCalibration::createFromTree(const char *fileName,
0511                                                                         const char *treeName) const {
0512   // Check for file existence on your own to work around
0513   // https://hypernews.cern.ch/HyperNews/CMS/get/swDevelopment/2715.html:
0514   TFile *file = nullptr;
0515   FILE *testFile = fopen(fileName, "r");
0516   if (testFile) {
0517     fclose(testFile);
0518     file = TFile::Open(fileName, "READ");
0519   }  // else not existing, see error below
0520
0521   TTree *tree = nullptr;
0522   if (file)
0523     file->GetObject(treeName, tree);
0524
0525   SiStripBackPlaneCorrection *result = nullptr;
0526   if (tree) {
0527     unsigned int id = 0;
0528     float value = 0.;
0529     tree->SetBranchAddress("detId", &id);
0530     tree->SetBranchAddress("value", &value);
0531
0532     result = new SiStripBackPlaneCorrection;
0533     const Long64_t nEntries = tree->GetEntries();
0534     for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
0535       tree->GetEntry(iEntry);
0536       result->putBackPlaneCorrection(id, value);
0537     }
0538   } else {  // Warning only since could be parallel job on no events.
0539     edm::LogWarning("Alignment") << "@SUB=SiStripBackplaneCalibration::createFromTree"
0540                                  << "Could not get TTree '" << treeName << "' from file '" << fileName
0541                                  << (file ? "'." : "' (file does not exist).");
0542   }
0543
0544   delete file;  // tree will vanish with file
0545   return result;
0546 }
0547
0548 //======================================================================
0549 //======================================================================
0550 // Plugin definition
0551
0552 #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationPluginFactory.h"
0553
0554 DEFINE_EDM_PLUGIN(IntegratedCalibrationPluginFactory, SiStripBackplaneCalibration, "SiStripBackplaneCalibration");