ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/CMSSW/Alignment/CommonAlignmentAlgorithm/plugins/SiStripLorentzAngleCalibration.cc
Revision: 1.7
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.6: +128 -176 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.1 /// \class SiStripLorentzAngleCalibration
2     ///
3 flucke 1.7 /// Calibration of Lorentz angle for the strip tracker, integrated in the
4     /// alignment algorithms. Note that not all algorithms support this...
5 flucke 1.1 ///
6 flucke 1.7 /// Use one instance for peak and/or one instance for deco mode data.
7 flucke 1.1 ///
8     /// \author : Gero Flucke
9     /// date : August 2012
10 flucke 1.7 /// $Revision: 1.6.2.14 $
11     /// $Date: 2013/05/31 08:37:12 $
12 flucke 1.1 /// (last update by $Author: flucke $)
13    
14     #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
15 flucke 1.7 #include "Alignment/CommonAlignmentAlgorithm/interface/TkModuleGroupSelector.h"
16     // include 'locally':
17     #include "SiStripReadoutModeEnums.h"
18     #include "TreeStruct.h"
19 flucke 1.1
20     #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
21     #include "CondFormats/DataRecord/interface/SiStripCondDataRecords.h"
22     #include "CondFormats/SiStripObjects/interface/SiStripLatency.h"
23     #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
24 flucke 1.7 //#include "CalibTracker/Records/interface/SiStripDependentRecords.h"
25     #include "CondFormats/SiStripObjects/interface/SiStripBackPlaneCorrection.h"
26 flucke 1.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     #include "FWCore/ServiceRegistry/interface/Service.h"
34 flucke 1.7 #include "FWCore/ParameterSet/interface/ParameterSet.h"
35 flucke 1.1
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     #include "TTree.h"
42     #include "TFile.h"
43     #include "TString.h"
44    
45     // #include <iostream>
46 flucke 1.7 #include <boost/assign/list_of.hpp>
47 flucke 1.1 #include <vector>
48 flucke 1.6 #include <map>
49 flucke 1.1 #include <sstream>
50 flucke 1.3 #include <cstdio>
51 flucke 1.7 #include <functional>
52 flucke 1.1
53     class SiStripLorentzAngleCalibration : public IntegratedCalibrationBase
54     {
55     public:
56     /// Constructor
57     explicit SiStripLorentzAngleCalibration(const edm::ParameterSet &cfg);
58    
59     /// Destructor
60     virtual ~SiStripLorentzAngleCalibration();
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 flucke 1.7 /// returns false if out-of-bounds, true otherwise.
83 flucke 1.1 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 flucke 1.7 /// Returns 0. if index out-of-bounds.
91 flucke 1.1 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 flucke 1.7 virtual void beginOfJob(AlignableTracker *tracker,
99     AlignableMuon *muon,
100     AlignableExtras *extras);
101    
102 flucke 1.1
103     /// Called at end of a the job of the AlignmentProducer.
104     /// Write out determined parameters.
105     virtual void endOfJob();
106    
107     private:
108     /// If called the first time, fill 'siStripLorentzAngleInput_',
109     /// later check that LorentzAngle has not changed.
110     bool checkLorentzAngleInput(const edm::EventSetup &setup, const EventInfo &eventInfo);
111     /// Input LorentzAngle values:
112     /// - either from EventSetup of first call to derivatives(..)
113     /// - or created from files of passed by configuration (i.e. from parallel processing)
114     const SiStripLorentzAngle* getLorentzAnglesInput();
115 flucke 1.6 /// in non-peak mode the effective thickness is reduced...
116     double effectiveThickness(const GeomDet *det, int16_t mode, const edm::EventSetup &setup) const;
117 flucke 1.1
118     /// Determined parameter value for this detId (detId not treated => 0.)
119     /// and the given run.
120     double getParameterForDetId(unsigned int detId, edm::RunNumber_t run) const;
121    
122 flucke 1.7 void writeTree(const SiStripLorentzAngle *lorentzAngle,
123     const std::map<unsigned int,TreeStruct> &treeInfo, const char *treeName) const;
124 flucke 1.1 SiStripLorentzAngle* createFromTree(const char *fileName, const char *treeName) const;
125 flucke 1.7
126 flucke 1.1 const std::string readoutModeName_;
127     int16_t readoutMode_;
128     const bool saveToDB_;
129     const std::string recordNameDBwrite_;
130     const std::string outFileName_;
131     const std::vector<std::string> mergeFileNames_;
132    
133     edm::ESWatcher<SiStripLorentzAngleRcd> watchLorentzAngleRcd_;
134    
135     SiStripLorentzAngle *siStripLorentzAngleInput_;
136 flucke 1.6
137 flucke 1.1 std::vector<double> parameters_;
138     std::vector<double> paramUncertainties_;
139 flucke 1.7
140     TkModuleGroupSelector *moduleGroupSelector_;
141     const edm::ParameterSet moduleGroupSelCfg_;
142 flucke 1.1 };
143    
144     //======================================================================
145     //======================================================================
146     //======================================================================
147    
148     SiStripLorentzAngleCalibration::SiStripLorentzAngleCalibration(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 flucke 1.7 siStripLorentzAngleInput_(0),
156     moduleGroupSelector_(0),
157     moduleGroupSelCfg_(cfg.getParameter<edm::ParameterSet>("LorentzAngleModuleGroups"))
158 flucke 1.1 {
159 flucke 1.7
160 flucke 1.1 // 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 flucke 1.7 readoutMode_ = kPeakMode;
165 flucke 1.1 } else if (readoutModeName_ == "deconvolution") {
166 flucke 1.7 readoutMode_ = kDeconvolutionMode;
167 flucke 1.1 } else {
168     throw cms::Exception("BadConfig")
169     << "SiStripLorentzAngleCalibration:\n" << "Unknown mode '"
170     << readoutModeName_ << "', should be 'peak' or 'deconvolution' .\n";
171     }
172    
173     }
174    
175     //======================================================================
176     SiStripLorentzAngleCalibration::~SiStripLorentzAngleCalibration()
177     {
178 flucke 1.7 delete moduleGroupSelector_;
179 flucke 1.1 // std::cout << "Destroy SiStripLorentzAngleCalibration named " << this->name() << std::endl;
180     delete siStripLorentzAngleInput_;
181     }
182    
183     //======================================================================
184     unsigned int SiStripLorentzAngleCalibration::numParameters() const
185     {
186     return parameters_.size();
187     }
188    
189     //======================================================================
190     unsigned int
191     SiStripLorentzAngleCalibration::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<SiStripLorentzAngleCalibration*>(this)->checkLorentzAngleInput(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 flucke 1.7 if (mode == readoutMode_) {
207 flucke 1.1 if (hit.det()) { // otherwise 'constraint hit' or whatever
208    
209 flucke 1.7 const int index = moduleGroupSelector_->getParameterIndexFromDetId(hit.det()->geographicalId(),
210     eventInfo.eventId_.run());
211 flucke 1.1 if (index >= 0) { // otherwise not treated
212     edm::ESHandle<MagneticField> magneticField;
213     setup.get<IdealMagneticFieldRecord>().get(magneticField);
214     const GlobalVector bField(magneticField->inTesla(hit.det()->surface().position()));
215     const LocalVector bFieldLocal(hit.det()->surface().toLocal(bField));
216 flucke 1.6 //std::cout << "SiStripLorentzAngleCalibration derivatives " << readoutModeName_ << std::endl;
217     const double dZ = this->effectiveThickness(hit.det(), mode, setup);
218 flucke 1.5 // shift due to LA: dx = tan(LA) * dz/2 = mobility * B_y * dz/2,
219 flucke 1.7 // '-' since we have derivative of the residual r = hit - trk and mu is part of trk model
220     // (see GF's presentation in alignment meeting 25.10.2012,
221     // https://indico.cern.ch/conferenceDisplay.py?confId=174266#2012-10-25)
222     // Hmm! StripCPE::fillParams() defines, together with
223     // StripCPE::driftDirection(...):
224     // drift.x = -mobility * by * thickness (full drift from backside)
225     // So '-' already comes from that, not from mobility being part of
226     // track model...
227 flucke 1.5 const double xDerivative = bFieldLocal.y() * dZ * -0.5; // parameter is mobility!
228 flucke 1.1 if (xDerivative) { // If field is zero, this is zero: do not return it
229     const Values derivs(xDerivative, 0.); // yDerivative = 0.
230     outDerivInds.push_back(ValuesIndexPair(derivs, index));
231     }
232     }
233     } else {
234 flucke 1.7 edm::LogWarning("Alignment") << "@SUB=SiStripLorentzAngleCalibration::derivatives1"
235 flucke 1.1 << "Hit without GeomDet, skip!";
236     }
237 flucke 1.7 } else if (mode != kDeconvolutionMode && mode != kPeakMode) {
238     // warn only if unknown/mixed mode
239     edm::LogWarning("Alignment") << "@SUB=SiStripLorentzAngleCalibration::derivatives2"
240 flucke 1.1 << "Readout mode is " << mode << ", but looking for "
241     << readoutMode_ << " (" << readoutModeName_ << ").";
242     }
243    
244     return outDerivInds.size();
245     }
246    
247     //======================================================================
248     bool SiStripLorentzAngleCalibration::setParameter(unsigned int index, double value)
249     {
250     if (index >= parameters_.size()) {
251     return false;
252     } else {
253     parameters_[index] = value;
254     return true;
255     }
256     }
257    
258     //======================================================================
259     bool SiStripLorentzAngleCalibration::setParameterError(unsigned int index, double error)
260     {
261     if (index >= paramUncertainties_.size()) {
262     return false;
263     } else {
264     paramUncertainties_[index] = error;
265     return true;
266     }
267     }
268    
269     //======================================================================
270     double SiStripLorentzAngleCalibration::getParameter(unsigned int index) const
271     {
272     return (index >= parameters_.size() ? 0. : parameters_[index]);
273     }
274    
275     //======================================================================
276     double SiStripLorentzAngleCalibration::getParameterError(unsigned int index) const
277     {
278     return (index >= paramUncertainties_.size() ? 0. : paramUncertainties_[index]);
279     }
280    
281 flucke 1.7 //======================================================================
282     void SiStripLorentzAngleCalibration::beginOfJob(AlignableTracker *aliTracker,
283     AlignableMuon * /*aliMuon*/,
284     AlignableExtras * /*aliExtras*/)
285     {
286     //specify the sub-detectors for which the LA is determined
287     const std::vector<int> sdets = boost::assign::list_of(SiStripDetId::TIB)(SiStripDetId::TOB); //no TEC,TID
288     moduleGroupSelector_ = new TkModuleGroupSelector(aliTracker, moduleGroupSelCfg_, sdets);
289    
290     parameters_.resize(moduleGroupSelector_->getNumberOfParameters(), 0.);
291     paramUncertainties_.resize(moduleGroupSelector_->getNumberOfParameters(), 0.);
292    
293     edm::LogInfo("Alignment") << "@SUB=SiStripLorentzAngleCalibration" << "Created with name "
294     << this->name() << " for readout mode '" << readoutModeName_
295     << "',\n" << this->numParameters() << " parameters to be determined."
296     << "\nsaveToDB = " << saveToDB_
297     << "\n outFileName = " << outFileName_
298     << "\n N(merge files) = " << mergeFileNames_.size()
299     << "\n number of IOVs = " << moduleGroupSelector_->numIovs();
300    
301     if (mergeFileNames_.size()) {
302     edm::LogInfo("Alignment") << "@SUB=SiStripLorentzAngleCalibration"
303     << "First file to merge: " << mergeFileNames_[0];
304     }
305     }
306 flucke 1.1
307    
308     //======================================================================
309     void SiStripLorentzAngleCalibration::endOfJob()
310     {
311     // loginfo output
312     std::ostringstream out;
313     out << "Parameter results for readout mode '" << readoutModeName_ << "'\n";
314     for (unsigned int iPar = 0; iPar < parameters_.size(); ++iPar) {
315     out << iPar << ": " << parameters_[iPar] << " +- " << paramUncertainties_[iPar] << "\n";
316     }
317     edm::LogInfo("Alignment") << "@SUB=SiStripLorentzAngleCalibration::endOfJob" << out.str();
318    
319 flucke 1.7 std::map<unsigned int, TreeStruct> treeInfo; // map of TreeStruct for each detId
320    
321 flucke 1.1 // now write 'input' tree
322     const SiStripLorentzAngle *input = this->getLorentzAnglesInput(); // never NULL
323     const std::string treeName(this->name() + '_' + readoutModeName_ + '_');
324 flucke 1.7 this->writeTree(input, treeInfo, (treeName + "input").c_str()); // empty treeInfo for input...
325 flucke 1.1
326 flucke 1.3 if (input->getLorentzAngles().empty()) {
327     edm::LogError("Alignment") << "@SUB=SiStripLorentzAngleCalibration::endOfJob"
328     << "Input Lorentz angle map is empty ('"
329     << readoutModeName_ << "' mode), skip writing output!";
330     return;
331     }
332    
333 flucke 1.7 const unsigned int nonZeroParamsOrErrors = // Any determined value?
334     count_if (parameters_.begin(), parameters_.end(), std::bind2nd(std::not_equal_to<double>(),0.))
335     + count_if(paramUncertainties_.begin(), paramUncertainties_.end(),
336     std::bind2nd(std::not_equal_to<double>(), 0.));
337    
338     for (unsigned int iIOV = 0; iIOV < moduleGroupSelector_->numIovs(); ++iIOV) {
339     cond::Time_t firstRunOfIOV = moduleGroupSelector_->firstRunOfIOV(iIOV);
340 flucke 1.1 SiStripLorentzAngle *output = new SiStripLorentzAngle;
341     // Loop on map of values from input and add (possible) parameter results
342     for (auto iterIdValue = input->getLorentzAngles().begin();
343     iterIdValue != input->getLorentzAngles().end(); ++iterIdValue) {
344     // type of (*iterIdValue) is pair<unsigned int, float>
345     const unsigned int detId = iterIdValue->first; // key of map is DetId
346 flucke 1.7 // Some code one could use to miscalibrate wrt input:
347     // double param = 0.;
348     // const DetId id(detId);
349     // if (id.subdetId() == 3) { // TIB
350     // param = (readoutMode_ == kPeakMode ? -0.003 : -0.002);
351     // } else if (id.subdetId() == 5) { // TOB
352     // param = (readoutMode_ == kPeakMode ? 0.005 : 0.004);
353     // }
354     const double param = this->getParameterForDetId(detId, firstRunOfIOV);
355     // put result in output, i.e. sum of input and determined parameter:
356     output->putLorentzAngle(detId, iterIdValue->second + param);
357     const int paramIndex = moduleGroupSelector_->getParameterIndexFromDetId(detId,firstRunOfIOV);
358     treeInfo[detId] = TreeStruct(param, this->getParameterError(paramIndex), paramIndex);
359 flucke 1.1 }
360    
361 flucke 1.7 if (saveToDB_ || nonZeroParamsOrErrors != 0) { // Skip writing mille jobs...
362     this->writeTree(output, treeInfo, (treeName + Form("result_%lld", firstRunOfIOV)).c_str());
363     }
364 flucke 1.1
365     if (saveToDB_) { // If requested, write out to DB
366     edm::Service<cond::service::PoolDBOutputService> dbService;
367     if (dbService.isAvailable()) {
368     dbService->writeOne(output, firstRunOfIOV, recordNameDBwrite_.c_str());
369     // no 'delete output;': writeOne(..) took over ownership
370     } else {
371     delete output;
372     edm::LogError("BadConfig") << "@SUB=SiStripLorentzAngleCalibration::endOfJob"
373     << "No PoolDBOutputService available, but saveToDB true!";
374     }
375     } else {
376     delete output;
377     }
378     } // end loop on IOVs
379     }
380    
381     //======================================================================
382     bool SiStripLorentzAngleCalibration::checkLorentzAngleInput(const edm::EventSetup &setup,
383     const EventInfo &eventInfo)
384     {
385     edm::ESHandle<SiStripLorentzAngle> lorentzAngleHandle;
386     if (!siStripLorentzAngleInput_) {
387     setup.get<SiStripLorentzAngleRcd>().get(readoutModeName_, lorentzAngleHandle);
388     siStripLorentzAngleInput_ = new SiStripLorentzAngle(*lorentzAngleHandle);
389 flucke 1.7 // FIXME: Should we call 'watchLorentzAngleRcd_.check(setup)' as well?
390     // Otherwise could be that next check has to check via following 'else', though
391     // no new IOV has started... (to be checked)
392 flucke 1.1 } else {
393     if (watchLorentzAngleRcd_.check(setup)) { // new IOV of input - but how to check peak vs deco?
394     setup.get<SiStripLorentzAngleRcd>().get(readoutModeName_, lorentzAngleHandle);
395     if (lorentzAngleHandle->getLorentzAngles() // but only bad if non-identical values
396     != siStripLorentzAngleInput_->getLorentzAngles()) { // (comparing maps)
397 flucke 1.7 // Maps are containers sorted by key, but comparison problems may arise from
398     // 'floating point comparison' problems (FIXME?)
399 flucke 1.1 throw cms::Exception("BadInput")
400     << "SiStripLorentzAngleCalibration::checkLorentzAngleInput:\n"
401     << "Content of SiStripLorentzAngle changed at run " << eventInfo.eventId_.run()
402     << ", but algorithm expects constant input!\n";
403     return false; // not reached...
404     }
405     }
406     }
407    
408     return true;
409     }
410    
411     //======================================================================
412 flucke 1.6 double SiStripLorentzAngleCalibration::effectiveThickness(const GeomDet *det,
413     int16_t mode,
414     const edm::EventSetup &setup) const
415     {
416     if (!det) return 0.;
417     double dZ = det->surface().bounds().thickness(); // it is a float only...
418 flucke 1.7 const SiStripDetId id(det->geographicalId());
419     edm::ESHandle<SiStripBackPlaneCorrection> backPlaneHandle;
420     // FIXME: which one? DepRcd->get(handle) or Rcd->get(readoutModeName_, handle)??
421     // setup.get<SiStripBackPlaneCorrectionDepRcd>().get(backPlaneHandle); // get correct mode
422     setup.get<SiStripBackPlaneCorrectionRcd>().get(readoutModeName_, backPlaneHandle);
423     const double bpCor = backPlaneHandle->getBackPlaneCorrection(id); // it's a float...
424     // std::cout << "bpCor " << bpCor << " in subdet " << id.subdetId() << std::endl;
425     dZ *= (1. - bpCor);
426    
427 flucke 1.6 return dZ;
428     }
429    
430     //======================================================================
431 flucke 1.1 const SiStripLorentzAngle* SiStripLorentzAngleCalibration::getLorentzAnglesInput()
432     {
433     // For parallel processing in Millepede II, create SiStripLorentzAngle
434     // from info stored in files of parallel jobs and check that they are identical.
435 flucke 1.7 // If this job has run on events, still check that LA is identical to the ones
436 flucke 1.1 // from mergeFileNames_.
437     const std::string treeName(((this->name() + '_') += readoutModeName_) += "_input");
438     for (auto iFile = mergeFileNames_.begin(); iFile != mergeFileNames_.end(); ++iFile) {
439     SiStripLorentzAngle* la = this->createFromTree(iFile->c_str(), treeName.c_str());
440     // siStripLorentzAngleInput_ could be non-null from previous file of this loop
441     // or from checkLorentzAngleInput(..) when running on data in this job as well
442     if (!siStripLorentzAngleInput_ || siStripLorentzAngleInput_->getLorentzAngles().empty()) {
443     delete siStripLorentzAngleInput_; // NULL or empty
444     siStripLorentzAngleInput_ = la;
445     } else {
446     // FIXME: about comparison of maps see comments in checkLorentzAngleInput
447 flucke 1.4 if (la && !la->getLorentzAngles().empty() && // single job might not have got events
448 flucke 1.1 la->getLorentzAngles() != siStripLorentzAngleInput_->getLorentzAngles()) {
449     // Throw exception instead of error?
450     edm::LogError("NoInput") << "@SUB=SiStripLorentzAngleCalibration::getLorentzAnglesInput"
451     << "Different input values from tree " << treeName
452     << " in file " << *iFile << ".";
453    
454     }
455     delete la;
456     }
457     }
458    
459     if (!siStripLorentzAngleInput_) { // no files nor ran on events
460     siStripLorentzAngleInput_ = new SiStripLorentzAngle;
461     edm::LogError("NoInput") << "@SUB=SiStripLorentzAngleCalibration::getLorentzAnglesInput"
462     << "No input, create an empty one ('" << readoutModeName_ << "' mode)!";
463     } else if (siStripLorentzAngleInput_->getLorentzAngles().empty()) {
464     edm::LogError("NoInput") << "@SUB=SiStripLorentzAngleCalibration::getLorentzAnglesInput"
465     << "Empty result ('" << readoutModeName_ << "' mode)!";
466     }
467    
468     return siStripLorentzAngleInput_;
469     }
470    
471     //======================================================================
472     double SiStripLorentzAngleCalibration::getParameterForDetId(unsigned int detId,
473     edm::RunNumber_t run) const
474     {
475 flucke 1.7 const int index = moduleGroupSelector_->getParameterIndexFromDetId(detId, run);
476 flucke 1.1
477     return (index < 0 ? 0. : parameters_[index]);
478     }
479    
480     //======================================================================
481     void SiStripLorentzAngleCalibration::writeTree(const SiStripLorentzAngle *lorentzAngle,
482 flucke 1.7 const std::map<unsigned int, TreeStruct> &treeInfo,
483 flucke 1.1 const char *treeName) const
484     {
485     if (!lorentzAngle) return;
486    
487     TFile* file = TFile::Open(outFileName_.c_str(), "UPDATE");
488     if (!file) {
489     edm::LogError("BadConfig") << "@SUB=SiStripLorentzAngleCalibration::writeTree"
490     << "Could not open file '" << outFileName_ << "'.";
491     return;
492     }
493    
494     TTree *tree = new TTree(treeName, treeName);
495     unsigned int id = 0;
496 flucke 1.2 float value = 0.;
497 flucke 1.7 TreeStruct treeStruct;
498 flucke 1.1 tree->Branch("detId", &id, "detId/i");
499 flucke 1.2 tree->Branch("value", &value, "value/F");
500 flucke 1.7 tree->Branch("treeStruct", &treeStruct, TreeStruct::LeafList());
501 flucke 1.1
502     for (auto iterIdValue = lorentzAngle->getLorentzAngles().begin();
503     iterIdValue != lorentzAngle->getLorentzAngles().end(); ++iterIdValue) {
504     // type of (*iterIdValue) is pair<unsigned int, float>
505     id = iterIdValue->first; // key of map is DetId
506     value = iterIdValue->second;
507 flucke 1.7 // type of (*treeStructIter) is pair<unsigned int, TreeStruct>
508     auto treeStructIter = treeInfo.find(id); // find info for this id
509     if (treeStructIter != treeInfo.end()) {
510     treeStruct = treeStructIter->second; // info from input map
511     } else { // if none found, fill at least parameter index (using 1st IOV...)
512     const cond::Time_t run1of1stIov = moduleGroupSelector_->firstRunOfIOV(0);
513     const int ind = moduleGroupSelector_->getParameterIndexFromDetId(id, run1of1stIov);
514     treeStruct = TreeStruct(ind);
515     }
516 flucke 1.1 tree->Fill();
517     }
518     tree->Write();
519 flucke 1.7 delete file; // tree vanishes with the file...
520 flucke 1.1 }
521    
522     //======================================================================
523     SiStripLorentzAngle*
524     SiStripLorentzAngleCalibration::createFromTree(const char *fileName, const char *treeName) const
525     {
526 flucke 1.3 // 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 flucke 1.1 TTree *tree = 0;
536     if (file) file->GetObject(treeName, tree);
537    
538     SiStripLorentzAngle *result = 0;
539     if (tree) {
540     unsigned int id = 0;
541 flucke 1.2 float value = 0.;
542 flucke 1.1 tree->SetBranchAddress("detId", &id);
543     tree->SetBranchAddress("value", &value);
544    
545     result = new SiStripLorentzAngle;
546     const Long64_t nEntries = tree->GetEntries();
547     for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
548     tree->GetEntry(iEntry);
549 flucke 1.2 result->putLorentzAngle(id, value);
550 flucke 1.1 }
551 flucke 1.3 } else { // Warning only since could be parallel job on no events.
552     edm::LogWarning("Alignment") << "@SUB=SiStripLorentzAngleCalibration::createFromTree"
553     << "Could not get TTree '" << treeName << "' from file '"
554     << fileName << (file ? "'." : "' (file does not exist).");
555 flucke 1.1 }
556    
557     delete file; // tree will vanish with file
558     return result;
559     }
560    
561    
562     //======================================================================
563     //======================================================================
564     // Plugin definition
565    
566     #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationPluginFactory.h"
567    
568     DEFINE_EDM_PLUGIN(IntegratedCalibrationPluginFactory,
569     SiStripLorentzAngleCalibration, "SiStripLorentzAngleCalibration");