22 #include <framework/database/DBImportObjPtr.h>
23 #include <framework/database/DBObjPtr.h>
24 #include <framework/database/DBStore.h>
25 #include <framework/datastore/StoreObjPtr.h>
26 #include <framework/datastore/DataStore.h>
27 #include <framework/dataobjects/EventMetaData.h>
28 #include <framework/database/Configuration.h>
29 #include <ecl/dbobjects/ECLCrystalCalib.h>
38 void setupDatabase(
int exp,
int run,
int eventNr = 1)
43 evtPtr.registerInDataStore();
57 int main(
int argc,
char** argv)
61 if (argc < 3 || argc > 5) {
62 std::cout <<
"incorrect number of arguments for eclCrystalEnergy" << std::endl;
65 int experiment = std::stoi(argv[1]);
66 int run = std::stoi(argv[2]);
67 std::string gtName =
"ECL_energy_calibrations";
68 if (argc >= 4) {gtName = argv[3];}
69 bool writeOutput =
true;
70 if (argc == 5) {writeOutput = std::stoi(argv[4]);}
72 std::cout <<
"eclCrystalEnergy called with arguments exp = " << experiment <<
" run = " << run <<
" GT = " << gtName <<
73 " writeOutput = " << writeOutput << std::endl;
79 conf.prependGlobalTag(gtName);
80 conf.prependTestingPayloadLocation(
"localdb/database.txt");
87 setupDatabase(experiment, run);
91 std::vector<float> GammaGammaCalib;
92 std::vector<float> GammaGammaCalibUnc;
93 GammaGammaCalib = GammaGamma->getCalibVector();
94 GammaGammaCalibUnc = GammaGamma->getCalibUncVector();
96 std::vector<float> ExistingCalib;
97 std::vector<float> ExistingCalibUnc;
98 ExistingCalib = Existing->getCalibVector();
99 ExistingCalibUnc = Existing->getCalibUncVector();
104 std::vector<float> NewCalib;
105 std::vector<float> NewCalibUnc;
106 NewCalib.resize(8736);
107 NewCalibUnc.resize(8736);
108 for (
int ic = 0; ic < 8736; ic++) {
109 if (GammaGammaCalib[ic] > 0.) {
110 NewCalib[ic] = GammaGammaCalib[ic];
111 NewCalibUnc[ic] = GammaGammaCalibUnc[ic];
113 NewCalib[ic] = ExistingCalib[ic];
114 NewCalibUnc[ic] = ExistingCalibUnc[ic];
120 for (
int ic = 0; ic < 9000; ic += 1000) {
121 std::cout << ic + 1 <<
" " << ExistingCalib[ic] <<
" " << ExistingCalibUnc[ic] <<
" " << GammaGammaCalib[ic] <<
" " <<
122 GammaGammaCalibUnc[ic] <<
" " << NewCalib[ic] <<
" " << NewCalibUnc[ic] << std::endl;
127 TString payloadTitle =
"ECLCrystalEnergy";
129 payloadTitle += experiment;
132 TString fname = payloadTitle;
134 TFile hfile(fname,
"recreate");
136 TString htitle = payloadTitle;
137 htitle +=
" existing values;cellID";
138 TH1F* existingPayload =
new TH1F(
"existingPayload", htitle, 8736, 1, 8737);
140 htitle = payloadTitle;
141 htitle +=
" new values;cellID";
142 TH1F* newPayload =
new TH1F(
"newPayload", htitle, 8736, 1, 8737);
144 htitle = payloadTitle;
145 htitle +=
" ratio new/old;cellID";
146 TH1F* payloadRatioVsCellID =
new TH1F(
"payloadRatioVsCellID", htitle, 8736, 1, 8737);
148 htitle = payloadTitle;
149 htitle +=
" ratio new/old";
150 TH1F* payloadRatio =
new TH1F(
"payloadRatio", htitle, 200, 0.95, 1.05);
152 for (
int cellID = 1; cellID <= 8736; cellID++) {
153 existingPayload->SetBinContent(cellID, ExistingCalib[cellID - 1]);
154 existingPayload->SetBinError(cellID, ExistingCalibUnc[cellID - 1]);
156 newPayload->SetBinContent(cellID, NewCalib[cellID - 1]);
157 newPayload->SetBinError(cellID, NewCalibUnc[cellID - 1]);
161 if (abs(ExistingCalib[cellID - 1]) > 1.0e-12) {
162 ratio = NewCalib[cellID - 1] / ExistingCalib[cellID - 1];
163 float rUnc0 = ExistingCalibUnc[cellID - 1] / ExistingCalib[cellID - 1];
165 if (abs(NewCalib[cellID - 1]) > 1.0e-12) {rUnc1 = NewCalibUnc[cellID - 1] / NewCalib[cellID - 1];}
166 ratioUnc = ratio * sqrt(rUnc0 * rUnc0 + rUnc1 * rUnc1);
169 payloadRatioVsCellID->SetBinContent(cellID, ratio);
170 payloadRatioVsCellID->SetBinError(cellID, ratioUnc);
172 payloadRatio->Fill(ratio);
178 std::cout << std::endl <<
"Values written to " << fname << std::endl;
183 std::cout <<
"Creating importer" << std::endl;
185 importer.construct();
186 importer->setCalibVector(NewCalib, NewCalibUnc);
188 std::cout <<
"Successfully wrote payload ECLCrystalEnergy with iov " << experiment <<
"," << run <<
",-1,-1" << std::endl;
190 std::cout <<
"Payload not written to database" << std::endl;