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");
|