CommonAlignmentAlgorithm/plugins/SiStripLorentzAngleCalibration.cc

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