1 |
/// \class SiStripLorentzAngleCalibration
|
2 |
///
|
3 |
/// Calibration of Lorentz angle for the strip tracker, integrated in the
|
4 |
/// alignment algorithms. Note that not all algorithms support this...
|
5 |
///
|
6 |
/// Use one instance for peak and/or one instance for deco mode data.
|
7 |
///
|
8 |
/// \author : Gero Flucke
|
9 |
/// date : August 2012
|
10 |
/// $Revision: 1.6.2.14 $
|
11 |
/// $Date: 2013/05/31 08:37:12 $
|
12 |
/// (last update by $Author: flucke $)
|
13 |
|
14 |
#include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
|
15 |
#include "Alignment/CommonAlignmentAlgorithm/interface/TkModuleGroupSelector.h"
|
16 |
// include 'locally':
|
17 |
#include "SiStripReadoutModeEnums.h"
|
18 |
#include "TreeStruct.h"
|
19 |
|
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 |
//#include "CalibTracker/Records/interface/SiStripDependentRecords.h"
|
25 |
#include "CondFormats/SiStripObjects/interface/SiStripBackPlaneCorrection.h"
|
26 |
|
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 |
#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 |
#include "TTree.h"
|
42 |
#include "TFile.h"
|
43 |
#include "TString.h"
|
44 |
|
45 |
// #include <iostream>
|
46 |
#include <boost/assign/list_of.hpp>
|
47 |
#include <vector>
|
48 |
#include <map>
|
49 |
#include <sstream>
|
50 |
#include <cstdio>
|
51 |
#include <functional>
|
52 |
|
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 |
/// 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 |
virtual void beginOfJob(AlignableTracker *tracker,
|
99 |
AlignableMuon *muon,
|
100 |
AlignableExtras *extras);
|
101 |
|
102 |
|
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 |
/// in non-peak mode the effective thickness is reduced...
|
116 |
double effectiveThickness(const GeomDet *det, int16_t mode, const edm::EventSetup &setup) const;
|
117 |
|
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 |
void writeTree(const SiStripLorentzAngle *lorentzAngle,
|
123 |
const std::map<unsigned int,TreeStruct> &treeInfo, const char *treeName) const;
|
124 |
SiStripLorentzAngle* createFromTree(const char *fileName, const char *treeName) const;
|
125 |
|
126 |
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 |
|
137 |
std::vector<double> parameters_;
|
138 |
std::vector<double> paramUncertainties_;
|
139 |
|
140 |
TkModuleGroupSelector *moduleGroupSelector_;
|
141 |
const edm::ParameterSet moduleGroupSelCfg_;
|
142 |
};
|
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 |
siStripLorentzAngleInput_(0),
|
156 |
moduleGroupSelector_(0),
|
157 |
moduleGroupSelCfg_(cfg.getParameter<edm::ParameterSet>("LorentzAngleModuleGroups"))
|
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 |
<< "SiStripLorentzAngleCalibration:\n" << "Unknown mode '"
|
170 |
<< readoutModeName_ << "', should be 'peak' or 'deconvolution' .\n";
|
171 |
}
|
172 |
|
173 |
}
|
174 |
|
175 |
//======================================================================
|
176 |
SiStripLorentzAngleCalibration::~SiStripLorentzAngleCalibration()
|
177 |
{
|
178 |
delete moduleGroupSelector_;
|
179 |
// 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 |
if (mode == readoutMode_) {
|
207 |
if (hit.det()) { // otherwise 'constraint hit' or whatever
|
208 |
|
209 |
const int index = moduleGroupSelector_->getParameterIndexFromDetId(hit.det()->geographicalId(),
|
210 |
eventInfo.eventId_.run());
|
211 |
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 |
//std::cout << "SiStripLorentzAngleCalibration derivatives " << readoutModeName_ << std::endl;
|
217 |
const double dZ = this->effectiveThickness(hit.det(), mode, setup);
|
218 |
// shift due to LA: dx = tan(LA) * dz/2 = mobility * B_y * dz/2,
|
219 |
// '-' 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 |
const double xDerivative = bFieldLocal.y() * dZ * -0.5; // parameter is mobility!
|
228 |
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 |
edm::LogWarning("Alignment") << "@SUB=SiStripLorentzAngleCalibration::derivatives1"
|
235 |
<< "Hit without GeomDet, skip!";
|
236 |
}
|
237 |
} else if (mode != kDeconvolutionMode && mode != kPeakMode) {
|
238 |
// warn only if unknown/mixed mode
|
239 |
edm::LogWarning("Alignment") << "@SUB=SiStripLorentzAngleCalibration::derivatives2"
|
240 |
<< "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 |
//======================================================================
|
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 |
|
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 |
std::map<unsigned int, TreeStruct> treeInfo; // map of TreeStruct for each detId
|
320 |
|
321 |
// now write 'input' tree
|
322 |
const SiStripLorentzAngle *input = this->getLorentzAnglesInput(); // never NULL
|
323 |
const std::string treeName(this->name() + '_' + readoutModeName_ + '_');
|
324 |
this->writeTree(input, treeInfo, (treeName + "input").c_str()); // empty treeInfo for input...
|
325 |
|
326 |
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 |
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 |
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 |
// 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 |
}
|
360 |
|
361 |
if (saveToDB_ || nonZeroParamsOrErrors != 0) { // Skip writing mille jobs...
|
362 |
this->writeTree(output, treeInfo, (treeName + Form("result_%lld", firstRunOfIOV)).c_str());
|
363 |
}
|
364 |
|
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 |
// 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 |
} 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 |
// Maps are containers sorted by key, but comparison problems may arise from
|
398 |
// 'floating point comparison' problems (FIXME?)
|
399 |
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 |
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 |
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 |
return dZ;
|
428 |
}
|
429 |
|
430 |
//======================================================================
|
431 |
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 |
// If this job has run on events, still check that LA is identical to the ones
|
436 |
// 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 |
if (la && !la->getLorentzAngles().empty() && // single job might not have got events
|
448 |
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 |
const int index = moduleGroupSelector_->getParameterIndexFromDetId(detId, run);
|
476 |
|
477 |
return (index < 0 ? 0. : parameters_[index]);
|
478 |
}
|
479 |
|
480 |
//======================================================================
|
481 |
void SiStripLorentzAngleCalibration::writeTree(const SiStripLorentzAngle *lorentzAngle,
|
482 |
const std::map<unsigned int, TreeStruct> &treeInfo,
|
483 |
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 |
float value = 0.;
|
497 |
TreeStruct treeStruct;
|
498 |
tree->Branch("detId", &id, "detId/i");
|
499 |
tree->Branch("value", &value, "value/F");
|
500 |
tree->Branch("treeStruct", &treeStruct, TreeStruct::LeafList());
|
501 |
|
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 |
// 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 |
tree->Fill();
|
517 |
}
|
518 |
tree->Write();
|
519 |
delete file; // tree vanishes with the file...
|
520 |
}
|
521 |
|
522 |
//======================================================================
|
523 |
SiStripLorentzAngle*
|
524 |
SiStripLorentzAngleCalibration::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 |
SiStripLorentzAngle *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 SiStripLorentzAngle;
|
546 |
const Long64_t nEntries = tree->GetEntries();
|
547 |
for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
|
548 |
tree->GetEntry(iEntry);
|
549 |
result->putLorentzAngle(id, value);
|
550 |
}
|
551 |
} 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 |
}
|
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");
|