SiStripLorentzAngle/plugins/MeasureLA.cc

0001 #include "CalibTracker/SiStripLorentzAngle/plugins/MeasureLA.h"
0002 #include "CalibTracker/SiStripLorentzAngle/interface/LA_Filler_Fitter.h"
0003 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0004 #include "CalibTracker/StandaloneTrackerTopology/interface/StandaloneTrackerTopology.h"
0005
0006 #include <TChain.h>
0007 #include <TFile.h>
0008 #include <regex>
0009 #include <fstream>
0010
0011 namespace sistrip {
0012
0013   void MeasureLA::store_methods_and_granularity(const edm::VParameterSet& vpset) {
0014     for (auto const& p : vpset) {
0015       methods |= p.getParameter<int32_t>("Method");
0016       byModule = byModule || p.getParameter<int32_t>("Granularity");
0017       byLayer = byLayer || !p.getParameter<int32_t>("Granularity");
0018     }
0019   }
0020
0021   MeasureLA::MeasureLA(const edm::ParameterSet& conf)
0022       : inputFiles(conf.getParameter<std::vector<std::string> >("InputFiles")),
0023         inFileLocation(conf.getParameter<std::string>("InFileLocation")),
0024         fp_(conf.getParameter<edm::FileInPath>("SiStripDetInfo")),
0025         reports(conf.getParameter<edm::VParameterSet>("Reports")),
0026         measurementPreferences(conf.getParameter<edm::VParameterSet>("MeasurementPreferences")),
0027         calibrations(conf.getParameter<edm::VParameterSet>("Calibrations")),
0028         methods(0),
0029         byModule(false),
0030         byLayer(false),
0031         localybin(conf.getUntrackedParameter<double>("LocalYBin", 0.0)),
0032         stripsperbin(conf.getUntrackedParameter<unsigned>("StripsPerBin", 0)),
0033         maxEvents(conf.getUntrackedParameter<unsigned>("MaxEvents", 0)),
0034         tTopo_(StandaloneTrackerTopology::fromTrackerParametersXMLFile(
0035             conf.getParameter<edm::FileInPath>("TrackerParameters").fullPath())) {
0036     store_methods_and_granularity(reports);
0037     store_methods_and_granularity(measurementPreferences);
0038     store_calibrations();
0039
0040     TChain* const chain = new TChain("la_data");
0041     for (auto const& file : inputFiles)
0042       chain->Add((file + inFileLocation).c_str());
0043
0044     LA_Filler_Fitter laff(methods, byLayer, byModule, localybin, stripsperbin, maxEvents, &tTopo_);
0045     laff.fill(chain, book);
0046     laff.fit(book);
0047     summarize_module_muH_byLayer(laff);
0048     process_reports();
0049
0050     setWhatProduced(this, &MeasureLA::produce);
0051   }
0052
0053   std::unique_ptr<SiStripLorentzAngle> MeasureLA::produce(const SiStripLorentzAngleRcd&) {
0054     auto lorentzAngle = std::make_unique<SiStripLorentzAngle>();
0055     /*
0056   std::map<uint32_t,LA_Filler_Fitter::Result>
0057     module_results = LA_Filler_Fitter::module_results(book, LA_Filler_Fitter::SQRTVAR);
0058
0059   BOOST_FOREACH(const uint32_t& detid, SiStripDetInfoFileReader::read(fp_.fullPath()).getAllDetIds()) {
0060     float la = module_results[detid].measure / module_results[detid].field ;
0061     lorentzAngle->putLorentzAngle( detid, la );
0062   }
0063   */
0064     return lorentzAngle;
0065   }
0066
0067   void MeasureLA::summarize_module_muH_byLayer(const LA_Filler_Fitter& laff) {
0068     for (int m = LA_Filler_Fitter::FIRST_METHOD; m <= LA_Filler_Fitter::LAST_METHOD; m <<= 1) {
0069       const LA_Filler_Fitter::Method method = (LA_Filler_Fitter::Method)m;
0070       for (auto& result : LA_Filler_Fitter::module_results(book, method)) {
0071         calibrate(calibration_key(result.first, method), result.second);
0072         std::string label =
0073             laff.layerLabel(result.first) + granularity(MODULESUMMARY) + LA_Filler_Fitter::method(method);
0074         label = std::regex_replace(label, std::regex("layer"), "");
0075
0076         const double mu_H = -result.second.calMeasured.first / result.second.field;
0077         const double sigma_mu_H = result.second.calMeasured.second / result.second.field;
0078         const double weight = pow(1. / sigma_mu_H, 2);
0079
0080         book.fill(mu_H, label, 150, -0.05, 0.1, weight);
0081       }
0082       for (Book::iterator it = book.begin(".*" + granularity(MODULESUMMARY) + ".*"); it != book.end(); ++it) {
0083         if (it->second->GetEntries())
0084           it->second->Fit("gaus", "LLQ");
0085       }
0086     }
0087   }
0088
0089   void MeasureLA::process_reports() const {
0090     for (auto const& p : reports) {
0091       const GRANULARITY gran = (GRANULARITY)p.getParameter<int32_t>("Granularity");
0092       const std::string name = p.getParameter<std::string>("ReportName");
0093       const LA_Filler_Fitter::Method method = (LA_Filler_Fitter::Method)p.getParameter<int32_t>("Method");
0094
0095       write_report_plots(name, method, gran);
0096       switch (gran) {
0097         case LAYER:
0098           write_report_text(name, method, LA_Filler_Fitter::layer_results(book, method));
0099           break;
0100         case MODULE:
0101           write_report_text(name, method, LA_Filler_Fitter::module_results(book, method));
0102           break;
0103         case MODULESUMMARY:
0104           write_report_text_ms(name, method);
0105           break;
0106       }
0107     }
0108
0109     {
0110       TFile widthsFile("widths.root", "RECREATE");
0111       for (Book::const_iterator it = book.begin(".*_width"); it != book.end(); it++)
0112         if (it->second)
0113           it->second->Write();
0114       widthsFile.Close();
0115     }
0116   }
0117
0118   void MeasureLA::write_report_plots(std::string name, LA_Filler_Fitter::Method method, GRANULARITY gran) const {
0119     TFile file((name + ".root").c_str(), "RECREATE");
0120     const std::string key = ".*" + granularity(gran) + ".*(" + LA_Filler_Fitter::method(method) + "|" +
0121                             LA_Filler_Fitter::method(method, false) + ".*)";
0122     for (Book::const_iterator hist = book.begin(key); hist != book.end(); ++hist)
0123       if (hist->second)
0124         hist->second->Write();
0125     file.Close();
0126   }
0127
0128   template <class T>
0129   void MeasureLA::write_report_text(std::string name,
0130                                     const LA_Filler_Fitter::Method& _method,
0131                                     const std::map<T, LA_Filler_Fitter::Result>& _results) const {
0132     LA_Filler_Fitter::Method method = _method;
0133     std::map<T, LA_Filler_Fitter::Result> results = _results;
0134     std::fstream file((name + ".dat").c_str(), std::ios::out);
0135     for (auto& result : results) {
0136       calibrate(calibration_key(result.first, method), result.second);
0137       file << result.first << "\t" << result.second << std::endl;
0138     }
0139     file.close();
0140   }
0141
0142   void MeasureLA::write_report_text_ms(std::string name, LA_Filler_Fitter::Method method) const {
0143     std::fstream file((name + ".dat").c_str(), std::ios::out);
0144     const std::string key = ".*" + granularity(MODULESUMMARY) + LA_Filler_Fitter::method(method);
0145     for (Book::const_iterator it = book.begin(key); it != book.end(); ++it) {
0146       const TF1* const f = it->second->GetFunction("gaus");
0147       if (f) {
0148         file << it->first << "\t" << f->GetParameter(1) << "\t" << f->GetParError(1) << "\t" << f->GetParameter(2)
0149              << "\t" << f->GetParError(2) << std::endl;
0150       }
0151     }
0152     file.close();
0153   }
0154
0155   void MeasureLA::store_calibrations() {
0156     for (auto const& p : calibrations) {
0157       LA_Filler_Fitter::Method method = (LA_Filler_Fitter::Method)p.getParameter<int32_t>("Method");
0158       std::vector<double> slopes(p.getParameter<std::vector<double> >("Slopes"));
0159       assert(slopes.size() == 14);
0160       std::vector<double> offsets(p.getParameter<std::vector<double> >("Offsets"));
0161       assert(offsets.size() == 14);
0162       std::vector<double> pulls(p.getParameter<std::vector<double> >("Pulls"));
0163       assert(pulls.size() == 14);
0164
0165       for (unsigned i = 0; i < 14; i++) {
0166         const std::pair<unsigned, LA_Filler_Fitter::Method> key(i, method);
0167         offset[key] = offsets[i];
0168         slope[key] = slopes[i];
0169         error_scaling[key] = pulls[i];
0170       }
0171     }
0172   }
0173
0174   inline void MeasureLA::calibrate(const std::pair<unsigned, LA_Filler_Fitter::Method> key,
0175                                    LA_Filler_Fitter::Result& result) const {
0176     result.calMeasured = std::make_pair<float, float>(
0177         (result.measured.first - offset.find(key)->second) / slope.find(key)->second,
0178         result.measured.second * error_scaling.find(key)->second / slope.find(key)->second);
0179   }
0180
0181   std::pair<uint32_t, LA_Filler_Fitter::Method> MeasureLA::calibration_key(
0182       const std::string layer, const LA_Filler_Fitter::Method method) const {
0183     std::regex format(".*(T[IO]B)_layer(\\d)([as]).*");
0184     const bool isTIB = "TIB" == std::regex_replace(layer, format, "\\1");
0185     const bool stereo = "s" == std::regex_replace(layer, format, "\\3");
0186     const unsigned layerNum = std::stoul(std::regex_replace(layer, format, "\\2"));
0187     return std::make_pair(LA_Filler_Fitter::layer_index(isTIB, stereo, layerNum), method);
0188   }
0189
0190   std::pair<uint32_t, LA_Filler_Fitter::Method> MeasureLA::calibration_key(
0191       const uint32_t detid, const LA_Filler_Fitter::Method method) const {
0192     const bool isTIB = SiStripDetId(detid).subDetector() == SiStripDetId::TIB;
0193     const bool stereo = isTIB ? tTopo_.tibStereo(detid) : tTopo_.tobStereo(detid);
0194     const unsigned layer = isTIB ? tTopo_.tibLayer(detid) : tTopo_.tobStereo(detid);
0195
0196     return std::make_pair(LA_Filler_Fitter::layer_index(isTIB, stereo, layer), method);
0197   }
0198
0199 }  // namespace sistrip