11 #include <ecl/dbobjects/ECLDatabaseImporter.h>
19 #include <boost/filesystem.hpp>
22 #include <ecl/dbobjects/ECLDigitEnergyConstants.h>
23 #include <ecl/dbobjects/ECLDigitTimeConstants.h>
24 #include <ecl/modules/eclShowerShape/ECLShowerShapeModule.h>
25 #include <ecl/dbobjects/ECLShowerShapeSecondMomentCorrection.h>
26 #include <ecl/dbobjects/ECLShowerCorrectorLeakageCorrection.h>
27 #include <ecl/dbobjects/ECLShowerEnergyCorrectionTemporary.h>
28 #include <ecl/dbobjects/ECLTrackClusterMatchingParameterizations.h>
29 #include <ecl/dbobjects/ECLTrackClusterMatchingThresholds.h>
30 #include <ecl/dataobjects/ECLShower.h>
33 #include <framework/database/IntervalOfValidity.h>
34 #include <framework/database/Database.h>
35 #include <framework/database/DBImportArray.h>
36 #include <framework/database/DBImportObjPtr.h>
41 #include <TClonesArray.h>
49 ECLDatabaseImporter::ECLDatabaseImporter(vector<string> inputFileNames,
const std::string& name)
52 for (
auto& inputFileName : inputFileNames)
53 m_inputFileNames.push_back(inputFileName);
59 void ECLDatabaseImporter::importDigitEnergyCalibration()
62 TClonesArray digitCalibrationConstants(
"Belle2::ECLDigitEnergyConstants");
68 for (
const string& inputFile : m_inputFileNames) {
70 TFile* f = TFile::Open(inputFile.c_str(),
"READ");
72 TIter next(f->GetListOfKeys());
74 while ((key = (TKey*) next())) {
76 string histconstants = key->GetName();
78 if (histconstants.compare(
"energy") == 0) {
79 energy = (TH1F*)f->Get(histconstants.c_str());
80 }
else if (histconstants.compare(
"amplitude") == 0) {
81 amplitude = (TH1F*)f->Get(histconstants.c_str());
84 else { B2FATAL(
"Key name does not match any of the following: energy, amplitude!"); }
90 if (nFiles != 1) { B2FATAL(
"Sorry, you must only import one file at a time for now!"); }
94 for (
int bin = 1; bin <= amplitude->GetNbinsX(); ++bin) {
95 float amplitudeval = amplitude->GetBinContent(bin);
96 float energyval = energy->GetBinContent(bin);
102 Database::Instance().storeData(m_name, &digitCalibrationConstants, iov);
105 void ECLDatabaseImporter::importDigitTimeCalibration()
108 TClonesArray digitCalibrationConstants(
"Belle2::ECLDigitTimeConstants");
113 for (
const string& inputFile : m_inputFileNames) {
115 TFile* f = TFile::Open(inputFile.c_str(),
"READ");
117 TIter next(f->GetListOfKeys());
119 while ((key = (TKey*) next())) {
121 string histconstants = key->GetName();
123 if (histconstants.compare(
"constantB") == 0) {
124 offset = (TH1F*)f->Get(histconstants.c_str());
125 }
else { B2FATAL(
"Key name does not match any of the following: constantC!"); }
131 if (nFiles != 1) { B2FATAL(
"Sorry, you must only import one file at a time for now!"); }
135 for (
int bin = 1; bin <= offset->GetNbinsX(); ++bin) {
136 float offsetval = offset->GetBinContent(bin);
142 Database::Instance().storeData(m_name, &digitCalibrationConstants, iov);
145 void ECLDatabaseImporter::importShowerCorrectorLeakageCorrections()
147 if (m_inputFileNames.size() > 1)
148 B2FATAL(
"Sorry, you must only import one file at a time for now!");
151 TFile* inputFile =
new TFile(m_inputFileNames[0].data(),
"READ");
153 if (!inputFile || inputFile->IsZombie())
154 B2FATAL(
"Could not open file " << m_inputFileNames[0]);
157 TTree* correctionTree = getRootObjectFromFile<TTree*>(inputFile,
"ParameterNtuple");
158 TTree* helperTree = getRootObjectFromFile<TTree*>(inputFile,
"ConstantNtuple");
164 int bgFractionBinNum;
169 float correctionFactor;
172 correctionTree->SetBranchAddress(m_bgFractionBinNumBranchName.c_str(), &bgFractionBinNum);
173 correctionTree->SetBranchAddress(m_regNumBranchName.c_str(), ®Num);
174 correctionTree->SetBranchAddress(m_phiBinNumBranchName.c_str(), &phiBinNum);
175 correctionTree->SetBranchAddress(m_thetaBinNumBranchName.c_str(), &thetaBinNum);
176 correctionTree->SetBranchAddress(m_energyBinNumBranchName.c_str(), &energyBinNum);
177 correctionTree->SetBranchAddress(m_correctionFactorBranchName.c_str(), &correctionFactor);
180 std::vector<int> m_bgFractionBinNum;
181 std::vector<int> m_regNum;
182 std::vector<int> m_phiBinNum;
183 std::vector<int> m_thetaBinNum;
184 std::vector<int> m_energyBinNum;
185 std::vector<float> m_correctionFactor;
187 for (
long iEntry = 0; iEntry < correctionTree->GetEntries(); ++iEntry) {
188 correctionTree->GetEntry(iEntry);
190 m_bgFractionBinNum.push_back(bgFractionBinNum);
191 m_regNum.push_back(regNum);
192 m_phiBinNum.push_back(phiBinNum);
193 m_thetaBinNum.push_back(thetaBinNum);
194 m_energyBinNum.push_back(energyBinNum);
195 m_correctionFactor.push_back(correctionFactor);
212 int numOfReg1ThetaBins;
213 int numOfReg2ThetaBins;
214 int numOfReg3ThetaBins;
221 std::vector<float> avgRecEns(m_numAvgRecEnEntries);
222 helperTree->SetBranchAddress(m_avgRecEnBranchName.c_str(), avgRecEns.data());
224 helperTree->SetBranchAddress(m_lReg1ThetaBranchName.c_str(), &lReg1Theta);
225 helperTree->SetBranchAddress(m_hReg1ThetaBranchName.c_str(), &hReg1Theta);
226 helperTree->SetBranchAddress(m_lReg2ThetaBranchName.c_str(), &lReg2Theta);
227 helperTree->SetBranchAddress(m_hReg2ThetaBranchName.c_str(), &hReg2Theta);
228 helperTree->SetBranchAddress(m_lReg3ThetaBranchName.c_str(), &lReg3Theta);
229 helperTree->SetBranchAddress(m_hReg3ThetaBranchName.c_str(), &hReg3Theta);
230 helperTree->SetBranchAddress(m_numOfBfBinsBranchName.c_str(), &numOfBfBins);
231 helperTree->SetBranchAddress(m_numOfEnergyBinsBranchName.c_str(), &numOfEnergyBins);
232 helperTree->SetBranchAddress(m_numOfPhiBinsBranchName.c_str(), &numOfPhiBins);
233 helperTree->SetBranchAddress(m_numOfReg1ThetaBinsBranchName.c_str(), &numOfReg1ThetaBins);
234 helperTree->SetBranchAddress(m_numOfReg2ThetaBinsBranchName.c_str(), &numOfReg2ThetaBins);
235 helperTree->SetBranchAddress(m_numOfReg3ThetaBinsBranchName.c_str(), &numOfReg3ThetaBins);
236 helperTree->SetBranchAddress(m_phiPeriodicityBranchName.c_str(), &phiPeriodicity);
239 std::vector<float> m_avgRecEn;
240 std::vector<float> m_lReg1Theta;
241 std::vector<float> m_hReg1Theta;
242 std::vector<float> m_lReg2Theta;
243 std::vector<float> m_hReg2Theta;
244 std::vector<float> m_lReg3Theta;
245 std::vector<float> m_hReg3Theta;
246 std::vector<int> m_numOfBfBins;
247 std::vector<int> m_numOfEnergyBins;
248 std::vector<int> m_numOfPhiBins;
249 std::vector<int> m_numOfReg1ThetaBins;
250 std::vector<int> m_numOfReg2ThetaBins;
251 std::vector<int> m_numOfReg3ThetaBins;
252 std::vector<int> m_phiPeriodicity;
254 for (
long iEntry = 0; iEntry < helperTree->GetEntries(); ++iEntry) {
255 helperTree->GetEntry(iEntry);
256 for (
unsigned int iIdx = 0; iIdx < avgRecEns.size(); ++iIdx) m_avgRecEn.push_back(avgRecEns[iIdx]);
258 m_lReg1Theta.push_back(lReg1Theta);
259 m_hReg1Theta.push_back(hReg1Theta);
260 m_lReg2Theta.push_back(lReg2Theta);
261 m_hReg2Theta.push_back(hReg2Theta);
262 m_lReg3Theta.push_back(lReg3Theta);
263 m_hReg3Theta.push_back(hReg3Theta);
264 m_numOfBfBins.push_back(numOfBfBins);
265 m_numOfEnergyBins.push_back(numOfEnergyBins);
266 m_numOfPhiBins.push_back(numOfPhiBins);
267 m_numOfReg1ThetaBins.push_back(numOfReg1ThetaBins);
268 m_numOfReg2ThetaBins.push_back(numOfReg2ThetaBins);
269 m_numOfReg3ThetaBins.push_back(numOfReg3ThetaBins);
270 m_phiPeriodicity.push_back(phiPeriodicity);
293 m_numOfReg1ThetaBins,
294 m_numOfReg2ThetaBins,
295 m_numOfReg3ThetaBins,
312 void ECLDatabaseImporter::importShowerShapesSecondMomentCorrections()
314 if (m_inputFileNames.size() > 1)
315 B2FATAL(
"Sorry, you must only import one file at a time for now!");
319 TFile* inputFile =
new TFile(m_inputFileNames[0].data(),
"READ");
321 if (!inputFile || inputFile->IsZombie())
322 B2FATAL(
"Could not open file " << m_inputFileNames[0]);
325 TGraph* theta_N1_graph = getRootObjectFromFile<TGraph*>(inputFile,
"SecondMomentCorrections_theta_N1");
326 dbArray.
appendNew(ECLShower::c_nPhotons, ECLShowerShapeModule::c_thetaType , *theta_N1_graph);
329 TGraph* phi_N1_graph = getRootObjectFromFile<TGraph*>(inputFile,
"SecondMomentCorrections_phi_N1");
330 dbArray.
appendNew(ECLShower::c_nPhotons, ECLShowerShapeModule::c_phiType , *phi_N1_graph);
333 TGraph* theta_N2_graph = getRootObjectFromFile<TGraph*>(inputFile,
"SecondMomentCorrections_theta_N2");
334 dbArray.
appendNew(ECLShower::c_neutralHadron, ECLShowerShapeModule::c_thetaType , *theta_N2_graph);
337 TGraph* phi_N2_graph = getRootObjectFromFile<TGraph*>(inputFile,
"SecondMomentCorrections_phi_N2");
338 dbArray.
appendNew(ECLShower::c_neutralHadron, ECLShowerShapeModule::c_phiType , *phi_N2_graph);
354 void ECLDatabaseImporter::importShowerEnergyCorrectionTemporary()
356 if (m_inputFileNames.size() > 1)
357 B2FATAL(
"Sorry, you must only import one file at a time for now!");
360 boost::filesystem::path path(m_inputFileNames[0]);
361 if (path.extension() !=
".root")
362 B2FATAL(
"Expecting a .root file. Aborting");
364 TFile* inputFile =
new TFile(m_inputFileNames[0].data(),
"READ");
366 TGraph2D* theta_geo_graph = getRootObjectFromFile<TGraph2D*>(inputFile,
"LeakageCorrections_theta_geometry");
367 TGraph2D* phi_geo_graph = getRootObjectFromFile<TGraph2D*>(inputFile,
"LeakageCorrections_phi_geometry");
368 TGraph2D* theta_en_graph = getRootObjectFromFile<TGraph2D*>(inputFile,
"LeakageCorrections_theta_energy");
369 TGraph2D* phi_en_graph = getRootObjectFromFile<TGraph2D*>(inputFile,
"LeakageCorrections_phi_energy");
370 TH1F* bg_histo = getRootObjectFromFile<TH1F*>(inputFile,
"LeakageCorrections_background_fraction");
372 double bkgFactor = bg_histo->GetBinContent(1);
374 double thetaMin = theta_en_graph->GetXmin();
375 double thetaMax = theta_en_graph->GetXmax();
376 double phiMin = phi_en_graph->GetXmin();
377 double phiMax = phi_en_graph->GetXmax();
379 double energyMin = theta_en_graph->GetYmin();
380 double energyMax = theta_en_graph->GetYmax();
383 B2DEBUG(28,
"Leakage DBobjects angle boundaries: thetaMin=" << thetaMin <<
" thetaMax=" << thetaMax <<
" phiMin= " << phiMin <<
384 " phiMax= " << phiMax <<
" enmin=" << energyMin <<
385 " enmax=" << energyMax);
394 if (std::abs(bkgFactor - 1.0) < 1e-9) {
397 dbPtr_theta_geo.
construct(*theta_geo_graph, thetaMin, thetaMax, energyMin, energyMax);
399 dbPtr_phi_geo.
construct(*phi_geo_graph, phiMin, phiMax, energyMin, energyMax);
401 dbPtr_theta_en.
construct(*theta_en_graph, thetaMin, thetaMax, energyMin, energyMax);
403 dbPtr_phi_en.
construct(*phi_en_graph, phiMin, phiMax, energyMin, energyMax);
406 dbPtr_theta_geo.
import(iov);
407 dbPtr_phi_geo.
import(iov);
408 dbPtr_theta_en.
import(iov);
411 if (std::abs(bkgFactor - 1.0) < 1e-9) {
414 dbPtr_theta_geo.
construct(*theta_geo_graph, thetaMin, thetaMax, energyMin, energyMax);
416 dbPtr_phi_geo.
construct(*phi_geo_graph, phiMin, phiMax, energyMin, energyMax);
418 dbPtr_theta_en.
construct(*theta_en_graph, thetaMin, thetaMax, energyMin, energyMax);
420 dbPtr_phi_en.
construct(*phi_en_graph, phiMin, phiMax, energyMin, energyMax);
423 dbPtr_theta_geo.
import(iov);
424 dbPtr_phi_geo.
import(iov);
425 dbPtr_theta_en.
import(iov);
432 void ECLDatabaseImporter::importTrackClusterMatchingThresholds()
434 if (m_inputFileNames.size() > 1)
435 B2FATAL(
"Sorry, you must only import one file at a time for now!");
438 boost::filesystem::path path(m_inputFileNames[0]);
439 if (path.extension() !=
".txt")
440 B2FATAL(
"Expecting a .txt file. Aborting");
442 vector<pair<double, double>> m_matchingThresholdPairsFWD;
443 vector<pair<double, double>> m_matchingThresholdPairsBWD;
444 vector<pair<double, pair<double, double>>> m_matchingThresholdPairsBRL;
445 pair<double, double> m_matchingThresholdPair;
446 pair<double, pair<double, double>> m_thetaMatchingThresholdPair;
447 double pt, threshold, thetalimit;
450 ifstream infile(m_inputFileNames[0]);
452 while (getline(infile, line)) {
453 istringstream iss(line);
455 if (eclregion ==
"FWD" || eclregion ==
"BWD") {
456 iss >> pt >> threshold;
457 m_matchingThresholdPair = make_pair(pt, threshold);
458 if (eclregion ==
"FWD") m_matchingThresholdPairsFWD.push_back(m_matchingThresholdPair);
459 else m_matchingThresholdPairsBWD.push_back(m_matchingThresholdPair);
460 }
else if (eclregion ==
"BRL") {
461 iss >> thetalimit >> pt >> threshold;
462 m_matchingThresholdPair = make_pair(pt, threshold);
463 m_thetaMatchingThresholdPair = make_pair(thetalimit, m_matchingThresholdPair);
464 m_matchingThresholdPairsBRL.push_back(m_thetaMatchingThresholdPair);
469 dbPtr.
construct(m_matchingThresholdPairsFWD, m_matchingThresholdPairsBWD, m_matchingThresholdPairsBRL);
477 void ECLDatabaseImporter::importTrackClusterMatchingParameterizations()
479 if (m_inputFileNames.size() > 1)
480 B2FATAL(
"Sorry, you must only import one file at a time for now!");
483 TFile* inputFile =
new TFile(m_inputFileNames[0].data(),
"READ");
485 if (!inputFile || inputFile->IsZombie())
486 B2FATAL(
"Could not open file " << m_inputFileNames[0]);
488 map<string, TF1> m_parametrizationFunctions;
489 vector<string> angles = {
"Theta",
"Phi"};
490 vector<string> regions = {
"BRL",
"BWD",
"FWD"};
491 vector<string> hittypes = {
"CROSS",
"DL",
"NEAR"};
493 for (
const auto& angle : angles) {
494 for (
const auto& region : regions) {
495 for (
const auto& hittype : hittypes) {
496 m_parametrizationFunctions.insert(make_pair(angle + region + hittype, *(getRootObjectFromFile<TF1*>(inputFile,
497 "RMSParameterization" + angle + region + hittype))));
503 dbPtr.
construct(m_parametrizationFunctions);