ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/CMSSW/Alignment/CommonAlignmentAlgorithm/plugins/SiStripBackplaneCalibration.cc
Revision: 1.1.2.7
Committed: Thu May 16 10:54:40 2013 UTC (11 years, 11 months ago) by jbehr
Content type: text/plain
Branch: branch53X_calibration
Changes since 1.1.2.6: +11 -3 lines
Log Message:
add configurability of optional parameter for reference run range.

File Contents

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