ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/CMSSW/Alignment/CommonAlignmentAlgorithm/plugins/SiStripBackplaneCalibration.cc
Revision: 1.2
Committed: Fri May 31 12:13:40 2013 UTC (11 years, 11 months ago) by flucke
Content type: text/plain
Branch: MAIN
CVS Tags: CMSSW_6_2_0, CMSSW_6_2_0_pre8, V04-00-14, V04-00-13, V04-00-09-53X-calib-V01, V04-00-12, HEAD
Changes since 1.1: +585 -0 lines
Log Message:
merging calibration development from branch 'branch53X_calibration' (tag
'br53-00-09') to HEAD:
- add backplane calibration
- configurable granularity in time (i.e. run) and space for lorentz angle
  calibrations
- more info in diagnostics tree

File Contents

# User Rev Content
1 flucke 1.2 /// \class SiStripBackplaneCalibration
2     ///
3     /// Calibration of back plane corrections for the strip tracker, integrated in the
4     /// alignment algorithms. Note that not all algorithms support this...
5     ///
6     /// Usally use one instance for deco mode data since peak mode should give the
7     // reference - but the strip local reco foresees also calibration parameters
8     // for peak data like for the Lorentz angle.
9     ///
10     /// \author : Gero Flucke
11     /// date : November 2012
12     /// $Revision: 1.1.2.13 $
13     /// $Date: 2013/05/31 08:37:12 $
14     /// (last update by $Author: flucke $)
15    
16     #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
17     #include "Alignment/CommonAlignmentAlgorithm/interface/TkModuleGroupSelector.h"
18     // include 'locally':
19     #include "SiStripReadoutModeEnums.h"
20     #include "TreeStruct.h"
21    
22     #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
23     #include "CondFormats/DataRecord/interface/SiStripCondDataRecords.h"
24     #include "CondFormats/SiStripObjects/interface/SiStripLatency.h"
25     #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
26     // #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
27     #include "CondFormats/SiStripObjects/interface/SiStripBackPlaneCorrection.h"
28    
29     #include "DataFormats/DetId/interface/DetId.h"
30     #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
31    
32     #include "FWCore/Framework/interface/ESHandle.h"
33     #include "FWCore/Framework/interface/ESWatcher.h"
34     #include "FWCore/MessageLogger/interface/MessageLogger.h"
35     #include "FWCore/ServiceRegistry/interface/Service.h"
36     #include "FWCore/ParameterSet/interface/ParameterSet.h"
37    
38     #include "Geometry/CommonDetUnit/interface/GeomDet.h"
39     #include "MagneticField/Engine/interface/MagneticField.h"
40     #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
41     #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
42    
43     #include "TTree.h"
44     #include "TFile.h"
45     #include "TString.h"
46    
47     // #include <iostream>
48     #include <boost/assign/list_of.hpp>
49     #include <vector>
50     #include <map>
51     #include <sstream>
52     #include <cstdio>
53     #include <functional>
54    
55     class SiStripBackplaneCalibration : public IntegratedCalibrationBase
56     {
57     public:
58     /// Constructor
59     explicit SiStripBackplaneCalibration(const edm::ParameterSet &cfg);
60    
61     /// Destructor
62     virtual ~SiStripBackplaneCalibration();
63    
64     /// How many parameters does this calibration define?
65     virtual unsigned int numParameters() const;
66    
67     // /// Return all derivatives,
68     // /// default implementation uses other derivatives(..) method,
69     // /// but can be overwritten in derived class for efficiency.
70     // virtual std::vector<double> derivatives(const TransientTrackingRecHit &hit,
71     // const TrajectoryStateOnSurface &tsos,
72     // const edm::EventSetup &setup,
73     // const EventInfo &eventInfo) const;
74    
75     /// Return non-zero derivatives for x- and y-measurements with their indices by reference.
76     /// Return value is their number.
77     virtual unsigned int derivatives(std::vector<ValuesIndexPair> &outDerivInds,
78     const TransientTrackingRecHit &hit,
79     const TrajectoryStateOnSurface &tsos,
80     const edm::EventSetup &setup,
81     const EventInfo &eventInfo) const;
82    
83     /// Setting the determined parameter identified by index,
84     /// returns false if out-of-bounds, true otherwise.
85     virtual bool setParameter(unsigned int index, double value);
86    
87     /// Setting the determined parameter uncertainty identified by index,
88     /// returns false if out-of-bounds, true otherwise.
89     virtual bool setParameterError(unsigned int index, double error);
90    
91     /// Return current value of parameter identified by index.
92     /// Returns 0. if index out-of-bounds.
93     virtual double getParameter(unsigned int index) const;
94    
95     /// Return current value of parameter identified by index.
96     /// Returns 0. if index out-of-bounds or if errors undetermined.
97     virtual double getParameterError(unsigned int index) const;
98    
99     // /// Call at beginning of job:
100     virtual void beginOfJob(AlignableTracker *tracker,
101     AlignableMuon *muon,
102     AlignableExtras *extras);
103    
104     /// Called at end of a the job of the AlignmentProducer.
105     /// Write out determined parameters.
106     virtual void endOfJob();
107    
108     private:
109     /// If called the first time, fill 'siStripBackPlaneCorrInput_',
110     /// later check that Backplane has not changed.
111     bool checkBackPlaneCorrectionInput(const edm::EventSetup &setup, const EventInfo &eventInfo);
112     /// Input BackPlaneCorrection values:
113     /// - either from EventSetup of first call to derivatives(..)
114     /// - or created from files of passed by configuration (i.e. from parallel processing)
115     const SiStripBackPlaneCorrection* getBackPlaneCorrectionInput();
116    
117     /// Determined parameter value for this detId (detId not treated => 0.)
118     /// and the given run.
119     double getParameterForDetId(unsigned int detId, edm::RunNumber_t run) const;
120    
121     void writeTree(const SiStripBackPlaneCorrection *backPlaneCorr,
122     const std::map<unsigned int,TreeStruct> &treeInfo, const char *treeName) const;
123     SiStripBackPlaneCorrection* createFromTree(const char *fileName, const char *treeName) const;
124    
125     const std::string readoutModeName_;
126     int16_t readoutMode_;
127     const bool saveToDB_;
128     const std::string recordNameDBwrite_;
129     const std::string outFileName_;
130     const std::vector<std::string> mergeFileNames_;
131    
132     edm::ESWatcher<SiStripBackPlaneCorrectionRcd> watchBackPlaneCorrRcd_;
133    
134     SiStripBackPlaneCorrection *siStripBackPlaneCorrInput_;
135    
136     std::vector<double> parameters_;
137     std::vector<double> paramUncertainties_;
138    
139     TkModuleGroupSelector *moduleGroupSelector_;
140     const edm::ParameterSet moduleGroupSelCfg_;
141    
142     };
143    
144     //======================================================================
145     //======================================================================
146     //======================================================================
147    
148     SiStripBackplaneCalibration::SiStripBackplaneCalibration(const edm::ParameterSet &cfg)
149     : IntegratedCalibrationBase(cfg),
150     readoutModeName_(cfg.getParameter<std::string>("readoutMode")),
151     saveToDB_(cfg.getParameter<bool>("saveToDB")),
152     recordNameDBwrite_(cfg.getParameter<std::string>("recordNameDBwrite")),
153     outFileName_(cfg.getParameter<std::string>("treeFile")),
154     mergeFileNames_(cfg.getParameter<std::vector<std::string> >("mergeTreeFiles")),
155     siStripBackPlaneCorrInput_(0),
156     moduleGroupSelector_(0),
157     moduleGroupSelCfg_(cfg.getParameter<edm::ParameterSet>("BackplaneModuleGroups"))
158     {
159    
160     // SiStripLatency::singleReadOutMode() returns
161     // 1: all in peak, 0: all in deco, -1: mixed state
162     // (in principle one could treat even mixed state APV by APV...)
163     if (readoutModeName_ == "peak") {
164     readoutMode_ = kPeakMode;
165     } else if (readoutModeName_ == "deconvolution") {
166     readoutMode_ = kDeconvolutionMode;
167     } else {
168     throw cms::Exception("BadConfig")
169     << "SiStripBackplaneCalibration:\n" << "Unknown mode '"
170     << readoutModeName_ << "', should be 'peak' or 'deconvolution' .\n";
171     }
172    
173     }
174    
175     //======================================================================
176     SiStripBackplaneCalibration::~SiStripBackplaneCalibration()
177     {
178     delete moduleGroupSelector_;
179     // std::cout << "Destroy SiStripBackplaneCalibration named " << this->name() << std::endl;
180     delete siStripBackPlaneCorrInput_;
181     }
182    
183     //======================================================================
184     unsigned int SiStripBackplaneCalibration::numParameters() const
185     {
186     return parameters_.size();
187     }
188    
189     //======================================================================
190     unsigned int
191     SiStripBackplaneCalibration::derivatives(std::vector<ValuesIndexPair> &outDerivInds,
192     const TransientTrackingRecHit &hit,
193     const TrajectoryStateOnSurface &tsos,
194     const edm::EventSetup &setup,
195     const EventInfo &eventInfo) const
196     {
197     // ugly const-cast:
198     // But it is either only first initialisation or throwing an exception...
199     const_cast<SiStripBackplaneCalibration*>(this)->checkBackPlaneCorrectionInput(setup, eventInfo);
200    
201     outDerivInds.clear();
202    
203     edm::ESHandle<SiStripLatency> latency;
204     setup.get<SiStripLatencyRcd>().get(latency);
205     const int16_t mode = latency->singleReadOutMode();
206    
207     if(mode == readoutMode_) {
208     if (hit.det()) { // otherwise 'constraint hit' or whatever
209    
210     const int index = moduleGroupSelector_->getParameterIndexFromDetId(hit.det()->geographicalId(),
211     eventInfo.eventId_.run());
212     if (index >= 0) { // otherwise not treated
213     edm::ESHandle<MagneticField> magneticField;
214     setup.get<IdealMagneticFieldRecord>().get(magneticField);
215     const GlobalVector bField(magneticField->inTesla(hit.det()->surface().position()));
216     const LocalVector bFieldLocal(hit.det()->surface().toLocal(bField));
217     //std::cout << "SiStripBackplaneCalibration derivatives " << readoutModeName_ << std::endl;
218     const double dZ = hit.det()->surface().bounds().thickness(); // it's a float only...
219     const double tanPsi = tsos.localParameters().mixedFormatVector()[1]; //float...
220    
221     edm::ESHandle<SiStripLorentzAngle> lorentzAngleHandle;
222     setup.get<SiStripLorentzAngleRcd>().get(readoutModeName_, lorentzAngleHandle);
223     // Yes, mobility (= LA/By) stored in object called LA...
224     const double mobility = lorentzAngleHandle->getLorentzAngle(hit.det()->geographicalId());
225     // shift due to dead back plane has two parts:
226     // 1) Lorentz Angle correction formula gets reduced thickness dz*(1-bp)
227     // 2) 'Direct' effect is shift of effective module position in local z by bp*dz/2
228     // (see GF's presentation in alignment meeting 25.10.2012,
229     // https://indico.cern.ch/conferenceDisplay.py?confId=174266#2012-10-25)
230     // const double xDerivative = 0.5 * dZ * (mobility * bFieldLocal.y() - tanPsi);
231     //FIXME: +tanPsi? At least that fits the sign of the dr/dw residual
232     // in KarimakiDerivatives...
233     const double xDerivative = 0.5 * dZ * (mobility * bFieldLocal.y() + tanPsi);
234     // std::cout << "derivative is " << xDerivative << " for index " << index
235     // << std::endl;
236     // if (index >= 10) {
237     // std::cout << "mobility * y-field " << mobility * bFieldLocal.y()
238     // << ", \n tanPsi " << tanPsi /*<< ", dZ " << dZ */ << std::endl;
239     // std::cout << "|tanPsi| - |mobility * y-field| "
240     // << fabs(tanPsi) - fabs(mobility * bFieldLocal.y())
241     // << std::endl;
242     // }
243     const Values derivs(xDerivative, 0.); // yDerivative = 0.
244     outDerivInds.push_back(ValuesIndexPair(derivs, index));
245     }
246     } else {
247     edm::LogWarning("Alignment") << "@SUB=SiStripBackplaneCalibration::derivatives1"
248     << "Hit without GeomDet, skip!";
249     }
250     } else if (mode != kDeconvolutionMode && mode != kPeakMode) {
251     // warn only if unknown/mixed mode
252     edm::LogWarning("Alignment") << "@SUB=SiStripBackplaneCalibration::derivatives2"
253     << "Readout mode is " << mode << ", but looking for "
254     << readoutMode_ << " (" << readoutModeName_ << ").";
255     }
256    
257     return outDerivInds.size();
258     }
259    
260     //======================================================================
261     bool SiStripBackplaneCalibration::setParameter(unsigned int index, double value)
262     {
263     if (index >= parameters_.size()) {
264     return false;
265     } else {
266     parameters_[index] = value;
267     return true;
268     }
269     }
270    
271     //======================================================================
272     bool SiStripBackplaneCalibration::setParameterError(unsigned int index, double error)
273     {
274     if (index >= paramUncertainties_.size()) {
275     return false;
276     } else {
277     paramUncertainties_[index] = error;
278     return true;
279     }
280     }
281    
282     //======================================================================
283     double SiStripBackplaneCalibration::getParameter(unsigned int index) const
284     {
285     return (index >= parameters_.size() ? 0. : parameters_[index]);
286     }
287    
288     //======================================================================
289     double SiStripBackplaneCalibration::getParameterError(unsigned int index) const
290     {
291     return (index >= paramUncertainties_.size() ? 0. : paramUncertainties_[index]);
292     }
293    
294     //======================================================================
295     void SiStripBackplaneCalibration::beginOfJob(AlignableTracker *aliTracker,
296     AlignableMuon * /*aliMuon*/,
297     AlignableExtras */*aliExtras*/)
298     {
299     //specify the sub-detectors for which the back plane correction is determined: all strips
300     const std::vector<int> sdets = boost::assign::list_of(SiStripDetId::TIB)(SiStripDetId::TID)
301     (SiStripDetId::TOB)(SiStripDetId::TEC);
302    
303     moduleGroupSelector_ = new TkModuleGroupSelector(aliTracker, moduleGroupSelCfg_, sdets);
304    
305     parameters_.resize(moduleGroupSelector_->getNumberOfParameters(), 0.);
306     paramUncertainties_.resize(moduleGroupSelector_->getNumberOfParameters(), 0.);
307    
308     edm::LogInfo("Alignment") << "@SUB=SiStripBackplaneCalibration" << "Created with name "
309     << this->name() << " for readout mode '" << readoutModeName_
310     << "',\n" << this->numParameters() << " parameters to be determined."
311     << "\nsaveToDB = " << saveToDB_
312     << "\n outFileName = " << outFileName_
313     << "\n N(merge files) = " << mergeFileNames_.size()
314     << "\n number of IOVs = " << moduleGroupSelector_->numIovs();
315    
316     if (mergeFileNames_.size()) {
317     edm::LogInfo("Alignment") << "@SUB=SiStripBackplaneCalibration"
318     << "First file to merge: " << mergeFileNames_[0];
319     }
320    
321     }
322    
323    
324     //======================================================================
325     void SiStripBackplaneCalibration::endOfJob()
326     {
327     // loginfo output
328     std::ostringstream out;
329     out << "Parameter results for readout mode '" << readoutModeName_ << "'\n";
330     for (unsigned int iPar = 0; iPar < parameters_.size(); ++iPar) {
331     out << iPar << ": " << parameters_[iPar] << " +- " << paramUncertainties_[iPar] << "\n";
332     }
333     edm::LogInfo("Alignment") << "@SUB=SiStripBackplaneCalibration::endOfJob" << out.str();
334    
335     std::map<unsigned int, TreeStruct> treeInfo; // map of TreeStruct for each detId
336    
337     // now write 'input' tree
338     const SiStripBackPlaneCorrection *input = this->getBackPlaneCorrectionInput(); // never NULL
339     const std::string treeName(this->name() + '_' + readoutModeName_ + '_');
340     this->writeTree(input, treeInfo, (treeName + "input").c_str()); // empty treeInfo for input...
341    
342     if (input->getBackPlaneCorrections().empty()) {
343     edm::LogError("Alignment") << "@SUB=SiStripBackplaneCalibration::endOfJob"
344     << "Input back plane correction map is empty ('"
345     << readoutModeName_ << "' mode), skip writing output!";
346     return;
347     }
348    
349     const unsigned int nonZeroParamsOrErrors = // Any determined value?
350     count_if (parameters_.begin(), parameters_.end(), std::bind2nd(std::not_equal_to<double>(),0.))
351     + count_if(paramUncertainties_.begin(), paramUncertainties_.end(),
352     std::bind2nd(std::not_equal_to<double>(), 0.));
353    
354     for (unsigned int iIOV = 0; iIOV < moduleGroupSelector_->numIovs(); ++iIOV) {
355     cond::Time_t firstRunOfIOV = moduleGroupSelector_->firstRunOfIOV(iIOV);
356     SiStripBackPlaneCorrection *output = new SiStripBackPlaneCorrection;
357     // Loop on map of values from input and add (possible) parameter results
358     for (auto iterIdValue = input->getBackPlaneCorrections().begin();
359     iterIdValue != input->getBackPlaneCorrections().end(); ++iterIdValue) {
360     // type of (*iterIdValue) is pair<unsigned int, float>
361     const unsigned int detId = iterIdValue->first; // key of map is DetId
362     // Some code one could use to miscalibrate wrt input:
363     // double param = 0.;
364     // const DetId id(detId);
365     // switch (id.subdetId()) {
366     // case 3:
367     // // delta = 0.025;
368     // {
369     // TIBDetId tibId(detId);
370     // param = tibId.layer() * 0.01;
371     // }
372     // break;
373     // case 5:
374     // //delta = 0.050;
375     // {
376     // TOBDetId tobId(detId);
377     // param = tobId.layer() * 0.015;
378     // }
379     // break;
380     // case 4: // TID
381     // param = 0.03; break;
382     // case 6: // TEC
383     // param = 0.05; break;
384     // break;
385     // default:
386     // std::cout << "unknown subdet " << id.subdetId() << std::endl;
387     // }
388     const double param = this->getParameterForDetId(detId, firstRunOfIOV);
389     // put result in output, i.e. sum of input and determined parameter:
390     output->putBackPlaneCorrection(detId, iterIdValue->second + param);
391     const int paramIndex = moduleGroupSelector_->getParameterIndexFromDetId(detId,firstRunOfIOV);
392     treeInfo[detId] = TreeStruct(param, this->getParameterError(paramIndex), paramIndex);
393     }
394    
395     if (saveToDB_ || nonZeroParamsOrErrors != 0) { // Skip writing mille jobs...
396     this->writeTree(output, treeInfo, (treeName + Form("result_%lld", firstRunOfIOV)).c_str());
397     }
398    
399     if (saveToDB_) { // If requested, write out to DB
400     edm::Service<cond::service::PoolDBOutputService> dbService;
401     if (dbService.isAvailable()) {
402     dbService->writeOne(output, firstRunOfIOV, recordNameDBwrite_.c_str());
403     // no 'delete output;': writeOne(..) took over ownership
404     } else {
405     delete output;
406     edm::LogError("BadConfig") << "@SUB=SiStripBackplaneCalibration::endOfJob"
407     << "No PoolDBOutputService available, but saveToDB true!";
408     }
409     } else {
410     delete output;
411     }
412     } // end loop on IOVs
413     }
414    
415     //======================================================================
416     bool SiStripBackplaneCalibration::checkBackPlaneCorrectionInput(const edm::EventSetup &setup,
417     const EventInfo &eventInfo)
418     {
419     edm::ESHandle<SiStripBackPlaneCorrection> backPlaneCorrHandle;
420     if (!siStripBackPlaneCorrInput_) {
421     setup.get<SiStripBackPlaneCorrectionRcd>().get(readoutModeName_, backPlaneCorrHandle);
422     siStripBackPlaneCorrInput_ = new SiStripBackPlaneCorrection(*backPlaneCorrHandle);
423     // FIXME: Should we call 'watchBackPlaneCorrRcd_.check(setup)' as well?
424     // Otherwise could be that next check has to check via following 'else', though
425     // no new IOV has started... (to be checked)
426     } else {
427     if (watchBackPlaneCorrRcd_.check(setup)) { // new IOV of input - but how to check peak vs deco?
428     setup.get<SiStripBackPlaneCorrectionRcd>().get(readoutModeName_, backPlaneCorrHandle);
429     if (backPlaneCorrHandle->getBackPlaneCorrections() // but only bad if non-identical values
430     != siStripBackPlaneCorrInput_->getBackPlaneCorrections()) { // (comparing maps)
431     // Maps are containers sorted by key, but comparison problems may arise from
432     // 'floating point comparison' problems (FIXME?)
433     throw cms::Exception("BadInput")
434     << "SiStripBackplaneCalibration::checkBackPlaneCorrectionInput:\n"
435     << "Content of SiStripBackPlaneCorrection changed at run " << eventInfo.eventId_.run()
436     << ", but algorithm expects constant input!\n";
437     return false; // not reached...
438     }
439     }
440     }
441    
442     return true;
443     }
444    
445     //======================================================================
446     const SiStripBackPlaneCorrection* SiStripBackplaneCalibration::getBackPlaneCorrectionInput()
447     {
448     // For parallel processing in Millepede II, create SiStripBackPlaneCorrection
449     // from info stored in files of parallel jobs and check that they are identical.
450     // If this job has run on events, still check that back plane corrections are
451     // identical to the ones from mergeFileNames_.
452     const std::string treeName(((this->name() + '_') += readoutModeName_) += "_input");
453     for (auto iFile = mergeFileNames_.begin(); iFile != mergeFileNames_.end(); ++iFile) {
454     SiStripBackPlaneCorrection* bpCorr = this->createFromTree(iFile->c_str(), treeName.c_str());
455     // siStripBackPlaneCorrInput_ could be non-null from previous file of this loop
456     // or from checkBackPlaneCorrectionInput(..) when running on data in this job as well
457     if (!siStripBackPlaneCorrInput_ || siStripBackPlaneCorrInput_->getBackPlaneCorrections().empty()) {
458     delete siStripBackPlaneCorrInput_; // NULL or empty
459     siStripBackPlaneCorrInput_ = bpCorr;
460     } else {
461     // about comparison of maps see comments in checkBackPlaneCorrectionInput
462     if (bpCorr && !bpCorr->getBackPlaneCorrections().empty() && // single job might not have got events
463     bpCorr->getBackPlaneCorrections() != siStripBackPlaneCorrInput_->getBackPlaneCorrections()) {
464     // Throw exception instead of error?
465     edm::LogError("NoInput") << "@SUB=SiStripBackplaneCalibration::getBackPlaneCorrectionInput"
466     << "Different input values from tree " << treeName
467     << " in file " << *iFile << ".";
468    
469     }
470     delete bpCorr;
471     }
472     }
473    
474     if (!siStripBackPlaneCorrInput_) { // no files nor ran on events
475     siStripBackPlaneCorrInput_ = new SiStripBackPlaneCorrection;
476     edm::LogError("NoInput") << "@SUB=SiStripBackplaneCalibration::getBackPlaneCorrectionInput"
477     << "No input, create an empty one ('" << readoutModeName_ << "' mode)!";
478     } else if (siStripBackPlaneCorrInput_->getBackPlaneCorrections().empty()) {
479     edm::LogError("NoInput") << "@SUB=SiStripBackplaneCalibration::getBackPlaneCorrectionInput"
480     << "Empty result ('" << readoutModeName_ << "' mode)!";
481     }
482    
483     return siStripBackPlaneCorrInput_;
484     }
485    
486     //======================================================================
487     double SiStripBackplaneCalibration::getParameterForDetId(unsigned int detId,
488     edm::RunNumber_t run) const
489     {
490     const int index = moduleGroupSelector_->getParameterIndexFromDetId(detId, run);
491    
492     return (index < 0 ? 0. : parameters_[index]);
493     }
494    
495     //======================================================================
496     void SiStripBackplaneCalibration::writeTree(const SiStripBackPlaneCorrection *backPlaneCorrection,
497     const std::map<unsigned int, TreeStruct> &treeInfo,
498     const char *treeName) const
499     {
500     if (!backPlaneCorrection) return;
501    
502     TFile* file = TFile::Open(outFileName_.c_str(), "UPDATE");
503     if (!file) {
504     edm::LogError("BadConfig") << "@SUB=SiStripBackplaneCalibration::writeTree"
505     << "Could not open file '" << outFileName_ << "'.";
506     return;
507     }
508    
509     TTree *tree = new TTree(treeName, treeName);
510     unsigned int id = 0;
511     float value = 0.;
512     TreeStruct treeStruct;
513     tree->Branch("detId", &id, "detId/i");
514     tree->Branch("value", &value, "value/F");
515     tree->Branch("treeStruct", &treeStruct, TreeStruct::LeafList());
516    
517     for (auto iterIdValue = backPlaneCorrection->getBackPlaneCorrections().begin();
518     iterIdValue != backPlaneCorrection->getBackPlaneCorrections().end(); ++iterIdValue) {
519     // type of (*iterIdValue) is pair<unsigned int, float>
520     id = iterIdValue->first; // key of map is DetId
521     value = iterIdValue->second;
522     // type of (*treeStructIter) is pair<unsigned int, TreeStruct>
523     auto treeStructIter = treeInfo.find(id); // find info for this id
524     if (treeStructIter != treeInfo.end()) {
525     treeStruct = treeStructIter->second; // info from input map
526     } else { // if none found, fill at least parameter index (using 1st IOV...)
527     const cond::Time_t run1of1stIov = moduleGroupSelector_->firstRunOfIOV(0);
528     const int ind = moduleGroupSelector_->getParameterIndexFromDetId(id, run1of1stIov);
529     treeStruct = TreeStruct(ind);
530     }
531     tree->Fill();
532     }
533     tree->Write();
534     delete file; // tree vanishes with the file...
535    
536     }
537    
538     //======================================================================
539     SiStripBackPlaneCorrection*
540     SiStripBackplaneCalibration::createFromTree(const char *fileName, const char *treeName) const
541     {
542     // Check for file existence on your own to work around
543     // https://hypernews.cern.ch/HyperNews/CMS/get/swDevelopment/2715.html:
544     TFile* file = 0;
545     FILE* testFile = fopen(fileName,"r");
546     if (testFile) {
547     fclose(testFile);
548     file = TFile::Open(fileName, "READ");
549     } // else not existing, see error below
550    
551     TTree *tree = 0;
552     if (file) file->GetObject(treeName, tree);
553    
554     SiStripBackPlaneCorrection *result = 0;
555     if (tree) {
556     unsigned int id = 0;
557     float value = 0.;
558     tree->SetBranchAddress("detId", &id);
559     tree->SetBranchAddress("value", &value);
560    
561     result = new SiStripBackPlaneCorrection;
562     const Long64_t nEntries = tree->GetEntries();
563     for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
564     tree->GetEntry(iEntry);
565     result->putBackPlaneCorrection(id, value);
566     }
567     } else { // Warning only since could be parallel job on no events.
568     edm::LogWarning("Alignment") << "@SUB=SiStripBackplaneCalibration::createFromTree"
569     << "Could not get TTree '" << treeName << "' from file '"
570     << fileName << (file ? "'." : "' (file does not exist).");
571     }
572    
573     delete file; // tree will vanish with file
574     return result;
575     }
576    
577    
578     //======================================================================
579     //======================================================================
580     // Plugin definition
581    
582     #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationPluginFactory.h"
583    
584     DEFINE_EDM_PLUGIN(IntegratedCalibrationPluginFactory,
585     SiStripBackplaneCalibration, "SiStripBackplaneCalibration");